Nonequilibrium Statistical Physics in Ecology: Vegetation Patterns, Animal Mobility and Temporal Fluctuations

TL;DR

This thesis applies nonequilibrium statistical physics tools to ecological modeling, exploring vegetation patterns, animal movement, and environmental fluctuations, revealing insights into ecosystem dynamics and collective behaviors.

Contribution

It introduces novel mathematical models integrating nonequilibrium physics concepts into ecology, addressing vegetation pattern formation, animal group behavior, and environmental stochasticity.

Findings

Vegetation patterns are explained by nonequilibrium models under scarce resources.

Animal communication influences group formation and search efficiency.

Environmental fluctuations significantly impact ecosystem stability and interactions.

Abstract

This thesis focuses on the applications of mathematical tools and concepts brought from nonequilibrium statistical physics to the modeling of ecological problems. The first part provides a short introduction where the theoretical concepts and mathematical tools that are going to be used in subsequent chapters are presented. Firstly, the different levels of description usually employed in the models are explained. Secondly, the mathematical relationships among them are presented. Finally, the notation and terminology that will be used later on are explained. The second part is devoted to studying vegetation pattern formation in regions where precipitations are not frequent and resources for plant growth are scarce. This part comprises two chapters. The third part of the thesis develops a series of mathematical models describing the collective movement and behavior of some animal species. Its primary objective is to investigate the effect that communication among foragers has on searching times and the formation of groups. It consists of two chapters. The fourth part covers the effect of stochastic temporal disorder, mimicking climate and environmental variability, on systems formed by many interacting particles. These models may serve as an example of ecosystems. The thesis ends with a summary and devising future research lines.

Equations

$$P_c(t+\Delta t)=\left(1-\sum _{c^{\prime }}{\omega }_{c\to c^{\prime }}\Delta t\right)P_c(t)+\sum _{c^{\prime }}{\omega }_{c^{\prime }\to c}\Delta t{P}_{c^{\prime }}(t),$$

$$\frac{\partial P_c(t)}{\partial t}=\sum _{c^{\prime }}{\omega }_{c^{\prime }\to c}{P}_{c^{\prime }}(t)-\sum _{c^{\prime }}{\omega }_{c\to c^{\prime }}{P}_c(t).$$

$$P(n,p,t+dt)$$

$$\frac{\partial P(n,p;t)}{\partial t}$$

$${\Omega }_i=\sum _{j\ne i}{\omega }_{i\to j},i=1,\dots ,M$$

$$dt=\frac{-\ln u_0}{{\Omega }_i}.$$

$$p_{i\to j}=\frac{{\omega }_{i\to j}}{{\Omega }_i}$$

$$⟨n(t)⟩=\sum _{p,n=0}^{\infty }nP(n,p;t),⟨p(t)⟩=\sum _{p,n=0}^{\infty }pP(n,p;t).$$

$$\frac{d}{dt}⟨n(t)⟩=k_b⟨n(t)⟩-k_p⟨n(t)p(t)⟩$$

$$\frac{d}{dt}⟨p(t)⟩=k_p⟨n(t)p(t)⟩-k_d⟨p(t)⟩.$$

$$\frac{dN}{dt}$$

$$\frac{dP}{dt}$$

$$u(\tau )=\frac{k_{p}N}{k_d},v(\tau )=\frac{k_{p}P}{k_b},\tau =k_{b}t,\alpha =\frac{k_d}{k_b},$$

$$\frac{du}{dt}=u(1-v),$$

$$\frac{dv}{dt}=\alpha v(u-1).$$

$$\frac{\partial P_c(t)}{\partial t}=\int \left[{\omega }_{c^{\prime }\to c}{P}_{c^{\prime }}(t)-{\omega }_{c\to c^{\prime }}{P}_c(t)\right]d{c}^{\prime }.$$

$${\omega }_{c\to c^{\prime }}=\omega (c;r),$$

$$\frac{\partial P_c(t)}{\partial t}=\int \omega (c-r;r)P_{c-r}(t)dr-P_c(t)\int \omega (c;-r)dr.$$

$$\displaystyle\omega(c^{\prime};r)\approx 0\ \ \ \$$

$$\displaystyle\omega(c^{\prime}+\Delta c;r)\approx\omega(c^{\prime};r)\ \ \ \$$

$$\frac{\partial P_c(t)}{\partial t}$$

$${\alpha }_{\nu }(c)={\int }_{-\infty }^{+\infty }{r}^{\nu }\omega (c;r)dr,$$

$$\frac{\partial P_c(t)}{\partial t}=-\frac{\partial }{\partial c}[{\alpha }_1(c)P_c(t)]+\frac{1}{2}\frac{{\partial }^2}{\partial c^2}[{\alpha }_2(c)P_c(t)].$$

$$\frac{dc}{dt}=f(c,t)+g(c,t)\eta (t),$$

$$⟨\eta (t)\eta (t^{\prime })⟩=\Gamma \delta (t-t^{\prime }),$$

$$c(t)-c(0)={\int }_0^{t}f[c(s),s]ds+{\int }_0^{t}g[c(s),s]\eta (s)ds,$$

$$dW(t)\equiv W(t+dt)-W(t)=\eta (t)dt.$$

$$c(t)-c(0)={\int }_0^{t}f[c(s),s]ds+{\int }_0^{t}g[c(s),s]dW(s),$$

$$\frac{\partial \varphi (\mathbf{r},t)}{\partial t}=f(\varphi (\mathbf{r},t),t)+{\nabla }^2\varphi (\mathbf{r},t)+g(\varphi (\mathbf{r},t),t)\eta (\mathbf{r},t).$$

$$\frac{\partial P(c,t)}{\partial t}=-\frac{\partial }{\partial c}{\alpha }_1(c)P(c,t)+\frac{1}{2}\frac{{\partial }^2}{\partial c^2}{\alpha }_2(c)P(c,t),$$

$$\frac{dc}{dt}=f(c,t)+g(c,t)\eta (t),$$

$$f(c,t)$$

Full text

Nonequilibrium Statistical Physics in Ecology: Vegetation Patterns, Animal Mobility and Temporal Fluctuations

Ricardo Martínez García

**Nonequilibrium Statistical Physics in Ecology: Vegetation Patterns, Animal Mobility and Temporal Fluctuations

**Ricardo Martínez García

Tesis realizada en el Instituto de Física Interdisciplinar y Sistemas Complejos, IFISC (CSIC-UIB).

Presentada en el Departamento de Física de la Universitat de les Illes Balears.

**For an updated version of this thesis please contact:

ricardo@ifisc.uib-csic.es or rmtzgarcia@gmail.com

**PhD Thesis

Director: Dr. Cristóbal López

Palma de Mallorca, 2 de Mayo de 2014.

The cover shows two aerial photographies of vegetation patterns taken from von Hardenberg et al., [2001], and a Mongolian gazelle offspring taken from www.lhnet.org.

Cristóbal López Sánchez, Profesor Titular de Universidad

CERTIFICA:

que esta tesis doctoral ha sido realizada por el doctorando Sr. Ricardo Martínez García bajo su dirección en el Instituto de Física Interdisciplinar y Sistemas Complejos y, para que conste, firma la presente

Director:

Dr. Cristóbal López Sánchez

Doctorando:

Ricardo Martínez García

Palma de Mallorca,                                                                            2 de Mayo de 2014.

A mis padres y hermano.

*Look up to the sky.

You will never find rainbows

if you are looking down.

*Charles Chaplin. *

Agradecimientos

En primer lugar, quisiera agradecer a Cristóbal el haberme dado la oportunidad de hacer esta tesis doctoral. Gracias por todo lo aprendido, pero sobre todo por tu buen humor y continuo apoyo. Gracias también por darme la libertad necesaria para equivocarme y tener mis propias ideas, o al menos intentarlo. I also want to thank Justin Calabrese who I have considered my second PhD supervisor since my stay in Front Royal. Thanks for your friendliness and for teaching me so many things. No quisiera olvidarme de Emilio Hernández-García y Federico Vázquez, siempre habéis tenido un rato para sentaros conmigo. Gracias especialmente a tí, Fede, por tantas horas juntos delante del ordenador buscando power laws. Finally, thank you Miguel Ángel Muñoz, Chris Fleming, Thomas Mueller and Kirk Olson for spending part of your time on me.

Empezando por el principio, quiero dedicar unas palabras a Juan Francisco, por dedicarme tanto tiempo en el colegio y por despertar el interés por la física en mí. También a Daniel Alonso, por ayudarme a recuperarlo, y a Santiago Brouard, por iniciarme en el mundo de la investigación en La Laguna.

En el IFISC he encontrado el ambiente ideal para disfrutar de esta tesis. Gracias a todos por hacerlo posible. En particular a Marta, Inma, Rosa Campomar y Rosa Rodríguez por hacer de la burocracia algo más sencillo. A Edu, Rubén, Antònia y David, porque sin su ayuda seguramente esta trabajo no habría ido para delante. Espero no haber dado mucho la lata, y si lo he hecho, os debo un desayuno. Siempre estaré agradecido a todos mis compañeros en la profundidad de la S07. Al Dr. Luis F. Lafuerza, por su acogida durante mi primer verano y por estar siempre dispuesto a echarme una mano. A Luca: Calviá y Milán nunca estuvieron tan cerca. A Pablo y Víctor, por nuestras conversaciones sobre arte precolombino y vuestro sentido del humor, y a Enrico y Simone, por las clases gratis de italiano. Muchas gracias también a todos los demás: Przemyslaw (espero haberlo escrito bien), Miguel Ángel, Pedro, Adrián, Julián, Neus, Juan, Xavi, Marie, Ismael, Leo, Alejandro, Toni Pérez…

Part of this Thesis has been done in other institutions. I have really enjoyed these experiences, but they would not have been the same without all the nice people that I could meet. Thanks to Nat, Ben, Fan, Meng, Tuya, Leah, Caroline, Bettina, Christian, Jan… for the time in Front Royal, y gracias a Pablo por alguna escapada por Dresden.

A nivel institucional, gracias al CSIC y a la Univeristat de les Illes Balears, pues sin sus fondos nunca podría haber completado este doctorado, y a los proyectos FISICOS e INTENSE@COSYP.

Estos años en Mallorca han sido mucho más que años de trabajo. A través del fútbol sala he podido conocer a muchísima gente. Gracias a todos, a los del Bahia’s y Ciutat de Mallorca, pero especialmente a mis compañeros de La Salle Pont d’Inca. Habéis sido una familia desde el primer día. Especialmente, me gustaría acordarme de Alberto, por abrirme las puertas y por tener siempre un rato para liarse. A los demás, ¿qué pasa …? Nombraros a todos ocuparía mucho espacio, pero lo voy a hacer: Bily, Biel, Raúl, Xevi, Héctor, Colo, Edu, Toni, Roberto, Berni, Baia (gran futuro delante de las cámaras), Toni Mir, Pitu, Majoni, Rafa, Perillas y Bernat. Muchísimas gracias también a todos los niños que he tenido la inmensa fortuna de entrenar. Sobre todo a los más pequeños, de quienes tenemos muchísimo que aprender.

He tenido la gran suerte de vivir con grandes compañeros. Miguel, uno de los mejores amigos que se pueden tener y probablemente la persona con mayor visión de futuro. Mario, importador a Mallorca de la última tecnología riojana: el motorabo. ¡Dale qué suene! Siempre. Juntos hemos compartido muchos de los mejores momentos. Gracias. Gracias Àngel por enseñarme la isla y por nuestros momentos de pesca. No me olvido de vosotros, Luis, Gloria, Kike, Laura y Adrián porque por vuestra culpa llegar a Mallorca fue fácil e irse será difícil.

Estas últimas líneas van para la gente que más quiero. A mis amigos de Tenerife, los de siempre. A Paula, por ser como eres, por tu sonrisa, y por acompañarme. A Guillermo, Bárbara y Miki, por preocuparos por mi. Pero sobre todo dedico este trabajo a mis padres, por ser el mejor ejemplo y por creer y confiar en mi. Ya que no os lo recuerdo muy a menudo, aprovecho estas palabras para deciros lo mucho que os quiero. A mi hermano, el mejor compañero y un ejemplo como científico. Gracias por mis primeros trabajos de campo a los 2 años, poca gente tiene la suerte de saber qué es un coleóptero antes de ir a la guardería111Orden de insectos masticadores que poseen un caparazón duro y dos alas, también duras, que cubren a su vez dos alas membranosas. Por último, a mis abuelos, a los que ya no están, Alejandro, Rosa e Isidro y a Cuquita, por poder disfrutarte cada día.

Resumen

Esta tesis doctoral se centra en la aplicación de técnicas propias de la física estadística del no equilibrio al estudio de problemas con trasfondo ecológico.

En la primera parte se presenta una breve introducción con el fin de contextualizar el uso de modelos cuantitativos en el estudio de problemas ecológicos. Para ello, se revisan los fundamentos teóricos y las herramientas matemáticas utilizadas en los trabajos que ocupan los capítulos siguientes. En primer lugar, se explican las distintas maneras de describir matemáticamente este tipo de sistemas, estableciendo relaciones entre ellas y explicando las ventajas e incovenientes que presenta cada una. En esta sección también se introducen la terminología y la notación que se emplearán más adelante.

En la segunda parte se comienzan a presentar resultados originales. Se estudia la formación de patrones de vegetación en sistemas en los que el agua es un factor que limita la aparición de nuevas plantas. Esta parte se divide en dos capítulos.

• •

El primero se centra en el caso particular de sabanas mésicas, con una precipitación media anual intermedia, y en las que los árboles coexisten con otros tipos de vegetación más baja (arbustos y hierbas). Se presenta un modelo en el que se incluyen los efectos de la competición por recursos y la presencia de incendios. En este último caso, la protección que los árboles adultos proporcionan a los jóvenes contra el fuego supone una interacción de facilitación a muy corto alcance entre la vegetación. El principal resultado de este estudio concluye que, incluso en el límite en el que los mecanismos facilitativos tienen un alcance muy corto (local), aparecen patrones en el sistema. Finalmente, incluyendo la naturaleza estocástica de la dinámica de nacimiento y muerte de los árboles se recuperan estructuras con formas más parecidas a las observadas en sabanas reales.

• •

El segundo capítulo de esta parte estudia la formación de patrones en sistemas áridos, cuyas formas son mucho más regulares que en las sabanas mésicas. Además, las precipitaciones también son más escasas. El origen de estas estructuras se atribuye tradicionalmente a la presencia de diferentes interacciones entre las plantas que actúan en distintas escalas espaciales. En particular, muchos de los trabajos previos defienden que se deben a la combinación de mecanismos que facilitan el crecimiento de vegetación a corto alcance (facilitación) con otros, de mayor alcance, que lo inhiben (competición). En este capítulo se presentan modelos en los que únicamente se incluyen interacciones competitivas, a pesar de lo cual se recupera la secuencia típica de patrones obtenida en modelos previos. Se introduce el concepto de zonas de exclusión como mecanismo biológico responsable de la formación de patrones.

En la tercera parte de la tesis se presentan modelos para el estudio del movimiento y comportamiento colectivo de animales. En concreto, se investiga la influencia que tiene la comunicación entre individuos en los procesos de búsqueda que estos llevan a cabo, con especial énfasis en la búsqueda de recursos. Consta de dos capítulos.

• •

En primer lugar, se analiza desde un punto de vista teórico la influencia de la comunicación en los tiempos de búsqueda. En general, comunicaciones a escalas intermedias resultan en tiempos de búsqueda menores, mientras que alcances más cortos o más largos proporcionan una cantidad de información insuficiente o excesiva al resto de la población. Esto impide a los individuos decidir correctamente en qué dirección moverse, lo cual da lugar a tiempos de búsqueda mayores. El capítulo se completa estudiando la influencia que tiene el tipo de movimiento de los individuos (browniano o Lévy) en los resultados del modelo.

• •

Esta parte finaliza presentando una aplicación del modelo desarrollado en el capítulo anterior al caso de las gacelas que habitan las estepas centroasiáticas (Procapra gutturosa). En los últimos años, se ha observado un gran decrecimiento en la población de esta especie. Esto se debe a la caza masiva de estos animales y a una pérdida y fragmentación de su hábitat provocada por la acción del hombre. Conocer sus hábitos migratorios y comportamiento resulta, por tanto, fundamental para desarrollar estrategias de conservación eficientes. En particular, en este capítulo se estudia la búsqueda de pastos por parte de estas gacelas, utilizando mapas reales de vegetación y medidas GPS del posicionamiento de un grupo de individuos. Se presta especial atención al efecto de la comunicación vocal entre animales, midiendo la eficiencia de la búsqueda en términos de su duración y de la formación de grupos en las zonas más ricas en recursos. Las gacelas encuentran buenos pastos de una manera óptima cuando se comunican emitiendo sonidos cuyas frecuencias coinciden con las obtenidas en medidas reales hechas en grabaciones de estos animales. Este resultado sugiere la posibilidad de que a lo largo de su evolución la gacela Procapra gutturosa haya optimizado su tracto vocal para facilitar la comunicación en la estepa.

En la cuarta parte, que consta de un único capítulo, se analiza el efecto que tiene un medio externo cuyas propiedades cambian estocásticamente en el tiempo sobre diferentes propiedades de un sistema compuesto por muchas partículas que interaccionan entre sí. Se estudian los tiempos de paso cuando el parámetro de control del problema fluctúa en torno a un valor medio. Se encuentra una región finita del diagrama de fases en la cual los tiempos escalan como una ley de potencia con el tamaño del sistema. Este resultado es contrario al caso puro, en el que el parámetro de control es constante y esto únicamente ocurre en el punto critico. Con estos resultados se extiende el concepto de Fases Temporales de Griffiths a un mayor número de sistemas.

La tesis termina con las conclusiones del trabajo y señalando posibles líneas de investigación que toman como punto de partida los resultados obtenidos.

Abstract

This thesis focuses on the applications of mathematical tools and concepts brought from nonequilibrium statistical physics to the modeling of ecological problems.

The first part provides a short introduction where the theoretical concepts and mathematical tools that are going to be used in subsequent chapters are presented. Firstly, the different levels of description usually employed in the models are explained. Secondly, the mathematical relationships among them are presented. Finally, the notation and terminology that will be used later on are explained.

The second part is devoted to studying vegetation pattern formation in regions where precipitations are not frequent and resources for plant growth are scarce. This part comprises two chapters.

• •

The first one studies the case of mesic savannas. These systems are characterized by receiving an intermediate amount of water and by a long term coexistence of layer made of grass and shrubs interspersed with irregular clusters of trees. A minimalistic model considering only long range competition among plants and the effect of possible fires is presented. In this later case, adult trees protect the growth of juvenile individuals against the fires by surrounding them and creating an antifire shell. This introduces a local facilitation effect for the establishment of new trees. Despite the range of facilitative interactions is taken to its infinitesimally short limit, the spectrum of patterns obtained in models with competitive and facilitative nonlocal interactions is recovered. Finally, considering the stochasticity in the birth and death dynamics of trees, the shapes of the structures reproduce the irregularity observed in aerial photographs of mesic savannas.

• •

The second chapter investigates the formation of patterns in arid regions, that are typically more regular than in mesic savannas. Previous studies attribute the origin of these structures to the existence of competitive and facilitative interactions among plants acting simultaneously but at different spatial scales. More precisely, to the combination of a short-range facilitation and a long-range competition (scale-depedent feedback). The findings of this chapter are based on the study of a theoretical model that assumes only long-range competitive interactions and shows the existence of vegetation patterns even under these conditions. This result suggests that the role of facilitative interactions could be superfluous in the development of these spatial structures. The biological concept of exclusion areas is proposed as an alternative to conventional scale-dependent feedback.

The third part of the thesis develops a series of mathematical models describing the collective movement and behavior of some animal species. Its primary objective is to investigate the effect that communication among foragers has on searching times and the formation of groups. It consists of two chapters:

• •

In the first one, the model is established and its properties studied from a theoretical point of view. The main novelty of this work is the inclusion of communication among searchers to share information about the location of the targets. Communication and amount of shared information are directly connected through the range of the signals emitted by successful searchers. In this context, searching processes are optimized in terms of duration when the individuals share intermediate amounts of information, corresponding to mid-range communication. Both a lack and an excess of information may worsen the search. The first implies an almost uninformed search, while the latter causes a loss in the directionality of the movement since individuals are overwhelmed with information coming from many targets. Finally, the influence of the type of movement on the search efficiency is investigated, comparing the Brownian and Lévy cases. Some analytical approximations and a continuum description of the model are also presented.

• •

This part ends with an application of the previous model to the foraging behavior of Mongolian gazelles (Procapra gutturosa). The population of this species has decreased in the last century because of massive hunting and a progressive habitat degradation and fragmentation caused by human disturbances in the Eastern steppe of Mongolia. Studying their mobility patterns and social behavior improves the development of conservation strategies. This chapter suggests possible searching strategies used by these animals to increase their forage encounters rate. The study is supported by the use of real vegetation maps based on satellite imagery and GPS data tracking the position of a group of gazelles. The main focus is on the effect that nonlocal vocal communication among individuals has on foraging times and group formation in the areas with better resources. According to the results of the model, the searching time is minimized when the communication takes place at a frequency that agrees with measurements made in gazelle’s acoustic signals. This suggests that, through its evolution, Procapra gutturosa may have optimized its vocal tract in order to facilitate the communication in the steppe.

The fourth part covers the effect of stochastic temporal disorder, mimicking climate and environmental variability, on systems formed by many interacting particles. These models may serve as an example of ecosystems. The temporal disorder is implemented making the control parameter fluctuating around a mean value close to the critical point. The effect of this external variability is quantified using passage times. The results show a change in the behavior of this magnitude compared with the pure case, that is, in the absence of external fluctuations. Within a finite region of the phase diagram, close to the critical point, the passage times scale as a power law with continuously varying exponent. In the pure model this behavior is only observed at the critical point. After these results, the concept of Temporal Griffiths Phases, introduced in the spreading of epidemics, is extended to a vast range of models.

The thesis ends with a summary and devising future research lines.

Contents

1. Resumen

2. Abstract

3. I Introduction

4. 1 Methods and tools

5. 1.1 From Individual Based to Population Level Models

6. 1.1.1 The Master equation

7. 1.1.2 The Fokker-Planck equation

8. 1.1.3 The Langevin equation

9. 1.2 Linear stability analysis

10. 1.3 First-passage times processes

11. II Vegetation Patterns

12. 2 Mesic savannas

13. 2.1 Introduction

14. 2.2 The deterministic description

15. 2.2.1 The nonlocal savanna model

16. 2.2.2 Linear stability analysis

17. 2.2.3 Numerical simulations

18. 2.3 Stochastic model

19. 2.4 Discussion

20. 2.5 Summary

21. 3 Semiarid systems

22. 3.1 Introduction

23. 3.2 Competition in a nonlocal nonlinear birth term

24. 3.3 Competition in a nonlocal nonlinear death term

25. 3.4 Competition in a nonlocal linear death term

26. 3.5 Summary and conclusions

27. A Linear stability analysis

28. B Numerical integration of Eq. (2.19)

29. C Derivation of the effective nonlocal description from tree-water dynamics

30. III Animal mobility

31. 4 Optimal search in interacting populations

32. 4.1 Introduction

33. 4.2 The Individual Based Model for Brownian searchers

34. 4.3 Lévy flights

35. 4.4 One-dimensional analytical approximations

36. 4.4.1 Brownian motion

37. 4.4.2 Lévy flights.

38. 4.5 Continuum approximation

39. 4.6 Brownian jumps vs. Lévy flights

40. 4.7 Summary and conclusions

41. 5 Foraging in * Procapra gutturosa*

42. 5.1 Introduction

43. 5.2 The model for acoustic communication

44. 5.3 Results and discussion

45. A Derivation of the macroscopic Eq. (4.3)

46. B Voronoi diagrams of the model.

47. IV Temporal fluctuations

48. 3 Temporal disorder in up-down symmetric systems

49. 3.1 Introduction

50. 3.2 Spatial disorder. Rare regions and the Griffiths Phase.

51. 3.3 Mean-field theory of $Z_{2}$-symmetric models with temporal disorder

52. 3.4 Ising transition with temporal disorder

53. 3.4.1 The Langevin equation

54. 3.4.2 Numerical results

55. 3.4.3 Analytical results

56. 3.5 Generalized Voter transition with temporal disorder

57. 3.5.1 The Langevin equation

58. 3.5.2 Numerical Results

59. 3.5.3 Analytical results

60. 3.6 Summary and conclusions

61. A Itô-Stratonovich discussion.

62. A.1 Stochastic integration.

63. A.2 Itô’s formula.

64. A.3 From Stratonovich to Itô.

65. A.4 Stratonovich / Itô dilemma.

66. B Analytical calculations on the escape time for the Ising Model

67. B.1 Case $\alpha\neq 1$

68. B.1.1 $\alpha<1$

69. B.1.2 $\alpha>1$

70. B.2 Case $\alpha=1$ Critical point.

71. V Conclusions and outlook

72. 3 Conclusions and outlook

Part I  **Introduction

**

Statistical physics focuses on the study of those systems that comprise a large number of simple components. Regardless of the particular nature of these fundamental entities, it describes the interactions among them and the global properties that appear at a macroscopic scale. These emergent phenomena are the hallmark of complex systems. Such systems are used to model processes in several disciplines, most of the times, far from the physical sciences. That’s why, during the last few years, statistical physics has become a powerful cross disciplinary tool, supplying the theoretical framework and the mathematical techniques that allow the study of many different problems in biology, economics or sociology. It provides a scenario that makes possible to encapsulate the huge number of microscopic degrees of freedom of a complex system into just a few collective variables.

On the other hand, ecology is concerned with the study of the relationships between organisms and their environment. In terms of this thesis, it is a paradigmatic example of complexity science. Ecological systems are formed by a huge number of heterogeneous constituents that interact and evolve stochastically in time. In addition, they are subject to changes and fluctuations in the surroundings that apart them from equilibrium222Here equilibrium refers to the thermodynamic equilibrium. It is a state of balance characterized by the absence of fluxes and currents in the system¡..

Because of this complex nature, ecology was originally an empirical science with purely descriptive purposes. Ancient Greek philosophers such as Hippocrates and Aristotles laid the foundations of ecology in their studies on natural history. However, over the years, the need for a mathematical formalism to tighten all the observations increased, until ecology adopted a more analytical approach in the late 19th century. The first models attracted the attention of many physicist and mathematicians that started developing new techniques and tools. Nowadays, theoretical ecology is a well established discipline that deals with several topics related not only with environmental conservation but also with evolutionary biology, ethology and genetics. It constitutes, together with recent technological advances, a potent instrument to better understand the natural environment.

Ecological systems show characteristic variability on a range of spatial, temporal and organizational scales [Levin, , 1992]. However, when we observe them, we do it in a limited range. Theoretical studies aim to comprehend how information is transferred from one level to other. They permit understanding natural phenomena in terms of the processes that govern them, and consequently develop management strategies. Without this knowledge, each stress must be evaluated separately in every system, and it would not be possible to extrapolate the knowledge obtained from one situation to another. But, what is the role of statistical physics in this task? On the one hand, most ecological systems can exhibit multistability, abrupt transitions, patterns or self-organization when a control parameter is varied. These concepts are characteristic of nonlinear systems, that have been traditionally studied by statistical physicists. Particularly interesting are those cases in which the dynamics at one level of organization can be understood as a consequence of the collective behavior of multiple similar identities. This reminds the definition of the systems that are the focus of statistical physics, which serves for developing simple models that retain and condense the essential information, omitting unnecessary details.

There is a large list of recent developments that may serve as examples of this relationship [Fort, , 2013]: collective animal movement [Cavagna et al., , 2010], demographic stochasticity in multiple species systems [McKane and Newman, , 2005; Butler and Goldenfeld, , 2009], evolutionary theory [Chia and Goldenfeld, , 2011], population genetics [Vladar and Barton, , 2011], species distribution [Harte et al., , 2008; Volkov et al., , 2003], complex ecological networks [Montoya et al., , 2006; Bastolla et al., , 2009], animal foraging [Méndez et al., , 2014; Viswanathan et al., , 2011], or species invasion [Seebens et al., , 2013]. In this thesis, I will abroad different problems within the framework of statistical physics, in particular vegetation pattern formation, animal behavior and ecosystem’s robustness. It is important to remark the diverse nature of each of these systems. Plants are inert, and so the development of patterns is a consequence of the interaction with the environment and the birth-death dynamics. On the other hand, animals usually show large migratory displacements and tend to form groups of individuals by coming together. Gathering these problems, the objective of this dissertation is to emphasize the connection between statistical physics and environmental sciences and its role in the establishment of ecological models.

The powerful of statistical physics as a cross disciplinary tool allows to tackle different questions depending on the particularities of each system. Here we wonder how external variability affects robustness and evolution of ecosystems and the mean lifetime of the species. We are also interested in disentangling the different facilitative and competitive interaction among plants in vegetation systems to unveil its role in the formation of patterns. Are both needed to maintain these regular structures? How efficient are inhomogeneous distributions of vegetation to avoid desertification in water-limited systems? Finally, we will try to shed light on the relationship between communication and foraging efficiency. This is one of the less investigated topics in the study of searching strategies. How can different communication mechanisms affect searching processes? Is the mean searching time a good metrics to quantify search efficiency? Does it exist an optimal communication range that accelerates the search? How does sharing information affect the collective use of a heterogeneous landscape? Answering these and other issues will be the goal of this thesis.

The results of each chapter can be found in the following publications:

• •

Chapter 2:

• –

R. Martínez-García, J.M. Calabrese, and C. López, (2013), Spatial patterns in mesic savannas: the local facilitation limit and the role of demographic stochasticity, Journal of Theoretical Biology, 333, 156-165.

• •

Chapter 3:

• –

R. Martínez-García, J.M. Calabrese, E. Hernández-García and C. López, (2013), Vegetation pattern formation in semiarid systems without facilitative mechanisms, Geophysical Research Letters, 40, 6143-6147.

• –

R.Martínez-García, J.M. Calabrese, E. Hernández-García and C. López, (2014), Minimal mechanisms for vegetation patterns in semiarid regions, Reviewed and resubmitted to Philosophical Transactions of the Royal Society A.

• •

Chapter 4:

• –

R. Martínez-García, J.M. Calabrese, T. Muller, K.A. Olson, and C. López, (2013), Optimizing the Search for Resources by Sharing Information: Mongolian Gazelles as a Case Study, Physical Review Letters, 110, 248106.

• –

R. Martínez-García, J.M. Calabrese, and C. López, (2014), Optimal search in interacting populations: Gaussian jumps versus Lévy flights, Physical Review E, 89, 032718,

• •

Chapter 5:

• –

R. Martínez-García, J.M. Calabrese, T. Muller, K.A. Olson, and C. López, (2013), Optimizing the Search for Resources by Sharing Information: Mongolian Gazelles as a Case Study, Physical Review Letters, 110, 248106.

• •

Chapter 6:

• –

R. Martínez-García, F. Vázquez, C. López, and M.A. Muñoz, (2012) Temporal disorder in up-down symmetric systems, Physical Review E, 85, 051125.

CHAPTER   1    Methods and tools

1.1 From Individual Based to Population Level Models

1.1.1 The Master equation

The master equation provides a complete description of a stochastic dynamics. It encapsulates, in the evolution of the probability of finding the system in a particular state, all the processes that occur with given transition rates. Let us consider an arbitrary system with $N$ possible states jumping from to other with exponentially distributed waiting times. In addition, let us consider that the state of the system at a given time only depends on the previous state, which is called the Markovian assumption. The probability of finding it in a particular one, $c$, at a time $t+\Delta t$ is

[TABLE]

where $c^{\prime}$ in the first term denotes the set of states that can be reached from $c$ while in the second one it refers to the states from which $c$ can be reached. The first term in Eq. (1.1) is the probability of having the system in the state $c$ at time $t$ and still remaining there at time $t+\Delta t$ (no transitions occur in the time interval $\Delta t$). The second one gives the probability of finding the system at any state $c^{\prime}$ at time $t$ and then jumping to $c$ in a time interval $\Delta t$.

In the limit of infinitely short time steps, $\Delta t\rightarrow dt$, Eq. (1.1) becomes an evolution equation for the probability of finding the system at each state $c$. This is the master equation:

[TABLE]

Gain and loss terms in Eq. (1.2) balance each other, so the probability distribution remains normalized. In addition, the coefficients $\omega_{c\rightarrow c^{\prime}}$ are rates rather than probabilities, so they have units of omverse of time and may be greater than one.

Master equations are often hard to solve because they involve a set of several, many times infinite, coupled first order ordinary differential equations. The most common techniques to obtain analytical solutions are based on the use of integral transformations such as the generating function, the Fourier or the Laplace transform [Redner, , 2001]. However, only in few simple cases the general time dependent solution $P_{c}(t)$ can be found, and most of the times numerical simulations of the underlying stochastic dynamics are done [Gillespie, , 1977].

To illustrate all the derivarions shown in this chapter, we will use a Lotka-Volterra model as a paradigmatic case of a stochastic dynamics that can be modelled at different levels. As it is shown in Fig. 1.1, several events can take place with given rates: a birth of a prey with rate $k_{b}$, a death of a predator with rate $k_{d}$ and a predation and birth of a predator with rate $k_{p}$.

The elementary processes occuring in the time interval $(t,t+dt)$ that contribute to $P(n,p;t+dt)$ are the following:

1. The population was $(n,p)$ at time $t$ and nothing happened.

2. The population was $(n-1,p)$ at time $t$ and a rabbit reproduced.

3. The population was $(n,p+1)$ at time $t$ and a fox died.

4. The population was $(n+1,p-1)$ at time $t$ and a fox ate a rabbit and reproduced.

These contributions imply a probability of having $n$ prey and $p$ predators at time $t+dt$ given by

[TABLE]

that in the limit $dt\rightarrow 0$, and retaining linear terms in $dt$, gives

[TABLE]

The master equation contains all the information about the stochastic dynamics, so it is possible to know the probability of finding the system in a particular state as a function of time. However, due to the difficulties that one usually finds to obtain its complete solution, many numerical techniques and analytical approximations have been developed to deal with it. This is the case of the Gillespie algorithm and the mean-field approximation, that will be explained next.

The Gillespie algorithm

The Gillespie algorithm [Gillespie, , 1977] is a Monte Carlo method used to simulate Poissonian 111Exponentially distributed waiting times between events stochastic processes where transitions from one state to another take place with different rates. The main objective of the algorithm is to calculate the time until the next transition takes place and the state where the system will move to. In principle, one should obtain the time at which every transition occurs, then select the one that happens first, and execute it. The advantage of Gillespie method is that it avoids simulating all the transitions and, instead, only the one that takes place first has to be reproduced.

The algorithm can be explained in four steps:

Considering that the system is initially in one of the possible $M$ states,$i$, we obtain the total escape rate from it

[TABLE]

where $j$ is the set of accesible states from $i$ and $\omega_{i\rightarrow j}$ are the individual transition rates from $i$ to each of the states labelled by $j$. 2. 2.

The time until the next jump, $dt$, is computed. It is drawn from an exponential distribution of mean $1/\Omega_{i}$. To this aim one generates a random number uniformly distributed, $u_{0}$, and computes $dt$ as

[TABLE] 3. 3.

The final state has to be determined. Each of the possible transitions takes place with a probability $p_{i\rightarrow j}$ that is proportional to the corresponding rate $\omega_{i\rightarrow j}$,

[TABLE] 4. 4.

The time is updated $t\rightarrow t+dt$

When simulated, a Gillespie realization represents a random walk trajectory for the stochastic variables that exactly represents the distribution of the master equation. It can be used, for instance, to reproduce the dynamics of the individual based Lotka-Volterra model of the Fig. 1.1, where birth, death or predation can be interpreted as a transition from a state with $n$ prey and $p$ predators to a new one with different population sizes depending on which event has occured. In Fig. 1.2 a simulation of the stochastic Lotka-Volterra dynamics using the Gillespie algorithm is shown.

Mean-field approximation

It is the simplest analytical approximation to deal with a master equation. It allows the derivation of deterministic differential equations for the mean values of the stochastic variables and establishes the simplest class of population-level models. Referring to the Lotka-Volterra model as a guiding example, we will derive the equations for the evolution of the mean number of preys, $n$, and predators, $p$. Given a multivariate probability density function with discrete variables, as it is $P(n,p;t)$, the expected values are defined as

[TABLE]

Multiplying the master equartion, Eq. (1.1.1), by $n$ and $p$ respectively and making the summation over both variables, one gets the equations for the temporal evolution of the mean values coupled to the higher moments $\langle n(t)p(t)>$

[TABLE]

It is possible to obtain the equation for the temporal evolution of $\langle n(t)p(t)\rangle$, but it would be again coupled to higher moments, leading to an infinite system of coupled differential equations. The main assumption of the mean-field approximation is to consider that both populations are independent, $\langle n(t)p(t)\rangle=\langle n(t)\rangle\langle p(t)\rangle$, so it is possible to write a closed system of deterministic differential equations for the mean value of preys and predators

[TABLE]

where $N(t)\equiv\langle n(t)\rangle$ and $P(t)\equiv\langle p(t)\rangle$.

For simplicity, the set of equations (1.1.1) can be nondimensionalised by writing [Murray, , 2002]

[TABLE]

and it becomes,

[TABLE]

The nondimensional system (1.1.1) can be solved analytically, although this is not the general case for nonlinear systems. Most of the times one has to use linear approximations and other techniques developed in the study of dynamical systems. Additionally, it is always possible to numerically integrate the equations. This has been done for equations (1.1.1) and the results are shown in Fig. 1.3.

The mean-field equations are a simplified version of the complete stochastic dynamics, but still contain most of the relevant information of the system. For instance, the oscillations in the populations are preserved for the Lotka-Volterra model. However, there are many other approximations that, although more complicated, are able to keep the inherent stochasticity of the system. The Fokker-Planck and the Langevin equations are two of them.

1.1.2 The Fokker-Planck equation

The master equation describes the dynamics of a physical system as a sequence of jumps from one state to another. We present in this section an approximation that considers the limit where these jumps are very short and the evolution of the system can be seen as a diffusive process. This leads to a simpler description in terms of the Fokker-Planck equation. The accuracy of this method is better the smaller are the jumps, so the master equation becomes a Fokker-Planck in the limit of infinitely small jumps.

There are many ways of deriving the Fokker-Planck equation. In this section, we focus on the Kramers-Moyal expansion [van Kampen, , 2007; Gardiner, , 1985]. This is not a completely rigourous derivation from a mathematical point of view, in fact many alternatives have been used in this thesis, but it is still one of the most common and intuitive approaches.

To begin with, we consider a system with several possible states. To ensure the accuracy of the expansion, we assume that all the jumps between two states are small enough, so that the set of possible states of the system can be considered as a continuous in the master equation,

[TABLE]

Next, we write, in Eq. (1.13), the transition rates as a function of the size of the jump, $r$, and of the starting point, $c$,

[TABLE]

with $r=|c^{\prime}-c|$. Then, the master equation, Eq. (1.13), becomes

[TABLE]

At this point, two assumptions have to be introduced to allow the expansion of the transition rates:

Only small jumps occur. That is, $\omega(c^{\prime};r)$ is a sharply peaked function of $r$ but varies smoothly with $c^{\prime}$. Mathematically, it means that

[TABLE] 2. 2.

The solution, $P_{c}(t)$, varies slowly with $c$ as it is expressed by Eq. (1.17).

Therefore, one can do a Taylor expansion up to second order in Eq. (1.13) to deal with the shift from $c$ to $c-r$:

[TABLE]

The first and fourth term in the right-hand side of Eq. (1.18) cancel each other, and defining the jump moments

[TABLE]

the final result can be written as

[TABLE]

This is the Fokker-Planck equation. It is important to remark that we have not shown a completely rigurous derivation. The election of the small parameter to perform the Taylor expansion has not been justified and there are many processes in which this expansion fails. This is the case of systems with jump size $\pm 1$ or some small integer, whereas typical sizes of the variable may be large, e.g., the number of molecules in a chemical reaction or the position of a random walker on a long lattice. In those cases expansions where the small parameter is explicitly taken are much more appropiate (See Chapter 3 for a rigorous derivation of the Fokker-Planck equation). Nevertheless, this description provides a good first contact with the Fokker-Planck equation, that allows the development of a large variety of population level spatial models.

On the other hand, many ecological systems, such as groups of animals and vegetation landscapes that will be studied in this thesis, are formed by many particles. Let us now suppose that we have a suspension of a very large number of identical individuals, and denote its local density by $\rho({\bf x},t)$. If the suspension is sufficiently diluted, to the extent that particles can be considered independent, then $\rho({\bf x},t)$ will obey the same Eq. (1.20) [Peliti, , 2011]. This family of models based on the density of individuals is the basis of the studies on vegetation patterns shown in the Part II of this thesis.

In either case, and independently of the way used to write it, the Fokker-Planck equation describes a large class of stochastic dynamics in which the system has a continuous sample path. The state of the system can be written as a stochastic and continuous function of time. From this picture, it seems obvious to seek a description in some direct probabilistic way and in terms of stochastic differential equations for the path of the system. This procedure is discused next.

1.1.3 The Langevin equation

In some cases it is useful to describe a system in terms of a differential equation, that gives the stochastic evolution of its state as a trajectory in the phase space. This is the Langevin equation, that has the general form

[TABLE]

where $c$ is a stochastic variable that gives the state of the system at every time. $f(c,t)$ and $g(c,t)$ are known functions and $\eta(t)$ is a rapidly fluctuating term whose average over single realizations is equal to zero, $\langle\eta(t)\rangle=0$. Any nonzero mean can be absorbed into the definition of $f(c,t)$. An idealization of a term like $\eta(t)$ must be that in which if $t\neq t^{\prime}$, $\eta(t)$ and $\eta(t^{\prime})$ are statistically independent (white noise), so

[TABLE]

where $\Gamma$ gives the strength of the random function.

To be rigorous, the differential equation (1.21) is not properly defined, although the corresponding integral equation,

[TABLE]

can be consistently defined understanding the integral of the white noise as a Wiener process $W(t)$ [van Kampen, , 2007; Gardiner, , 1985]:

[TABLE]

Hence

[TABLE]

where the second integral can be seen like a kind of Riemann integral with respect to a sample function $W(t)$.

The definition of the Langevin equation (1.21), requires a careful interpretation due to this lack of mathematical rigor. When the noise term appears multiplicatively, that is, $g(c,t)$ is not a constant, ambiguities appear in some mathematical expressions. Giving a sense to the undefined expressions constitutes one of the main goals when integrating a Langevin equation. The most widely used interpretations are those of Itô and Stratonovich (Appendix A). The Itô integral is preferred by mathematicians [van Kampen, , 2007], but it is not always the most natural choice from a physical point of view. The Stratonovich integral is more suitable, for instance, when $\eta(t)$ is a real noise with finite correlation time where the vanishing correlation time limit wants to be taken. (In the Appendix A we show a more detailed discussion). The matter is not what is the right definition of the stochastic integral, but how stochastic processes can model real systems. That is, in what situations either Itô or Stratonovich choice is the most suitable.

Langevin equations are also valid to go beyond a mean-field description. In these cases a new term enters in the equation to include diffusion, besides other spatial couplings and degrees of freedom. The variable $c(t)$ becomes a continuous field $\phi({\bf r},t)$ that depends on space and time. The Langevin equation becomes a stochastic partial differential equation of the type

[TABLE]

This approach is quite useful for spatially extended systems or to study the formation of patterns.

From the Fokker-Planck to Langevin equation and vice versa.

To close this overview on the modeling of stochastic systems, we will show the relationship between Fokker-Planck and Langevin equations. Starting from a Fokker-Planck equation for the probability distribution of the variable $c$

[TABLE]

it is easy to write down a Langevin equation of the type (1.21) [Gardiner, , 1985; van Kampen, , 2007]

[TABLE]

where $\eta(t)$ is a white, Gaussian and zero mean noise.

The coefficients of the equations are related according to

[TABLE]

provided that the Itô interpretation is chosen.

The first term in Eq. (1.27) is called drift, because it leads to the deterministic part of the Langevin equation, and the second one, the diffusion term, since it determines the stochastic part of the Langevin equation.

In the Stratonovich scheme an additional drift appears,

[TABLE]

On the other hand, if the starting point is a Langevin equation

[TABLE]

to obtain the Fokker-Planck equation one has to specify if the Itô or the Stratonovich calculus will be used. In the Stratonovich interpretation the Fokker-Planck is

[TABLE]

while in the Itô case it is

[TABLE]

The diffusion term vanishes typically with the number of components as $N^{-1/2}$, so it is negligible if the system is large enough. Therefore, in the thermodynamic limit where $N$ and the volume $V$ tend to infinity keeping $N/V$ finite, a deterministic mean-field approximation gives an accurate description. Sometimes, this way is walked on the inverse sense. One may start with a deterministic equation and, using heuristic arguments, add noise to obtain the Langevin equation. Then, following the steps that have been explained in this section it is possible to get a Fokker-Planck equation.

1.2 Linear stability analysis

Linear stability analysis is the simplest analytical tool used to study the formation of patterns in deterministic spatially extended systems. It assumes an ideal infinite system and uses Fourier analysis to investigate the stability of its homogeneous state. We will consider in this section the two dimensional case. The starting point is the equation for the evolution of a field $\phi$

[TABLE]

where $R$ is the control parameter. The linear stability analysis assumes that the system is at the homogeneous (spatially independent) stationary state $\phi({\bf x},t)=\phi_{0}$ and studies its stability against small perturbations that will be denoted by $\psi({\bf x},t)$, with $|\psi|\ll 1$. The technique is applied in the Appendix A to one particular case and the calculations explained in detail. In this section we will introduce and discuss the theoretical basis and the main results that can be obtained. Plugging the ansatz $\phi({\bf x},t)=\phi_{0}+\psi({\bf x},t)$ into the model Eq. (1.35) and retaining only linear terms in the perturbation, one obtains a linear equation for the evolution of the perturbation at short times that can be solved using the Fourier transform. Then, the final task is to solve the transformed equation for the perturbation, $\hat{\psi}({\bf k},t)$. Assuming that at short time scales the temporal dependence is $\hat{\psi}({\bf k},t)\propto\exp(\lambda({\bf k})t)$, where $\lambda$ is the growth rate, then $\hat{\psi}({\bf k},t)=\lambda({\bf k})\hat{\psi}({\bf k},t)$. Finally an expression for $\lambda({\bf k})$ can be obtained. It is called the dispersion relation and contains all the information about the evolution of the Fourier modes of $\hat{\psi}({\bf k},t)$. The modes ${\bf k}$ with a negative growth rate will be stable while those corresponding to $\lambda\geq 0$ are unstable and lead to perturbations growing in time and, therefore, to spatial patterns in the system. The dispersion relation also allows to obtain the characteristic wavelength of the pattern through the value of the most unstable Fourier mode, ${\bf k}_{c}$, that most of the times corresponds with the one with the highest growth rate.

Depending on the functional form of the dispersion relation, it is possible to establish a classification of the different types of linear instabilities appearing in natural systems [Hohenberg and Halperin, , 1977]. These classes are shown in Fig. 1.4, where the real part of $\lambda$ is sketched as a function of the wave number, ${\bf k}$. They are:

• •

Type I. For $R<R_{c}$ the homogeneous state is stable and ${\rm Re}\lambda<0$, whereas for $R=R_{c}$ the instability sets in (${\rm Re}\lambda=0$) at a wave vector $k_{c}$. For $R>R_{c}$ there is a band of wave vectors $k_{-}<k<k_{+}$ for which the uniform state is unstable. The patterns observed in these system will be dominated by a wavelength given by one of this unstable modes, typically by that with the highest growth rate, $k_{c}$. This case is represented in the left panel of Fig. 1.4.

• •

Type II. This is a different type of instability appearing when, for some reason (usually a conservation law), ${\rm Re}\lambda(k=0)=0$ independently of the value of the control parameter $R$. This corresponds with the central panel of Fig. 1.4. The critical wave vector, the one that becomes unstable by the first time, is now $k_{c}=0$, and a band of unstable modes appears between [math] and $k_{+}$ for $R>R_{c}$. The pattern occurs on a long length scale. This case is remarkable because the critical wave vector is different from that with the highest growth rate.

• •

Type III. In this case both the instability and the maximum growth rate occur at $k_{c}=0$. There is not an intrinsic length scale, and patterns will occur over a length scale defined by the system size or the dynamics. This situation is depicted in the right panel of Fig. 1.4.

Finally, there are two subtypes for each type of instability depending on the temporal instability: stationary if ${\rm Im}\lambda=0$, and oscillatory if ${\rm Im}\lambda\neq 0$.

Linear stability analysis provides analytical results about the formation of patterns in spatially extended systems, such as the dominant wavelength and the type of instability leading the structure. However, it is important to remark that the analysis assumes that the perturbations of the uniform state are small. This assumption is good at short times and for an initial condition that has a small magnitude, but at long times the nonlinear terms left out in the linear approximation become important [Cross and Greenside, , 2009]. One effect of nonlinearity is to quench the assumed exponential growth. Further analysis, such as weakly nonlinear stability analysis [Cross and Hohenberg, , 1993], must be used in these cases.

1.3 First-passage times processes

First-passage phenomena are of high relevance in stochastic processes that are triggered by a first-passage event [Redner, , 2001] and play a fundamental role quantifying and limiting the success of different processes that can be mapped into random walks. Ecology and biology offer some examples such as the lifetime of a population or the duration of a search or a biochemical reaction.

In this section we will present some results on first-passage times in the simple case of a discrete symmetric random walk moving in a finite interval $[x_{-},x_{+}]$ [Redner, , 2001]. The extension to higher dimensions is straightforward. Let us denote the mean time to exit the interval starting at $x$ by $T(x)$. This quantity is equal to the exit time of a given trajectory times the probability of that path, averaged over all the trajectories,

[TABLE]

where $t_{p}$ is the exit time of the trajectory $p$ that starts at $x$ and $\mathcal{P}_{p}$ the probability of the path. Because of the definition of a symmetric random walk on a discrete space with jump length $\delta x$ and time step $\delta t$, the mean exit time also obeys

[TABLE]

with boundary conditions $T(x_{-})=T(x_{+})=0$ which correspond to a mean exit time equal to zero if the particle starts at either border of the interval. $\delta x$ is the jumping length. This recursion relation expresses the mean exit time starting at $x$ in terms of the outcome one step in the future, for which the initial walk can be seen as restarting in $x\pm\delta x$ (each with probability $1/2$) but also with the time incremented by $\delta t$.

Doing a Taylor expansion to the lowest nonvanishing order in Eq. (1.37), and considering the limit of continuous time and space, it yields

[TABLE]

where $D=\delta x^{2}/2\delta t$ is the difussion constant. In the case of a two dimensional domain Eq. (1.38) is

[TABLE]

These results can be extended to the case of general jumping processes with a single-step jumping probability given by $p_{\mathbf{x}\rightarrow\mathbf{x}^{\prime}}$. The equivalent of Eq. (1.37) is

[TABLE]

that provides an analog of Eq. (1.39) that is

[TABLE]

where ${\bf v}({\bf x})$ is a local velocity that gives the mean displacement after a single step when starting from ${\bf x}$ in the hopping process. This equation can be solved in each particular case. We have used it in this thesis as an starting point of many of the calculations in the Part IV. See Appendix B for a detailed calculation.

Part II  **Vegetation Patterns

**

CHAPTER   2    Mesic savannas

In this chapter we propose a continuum description for the dynamics of tree density in mesic savannas Martínez-García et al., 2013a inspired on the individual based model introduced in Calabrese et al., [2010]. It considers only long-range competition among trees and the effect of fires resulting in a local facilitation mechanism. Despite short-range facilitation is taken to the local-range limit, the standard full spectrum of spatial structures obtained in general vegetation models is recovered. Long-range competition is thus the key ingredient for the development of patterns. This result opens new questions on the role that facilitative interactions play in the maintenance of vegetation patterns. The long time coexistence between trees and grass, the effect of fires on the survival of trees as well as the maintenance of the patterns are also studied. The influence of demographic noise is analyzed. The stochastic system, under parameter constraints typical of more humid landscapes, shows irregular patterns characteristic of realistic situations. The coexistence of trees and grass still remains at reasonable noise intensities.

2.1 Introduction

Savanna ecosystems are characterized by the long-term coexistence between a continuous grass layer and scattered or clustered trees [Sarmiento, , 1984]. Occurring in many regions of the world, in areas with very different climatic and ecological conditions, the spatial structure, persistence, and resilience of savannas have long intrigued ecologists [Scholes and Archer, , 1997; Sankaran et al., , 2005; Borgogno et al., , 2009; Belsky, , 1994]. However, despite substantial research, the origin and nature of savannas have not yet been fully resolved and much remains to be learned.

Savanna tree populations often exhibit pronounced, non-random spatial structures [Skarpe, , 1991; Barot et al., , 1999; Jeltsch et al., , 1999; Caylor et al., , 2003; Scanlon et al., , 2007]. Much research has therefore focused on explaining how spatial patterning in savannas arises [Jeltsch et al., , 1996, 1999; Scanlon et al., , 2007; Skarpe, , 1991; Calabrese et al., , 2010; Vazquez et al., , 2010]. In most natural plant systems both facilitative and competitive processes are simultaneously present [Scholes and Archer, , 1997; Vetaas, , 1992] and hard to disentangle [Veblen, , 2008; Barbier et al., , 2008]. Some studies have pointed toward the existence of short-distance facilitation [Caylor et al., , 2003; Scanlon et al., , 2007], while others have demonstrated evidence of competition [Skarpe, , 1991; Jeltsch et al., , 1999; Barot et al., , 1999], with conflicting reports sometimes arriving from the same regions.

Different classes of savannas, which can be characterized by how much rainfall they typically receive, should be affected by different sets of processes. For example, in semiarid savannas water is extremely limited (low mean annual precipitation) and competition among trees is expected to be strong, but fire plays little role because there is typically not enough grass biomass to serve as fuel. In contrast, humid savannas should be characterized by weaker competition among trees, but also by frequent and intense fires. In-between these extremes, in mesic savannas, trees likely have to contend with intermediate levels of both competition for water and fire [Calabrese et al., , 2010; Sankaran et al., , 2005, 2008; Bond et al., , 2003; Bond, , 2008; Bucini and Hanan, , 2007].

Competition among trees is mediated by roots that typically extend well beyond the crown [Borgogno et al., , 2009; Barbier et al., , 2008]. Additionally, fire can lead to local facilitation due to a protection effect, whereby vulnerable juvenile trees placed near adults are protected from fire by them [Holdo, , 2005]. We are particularly interested in how the interplay between these mechanisms governs the spatial arrangement of trees in mesic savannas, where both mechanisms may operate. On the other side, it has frequently been claimed that pattern formation in arid systems can be explained by a combination of long-distance competition and short-distance facilitation [Klausmeier, , 1999; Lefever and Lejeune, , 1997; Lefever et al., , 2009; Lefever and Turner, , 2012; Rietkerk et al., , 2002; von Hardenberg et al., , 2001; D’Odorico et al., 2006b, ]. This combination of mechanisms is also known to produce spatial structures in many other natural systems [Cross and Hohenberg, , 1993]. Although mesic savannas do not display the same range of highly regular spatial patterns that arise in arid systems (e.g., tigerbush), similar mechanisms might be at work. Specifically, the interaction between long-range competition and short-range facilitation might still play a role in pattern formation in savanna tree populations, but only for a limited range of parameter values and possibly modified by demographic stochasticity.

Although the facilitation component has often been thought to be a key component in previous vegetation models [D’Odorico et al., 2006b, ; D’Odorico et al., 2006c, ; Rietkerk et al., , 2002; Scanlon et al., , 2007], Rietkerk and Van de Koppel [Rietkerk and van de Koppel, , 2008], speculated, but did not show, that pattern formation could occur without short-range facilitation in the particular example of tidal freshwater marsh. In the case of savannas, as stated before, the presence of adult trees favor the establishment of new trees in the area, protecting the juveniles against fires. Considering this effect, we take the facilitation component to its infinitesimally short spatial limit, and study its effect in the emergence of spatially periodic structures of trees. To our knowledge, this explanation, and the interrelation between long-range competition and local facilitation, has not been explored for a vegetation system.

To this aim, we develop a minimalistic model of savannas that considers two of the factors, as already mentioned, thought to be crucial to structure mesic savannas: tree-tree competition and fire, with a primary focus on spatially nonlocal competition. Employing standard tools used in the study of pattern formation phenomena in physics (stability analysis and the structure function) [Cross and Hohenberg, , 1993], we explore the conditions under which the model can produce non-homogeneous spatial distributions. A key strength of our approach is that we are able to provide a complete and rigorous analysis of the patterns the model is capable of producing, and we identify which among these correspond to situations that are relevant for mesic savannas. We further examine the role of demographic stochasticity in modifying both spatial patterns and the conditions under which trees persist in the system in the presence of fire, and discuss the implications of these results for the debate on whether the balance of processes affecting savanna trees is positive, negative, or is variable among systems. This is the framework of our study: the role of long-range competition, local facilitation and demographic fluctuations in the spatial structures of mesic savannas.

2.2 The deterministic description

In this section we derive the deterministic equation for the local density of trees, such that dynamics is of the logistic type and we only consider tree-tree competition and fire. We study the formation of patterns via stability analysis and provide numerical simulations, showing the emergence of spatial structures.

2.2.1 The nonlocal savanna model

Calabrese et al., [2010] introduced a simple discrete-particle lattice savanna model that considers the birth-death dynamics of trees, and where tree-tree competition and fire are the principal ingredients. These mechanisms act on the probability of establishment of a tree once a seed lands at a particular point on the lattice. In the discrete model, seeds land in the neighborhood of a parent tree with a rate $b$, and establish as adult trees if they are able to survive both competition neighboring trees and fire. As these two phenomena are independent, the probability of establishment is $P_{E}=P_{C}P_{F}$, where $P_{C}$ is the probability of surviving the competition, and $P_{F}$ is the probability of surviving a fire event. From this dynamics, we write a deterministic differential equation describing the time evolution of the global density of trees (mean field), $\rho(t)$, where the population has logistic growth at rate $b$, and an exponential death term at rate $\alpha$. It reads:

[TABLE]

Generalizing Eq. (2.1), we propose an evolution equation for the space-dependent (local) density of trees, $\rho({\bf x},t)$:

[TABLE]

We allow the probability of overcoming competition to depend on tree crowding in a local neighborhood, decaying exponentially with the density of surrounding trees as

[TABLE]

where $\delta$ is a parameter that modulates the strength of the competition, and $G({\bf x})$ is a positive kernel function that introduces a finite range of influence. This model is related to earlier one of pattern formation in arid systems [Lefever and Lejeune, , 1997], and subsequent works [Lefever et al., , 2009; Lefever and Turner, , 2012], but it differs from standard kernel-based models in that the kernel function accounts for the interaction neighborhood, and not for the type of interaction with the distance. Note also that the nonlocal term enters nonlinearly in the equation.

Following Calabrese et al., [2010], $P_{F}$ is assumed to be a saturating function of grass biomass, $1-\rho({\bf x},t)$, similar to the implementation of fire of Jeltsch et al. in [Jeltsch et al., , 1996]

[TABLE]

where $\sigma$ governs the resistance to fire, so $\sigma=0$ means no resistance to fires. Notice how our model is close to the one in [Calabrese et al., , 2010] through the definitions of $P_{C}$ and $P_{F}$, although we consider the probability of surviving a fire depending on the local density of trees, and in [Calabrese et al., , 2010] it depends on the global density. The final deterministic differential equation that considers tree-tree competition and fire for the spatial tree density is

[TABLE]

where

[TABLE]

Thus, we have a logistic-type equation with an effective growth rate that depends nonlocally on the density itself, and which is a combination of long-range competition and local facilitation mechanisms (fire). The probability of surviving a fire is higher when the local density of trees is higher, as can be seen from the definition in Eq. (2.4).

In Fig. 2.1 we show numerical solutions for the mean field Eq. (2.1) (lines) and the spatially explicit model (equation 2.5) (dots) in the stationary state $(t\rightarrow\infty)$ using different values of the competition. We have used a top-hat function as the competition kernel, $G({\bf x})$ (See Sec. 2.2.2 for more details on the kernel choice). We observe a very good agreement of both descriptions which becomes worse when we get closer to the critical point $\sigma^{*}$, where the model presents a phase transition from a tree-grass coexistence to a grassland state. This disagreement appears because while the mean field equation describes an infinite system, the Eq. (2.5) description forces us to choose a size for the system.

The model reproduces the long-term coexistence between grass and trees that is characteristic of savannas. To explore this coexistence, we study the long-time behavior of the system and analyze the homogeneous stationary solutions of Eq. (2.5), which has two fixed points. The first one is the absorbing state representing the absence of trees, $\rho_{0}=0$, and the other can be obtained, in the general case, by numerically solving

[TABLE]

In the regime where $\rho_{0}$ is small (near the critical point), if competition intensity, $\delta$, is also small, it is possible to obtain an analytical expression for the critical value of the probability of surviving a fire, $\sigma^{*}$,

[TABLE]

Outside of the limit where $\delta\ll 1$, we can solve Eq. (2.7) numerically in $\rho_{0}$ to show that the critical value of the fire resistance parameter, $\sigma^{*}$, does not depend on competition. A steady state with trees is stable for higher fire survival probability (Fig. 2.1).

There is, then, a transition from a state where grass is the only form of vegetation to another state where trees and grass coexist at $\sigma^{*}$. In what follows, we fix $\alpha=1$, so we choose our temporal scale in such a way that time is measured in units of $\alpha$. This choice does not qualitatively affect our results.

2.2.2 Linear stability analysis

The spatial patterns can be studied by performing a linear stability analysis [Cross and Hohenberg, , 1993] of the stationary homogeneous solutions of Eq. (2.5), $\rho_{0}=\rho_{0}(\sigma,\delta)$. The stability analysis is performed by considering small harmonic perturbations around $\rho_{0}$, $\rho({\bf x},t)=\rho_{0}+\epsilon{\rm e}^{\lambda t-i{\bf k}\cdot{\bf x}}$, $\epsilon\ll 1$. After some calculations111A linear stability analysis in a similar equation modeling vegetation in arid systems is shown in detail in Appendix A., one arrives at a perturbation growth rate given by

[TABLE]

where ${\hat{G}}(k)$, $k=|{\bf k}|$, is the Fourier transform of the kernel,

[TABLE]

The critical values of the parameters of the transition to pattern, $\delta_{c}$ and $\sigma_{c}$, and the fastest growing wavenumber $k_{c}$, are obtained from the simultaneous solution of

[TABLE]

Note that $k_{c}$ represents the most unstable mode of the system, which means that it grows faster than the others and eventually dominates the state of the system. Therefore, it determines the length scale of the spatial pattern. These two equations yield the values of the parameters $\delta$ and $\sigma$ at which the maximum of the curve $\lambda(k)$, right at $k_{c}$, starts becoming positive. This signals the formation of patterns in the solutions of Eq. (2.5). As Eq. (2.12) is explicitly written as

[TABLE]

the most unstable wavenumber $k_{c}$ can be obtained by evaluating the zeros of the derivative of the Fourier transform of the kernel.

Eq. (2.9) shows that competition, through the kernel function, fully determines the formation of patterns in the system. The local facilitation appears in $b_{eff}(\rho_{0})$ and it is not relevant in the formation of spatial structures. If the Fourier transform of $G$ never takes positive values, then $\lambda(k;\sigma,\delta)$ is always negative and only the homogeneous solution is stable. However, when $\hat{G}$ can take negative solutions then patterns may appear in the system. What does this mean in biological terms? Imagine that we have a family of kernels described by a parameter $p$: $G({\bf x})=\exp(-|({\bf x})/R|^{p})$ ($R$ gives the range of competition). The kernels are more peaked around ${\bf x}=0$ for $p<2$ and more box-like when $p>2$. It turns out that this family of functions has non-negative Fourier transform for $0\leq p<2$, so that no patterns appear in this case. A lengthy discussion of this property in the context of competition of species can be found in Pigolotti et al., [2007]. Thus, the shape of the competition kernel dictates whether or not patterns will appear in the system. If pattern formation is possible, then the values of the fire and competition parameters govern the type of solution (see Sec. 2.2.3).

Our central result for nonlocal competition is that, contrary to conventional wisdom, it can, in the limit of infinitesimally short (purely local) facilitation, promote the clustering of trees. Whether or not this occurs depends entirely on the shape of the competition kernel. For large $p$ we have a box-like shape, and in these cases trees compete strongly with other trees, roughly within a distance $R$ from their position. The mechanism behind this counterintuitive result is that trees farther than $R$ away from a resident tree area are not able to invade the zone defined by the radius R around the established tree (their seeds do not establish there), so that an exclusion zone develops around it. For smaller $p$ there is less competition and the exclusion zones disappear. We will develop longer this concept in the next chapter.

For a more detailed analysis, one must choose an explicit form for the kernel function. Our choice is determined by the original $P_{C}$ taken in [Calabrese et al., , 2010], so that it decays exponentially with the number of trees in a neighborhood of radius $R$ around a given tree. Thus, for $G$ we take the step function (limit $p\rightarrow\infty$)

[TABLE]

As noticed before, the idea behind the nonlocal competition is to capture the effect of the long roots of a tree. The kernel function defines the area of influence of the roots, and it can be modeled at first order with the constant function of Eq. (2.17). Thus the parameter R, which fixes the nonlocal interaction scale, must be of the order of the length of the roots [Borgogno et al., , 2009]. Since the roots are the responsible for the adsorption of resources (water and soil nutrients), a strong long-range competition term implies strong resource depletion. For this kernel the Fourier transform is [López and Hernández-García, , 2004] ${\hat{G}}(k)=2\pi R^{2}J_{1}(kR)/kR$ and its derivative is ${\hat{G}}^{\prime}(k)=-2\pi R^{2}J_{2}(kR)/k$, where $k\equiv|{\bf k}|$, and $J_{i}$ is the $i^{th}$-order Bessel function. Since ${\hat{G}}(k)$ can take positive and negative values, pattern solutions may arise in the system, that will in turn depend on the values of $\delta$ and $\sigma$. The most unstable mode is numerically obtained as the first zero of $\lambda^{\prime}(k)$, Eq. (2.13), which means the first zero of the Bessel function $J_{2}(kR)$. This value only depends on $R$, being independent of the resistance to fires and competition, and it is $k_{c}=5.136/R$. Because a pattern of $n$ cells is characterized by a wavenumber $k_{c}=2\pi n/L$, where $L$ is the system size, the typical distance between clusters, $d_{t}=L/n$, using the definition of the critical wavenumber is given by $d_{t}\approx 1.22R$. In other words, it is approximately the range of interaction $R$. This result is also independent of the other parameters of the system.

Since we are interested in the effect of competition and fire on the distribution of savanna trees, we will try to fix all the parameters but $\sigma$ and $\delta$. We will explore the effect of different values of these parameters on the results. First, we have chosen, as in Calabrese et al., [2010], the death rate $\alpha=1$, and solving Eq. (2.7) we will roughly estimate the birth rate, $b$. We will work in the limit of intermediate to high mean annual precipitation, so water is non-limiting and thus we can neglect the effects of competition ($\delta=0$). At this intermediate to high mean annual precipitation the empirically observed upper limit of savanna tree cover is approximately $\rho_{0}=0.8$ [Sankaran et al., , 2005; Bucini and Hanan, , 2007]. To reach this upper limit in the tree cover, disturbances must also be absent, implying no fire ($\sigma\rightarrow\infty$). In this limit, the mean field Eq. (2.1) is quantitatively accurate, as it is shown in Fig. 2.1, and the stationary mean field solution of the model depends only on the birth rate

[TABLE]

It can be solved for $b$ for a fixed $\rho_{0}=0.8$, and it yields $b=5$ [Calabrese et al., , 2010]. In the following we just consider the dependence of our results on $\delta$ and $\sigma$. In particular, $\rho_{0}=\rho_{0}(\sigma,\delta)$.

The phase diagram of the model, computed numerically, is shown in Fig. 2.2, where we plot the spatial character of the steady solution (homogeneous or inhomogeneous) as a function of $\delta$ and $\sigma$. Note that increasing competition enhances the inhomogeneous or pattern solution. This is because, as we are now in the case of a kernel giving rise to clusters, increasing $\delta$ makes it more difficult to enter the exclusion zones in-between the clusters. For very strong competition (high, unrealistic, $\delta$), fire has no influence on the pattern.

The critical line separating these two solutions (pattern and homogeneous) can be obtained analytically as a function of the parameters $\delta$, $\sigma$, $\rho_{0}$ and $\hat{G}(k_{c})$. Taking $b=5$ and $\alpha=1$, it is

[TABLE]

This complicated expression must be evaluated numerically together with the solution of Eq. (2.7) for the stationary density of trees, which is also a function of the competition and fire parameters. We show the results in Fig. 2.2, where the curve, represented with the black crosses, fits perfectly with the numerical results from the linear stability analysis.

With $b=5$, in the absence of fire ($\sigma\rightarrow\infty$), and for weak competition, we can take the limits $\delta\rightarrow 0$ and $\sigma\rightarrow\infty$ of the dispersion relation Eq. (2.9), leading to

[TABLE]

In Fig. 2.1, for large $\sigma$, it can be seen that typically $\rho_{0}>0.4$, so Eq. (2.20) becomes negative. This result means that in this limit, trees are uniformly distributed in the system as there is no competition, and space does not play a relevant role in the establishment of new trees. Such situation could be interpreted as favorable to forest leading to a fairly homogeneous density of trees. This result agrees with the phase plane plotted in Fig. 2.2. In biological terms, there are no exclusion zones in the system because there is no competition.

2.2.3 Numerical simulations

The previous analysis provides information, depending on the competition and fire parameters, about when the solution is spatially homogenous and when trees arrange in clusters. However, the different shapes of the patterns have to be studied via numerical simulations [Ridolfi et al., , 2011] of the whole equation of the model. We have taken a finite square region of savanna with an area of $1$ ha., allowed competition to occur in a circular area of radius $R=8\ m$, and employed periodic boundary conditions and a finite differences algorithm to obtain the numerical solution. Similarly to what has been observed in studies of semiarid water limited systems [D’Odorico et al., 2006b, ; Rietkerk et al., , 2002], different structures, including gaps, stripes, and tree spots, are obtained in the stationary state as we increase the strength of competition for a fixed value of the fire parameter or, on the other hand, as we decrease the resistance to fires for a given competition intensity. In both equivalent cases, we observe this spectrum of patterns as far as we go to a more dry state of the system, where resources (mainly water) are more limited (see Figs. 2.3(a)-2.3(c)) and competition is consequently stronger. This same sequence of appearance of patterns has been already observed in the presence of different short-range facilitation mechanisms [Lejeune and Tlidi, , 1999; Rietkerk et al., , 2002]. It indicates that, when $\delta$ is increased (i.e. the probability of surviving competition is decreased), new trees cannot establish in the exclusion areas so clustering is enhanced.

On the other hand, in the case of fire-prones savannas, previous works had only shown either tree spot [Lejeune et al., , 2002] or grass spots [D’Odorico et al., , 2007]. Therefore, at some values of the parameter space (see Fig. 2.3b), the patterns in our deterministic approach are not observed in mesic savannas, and should correspond to semiarid systems. However, we will show in the following sections that under the parameter constraints of a mesic savanna, and considering the stochastic nature of the tree growth dynamics in the system (i.e. demographic noise), our model shows realistic spatial structures.

A much more quantitative analysis of the periodicity in the patterns can be performed via the structure function. This will be helpful to check the previous results and, especially, for the analysis of the data of the stochastic model of the next section, for which we will not present analytical results. The structure function is defined as the modulus of the spatial Fourier transform of the density of trees in the stationary state,

[TABLE]

where the average is a spherical average over the wavevectors with modulus $k$. The structure function is helpful to study spatial periodicities in the system, similar to the power spectrum of a temporal signal. Its maximum identifies dominant periodicities, which in our case are the distances between tree clusters. Note that the geometry of the different patterns cannot be uncovered with the structure function, since it involves a spherical average. In Fig. 2.4, we show the transition to patterns using the maximum of the structure function as a function of the competition parameter. A peak appears when there are spatial structures in the system, so $Max[S(k)]\neq 0$. However, we do not have information about the values where the shapes of the patterns change. Taking $R=8~{}m$, the peak is always at $\lambda_{c}=10m$ for our deterministic savanna model, independently of the competition and fire resistance parameters, provided that they take values that ensure the emergence of patterns in the system (see the line labeled by $\Gamma=0$ in Fig. 2.5; for the definition of $\Gamma$ see Sec. 2.3). This result is in good agreement with the theoretical one provided for the wavelength by the linear stability analysis $\lambda=2\pi/k_{max}=9.78~{}m$, which is also independent of competition and resistance to fires.

2.3 Stochastic model

The perfectly periodic patterns emerging in Fig. 2.3 seem to be far from the disordered ones usually observed in aerial photographs of mesic savannas and shown by individual based models [Calabrese et al., , 2010; Jeltsch et al., , 1999; Barot et al., , 1999; Caylor et al., , 2003]. We have so far described a savanna system in terms of the density of trees with a deterministic dynamics. The interpretation of the field $\rho(\mathbf{x},t)$ is the density of tree (active) sites in a small volume, $V$. If we think of trees as reacting particles which are born and die probabilistically, then to provide a reasonable description of the underlying individual-based birth and death dynamics, we have to add a noise term to the standard deterministic equation. It will take into account the intrinsic stochasticity present at the individual level in the system.

If we take a small volume, $V$, the number of reactions taking place is proportional to the number of particles therein, $N$, with small deviations. If $N$ is large enough, the central limit theorem applies to the sum of $N$ independent random variables and predicts that the amplitude of the deviation is of the order of $\sqrt{N}\propto\sqrt{\rho(\mathbf{x},t)}$ [Gardiner, , 1985]. This stochasticity is referred to as demographic noise. The macroscopic equation is now stochastic,

[TABLE]

where $\Gamma\propto\sqrt{b_{eff}}$ (but we take it as a constant, [Dickman, , 1994]) modulates the intensity of $\eta(\mathbf{x},t)$, a Gaussian white noise term with zero mean and correlations given by Dirac delta distributions

[TABLE]

The complete description of the dynamics in Eq.(2.22) should have the potential to describe more realistic patterns.

We first investigate the effect of demographic noise on the persistence of trees in the system. We show in (Fig. 2.6) that the critical point, $\sigma^{*}$, depends on the value of the competition parameter $\delta$. This effect is rather small, so that when $\delta$ increases the transition to the grassland state appears only for a slightly larger $\sigma$ (i.e, less frequent fire). The reason seems to be that fire frequency and intensity depend on grass biomass. Seasonally wet savannas support much more grass biomass that serves as fuel for fires during the dry season [D’Odorico et al., 2006a, ; Hanan et al., , 2008]. Dry savannas have much lower grass biomass, so they do not burn as often or as intensely. The shift of the critical value of $\sigma$ when competition is stronger is consistent with the one showed in [Calabrese et al., , 2010], as can be seen comparing Fig. 2 therein with Fig. 2.6 here. Besides, the values obtained for $\sigma^{*}$ are larger when we consider the demographic stochasticity [Stanley, , 1987] neglected in the deterministic field approach.

We explore numerically the stochastic savanna model using an algorithm developed in Dickman, [1994] (See B). Note that the noise makes the transition to pattern smoother so the change from homogeneous to inhomogeneous spatial distributions is not as clear as it is in the limit where the demographic noise vanishes (See Fig. 2.4). The presence of demographic noise in the model, as shown in Fig. 2.2 (red stars), also decreases the value of the competition strength at which patterns appear in the system, as has been observed in other systems. Mathematically, these new patterns appear since demographic noise maintains Fourier modes of the solution which, due to the value of the parameters, would decay in a deterministic approach [Butler and Goldenfeld, , 2009]. Biologically, exclusion zones are promoted by demographic noise, since it does not affect regions where there are not trees. On the other hand in vegetated areas fluctuations may enhance tree density, leading to stronger competition. The presence of demographic noise in the model allows the existence of patterns under more humid conditions. This result is highly relevant for mesic savannas, as we expect competition to be of low to intermediate strength in such systems. We show two examples of these irregular patterns in Fig. 2.7(a) and Fig. 2.7(b). Unrealistic stripe-like patterns no longer appear in the stochastic model.

We have studied the dynamics of the system for some values of the fire and competition parameters. Demographic noise influences the spatial structures shown by the model. The deterministic approach shows a full spectrum of patterns which are not visually realistic for mesic savannas (but for arid systems). The role of the noise is to transform this spectrum of regular, unrealistic patterns into more irregular ones (Figures 2.7(a)-2.7(d)) that remind the observed in aerial photographs of real mesic savannas. On the other hand, these patterns are statistically equivalent to the deterministic ones, as it is shown with the structure function in Fig. 2.5. The dominant scale in the solution is given by the interaction radio, $R$, and it is independent of the amplitude of the noise (see the structure function in Fig. 2.5, peaked around $\lambda=10~{}m$ independently of the noise). Besides, over a certain treshold in the amplitude, demographic noise destroys the population of trees. Therefore, the model presents an active-absorbing transition with the noise strength, $\Gamma$, being the control parameter.

2.4 Discussion

Understanding the mechanisms that produce spatial patterns in savanna tree populations has long been an area of interest among savanna ecologists [Skarpe, , 1991; Jeltsch et al., , 1999; Barot et al., , 1999; Caylor et al., , 2003; Scanlon et al., , 2007]. A key step in such an analysis is defining the most parsimonious combination of mechanisms that will produce the pattern in question. In this chapter the combination of long-range competition for resources and the facilitation induced by fire are considered the responsible of the spatial structures, in the line of studies of vegetation pattern formation in arid systems, where also a combination of long-range inhibition and short-range facilitation is introduced [Klausmeier, , 1999; Lefever and Lejeune, , 1997; Rietkerk et al., , 2002; von Hardenberg et al., , 2001]. The main difference is that the facilitation provided by the protection effect of adult trees against fires in our savanna model takes the short-range facilitation to its infinitesimally short limit (i.e, local limit). Under this assumption we have studied the conditions under which our model could account for patterns. We have shown that nonlocal competition combined with local facilitation induces the full range of observed spatial patterns, provided the competition term enters nonlinearly in the equation for the density of trees, and that competition is strong enough.

The key technical requirement for this effect to occur is that the competition kernel must be an almost constant function in a given competition region, and decay abruptly out of the region. We verify this condition working with supergaussian kernel functions. In practice, this means that competition kernels whose Fourier transform takes negative values for some wavenumber values, will lead to competition driven clustering.

The other mechanism we have considered for a minimalistic but realistic savanna model, fire, has been shown to be relevant for the coexistence of trees and grass and for the shape of the patterns. However, competition is the main ingredient allowing pattern solutions to exist in the model. If the shape of the kernel allows these types of solutions, then the specific values of fire and competition parameters determine the kind of spatial structure that develops. It is also worth mentioning that one can observe the full spectrum of patterns in the limit where fires vanish ($\sigma\rightarrow\infty$), so there is no facilitation at all, provided competition is strong enough. However, when there is no competition, $\delta=0$, no patterns develop regardless of the value of the fire term. Therefore, we conclude that the nonlocal competition term is responsible for the emergence of clustered distributions of trees in the model, with the fire term playing a relevant role only to fix the value of the competition parameter at which patterns appear. In other words, for a given competition strength, patterns appear more readily when fire is combined with competition. A similar mechanism of competitive interactions between species has been shown to give rise to clusters of species in the context of classical ecological niche theory. Scheffer and van Nes [Scheffer and van Nes, , 2006] showed that species distribution in niche space was clustered, and Pigolotti et al. [Pigolotti et al., , 2007] showed that this arises as an instability of the nonlocal nonlinear equation describing the competition of species.

Long-distance competition for resources in combination with the local facilitation due to the protection effect of adult trees in the establishment of juvenile ones can explain the emergence of realistic structures of trees in mesic savannas. In these environmental conditions, competition is limited, so we should restrict to small to intermediate values of the parameter $\delta$, and the effect of fires is also worth to be taken into account. However, these two ingredients give a full range of patterns observed in vegetated systems, but not in the particular case of savannas. It is necessary to consider the role of demographic noise, which is present in the system through the stochastic nature of the birth and death processes of individual trees. In this complete framework our model shows irregular patterns of trees similar to the observed in real savannas.

The other important feature of savannas, the characteristic long-time coexistence of trees and grass is well captured with our model (Figures 2.1 and 2.6). Besides, the presence of demographic noise, as it is shown in Fig. 2.6, makes our approach much more realistic, since the persistence of trees in the face of fires is related to the water in the system. On the other hand, demographic stochasticity causes tree extinction at lower fire frequencies (larger $\sigma$) than in the deterministic case. This is because random fluctuations in tree density are of sufficient magnitude that this can hit zero even if the deterministic stationary tree density (for a given fire frequency) is greater than zero. This effect vanishes if we increase the system size. The demographic noise is proportional to the density of trees (proportional to $(L_{x}\times L_{y})^{-1}$), so fluctuations are smaller if we study bigger patches of savannas. As usually happens in the study of critical phenomena in Statistical Mechanics, the extinction times due to demographic noise increase exponentially with the size of the system for those intensities of competition and fire that allow the presence of trees in the stationary state. Over the critical line, this time will follow a power law scaling, and a logarithmic one when the stationary state of the deterministic model is already absorbing (without trees) [Marro and Dickman, , 2008].

2.5 Summary

We have shown the formation of patterns in a minimal savanna model, that considers the combination of long-range competition and local facilitation mechanisms as well as the transition from trees-grass coexistence to a grass only state.

The salient feature of the model is that it only considers nonlocal (and nonlinear) competition through a kernel function which defines the length of the interaction, while the facilitation is considered to have an infinitesimally short influence range. Our model thus differs from standard kernel-based savanna models that feature both short-range facilitation and long-range competition. The same sequence of spatial patterns appears in both approaches, confirming Rietkerk and van de Koppel’s [Rietkerk and van de Koppel, , 2008] suggestion that short-range facilitation does not induce spatial pattern formation by itself, and long-distance competition is also needed. It also suggests that long-range competition could be not only a necessary, but also a sufficient condition to the appearance of spatial structures of trees.

Inspired by [Calabrese et al., , 2010], we have proposed a nonlocal deterministic macroscopic equation for the evolution of the local density of trees where fire and tree-tree competition are the dominant mechanisms. If the kernel function falls off with distance very quickly (the Fourier transform is always positive) the system only has homogenous solutions. In the opposite case, patterns may appear depending on the value of the parameters ($\delta$ and $\sigma$), and in a sequence similar to the spatial structures appearing in standard kernel-based models. Under less favorable environmental conditions, trees tend to arrange in more robust structures to survive (Fig. 2.3(d)). Biologically, trees are lumped in dense groups, separated by empty regions. Entrance of new trees in these exclusion zones is impossible due to the intense competition they experience there.

A great strength of our approach is that our deterministic analysis is formal, and we have shown the different spatial distributions of the trees that occur as competition becomes more intense, concluding that self organization of trees is a good mechanism to promote tree survival under adverse conditions [Rietkerk et al., , 2002]. Trees tend to cluster in the high competition (low resources) limit (Fig. 2.3(d)), due to the formation of exclusion zones caused by nonlocal competition, and not as a result of facilitation. However, because we are dealing with a deterministic model, the patterns are too regular and the transition between the grass-only and a tree-populated states is independent of tree competition. We therefore considered stochasticity coming from the stochastic nature of individual birth and death events, to provide a more realistic description of savanna dynamics. Calabrese et al. [Calabrese et al., , 2010] also noted that savanna-to-grassland transition was independent of competition intensity in the mean field approach, but not when demographic noise was included. In the present model, both the grassland to savanna transition and the spatial structures that develop are influenced by demographic stochasticity. In the case of spatial structures, demographic noise is specially relevant, since it turns much of the unrealistic patterns of the deterministic model into more realistic ones, that remind the observed in real savannas. It also allows the existence of periodic arrangements of trees in more humid systems, which means environmental conditions closer to mesic savannas.

We have quantified the characteristic spacing of spatial patterns through the structure function. The irregular patterns produced by the stochastic model still have a dominant wavelength whose value is the same as in the deterministic model and depends only on the value of the range of the interaction, $R$, in the kernel function. The match between the typical spatial scale of the patterns and the characteristic distance over which nonlocal competition acts suggests that it could be responsible for the presence of clustered spatial structures. In the next chapter we will propose a competition model, neglecting every facilitative interaction, to confirm this hypothesis in arid to semiarid systems.

CHAPTER   3    Semiarid systems

Regular vegetation patterns in semiarid ecosystems have been traditionally believed to arise from the interplay between long-range competition and facilitation processes acting at smaller distances. In this chapter, it is shown that under rather general conditions, long-range competition alone may be enough to shape these patterns. To this end we propose three simple, general models for the dynamics of vegetation, that include only long-range competition between plants through a nonlocal term, where the kernel function quantifies the intensity and range of the interaction. Firstly, long-range competition is introduced influencing the growth of vegetation, secondly it is assumed to affect its death. Finally, it is considered as a term independent of the local birth-death vegetation dynamics, entering linearly in the equation. In all the situations, regardless of the way in which competition acts, we recover the vegetation spatial structures that account for facilitation in addition to competition. Models only have consider the finite range of the competition among plants, given by the length of the roots.

3.1 Introduction

Regular patterns and spatial organization of vegetation have been observed in many arid and semiarid ecosystems worldwide (Fig. 3.1), covering a diverse range of plant taxa and soil types [Klausmeier, , 1999; Rietkerk and van de Koppel, , 2008; Thompson et al., , 2009]. A key common ingredient in these systems is that plant growth is severely limited by water availability, and thus plants likely compete strongly for water [Rietkerk et al., , 2002]. The study of such patterns is especially interesting because their features may reveal much about the underlying physical and biological processes that generated them in addition to giving information on the characteristics of the ecosystem. It is possible, for instance, to infer their resilience against anthropogenic disturbances or climatic changes that could cause abrupt shifts in the system and lead it to a desert state [van de Koppel et al., , 2002; D’Odorico et al., 2006c, ; D’Odorico et al., 2006b, ].

Much research has therefore focused on identifying the mechanisms that can produce spatial patterning in water limited systems [Lefever and Lejeune, , 1997; Klausmeier, , 1999; Pueyo et al., , 2008]. An important class of deterministic vegetation models (i.e., those not considering noise associated with random disturbances) that can produce regular patterns are the kernel-based models [D’Odorico et al., 2006c, ]. These models produce patterns via a symmetry-breaking instability (i.e., a mechanism by which the symmetric-homogeneous state loses stability and a periodic pattern is created) that has its origins in the interplay between short-range facilitation and long-range competition [D’Odorico et al., 2006b, ; Rietkerk and van de Koppel, , 2008; Borgogno et al., , 2009], with field observations confirming this hypothesis in some landscapes [Dunkerley, , 2002]. Therefore it has been long assumed that both of these mechanisms must be present in semiarid systems to account for observed vegetation patterns, although quantifying the importance of each one has proven to be a difficult and contentious task [Barbier et al., , 2008; Veblen, , 2008]. A key role theory can play here is to identify the minimal requirements for pattern formation to occur. Some authors have speculated that pattern formation, under certain conditions, could occur without short-range facilitation [Rietkerk and van de Koppel, , 2008]. In this line, in the previous chapter we proposed a model for mesic savannas that includes fire and plant-plant competition as key ingredients. Fire introduces a positive feedback so that this model considers both competition and facilitation mechanisms. However, the model still produces regular patterns even when the facilitative interaction, fire, is considered at its very short-range (in fact, local) limit. These considerations suggest that local facilitation may be superfluous for pattern formation, and that a deeper exploration of the range of conditions under which pattern formation can occur in the absence of facilitation is therefore warranted. This will be the major concern of this chapter.

We will study simple, but quite general, single-variable models for vegetation density in water-limited regions. Only competitive interactions are considered, modeled in different ways to analyse if patterns depend on how competition enters in the dynamical equations Martínez-García et al., [2013]; Martínez-García et al., 2014a . The role of nonlinearities is also investigated. We show that when only a single broadly applicable condition is met, that competitive interactions have a finite range, the full set of regular patterns formerly attributed to the interaction between short-range facilitation and long-distance competition can be produced in the absence of facilitation.

3.2 Competition in a nonlocal nonlinear birth term

Arid and semiarid ecosystems are typified by patches of vegetation interspersed with bare ground. Water is a very limited resource for which juvenile plants must compete with those that have already established. Logistic-type population models have been used in a wide variety of applications including semiarid systems and, as shown in the previous chapter, savannas [Calabrese et al., , 2010]. They thus form a reasonable and very general starting point. Specifically, we consider the large-scale long-time description of the model in terms of a continuous-time evolution equation for the density of trees, $\rho({\bf x},t)$. Death occurs at a constant rate $\alpha$, whereas population growth occurs via a sequence of seed production, dispersal, and seed establishment processes. Seed production occurs at a rate $\beta_{0}$ per plant. For simplicity we consider dispersal to be purely local and then if all seeds would give rise to new plants the growth rate would be $\beta_{0}\rho({\bf x},t)$. But once a seed lands, it will have to overcome competition in order to establish as a new plant. We consider two different competition mechanisms. First, space availability alone limits density to a maximum value given by $\rho_{max}$. Thus, $0\leq\rho({\bf x},t)\leq\rho_{max}$. The proportion of available space at site ${\bf x}$ is $1-\rho({\bf x},t)/\rho_{max}$, so that the growth rate given by seed production should be reduced by this factor. Second, once the seed germinates, it has to overcome competition for resources with other plants. This is included in the model by an additional factor $r=r(\tilde{\rho},\delta)$, $0\leq r\leq 1$, which is the probability of overcoming competition. This probability decreases with increasing average vegetation density within a neighborhood $\tilde{\rho}$, and the strength of this decrease depends on the competition intensity parameter, $\delta$. Higher values of $\delta$ represent more arid lands, and thus stronger competition for water. In the following, we measure density in units so that $\rho_{max}=1$. Combining all processes, the evolution equation for the density then takes the form:

[TABLE]

$\tilde{\rho}=\tilde{\rho}({\bf x},t)$ is the nonlocal density of vegetation that is obtained by averaging (with a proper weighting function) the density of plants in a neighborhood:

[TABLE]

where $G({\bf x})$ is a normalized kernel function, which accounts for the weighted mean vegetation density, and defines the neighborhood of the plant. A Laplacian term could be included in the right-hand side of Eq. (3.1) as a way to model seed dispersal, but doing so would not qualitatively change our results, so we have left it out.

We have presented a phenomenological derivation of the model. An open problem is to infer this type of description from a mechanistic one where the explicit interactive dynamics of vegetation competing for water is considered. Preliminary results on this derivation are shown in the Appendix C.

In previous kernel-based vegetation models [Lefever and Lejeune, , 1997; D’Odorico et al., 2006c, ], the kernel function contained information on the class of interactions present in the system, that were both facilitative and competitive. That is, it can take positive and negative values either enhancing or inhibiting the vegetation growth. On the contrary, we introduce purely competitive interactions through the nonlocal function $r(\tilde{\rho},\delta)$, where the kernel is always positive and thus defines the area of influence of a focal plant, and how this influence decays with distance. These two kernels are compared in Fig. 3.2. Competition is included by assuming that the probability of establishment $r$ decreases with increasing vegetation density in the surroundings:

[TABLE]

As $\delta$ modulates the strength of the competition, it must be that $r(\tilde{\rho},\delta=0)=1$, and $r(\tilde{\rho},\delta\rightarrow\infty)=0$. This means that when water is abundant ($\delta=0$) competition for water is not important ($r=1$), whereas new plants cannot establish in the limit of extremely arid systems, $\delta\rightarrow\infty$.

Note the generality of this vegetation competition model: a spatially nonlocal population growth term of logistic type with rate fulfilling Eq. (3.3), and a linear death term. A complete description of our model should specify both the kernel function $G$ as well as $r$, but we can go further with the analysis in general terms.

The possible homogenous stationary values of the density for Eq. (3.1) are: a) no vegetation $\rho=0$, and b) the vegetated state $\rho=\rho_{0}$. The system will show either one or the other depending on the relationship between the seed production and death rates, $\beta_{0}$ and $\alpha$ [Calabrese et al., , 2010]. The non-trivial homogeneous stationary solution, $\rho_{0}$, can be obtained by solving

[TABLE]

that has only one solution in the interval $\rho_{0}\in[0,1]$ because of the conditions imposed on the function $r$ in equation(3.3). We now ask if this stationary solution gives rise to periodic structures via a symmetry-breaking instability as happens in other models that include not only competition but also facilitation mechanisms in the interactions [Borgogno et al., , 2009]. To explore this possibility in our model, we perform a linear stability analysis [Cross and Hohenberg, , 1993] introducing a small perturbation to the stationary solution, so $\rho({\bf x},t)=\rho_{0}+\epsilon\psi({\bf x},t)$, with $\epsilon\ll 1$. Technical details of this derivation may be found in Appendix A. The perturbation growth rate is

[TABLE]

where $\hat{G}({\bf k})$ is the Fourier transform of the kernel, $\hat{G}({\bf k})=\int G({\bf x})\exp(i{\bf k}\cdot{\bf x})d{\bf x}$, and $r^{\prime}(\rho_{0},\delta)\equiv\left(\frac{\partial r}{\partial\tilde{\rho}}\right)_{\tilde{\rho}=\rho_{0}}$.

Patterns appear if the maximum of the growth rate (i.e., of the most unstable mode), $\lambda(k_{c})$, is positive, which means that the perturbation grows with time. From Eq. (3.5), this is only possible if the Fourier transform of the kernel function, $\hat{G}({\bf k})$, takes negative values, since $r^{\prime}(\rho_{0},\delta)<0$. As has been in the previous chapter, this happens, for example, for all stretched exponentials $G(|{\bf x}|)\propto\exp{(-|{\bf x}/R|^{p}})$ with $p>2$, where $R$ is a typical interaction length [Pigolotti et al., , 2007, 2010]. Kernels satisfying this criterion have broader shoulders and shorter tails (i.e., are more platykurtic) than the Gaussian function, which is obtained for $p=2$. In reality, any competitive interaction among plants will have finite range because their roots, which mediate the interaction, have finite length. The interaction range $R$ between two plants will be twice the typical root length. Kernels with finite range can, in general, be modeled by considering a truncated function such that $G(|x|)=CF(|x|)\Pi(|x|)$, where $C$ is a normalization constant, $\Pi(x)$ is a unit-step function defined as $\Pi(x)=1$ if $|x|\leq R$ and $\Pi(x)=0$ if $|x|>R$, and $F(|x|)$ is a function of the distance that models the interactions among the plants. Because of the finite range in the kernel function, the Fourier transform will show oscillations and thus will always take negative values. The functional form of the probability of surviving the competition, $r(\tilde{\rho},\delta)$, changes only the parameter regime where patterns first develop, but they will appear in the system, regardless of its form, for $r^{\prime}(\rho_{0},\delta)/r(\rho_{0},\delta)$ large enough.

For the rest of our analysis, we will use $F(x)=1$, so the kernel is given by $G(x)=1/\pi R^{2}$ if $|x|\leq R$ and $G(x)=0$ if $|x|>R$, which defines an interaction area of radius $R$ (that is, roots of typical length $R/2$). Its Fourier transform (in two dimensions) is

[TABLE]

where $J_{1}(|{\bf k}|R)$ is the first-order Bessel function. We will further specify the model by assuming particular forms for the growth rates. Let us consider a probability of surviving competition given by

[TABLE]

with $q>0$. In the particular case of $q=1$, the homogeneous density, $\rho_{0}$, and the perturbation growth rate, $\lambda$, can be obtained analytically. Numerical evaluations must be done if $q\neq 1$. In the following, for simplicity, we consider the case $q=1$ and only briefly discuss other values. The nontrivial stationary solution, $\rho_{0}\neq 0$, can be obtained analytically

[TABLE]

where $\beta_{0}\geq\alpha$. Eq. (3.8) shows that the homogeneous density of trees in the stationary state decays as $\sim\delta^{-1}$ with increasing competition strength (i.e., large $\delta$). It can be analytically shown that the same dependence of $\rho_{0}$ on large $\delta$ occurs for any value of $q$.

From Eq. (3.5), the growth rate of the perturbations can also be calculated

[TABLE]

and is shown in Fig. 3.3 (Left) for different values of the competition strength. When the growth rate of the most unstable mode (i.e. the maximum of $\lambda(k)$), $k_{c}$, becomes positive, patterns emerge in the system [Borgogno et al., , 2009]. To obtain the critical value of the competition parameter at the transition to patterns, $\delta_{c}$, we have to calculate the most unstable mode as the first extreme of $\lambda(\bf k)$ at $k\neq 0$, i.e., the first zero of the derivative of $\hat{G}(\bf k)$. This was already done in the previous chapter, and gives $k_{c}=5.136/R$. It depends only on the the range defining $G(\bf r)$). This value changes depending on the kernel, but in the case of kernels with a finite range (i.e. truncated by a unit step function of radius $R$) it is always of this order. The critical wavenumber is determined mainly by the contribution of the unit step function to the Fourier transform, which is always the same. This result is also independent of the other parameters of the system, and shows that the nonlocal competition mechanism is responsible for the formation of patterns in the system.

To identify the parameter values for the transition to patterns, we solve $\lambda({\bf k}_{c})=0$ in Eq. (3.9), which shows that patterns emerge when competition strength exceeds $\delta_{c}=-\beta_{0}/\alpha\hat{G}(\mathbf{k_{c}})$, which is positive because $\hat{G}(\mathbf{k_{c}})<0$. Figure 3.3 (Right) shows the growth rate of the most unstable mode as a function of competition strength for different values of the exponent $q$ for fixed values $\beta_{0}=1$, and $\alpha=0.5$. Note that the critical value of the competition parameter depends on the functional form of $r$. This dependence could be used to tune the value of $q$ to have a realistic competition strength for the transition to patterns, provided that one has sufficient data.

We can also explain the separation length between clusters of plants using ecological arguments adn expanding the concept of exclusion areas that was mentioned in the previous chapter. Consider a random and inhomogeneous distribution of plants. Maxima of this distribution identify places with the highest plant density. Imagine that two such maxima occur at a distance larger than $R$ but smaller than $2R$ from each other. There will be no direct interaction between the roots of plants in these different patches because they are separated by a distance larger than the interaction range $R$ (twice the root extension as first order approximation). But there is an area in-between which is simultaneously within the range of both patches. Compared with plants occurring inside a cluster, which only have to compete with plants in their own cluster, those that occur in-between clusters will experience stronger competition and will therefore tend to disappear (Fig. 3.4). We call these regions exclusion areas [Hernández-García and López, , 2004; Pigolotti et al., , 2007, 2010]. The disappearance of plants in these exclusion areas in turn reduces competition on the two well-populated patches, so that a positive feedback appears reinforcing the establishment of plants in patches periodically separated with a distance between $R$ and $2R$. We stress again that competition alone is responsible for the symmetry breaking instability, and no facilitative interactions are needed for pattern formation.

This mechanism does not work only in the particular case of considering competition among plants modifying the birth rate. It also leads to the formation of patterns if nonlocal competition is introduced in the death rate or entering linearly in the dynamics. This will be shown in Sections 3.3 and 3.4 respectively.

Finally, we have numerically integrated Eq. (3.1) in a patch of $10^{4}$ m2 with periodic boundary conditions and a competition range of $R=8$ m. Time stepping is done with an Euler algorithm. The results (see Fig. 3.5) exhibit steady striped and spotted vegetation patterns typically arising from symmetry breaking.

Similar results can be obtained for different growth rates, for example considering a family of stretched exponentials in the probability of overcoming long-range competition $r(\tilde{\rho},\delta)=\exp\left(-\delta\tilde{\rho}^{s}\right)$.

This gives a perturbation growth rate of the form

[TABLE]

where the stationary density, $\rho_{0}$ has to be obtained numerically from Eq. (3.4). The value of the competition parameter at the transition to patterns $\delta_{c}$ also has to be obtained numerically for a given set of values of $s$, $\alpha$, and $\beta_{0}$. This critical value is shown in Fig. 3.6 using a top-hat kernel for different values of $s$. It is represented the value of the perturbation growth rate of the most unstable mode $\lambda({\bf k_{c}})$ as a function of the competition strength $\delta$.

This further confirms our result that competition is the only necessary ingredient for the formation of vegetation patterns in the present framework, and that this does not depend on the functional form of the probability of surviving competition (growth rate) provided it verifies the requirements given by Eq. (3.3).

3.3 Competition in a nonlocal nonlinear death term

In the previous section we have considered that plant death occurs at a constant rate and the birth process takes place following a sequence of seed production, dispersal and establishment. The establishment of new seeds has been studied as a process depending on the density of vegetation in a given neighborhood. However, it is possible to consider it affecting the death rate.

In this framework, the population growth has two stages, namely seed production at rate $\beta_{0}$ and local seed dispersal. This means that once a seed is produced by a plant, it still has to compete for available space, limiting the density of vegetation to a maximum value $\rho_{max}$. As it was done before, we choose units so that $\rho_{max}=1$. On the other hand, the death of plants is influenced by the availability of resources and thus include nonlocal interactions. When water is abundant the competition for it is not relevant and plants die at a typical constant rate $\alpha_{0}$. However, the scarcity of resources promotes the death of vegetation, and this is included in the model by an additional factor $h$ in the death term, the probability of dying because of competition, $0\leq h\leq 1$. That is, the probability of not being able to overcome competition $h(\tilde{\rho},\delta)=1-r(\tilde{\rho},\delta)$, where $r$ is the function introduced in Sec. 3.2. The probability of dying because of competition for resources has to increase with the density of vegetation

[TABLE]

and it must tend to its maximum value in the limit of extremely arid systems $(\delta\rightarrow\infty)$ and to vanish when water is not a constraint for vegetation growth $(\delta\rightarrow 0)$. These properties can be derived from the properties of the function $r$ and its relationship with $h$.

Under these conditions, the model equation for the density of vegetation is

[TABLE]

where $\tilde{\rho}({\bf x},t)$ is the nonlocal density of vegetation at ${\bf x}$, $\tilde{\rho}({\bf x},t)=\int\rho({\bf x^{\prime}},t)G(|{\bf x}-{\bf x^{\prime}}|)dx^{\prime}$, and $G({\bf x})$ is the kernel function that defines an interaction range and modulates its strength with the distance from the focal plant, as it was in Sec. 3.2. Also following the step of this previous section, we study the existence of patterns in this model.

First of all, the stationary solutions of Eq. (3.12), $\rho_{0}$, are obtained solving

[TABLE]

that has a trivial solution, $\rho_{0}=0$ referring to the desert state. The vegetated state must be obtained from

[TABLE]

once the function $h$ has been chosen.

Second of all, the formation of patterns in the system has to be studied through a linear stability analysis of Eq. (3.12), introducing a small perturbation to the stationary homogeneous state, $\rho_{0}$. Considering that $\rho({\bf x},t)=\rho_{0}+\epsilon\psi({\bf x},t)$, with $\epsilon\ll 1$, the perturbation evolves according to the following linear integro-differential equation

[TABLE]

that can be solved using the Fourier transform to obtain the growth rate of the perturbation,

[TABLE]

where, again, $\hat{G}({\bf k})$ is the Fourier transform of the Kernel function. Using Eq. (3.13), one finally gets

[TABLE]

that shows that, as it was in the case of models with nonlocal interactions in the birth rate, that the patterns appear only when the Fourier transform of the kernel takes negative values. To investigate the particular behavior of this model in a simple situation, one needs to choose a function $h$.

The simplest case that allows a complete analysis is to choose a probability of dying because of competition growing linearly with the nonlocal density of vegetation,

[TABLE]

The behavior of this function $h$ is shown in the left panel of Fig. 3.7 for different values of $\delta$. The choice of this functional form in $h$ restricts the domain of the parameter $\delta$. It has to be $\delta\in[0,1]$ since it must be $h\leq 1$ to represent a probability.

The nontrivial stationary solution, Eq. (3.14), is

[TABLE]

This solution has the proper behavior in the limits $\delta\rightarrow 0$ and $\delta\rightarrow\infty$. The growth rate of the perturbation is

[TABLE]

from where we obtain a transition to pattern at,

[TABLE]

where $k_{c}$ is the most unstable mode, i.e., that with the highest growing rate. Eq. (3.21) limits the possible value of the birth and the death rates for a given kernel. They have to take values such as $\delta_{c}\leq 1$. This is a new condition that was not present in the model presented in Sec. 3.2, and appears because of the constraints imposed on the values of the competition parameter by the linear function $h$. For $\alpha_{0}=1$, $\beta=0.1$ and a top-hat kernel of radius $R=8$, patterns emerge for a competition strength $\delta_{c}\approx 0.747$. As it can be observed in Fig. 3.8, the homogeneous distribution is the stationary solution when $\delta<\delta_{c}$ (Panel a), while patterns (stripes and spots) appear otherwise (Panels b and c respectively).

A different function for the probability of dying might be chosen 111It is important to justify that $\delta$ has no dimension, and there is a saturating density, $\kappa$, that we have set to $\kappa=1$ so that $r=\frac{\delta\tilde{\rho}}{\delta\kappa+\tilde{\rho}}$.,

[TABLE]

that is shown in the right panel of Fig. 3.7 for different values of $\delta$. Although this makes the analysis much more complicated, it allows the competition strength parameter $\delta$ to take any positive real value. We present some results showing that the functional form of $h$ is not relevant for pattern formation.

The nontrivial stationary solution is

[TABLE]

where we have chosen the positive sign in the square root because it must be $\rho_{0}>0$. The growth rate of the perturbation is

[TABLE]

Plugging Eq. (3.23) into Eq. (3.24) and considering the most unstable mode, ${\bf k_{c}}$, we can obtain an expression for the value of the competition parameter at the transition to patterns,

[TABLE]

Using $\alpha_{0}=1$, $\beta=0.1$ and a constant top-hat kernel of radius $R=8$ one obtains $\delta_{c}\approx 0.663$. This result can be checked integrating numerically Eq. (3.12) with these values of the rates and the interaction range. We do not show this patterns again because they are equivalent as the ones in Fig. 3.8.

3.4 Competition in a nonlocal linear death term

The study that we are going to present in this section is inspired by the kernel-based models reviewed in [Borgogno et al., , 2009], where it is assumed that both competitive and facilitative interactions between plants must come into play for pattern formation. These models are not well posed, since the density of vegetation can become negative at some spatial points because the long range interactions entering linearly. This incovenience is often overcome assuming that whenever the density becomes negative, it takes a value $\rho({\bf x},t)=0$. In this section, we extend these models to the situations where there is only competition among trees, and also study the effect of seed dispersal through a diffusive term. The vegetation density changes in time because of its local dynamics (logistic growth) and the spatial interactions (competition) with other points of the space,

[TABLE]

where $\kappa$ is the carrying capacity of the system (i.e, the maximum density of vegetation in the absence of competition) and $\Omega$ is the interaction strength parameter. To account for purely competitive interactions $\Omega$ and the Kernel function $G({\bf x})$ are positive. It is important to remark, once more, that the dynamics of this model is not bounded at $\rho(x,t)=0$, so negative values must be avoided by imposing that $\rho=0$ when it becomes negative in the numerical simulations.

The stationary solutions are, $\rho_{0}=0$ (no vegetated state), and

[TABLE]

that limits the values of the interaction strength parameter. It must be $\Omega\leq 1$ to verify that $\rho_{0}\geq 0$, while in the limit $\Omega=0$, where the are not long range interactions, the population of plants grows logistically until it reaches its carrying capacity.

The linear stability analysis of Eq. (3.4) gives a peturbation gorwth rate,

[TABLE]

that, in the absence of diffusion, $D=0$, allows the existence of patterns if

[TABLE]

where $k_{c}$ is the most unstable mode. Provided that we choose, again, a normalized constant kernel of width $2R$ the transition to patterns appears at $\Omega_{c}\approx 0.47$. Integrating the full Eq. (3.4) including long-range seed dispersal $(D\neq 0)$ laberynthic and spotted patterns emerge. This result is shown in Figures 3.9(a) and 3.9(b) for different values of $\Omega$.

3.5 Summary and conclusions

In this chapter we have presented several models of vegetation in water-limited regions where, despite the absence of facilitative interactions, regular structures may still appear: stripes and spots of vegetation interspersed on the bare soil forming a hexagonal lattice. Patterns consisting on spots of bare soil characteristic of models with competition and facilitation are not observed when positive interactions are neglected. In fact, in the previous chapter, where we studied a model for mesic savannas with local facilitation, the gapped distributions appear only when the system is on the transition line to patterns. This result is different from what is observed in models including nonlocal facilitation, that show the whole sequence of patterns (gaps, stripes and spots), with each type of structure in a wide region of the parameter space. This suggests that facilitative interactions, although not indispensable for the formation of patterns in arid systems, could be relevant to substain some of the structures that have been reported in field observations. Future work is planned to address this problem.

From a mathematical point of view, nonlocality enters through an influence function that determines the number of plants competing within a range with any given plant. A first-order approximation of this distance can be given by (twice) the typical length of the roots, but field measurements are needed in order to determine the resources uptaking range of a given plant species. A necessary condition for pattern transitions, common to all the presented models, is the existence of negative values of the Fourier transform of the kernel function, which always happens, but not only, for kernel functions with finite range.

From a biological point of view, competitive interactions alone may give rise to spatial structures because of the development of regions (typically located between maxima of the plant density) where competition is stronger, preventing the growth of more vegetation.

An unfortunate consequence of the universal character of these models is that the information that is possible to gain on the underlying biophysical mechanisms operating in the system just by studying the spatial structures of the vegetation is limited. Many different mechanisms lead to the same patterns. Although they are universal models should be specific to each system, emphasizing the importance of empirical studies. Field work may help theoretical efforts by placing biologically reasonable bounds on the shape and extent of the kernel functions used in the models, and also by approximations to the probablity of overcoming competition, $r(\tilde{\rho},\delta)$, or of dying because of it $h(\tilde{\rho},\delta)$.

Finally, with this work, we aim at showing that competitive interactions alone may be responsible of patterns in arid to semiarid systems, regardless of how they are introduced in the different modeling approaches. Under certain conditions, nonlocal competition alone may be responsible of the formation of patterns regardless of entering in the equations either linearly or nonlinearly. We hope that our results shed light on the task of understanding the fundamental biological mechanisms -and the possible absence of facilitation- that could be at the origin of pattern formation in semiarid systems.

Appendix A Linear stability analysis

In chapters 2 and 3, we have referred many times to the linear stability analysis of the equations in order to gain insight into the possible emergence of patterns. This is a broadly used technique. Its objective is to obtain the temporal evolution of small perturbations to the stationary homogeneous state of the system. To this aim we assume idealized conditions that can be summarized in three points: choosing a uniform geometry for the walls that confine the medium; replacing the lateral boundaries with ones consistent with translational symmetry; and choosing to study the linear instability towards time-independent patterns that are consistent with translational symmetry. We impose periodic boundary conditions and study the formation of stationary patterns.

Our starting point is Eq. (3.1),

[TABLE]

Considering that the system is in the homogeneous stationary state, $\rho({\bf x},t)=\rho_{0}$, how will perturbations to this state evolve with time? Will they vanish, taking the system back to the homogeneous state, or will they grow up allowing the formation of stationary spatial structures? The key quantity will be the growth rate of the perturbation, restricted to the linear regime, so that only small perturbations in the short time limit will be considered. If a small perturbation is done to the stationary homogeneous state, the density is

[TABLE]

with $\epsilon\ll 1$. Plugging Eq. (A.2) into Eq. (A.1), it is

[TABLE]

where the spatiotemporal dependence in the perturbation $\psi$ is omitted to simplify the notation. We have also used that $\rho_{0}$ is stationary and thus its time derivative equals to zero together with the normalization of the kernel to simplify the dependence on $\rho_{0}$ in the function $r$.

Eq. (A.3) can be simplified because we are interested in small perturbations. The function $r$ can be expanded in a Taylor series and higher order terms in $\epsilon$ neglected. This gives a linearized integro-differential equation for the evolution of $\psi$,

[TABLE]

that can be solved using Fourier transform. Notice that the stationary solution of Eq. (A.1),

[TABLE]

has been used to obtain Eq. (A). We apply the Fourier transform to Eq. (A),

[TABLE]

where $\hat{\psi}({\bf k},t)=\int{\rm e}^{ik\cdot x}\psi({\bf x},t)d{\bf x}$ is the Fourier transform of the perturbation, and equivalently, $\hat{G}(k)$ is the Fourier transform of the kernel.

Finally, Eq. (A.6) is solved considering $\hat{\psi}({\bf k},t)\propto\exp(\lambda({\bf k})t)$, so that a linear growth rate of the perturbation is obtained

[TABLE]

Appendix B Numerical integration of Eq. (2.19)

The integration of stochastic equations where the noise amplitude depends on the square root of the variable, $\rho$, and there are absorbing states (i.e, states where the system stays indefinitely), has awaken a great interest, specially in the study of critical phenomena (i.e, properties of the system that appear when it is close to the critical point). The amplitude of the fluctuations tends to zero there, and thus numerical instabilities may appear. Recently [Dornic et al., , 2005; Pechenik and Levine, , 1999] a very efficient method has been developed, but we have used in this work an older one, presented in [Dickman, , 1994], since its implementation is easier and it gives precise results working far from the transition point. It consists on discretizing the Langevin equation, taking a step size $\Delta\rho$ in the variable.

First of all, to integrate the Eq. (2.22) we discretize the space, whose area is $L_{x}\times L_{y}$, in a $N_{x}\times N_{y}$ regular grid and compute the integral in the exponential term. It is approximated by a sum of the field evaluated in the nodes

[TABLE]

where $\Delta x$ and $\Delta y$ are the spatial steps in both directions. They are given by $\Delta_{x,y}=L_{x,y}/(N_{x,y}-1)$.

Then, we integrate the temporal dependence. The key of the algorithm is to prevent $\rho+\Delta\rho$ to take negative values. From a general equation

[TABLE]

where $\eta(t)$ is a Gaussian white noise with zero mean and delta correlated, it is

[TABLE]

where $\Delta W=\sqrt{\Delta t}Y$. $Y$ is a Gaussian number with zero mean and unit variance. To avoid negative values of $\rho+\Delta\rho$, Dickman, [1994] proposes to dicretize the density setting $\rho=n\rho_{min}$ and to simetrically truncate the Gaussian distribution from where $Y$ is drawn, so that $|Y|\leq Y_{max}$. The negative values are avoided requiring

[TABLE]

This can be seen adding $\rho$ in both sides of Eq. (B.3), considering the discretization of the density and assuming that $Y$ takes its most negative value according to the constraint (B.4), so

[TABLE]

that is always positive for $n$ integer.

Eq. (B.4) can be verified in many ways but, following Dickman, [1994] again, we use

[TABLE]

Finally, rescaling the equation, we can achieve a discretized version in which positive and zero-mean noise are ensured at the cost of a “quantized” density.

Appendix C Derivation of the effective nonlocal description from tree-water dynamics

In this section we present a preliminar (and not fully satisfactory) attempt to derive the model presented in Section 3.2 of Chapter 3 (the derivation corresponding to the nonlocal death model in 3.3 is a straightforward extension of this calculation).

Let us consider a system involving dimensionless vegetation density, $\rho({\bf x},t)$, and soil-water $w({\bf x},t)$. The dynamics is purely local and competitive and takes the form:

[TABLE]

where the nondimensional positive parameters are: the seed production rate $\beta$; the vegetation death rate $\alpha$; the consumption rate of water by vegetation, $\mu$; the evaporation rate $\gamma$, and the rainfall, $I$. Water percolation in the ground is modeled by a diffusion constant $D_{w}$. Note that this model is a simplified version, which only includes competitive interactions, of the model presented in Gilad et al., [2004].

Since the characteristic time scale of the water is much faster than the one of the biomass we can do an adiabatic elimination of the variable $w$ (i.e. $\frac{\partial w}{\partial t}=0$) so that

[TABLE]

and thus

[TABLE]

whose formal solution can be obtained using Green’s functions, $G_{d}$,

[TABLE]

with the boundary conditions $w(x\rightarrow\pm\infty)=0$. For simplicity we now consider a one-dimensional situation, although analogous calculations can be done in two dimensions. The Green’s function is the solution of

[TABLE]

and it is given by

[TABLE]

Taking the nondimensional small number $\mu$ as the perturbative parameter, we can further obtain an approximate expression for $w$ from Eq. (C.5)

[TABLE]

where $G_{d0}=\int G_{d}({{\bf x}}-{{\bf s}})d{{\bf s}}<0$, since the Green’s function is always negative. Plugging this in the equation for the biomass density (C.1), we obtain the closed expression:

[TABLE]

Defining the positive nonlocal density $\tilde{\rho}=\int G_{c}({{\bf x}}-{{\bf s}})\rho({{\bf s}})d{{\bf s}}$, where $G_{c}=-G_{d}$, we can write Eq. (C.9) as

[TABLE]

where we have defined $\mathrm{bar}r(\tilde{\rho})=I|G_{d0}|\left(1-\mu\tilde{\rho}\right)$.

To have a good agreement with the effective nonlocal dynamics Eq. (3.1), $\mathrm{bar}r>0$ since it represents a probability. This is certainly the case for small $\mu$. Note that some additional conditions on the normalization of the Green’s function have to be imposed to limit $r$ to values less than $1$. Also $\mathrm{bar}r^{\prime}(\tilde{\rho})=-I\mu|G_{d0}|$ is always negative, as we expected.

In this particular example we obtained the exponential kernel of Eq. (C.7), which does not have the finite-range support that would be associated to the finite root extent. As a consequence, the Fourier transform of this kernel has no negative components and then does not lead to pattern formation. The simple modeling of water dispersion by means of a diffusion constant does not contain the additional spatial scale associated to root size, and should be replaced by some mechanism implementing root effects. On the other side, the finite-range of the kernel is a sufficient but not a necessary condition for its Fourier transform to have negative values. It is well-known the existence of infinite-range kernels whose Fourier transform has negative values. This is the case of all stretched exponentials $G(x)\propto\exp({-|{\bf x}|^{p}})$ with $p>2$ as has been already mentioned several times before. Kernels satisfying this are more platykurtic than the Gaussian function. Work is in progress along this line.

Part III  **Animal mobility

**

CHAPTER   4    Optimal search in interacting populations

In this chapter we investigate the relationships between search efficiency, movement strategy, and nonlocal communication in the biological context of animal foraging. We consider situations where the members of a population of foragers perform either Gaussian jumps or Lévy flights, and show that the search time is minimized when communication among individuals occurs at intermediate ranges, independently of the type of movement Martínez-García et al., 2014b ; Martínez-García et al., 2013b . While Brownian strategies are more strongly influenced by the communication mechanism, Lévy flights still result in shorter overall search durations.

4.1 Introduction

Situations where a single individual or a group of searchers must find an object (target) appear in many different fields including chemistry [Bénichou et al., , 2011], information theory [Pirolli and Card, , 1999], and animal foraging [Viswanathan et al., , 2011; Méndez et al., , 2014]. The study of these searching problems has generated an increasing amount of work in the last years, many of them oriented towards the identification of efficient strategies [Bénichou et al., , 2011; Viswanathan et al., , 2011; Vergassola et al., , 2007; Bénichou et al., , 2006]. Many remarkable examples can be found in the context of biological encounters, such as proteins searching for targets on DNA [Taylor and Halford, , 1989], or animals searching for a mate, shelter or food [Campos et al., , 2013; Bénichou et al., , 2011; Viswanathan et al., , 2008; Mejía-Monasterio et al., , 2011; Bartumeus et al., , 2003]. In these cases, the search time is generally limiting and minimizing it can increase individual fitness or reaction rates.

The optimality of a search strategy depends strongly on the nature of both the targets and the searchers [Bartumeus et al., , 2005, 2002]. In the context of animal foraging, which is our focus here, searchers may move randomly, may use memory and experience to locate dispersed targets or they may also combine random search with memory-based search. In highly social species, groups of searchers may share information when no single individual is sufficiently knowledgeable. This is based on the many wrong hypothesis [Hoare et al., , 2004; Torney et al., , 2009], that states that error in sensing of individuals can be reduced by interacting with the rest of the group, where all individuals can act as sensors.

It is well known that individual movement plays a central role in search efficiency, and many studies have focused on the comparative efficiency of Lévy and Brownian movement strategies [Bénichou et al., , 2011; Viswanathan et al., , 2008; Bartumeus et al., , 2005, 2002]. Lévy flights are more efficient in some random search scenarios [Viswanathan et al., , 1999, 2000], but whether or not they are used in real animal search strategies is still an open and contentious topic [Edwards, , 2011; Edwards et al., , 2007]. Much less effort, however, has been spent on trying to understand the long-range (i.e. nonlocal) interaction mechanisms among social searchers. While diverse observations suggest that such interactions occur in many taxa, including bacteria [Liu and Passino, , 2002], insects, and mammals [Zuberbühler et al., , 1997; McComb et al., , 2003], previous studies have focused almost exclusively on how the collective movements of a group of animals can emerge from local interactions among individuals [Mishra et al., , 2012; Kolpas et al., , 2013; Couzin et al., , 2002]. The distance at which communication can be maintained strongly depends on the species. A variety of mammalian species are known to communicate acoustically over distances of up to several kilometers [McComb et al., , 2003; Mathevon et al., , 2010; Estes, , 1991], but while group formation via vocalizations has been well studied [McComb et al., , 2003; Ford, , 1991; Smith et al., , 2011], incidental benefits such as increased foraging efficiency have received little research attention. Therefore many open questions remain, particularly on the interrelation between communication and optimal search for resources. How can communication facilitate group formation and identification of areas of high quality resources? Does a communication range exist that optimizes foraging efficiency? To what degree does search efficiency depend on the communication mechanism? Finally, how does communication affect individual space use in a heterogeneous environment?

Here, we present a theoretical model to study the effect of nonlocal communication on the search efficiency of a group of individuals. Communication and the amount of shared information are directly connected through the distance travelled by the signals emitted by successful searchers (communication range). To introduce the model, we start from an individual based description and compare the situations where the individuals employ either Lévy or Brownian random search strategies. In the first section we study many properties of the model such as the existence of an optimal communication range and the influence of the distance between targets and the movement strategy on its value. For tractability, we will next consider a simplified, one-dimensional version of the model and compute analytically the search time for both Brownian and Lévy searchers as a function of the communication length scale. This simplified model allows us to unveil the dependence of this time on both the parameters governing individual mobility, and on the distance between targets. Finally, a continuum description in terms of an equation for the density of searchers is derived. The chapter ends with a discussion comparing Brownian and Lévy strategies.

4.2 The Individual Based Model for Brownian searchers

We consider $N$ particles which undergo a $2D$ Brownian random walk. Correlated random walks, often more appropriate to model directional persistence in animal movement, reduce to Brownian motion for large spatiotemporal scales [Turchin, , 1998]. The movement is biased by the gradients of the landscape quality (local information), and by the interaction among individuals through a communication mechanism that is activated when good resources are found, thus providing information on habitat quality in other areas (nonlocal information). The dynamics of any of the particles $i=1,...,N$ is

[TABLE]

where $\mathbf{\eta}_{i}(t)$ is a Gaussian white noise term characterized by $\langle\eta_{i}(t)\rangle=0$, and $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2D\delta_{ij}\delta(t-t^{\prime})$, with $D$ the diffusion coefficient. The term $B_{g}\nabla g(\mathbf{r}_{i})$ refers to the local search, where $g(\mathbf{r})$ is the environmental quality function (amount of grass, prey, etc…) and $B_{g}$ is the local search bias parameter. $B_{C}\nabla S(\mathbf{r}_{i})$ is the nonlocal search term, with $B_{C}$ the nonlocal search bias parameter and $S(\mathbf{r}_{i})$ is the available information function of the individual $i$. It represents the information arriving at the spatial position of the animal $i$ as a result of the communication with the rest of the population. This term makes the individuals move along the gradients of the information received. It is a function of the superposition of pairwise interactions between the individual $i$ and each one of its conspecifics,

[TABLE]

$F$ is an arbitrary perception function that must be set in each application of the model, $V(\mathbf{r}_{i},\mathbf{r}_{j})$ is the interaction between the receptor particle $i$ depending on its position $\mathbf{r}_{i}$ and the emitting particle fixed at $\mathbf{r}_{j}$, and $A[g(\mathbf{r}_{j})]$ is the activation function (typically, a Heaviside function) that indicates that the individual at $\mathbf{r}_{j}$ calls the others if it is in a good habitat.

We begin with Monte Carlo simulations of the individual based dynamics in Eq. (4.1) using a square system ($L_{x}=L_{y}=1$) with periodic boundary conditions, and a population of $N=100$ individuals. We use a theoretical landscape quality function, $g(\mathbf{r})$, consisting of three non-normalized Gaussian functions, to ensure that $g(\mathbf{r})\in[0,1]$, centered at different spatial points (Fig. 4.1). As a first approach, we consider a Gaussian-like interaction kernel. Manipulating its typical range via the standard deviation, $\sigma$, the information sent by a successful searcher will travel farther and the number of individuals with which the communication is established will be larger. The available information function of the individual $i$ depending on its position is

[TABLE]

where, as mentioned before, $A[g(\mathbf{r})]$ is a theta Heaviside function that activates the interaction when the quality is over a certain threshold $\kappa$, $A[g(\mathbf{r})]=\Theta(g(\mathbf{r})-\kappa)$. The perception function has been chosen as the identity for simplicity $F=I$.

The question is how the typical communication distance affects the average efficiency of individuals searching for targets in space (areas of high-quality forage). We give an answer in terms of spatial distributions of individuals at long times starting from a random initial condition and the mean first arrival time to the targets, $T$, as it is done in related works [Bénichou et al., , 2005]. This quantity (Fig. 4.2) may be optimized with a communication range parameter, $\sigma$, of intermediate scale. As it was said before, the number of individuals from which a given animal receives a signal will typically increase with the interaction scale. When this scale is too small, individuals receive too little information (no information when $\sigma=0$), and thus exhibit low search efficiency (Fig. 4.2). Similarly, interaction scales that are too large lead to individuals being overwhelmed with information from all directions, also resulting in inefficient search (Fig. 4.2). In this case, the information received by any individual is constant over the whole space, so that it does not have gradients to follow. Only intermediate communication scales supply the receiving individual with an optimal amount of information with which to efficiently locate the callers and the high-quality habitat areas they occupy. The values of the threshold $\kappa$, as long as they fall within a reasonable range, only change the absolute time scales of the searching process.

4.3 Lévy flights

In the case of Lévy flights individual trajectories show sequences of short displacements interspersed with long straight displacements. Alternatively, Lévy flights have been shown as a good searching strategy that may be used by some species. However, empirical studies have generated controversy, since many of the statistical methods used to support the presence of Lévy flights in nature have been questioned, and the issue remains unresolved [Humphries et al., , 2010; Edwards et al., , 2007; Edwards, , 2011]. In this section the case of Lévy searchers is considered. The results will show that neither the behavior of the model, nor the existence of an intermediate optimal communication scale, depend on the characteristics of the motion of the individuals.

Lévy flights do not have a typical length scale and thus the searcher can, in principle, make jumps as larger as the size of the system. The lengths of the jumps, $l>0$, are sorted from a probability distribution with a long tail [Metzler and Klafter, , 2000; Klages et al., , 2008]

[TABLE]

with $l\gg\tilde{l}$, and $0<\mu<2$, where $\tilde{l}$ is a characteristic length scale of the system. This distribution is not defined for $\mu<0$, its mean and variance are unbounded for $0<\mu\leq 1$, and it has a mean but no variance for $1<\mu<2$. Finally, for $\mu\geq 2$, the two first moments exist and thus it obeys the central limit theorem. The Brownian motion limit is recovered in this latter case, while very long jumps are more frequent when $\mu\rightarrow 0$. This extreme is usually referred as the ballistic limit, with a high abundance of straight-line long displacements [Viswanathan et al., , 2011; Méndez et al., , 2014]. The cumulative distribution corresponding to Eq. (4.4) is

[TABLE]

As a simple normalizable cumulative distribution function, with the asymptotic behavior of Eq. (4.5), we will use [Heinsalu et al., , 2010]

[TABLE]

whose probability distribution, $P_{\mu}(l)=\Psi_{\mu}^{\prime}(l)$, is given by

[TABLE]

with $0<\mu<2$, and $b=[\Gamma(1-\mu/2)\Gamma(\mu/2)]/\Gamma(\mu)$. We fix $\tilde{l}=h=1$. In addition, individuals will stop if they find a target during a displacement of length $l$. This naturally introduces a cutoff in the length of the jumps, which becomes more important as target density increases [Viswanathan et al., , 2011]. However, as we will focus on a situation where target density is low, we introduce an exponential cutoff of the order of the system size in the jump length probability distribution to ensure that very long jumps without physical meaning (they imply very high velocities) do not occur

[TABLE]

where $C=\int_{0}^{\infty}\varphi_{\mu}(l)dl$ is the normalization constant, and $\tilde{l}=h$.

Finally for the pairwise interaction we choose a family of functions given by

[TABLE]

where $\sigma^{1/p}$ gives the typical communication scale. For simplicity, and without loss of generality, we will consider only the case $p=2$. Indeed, the choice of the function $V$ is not relevant for the behavior of the model, provided that it defines an interaction length scale through the parameter $\sigma$. This scale must tend to zero in the limit $\sigma\rightarrow 0$ and to infinity in the limit $\sigma\rightarrow\infty$. This assures that the gradient of the calling function vanishes in these limits.

Generally, the search is faster when long displacements occur more frequently. Fig. 4.3 shows search time for different values of $\mu$ and the mean searching time at the optimal communication range as a function of the Lévy exponent (inset). As the frequence of long displacements decreases (increasing $\mu$) the search is slower. Again, the effect of the communication mechanism is more important when we approach the Brownian limit ($\mu\rightarrow 2$), as will be explained in Sec. 4.6.

4.4 One-dimensional analytical approximations

To gain clearer insight and provide analytical arguments, we study a minimalistic version of the model. Consider a single individual in a one-dimensional space of length $L$, so that the highest quality areas are located beyond the limits of the system, i.e. at $x=-1$ and $x=L+1$ (see Fig. 4.4). Note that this would correspond to the ideal situation where all the members of the population but one -the searcher- have already reached one of the targets. A landscape quality function, $g(x)$ must also be defined. Provided it is a smooth, well-behaved function, its particular shape is not relevant. We therefore assume a Gaussian-like quality landscape,

[TABLE]

where $\sigma_{r}$ gives the characteristic width of a high quality region. Notice that $g(x)$ is defined so that highest quality areas are located, as mentioned, at $x=-1$ and $x=L+1$. This ensures that the gradient of the function does not vanish at the extremes of the system (Fig. 4.4), and it is equivalent to setting the value of the threshold $\kappa$ such that the targets start at $x=0$ and $x=L$. We assume that a foraging area is good enough when its quality is higher than $80\%$ of the ideal environment, which means $\kappa=0.8$. As we center the patches of resources at $x=-1$ and $x=L+1$, fixing a good quality treshold at $\kappa=0.8$ is equivalent to fix the width of the environmental quality function at $\sigma_{r}=4.5$, to ensure that $g(0)=g(L)\approx 0.80$. However, the qualitative behavior of the model is independent of this choice.

Finally, for the pairwise communication function we choose the functions defined in Eq. (4.9) with $p=2$ again. The combination of local and nonlocal information gives the total available information for the searcher, $R(x)=B_{g}g(x)+B_{c}S(x)$.

To obtain analytical results, we work in the following on a discrete space. The stochastic particle dynamics equivalent to Eq. (4.1) considers left and right jumping rates which are defined for every individual using the total information function,

[TABLE]

where $\alpha$ is a small positive constant to avoid negative rates that has been given an arbitrary value $(\alpha=10^{-4})$, and $h$ is the spatial discretization ($h=1$). Finally, $\tau_{0}$ is the jumping rate of an individual in absence of information, and it is related to the diffusion component of the dynamics of Eq. (4.1). Given the transition rates of Eq. (4.11), the movement with a higher gain of information has a higher rate, and therefore a larger probability of taking place.

The simplest situation, which allows an analytical treatment of the problem, is to consider only $N=3$ individuals. Two of them are located in the top quality areas just beyond the frontiers of the system limit, $x=-1$ and $x=L+1$, and the other one is still searching for a target. Under these considerations, using the environmental quality function defined in Eq. (4.10), and the pairwise potential of Eq. (4.9), the total available information for the searcher is

[TABLE]

As it was done in Sec. 4.2, the efficiency of the search process is measured in terms of the first arrival time at one of the high quality areas, either at $x=0$ or $x=L$, starting from $x_{0}=L/2$. From the definition of the transition rates in Eq. (4.11), $\tau_{+}(L-1)\gg\tau_{-}(L)$, and equivalently $\tau_{-}(1)\gg\tau_{+}(0)$. This means that at both extremes of the system, the rate at which particles arrive is much higher than the rate at which they leave, so particles do not move when they arrive in the top quality areas. This allows us to consider both extremes $x=0$ and $x=L$ of the system as absorbing, and the first arrival time may be obtained from the flux of presence probability of the searcher there [Redner, , 2001]

[TABLE]

This definition will be used in the following sections to investigate the influence of sharing information (i.e., of the interaction mechanism) on search times. The results will be compared with those obtained using a deterministic approximation of the movement of the searcher. We study two different random strategies: Brownian and Lévy.

4.4.1 Brownian motion

We start studying the Brownian case, where the searcher only jumps -with a given rate- to its nearest neighbors. Therefore the coupling of the set of differential equations describing the occupancy probability of every site of the system is (notice that lattice spacing $h=1$),

[TABLE]

with $i=1,...,L-1$. If the initial position of the particle is known, it is possible to solve Eq. (4.4.1) using the Laplace transform. Once the probability distribution of each point has been obtained, it is possible to obtain the mean first arrival time using Eq. (4.13). The thick line in Fig. 4.5 shows this result, indicating that the searching process is optimal (minimal time to arrive to one of the good quality areas) for intermediate values of $\sigma$. A particularly simple limit in Eq. (4.4.1) appears when $\tau_{+}>>\tau_{-}$ when $x>L/2$ (and the contrary on the other half of the system). The search time is $T(\sigma)=\frac{L}{2\tau_{+}}$. This is the expected result since the movement is mainly in one direction and at a constant rate.

In biological terms this means that the optimal situation for the individuals is to deal with intermediate amounts of information. Extreme situations, where too much ($\sigma\rightarrow\infty$) or too little ($\sigma\rightarrow 0$) information is provided by the population, have the same effect on the mean first arrival time, which tends to the same asymptotic value in both limits. In both cases, the search is driven only by the local perception of the environment.

This calculation gives exact results, but it implies fixing the system size, solving a set of equations of dimension $L$, and finally obtaining the inverse Laplace transform of the solutions. The main disadvantage of this approach is that it is not possible to study the influence of the distance between targets on the optimal communication length. To circumvent this we use a deterministic approach in the continuum limit $h\rightarrow 0$ and define, using the symmetry of the system, a mean drift velocity towards one of the high quality areas, $x=L$,

[TABLE]

Substituting the definition of the transition rates Eq. (4.11), the drift velocity is,

[TABLE]

and therefore the search time is

[TABLE]

We compute the searching time using Eq. (4.17) with the same values of the parameters used before ($\sigma_{r}=4.5$, $B_{g}=1$, and $B_{c}=1$, $L=9$) to compare it with the results given by Eq. (4.13) (inset of Fig. 4.5). The approach in Eq. (4.17) (dotted-dashed line) reproduces the qualitative behavior of the searching time although underestimates the value of the optimal communication range ($\sigma_{opt}=7.2$ while Eq. (4.13) produces $\sigma_{opt}=12.5$). This can be fixed excluding from the average in Eq. (4.15) the boundary of the system introducing a parameter $\epsilon$ in the limits of the integration. To estimate the value of $\epsilon$ it is useful to plot $\tau_{+}(x)-\tau_{-}(x)$ versus $x$ (not shown). The difference between rates, although depending on $\sigma$, starts increasing quickly when $x\geq L-2$, so one can estimate $\epsilon=2$. The inset of Fig. 4.5 shows the exit time as a function of the communication range computed with this approach (dashed line). Its optimal value is in good agreement with the result obtained using the definition of the search time (thick line), with $\sigma_{opt}\approx 12.5$ for both approaches. However the temporal scale of the problem (the absolute values of the times), although higher than with $\epsilon=0$, is still lower in this calculation. Results for $\epsilon=2$ correspond to the dashed line in the inset of Fig. 4.5.

Regardless of the value of $\epsilon$ used in the average, this approximation underestimates the temporal scale of the problem (the absolute values of the times). This is because it is assumed that the searcher follows a deterministic movement to the target neglecting any fluctuation that may slow the process.

Finally, increasing $\sigma$ beyond its optimal value, there is a maximum for the search time for any of the approaches. For these values of the communication range, the nonlocal information at the middle of the system coming from both targets is higher than in the extremes and thus there is a bias to the middle in the movement of the searcher. This small effect, that vanishes when $\sigma$ increases and the information tends to be constant in the whole system, seems to be an artifact of the particular arrangement of the simplified $1D$ system, and does not seem relevant for any real-world consideration of this kind of model. In addition, it does not substantially affect the dynamics because local perception of the environment pushes the individual towards one of the targets.

Finally, within this deterministic approximation, besides studying larger systems with no additional computational cost, it is possible to obtain the optimal value of the interaction range parameter, $\sigma_{opt}$:

[TABLE]

which has to be solved numerically for different sizes of the system. The typical optimal communication scale defined by $\sigma^{1/p}$, (i.e., by $\sigma^{1/2}$ since $p=2$) grows approximately linearly with the distance between targets in the asymptotic limit. Using a regression of the results obtained from the integration of Eq. (4.18) yields an exponent $\sigma_{opt}^{1/2}\propto L^{0.93}$ for $L\gg 1$ (Fig. 4.6).

4.4.2 Lévy flights.

Proceeding as in the case of Brownian motion, considering a Lévy searcher the set of equations for the probability of occupancy is

[TABLE]

with $i=1,\ldots,L-1$.

We assume that if a jump of length in between $j-1$ and $j$ takes place, the individual gets the position $j$. To this aim, the coefficients $A_{j}$ enter in the set of equations (4.4.2) and are defined as $A_{j}=\int_{j-1}^{j}\Psi_{\mu}(l)dl$. They give the probability of a jump of length between $j-1$ and $j$ to happen. The coefficients $B_{j}$ are defined as $B_{j}=\int_{j-1}^{\infty}\Psi_{\mu}(l)dl$, to take into account that the searcher stops if it arrives to a target. This introduces a cutoff in the jumping length distribution Eq. (4.7).

Given the size of the system, $L$, which fixes the dimension of the system of equations (4.4.2), it is possible to obtain an analytical solution for the occupancy probabilities and the mean arrival time to the targets using Eq. (4.13). This is shown in Fig. 4.5, where the Brownian limit is recovered when $\mu\rightarrow 2$. It is also observed that when long jumps are frequent the search is much faster, although the gain in search efficiency due to the communication mechanisms is lower close to the ballistic limit (i.e., $\mu\rightarrow 0$). This will be explained later in Sec. 4.6.

Similarly to the Brownian case, a particularly simple limit in Eq. (4.4.2) appears when $\tau_{+}>>\tau_{-}$ for $x>L/2$ (and the contrary on the other half of the system). The search time is

[TABLE]

where the proportionality constant is a combination of the coefficients $A_{i}$ that depends on the size of the system.

4.5 Continuum approximation

From the Langevin Eq. (4.1), and following the standard arguments presented in Dean, [1996] and Marconi and Tarazona, [1999] it is possible to write an equation for the evolution of the density of individuals, $\rho(\mathbf{r},t)$ 111In Appendix A we show this derivation in detail.. This approach will allow us to fix the parameters of the problem having a better understanding of the role they are playing in the dynamics through a dimensional analysis. However, in the case of the large grazing mammals we are going to study later, it is not very suitable to describe a population as a continuum since the number of individuals is not very high and the typical distances among them is large. Neglecting fluctuations the continuum equation for the density is

[TABLE]

which is quite similar to the one derived in Savel’ev et al., [2005] to study the transport of interacting particles on a substrate.

The same behavior is also shown by the macroscopic Eq. (4.5) (Fig. 4.7). Now $T$ is defined as the time that passes until half of the population has found a target, that is $\int_{g(\mathbf{r})\geq\kappa}\rho(\mathbf{r},t)d\mathbf{r}\geq N/2$. We have integrated the Eq. (4.5) in $1D$ system of length $L=1$, using a single Gaussian patch of resources centered at $L/2$ and periodic boundary conditions for a random initial condition. This is equivalent to the case of an infinite system with equidistant high quality areas. We have taken the calling bias as being much stronger than the resource bias to make the nonlocal mechanism much more important in the search process and thus easier to see how the communication range parameter affects the search time. The differences between the $2D$ individual-based and the $1D$ deterministic density equation description, coupled with the parameter choices (stronger bias in the density equation), explain the different observed time-scales in Figs. 4.2 and 4.7. Due to the simplicity of resources landscape, the stationary distribution of individuals is Dirac delta peaked at the maximum of resources (the center of the system). In Fig. 4.8 we show the density of individuals in the long time limit (not stationary). All the animals are in good habitats, i.e., in areas where the maxima of the $g$ function occur.

4.6 Brownian jumps vs. Lévy flights

As a general result of the model, searching is faster when individuals have intermediate amounts of information, regardless of the kind of movement strategy followed by the population (Brownian or Lévy). However, communication has a larger impact on Brownian motion, i.e., the depth of the well at $\sigma_{opt}$ is larger (Figs. 4.3 and 4.5).

A measure of the improvement in search performance at the optimal communication range is given by the ratio between the search time without communication and that at the optimal communication range, $Q=T_{\sigma\rightarrow 0}/T_{\sigma_{opt}}$. This quantity is plotted in Fig. 4.9 for different Lévy exponents. As previously mentioned, Brownian searchers that are not able to perform long displacements benefit more from communication than Lévy searchers. This is because introducing an additional source of information increases the directionality of the random motion and prevents the searcher from revisiting the same place many times, which is the key problem with Brownian search strategies [Viswanathan et al., , 2011]. A Brownian walker has no directionality in the movement, so provided with sources of information (communication together with the local quality of the landscape) it can search much more efficiently. This effect is less important for Lévy searchers due to the presence of long, straight-line moves that, by themselves, decrease the number of times that a particular area is revisited. In addition the return probability to a given point is much higher in $1D$ than in $2D$. Therefore the directionality introduced by the communication has a stronger effect in the simple $1D$ scenario that we have studied. It is also important to remark that in this case the walker only can move either to the right or to the left at each step. This will make the influence of the bias due to communication much stronger in the jumping probabilities.

In summary, the communication mechanism is less important in Lévy strategies, so that its effect is less noticeable as it is shown in Fig. 4.9 both in $1D$ and $2D$.

However, the value of the optimal interaction range changes with the kind of motion. This is shown in the $2D$ model by the dependence of the mean search time on the communication range for different Lévy exponents (Fig. 4.3). The value of $\sigma_{opt}$ increases with the Lévy exponent, so Brownian searchers $(\mu\rightarrow 2)$ need to spread the information farther (a larger value of $\sigma_{opt}$) than Lévy $(\mu=1)$ walkers to obtain the maximum benefit. In Fig. 4.10 we show the value of the optimal communication range, $\sigma_{opt}$, as a function of the Lévy exponent. Lévy trajectories show clusters of short displacements with frequent turns occasionally broken up by long linear displacements, which account for most of the target encounters. However, because these steps are often much longer than the average distance between targets they are not positively influenced by communication, so any benefit a Lévy strategy gains from communication occurs during the series of short displacements. The time that an individual spends doing short movements is limited by the interarrival time of the large steps, so unless an individual is already relatively close to a target, it will not have time to reach a target before the next big step comes and moves it far away from that original target. Therefore the optimal communication range decreases with decreasing Lévy exponent, $\mu$, as longer displacements become more frequent at lower $\mu$ values.

In addition, the value of $\sigma_{opt}$ depends on both the number of targets and their spatial distribution, as was shown in Sec. 4.4 for a simple $1D$ situation where $\sigma_{opt}\sim L$.

4.7 Summary and conclusions

In this chapter we compared Brownian and Lévy search strategies using a population of individuals that exchange information about the location of spatially distributed targets. Using a simple $1D$ model we have provided analytical results on both cases, concluding that frequent long jumps ($\mu\rightarrow 0$, ballistic limit) minimize the searching times.

However the effect of a communication mechanism is more pronounced in the limit of short jumps i.e., Brownian motion. This means that a population of individuals employing Brownian motion gains proportionally more benefit from communicating and sharing information than does a population of Lévy walkers, where long jumps are more or less frequent depending on the value of the Lévy exponent $\mu$. When messages are exchanged in a range that minimizes search duration, communication is the driving force in the Brownian limit, but occasional long jumps are still responsible for most of the encounters with targets in the case of long-tailed step-length distributions.

The main result of this work is rather general: independently of the kind of communication performed by the population, and of the spatial distribution of the targets, a population of individuals with the ability to communicate will find the targets in a shorter time if the information is spread at intermediate ranges. Both an excess and a lack of information increase the search time. However, the communication mechanism does not have the same quantitative effect on the different moving strategies (i.e., ballistic, Lévy or Brownian). Uninformed Brownian individuals perform a random movement revisiting the same position many times, so having an external source of information introduces directionality on the movement, decreasing the number of times that a point in the space is visited. In the case of Lévy and ballistic strategies $(\mu\rightarrow 0)$, communication is less noticeable because individuals are able to do long jumps. This is already a source of directionality that prevents individuals from revisiting the same points in space many times, and thus weakening the effect of the directionality introduced by communication.

CHAPTER   5    Foraging in * Procapra gutturosa*

In this chapter we show an application of the model presented in Chapter 4 to the particular case of acoustic communication among Mongolian gazelles (Procapra gutturosa) Martínez-García et al., 2013b , for which data are available, searching for good habitat areas. Using Monte Carlo simulations, our results point out that the search is optimal (i.e. the mean first hitting time among searchers is minimum) at intermediate scales of communication. We also present this result in terms of the frequency of the sounds, showing a good agreement with field measurements of the sounds emitted by these gazelles in the wild. The formation of groups in the populations is also studied.

5.1 Introduction

Many living organisms, including bacteria [Liu and Passino, , 2002], insects, and mammals [Zuberbühler et al., , 1997; McComb et al., , 2003] communicate for a variety of reasons including facilitation of social cohesion [Cap et al., , 2008; Pfefferle and Fischer, , 2006], defense against predators [Zuberbühler, , 2001], maintenance of territories [Slater et al., , 1994; Frey et al., , 2007], and to pool information on resource locations when no single individual is sufficiently knowledgeable [von Frisch, , 1967; Hoare et al., , 2004; Berdahl et al., , 2013; Simons, , 2004; Torney et al., , 2009]. Communication among individuals frequently leads to group formation [Eftimie et al., , 2007], which often has clear direct benefits such as reducing individual vulnerability to predators. Such strategies may, however, also have important incidental benefits. For example, an individual that has found a good foraging patch might try to attract conspecifics to reduce its risk of predation, but also provides its conspecifics with information on the location of good forage, thus increasing the foraging efficiency of those responding to the call.

We apply a specialized version of the model introduced in Chapter 4 to the particular case of acoustic communication among Mongolian gazelles, the dominant wild herbivore in the Eastern Steppe of Mongolia (Fig. 5.1a) among 14 species of ungulates. A population of about one milion of animals is estimated, but it is difficult to give a good measurement because of large fluctuations due to the extreme conditions in the steppe that cause periods with important population losses. In addition, the nomadism of the species, travelling long distances during the year, makes more difficult a demographic control. In any case, the species is still recognized as one of Asia’s largest remaining wildlife populations [Olson et al., , 2005], although it has experienced a major reduction in range during the past century, and is further threatened by excessive hunting and continued habitat loss and fragmentation (grassland steppes are increasingly being carved up by fences, roads, agriculture, and densely settled areas while oil fields and pipelines are being developed in the region) (Fig. 5.1 b-d). In fact, Mongolian gazelles were formerly distributed across the whole area of the Republic of Mongolia but the range of this species, between the 1940s and 1960s, was reduced by 70% owing to excessive official hunting and poaching. Nowadays, although individuals or small groups are found across a wider geographical range, higher concentrations of this gazelle species are now limited to the Eastern Steppe (Fig. 5.1a) where they avoid narrow valleys, forest, sand dunes or cultivated fields unless driven there by exceptional circumstances. The plant cover of the dry Steppes of Eastern Mongolia is extremely sparse, generally 5–20% (Fig. 5.1, right panel), rarely reaching 30–40% [Frey and Gebler, , 2003; Mueller et al., , 2008].

In summary, gazelles must find each other to form grazing groups less susceptible to preadator’s attacks (Fig. 5.1), and relatively small areas of good forage in a vast landscape where sound can travel substantial distances [Frey et al., , 2008](Fig. 5.2a) . We aim to explore whether acoustic communication in the steppe could lead to the formation of observed large aggregations of animals (Fig. 5.2b) [Olson et al., , 2009], and how search efficiency depends on the distances over which calls can be perceived. We wonder if the frequency of the voice of the gazelles is optimal to communicate in the steppe, and if the call length-scales that optimize search in real landscapes are biologically and physically plausible. To do this, we couple an individual-based representation of our model with remotely-sensed data on resource quality in the Eastern Steppe.

5.2 The model for acoustic communication

A detailed analysis of gazelle relocation data has shown that, over the spatiotemporal scales relevant to searching for resources (days to weeks), the movement of Mongolian gazelles can be closely approximated by simple Brownian motion for the spatiotemporal scales involved in foraging. Therefore, the starting point of the modeling is

[TABLE]

where $\mathbf{\eta}_{i}(t)$ is a Gaussian white noise term to characterize Brownian motion. The function $\mathbf{\eta}_{i}(t)$ is defined by its statistical properties: its mean value $\langle\eta_{i}(t)\rangle=0$, and the correlation $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2D\delta_{ij}\delta(t-t^{\prime})$. $D$ is the diffusion coefficient. The terms $B_{g}\nabla g(\mathbf{r}_{i})$ and $B_{C}\nabla S(\mathbf{r}_{i})$ are referred to the local and the nonlocal search, as it was defined in Chapter 4. These terms drive the movement to the best grazing areas.

The function $g(\mathbf{r})$ quantifies the habitat quality in the Eastern Steppe of Mongolia through the Normalized Difference Vegetation Index (NDVI), one of the most widely used vegetation quality estimators. It can be calculated from satellite imagery, and has been already applied to gazelle habitat associations in the Mongolian Steppe [Mueller et al., , 2008]. NDVI is characterized by the function $g_{d}({\bf r})$, a continuous function taking values between [math] (no vegetation) and $1$ (fully vegetated). As the vegetation at low NDVI is too sparse, and at high NDVI is too mature and indigestible, gazelles typically seek forage patches characterized by intermediate NDVI values [Mueller et al., , 2008]. To make gradients of resources drive the movement of the individuals to regions with intermediate NDVI values, we apply to the data a linear transformation:

[TABLE]

The new function $g({\bf r})$ defines a resources landscape with values between [math] and $0.5$, where [math] represents both fully vegetated and no vegetation (i.e., low quality forage). We study a subregion with an extension of around $23000$ km2 (Coordinates: $46.50\text{\textdegree}$N, $114.15\text{\textdegree}$E; $46.50\text{\textdegree}$N, $116.80\text{\textdegree}$E; $47.65\text{\textdegree}$N, $114.15\text{\textdegree}$E; $47.65\text{\textdegree}$N, $116.80\text{\textdegree}$E) and assume that the resources remain constant in time during the search. It is crucial now to properly choose the perception function in order to realistically model the case of gazelles performing acoustic communication. It is well known that the sensitivity of the response of the ear does not follow a linear scale, but approximately a logarithmic one. That is why the bel and the decibel are quite suitable to describe the acoustic perception of a listener. Therefore we choose an acoustic perception function of the form

[TABLE]

where $F$ is taken as a logarithmic function. The sound calling of $j$, $V(\mathbf{r},\mathbf{r}_{j})$, plays the role of a two body interaction potential, and $I_{0}$ is the low perception threshold. We take the value of a human ear, $I_{0}=10^{-12}~{}W~{}m^{-2}$, which is similar for most other mammals [Fletcher, , 1992], and in any case, is just a reference value on which our results will not depend. The interaction potential mimicking acoustic communication is

[TABLE]

considering that sound from an acoustic source attenuates in space mainly due to the atmospheric absorption (exponential term), and the spherical spreading of the intensity ($4\pi r^{-2}$ contribution), and neglecting secondary effects [Naguib and Wiley, , 2001]. $P_{0}$ may be understood as the power of the sound at a distance of $1~{}m$ from the source. According to Stoke’s law of sound attenuation [Fletcher, , 1992], the absorption coefficient, $\gamma$, is given by

[TABLE]

where $\eta$ is the viscosity of the air, $\rho$ its density, $v$ the propagation velocity of the acoustic signal (which depends on the temperature and the humidity), and $\nu$ its frequency. We work under environmental conditions of $T=20^{o}C$, and relative humidity of $HR=50\%$, which are quite close to the corresponding empirical values for the summer months from the Baruun-Urt (Mongolia) weather station, averaged over the last 4 years. These values give an absorption coefficient of $\gamma\approx 10^{-10}\nu^{2}~{}m^{-1}$. The inverse of the absorption coefficient, $\gamma^{-1}$, gives the typical length scale for the communication at each frequency, and thus plays the same role as the standard deviation, $\sigma$, did in the Gaussian interaction used in the general model. From its functional dependence, different values of the frequency will modify the value of the absorption coefficient, and consequently, will lead to different communication ranges. Therefore, we will use sound frequency, $\nu$, as the control parameter of the interaction range.

From a statistical analysis of GPS data tracking the positions of $36$ gazelles between 2007 and 2011, we estimate a diffusion constant of $D=74$ km2 day*-1* (Fig. 5.3). To give empirically-based values to the bias parameters, we define a drift velocity, and based on previous field work [Mueller et al., , 2008] we set $v_{drift}=B_{g}\nabla g(\mathbf{r})+B_{c}\nabla S(\mathbf{r})=10~{}{\rm km~{}day}^{-1}$. The local search mechanism is responsible for short-range slow movements, while nonlocal communication gives rise to long and faster movements, and thus we require $B_{g}\nabla g(\mathbf{r})\ll B_{c}\nabla S(\mathbf{r})$.

5.3 Results and discussion

We explore the dependence of this metric on the communication length, $\gamma^{-1}$, or equivalently the frequency, $\nu$ (Fig. 5.4). To this aim, we couple an individual-based model following the dynamics of Eq. (5.1), with a data-based resources landscape sampled every $500$ m (shown in Fig. 5.5), and quantify the efficiency of the search for areas of high quality resources in terms of the mean first arrival time of the population. Similarly to other species, such as lions [Estes, , 1991] or hyenas [Mathevon et al., , 2010], the optimal foraging time ($41$ hours) is obtained for $\gamma^{-1}$ of the order of kilometers (around $6$ km). This result cannot currently be checked with data. However, switching to frequencies, the optimal search is obtained when gazelles communicate at a frequency of $1.25$ kHz, which lies inside the measured interval of frequencies of the sounds emitted by gazelles, $[0.4,2.4]$ kHz [Frey and Gebler, , 2003; Frey et al., , 2008]. This means that the search is optimal when the receiving individual has an intermediate amount of information. A lack of information leads to a slow, inefficient search, while an overabundance of information makes the individual to get lost in the landscape. These different regimes are also observed in the long time spatial distributions (i.e. efficiency of the search in terms of quality) of the Fig. 5.5. For frequencies out of the optimal range, either smaller (Fig. 5.5 top) or larger (Fig. 5.5 bottom right), some animals are still in low-quality areas at the end of the simulation period. At intermediate communication scales, $\nu=1$ kHz, (Fig. 5.5 bottom left) all of the animals end up in regions with the best resources, regardless of where they started from.

In summary, communication over intermediate scales results in faster search, and all the individuals form groups in areas of good resources. While this has obvious advantages in terms of group defense and predator swamping, it will also lead to rapid degradation of the forage (and thus habitat quality) at those locations. This is the problem of foraging influencing the patterns of vegetation, which could be treated in future investigations. Shorter-scale communication implies an almost individual search, which helps preserve local forage quality, but has clear disadvantages in terms of group defense strategies. On the other hand, longer scales lead to the formation of big groups (faster degradation of foraging), and animals need more time to join a group, which has negative consequences against predation. Furthermore, acoustic communication scales significantly larger than the optimal scale for foraging efficiency identified here would be biologically implausible, even if ultimate group size (and not rate of group formation) was the most important aspect of an antipredation strategy.

Our study clarifies some questions on the relationship between communication and optimal search for resources. Our key result is that, in general, intermediate communication distances optimize search efficiency in terms of time and quality. Individuals are able to find the best quality resource patches regardless of where they start from, opening new questions about the distribution of individuals in heterogeneous landscapes. The existence of maximum search efficiency at intermediate communication ranges is robust to the choice of functional form of $V(\mathbf{r})$, allowing the model to be generalized to many different ways of sharing information. Also considering different species on the model (preys and predators or males and females to studying the case of mating) would be interesting extensions of this work. Finally, regarding to the formation of groups because of communication among individuals, exploring tradeoffs between group defense and individual foraging efficiency in highly dynamic landscapes may be a promising avenue for future research.

Appendix A Derivation of the macroscopic Eq. (4.3)

In this appendix we will show the derivation of the macroscopic equation (4.5) in Chapter 4, starting from the Langevin equation for the movement of a singles individual.

Considering a single individual, the Langevin equation is

[TABLE]

where $\eta_{i}(t)$ is a Gaussian white noise with zero mean and correlation delta-correlated in space and time. The available information function, $S_{i}$ is given by

[TABLE]

To obtain the equation for the density of individuals, we will derive a density equation for the case of a single particle and then extend the result to a population with $N$ individuals [Dean, , 1996].

In the case of a single particle the density is given by

[TABLE]

where $\mathbf{X_{i}}(t)$ is the stochastic trajectory of the particle. Then, using this equation A.3 for the density of searchers and the definition of the Dirac delta,

[TABLE]

where $f(\mathbf{r})$ is an arbitrary function. Its time derivative is

[TABLE]

Using the Itô’s formula Eq. (A.5) can be expanded,

[TABLE]

and then, rearranging and integrating by parts each term111It is straightforward just choosing $\nabla f(x)=dV$ in the integration by parts. In the case of the laplacian term, we must integrate twice.,

[TABLE]

On the other hand the time derivative of $f(\mathbf{r})$ can be written as

[TABLE]

As both Eq. (A.7) and (A.8) are true in the case of an arbitrary function, $f(\mathbf{r})$, it is possible to write:

[TABLE]

Finally, we neglect the last stochastic term to get a deterministic approximation as Eq. (4.5).

Appendix B Voronoi diagrams of the model.

The behavior of the model, resulting in optimal searches at intermediate communication ranges, can be explained in terms of Voronoi diagrams Okabe et al., [1992]. Consider every target as a seed that has associated a Voronoi cell formed by those points whose distance to that seed is less than or equal to its distance to any other one (See Figure B.1 (top) for a distribution of the space in $5$ Voronoi cells for an initial distribution of particles with five targets (crosses)). The searching time will be minimized when the information coming from the individuals located on one target covers the full associated Voronoi cell, but only that cell. In this situation, the searchers within that cell will receive information coming only from that target and move towards it. $\sigma_{opt}$ is the communication range that maximizes the gradient (approximately the smallest value of $\sigma$ that makes the calling function not vanishing) of the calling function at the frontiers of the Voronoi cells. Increasing the communication range provides individuals with information coming from different targets, and makes them get overwhelmed in the limit $\sigma\rightarrow\infty$. This Voronoi construction may also help to explain the improvement of the searching strategies because of sharing information. The difference between Brownian and Lévy strategies can be seen in Figure B.1 (Bottom). They show the origin of the individuals that are at each target at the end of a Lévy (Left) and a Brownian search (Right) (i.e., in which Voronoi cell they were at the beginning). In the case of Brownian individuals most of the particles at every target were initially in its Voronoi cell. For Lévy flights the long displacements mix the population in the stationary state (i.e., individuals at a target come from different cells). The communication mechanism is less important in Lévy strategies, so that its effect is less noticeable and the encounters of individuals with targets are caused mainly by the long displacements.

Part IV  **Temporal fluctuations

**

CHAPTER   3    Temporal disorder in up-down symmetric systems

In this chapter we study the effect of temporal fluctuations on systems with up-down symmetry through the behavior of the first-passage times Martínez-García et al., [2012]; Martínez-García, [2011]. This is a relevant question in the modeling of ecosystems, since they are subject to environmental changing conditions. Therefore, it is important to have models that include temporal disorder. We analyze two well-known families of phase transitions in statistical physics —the Ising and the generalized voter universality classes— and scrutinize the consequences of placing them under fluctuating global conditions. It is observed that the variability of the control parameter induces in both classes “temporal Griffiths phases” (TGPs), characterized by broad regions in the parameter space in which the mean first-passage times scale algebraically with system size. In an ecological context, first-passage times are related to typical extinction times, and studying how they are affected by the size of the system (e.g. habitat fragmentation) is a problem of outmost relevance.

3.1 Introduction

Systems with up-down $Z_{2}$ symmetry –including the Ising model– are paradigmatic in mathematical ecology. They allow to address a big variety of problems ranging from species competition [Clifford and Sudbury, , 1973] and neutral theories of biodiversity [Durrett and Levin, , 1996] to allele frequency in genetics [Baxter et al., , 2007]. Some of them, such as the voter model, exhibit absorbing states, a distinctive feature of nonequilibrium dynamics. Once these particular configurations are reached, the system cannot escape from them so they imply the presence of currents [Hinrichsen, , 2000; Ódor, , 2004; Grinstein and Muñoz, , 1996; Marro and Dickman, , 2008].

Phase transitions into absorbing states are quite universal and they depend on few general properties of the system, such as symmetries and its dimensionality, and are insensitive to the underlying microscopic properties. This universality makes possible to establish a classification of the phase transitions into different classes. Those systems exhibiting one absorbing state belong generically to the so called Directed Percolation (DP) universality class and share the same set of critical exponents and scaling functions. However, when there is some additional symmetry or conservation, the phase transition exhibits critical scaling differing from DP. This is the case of systems with two symmetric absorbing states [Hinrichsen, , 2000; Ódor, , 2004; Grinstein and Muñoz, , 1996; Marro and Dickman, , 2008], that show a phase transition usually referred to as Generalized voter (GV). In an ecological context, this is a relevant class of models that can be used in many situations with two equivalent species.

Analytical and numerical studies [Dornic et al., , 2001; Droz et al., , 2003; Al Hammal et al., , 2005; Vazquez and López, , 2008; Russell and Blythe, , 2011] have shown that, depending on some details such GV transition can split into two separate ones: an Ising-like transition in which the up-down symmetry is broken, and a second DP-like transition below which the broken-symmetry phase collapses into the corresponding absorbing state. In particular, a general stochastic theory, aimed at capturing the phenomenology of these systems, was proposed in Al Hammal et al., [2005]; depending on general features they may exhibit a DP, an Ising, or a GV transition.

When they try to mimick ecological systems, these models should not be isolated but, instead, affected by external conditions or by environmental fluctuations. The question of how external variability affects diversity, robustness, and evolution of complex systems, is of outmost relevance in ecology. Take, for instance, the example of the neutral theory of biodiversity: if there are two $Z_{2}$-symmetric (or neutral) species competing, what happens if depending on environmental conditions one of the two species is favored at each time step in a symmetric way? Does such environmental variability enhance species coexistence or does it hinder it? [Giles Leigh Jr., , 1981; Vazquez et al., , 2010; Borgogno et al., , 2009].

Motivated by these questions, we study how basic properties of up-down symmetric systems, such as response functions and first-passage times, are affected by the presence of temporal disorder.

Some previous works have explored from a theoretical point of view the effects of fluctuating global conditions in simple models that exhibit phase transitions [Jensen, , 1996; Alonso and Muñoz, , 2001; Kamenev et al., , 2008]. Recently, a modified version of the simplest representative of the DP class –i.e. the Contact Process– equipped with temporal disorder was studied in Vazquez et al., [2011]. In this model, the control parameter (birth probability) was taken to be a random variable, varying at each time unit. As the control parameter is allowed to take values above and below the transition point of the pure contact process, the system alternates between the tendencies to be active or absorbing. As shown in Vazquez et al., [2011] this dynamical frustration induces a logarithmic type of finite-size scaling at the transition point and generates a subregion in the active phase characterized by a generic algebraic scaling of the extinction times with system size. More strikingly, this subregion is also characterized by generic divergences in the system susceptibility, a property which is reserved for critical points in pure systems. This phenomenology is akin to the one in systems with quenched “spatial” disorder [Vojta, , 2006], which show algebraic relaxation of the order parameter, and singularities in thermodynamic potentials in broad regions of parameter space: the so-called, Griffiths Phases [Bray, , 1987].

In order to investigate whether the anomalous behavior that leads to TGPs around absorbing state (DP) phase transitions is a universal property of systems in other universality classes –and in particular, in up-down symmetric systems– we study the possibility of having TGPs around Ising and GV transitions. We scrutinize simple models in these two classes and assume that the corresponding control parameter changes randomly in time, fluctuating around the transition point of the corresponding pure model, and study mean-first passage times.

3.2 Spatial disorder. Rare regions and the Griffiths Phase.

The presence of noise is an intrinsic property of natural system and it may change its behavior when compared with an ideal situation. Knowing whether and how the critical behavior changes when introducing a small amount of impurities is important in order to apply criticality to real systems. This is the case of the brain, where Griffiths Phases and Lifshitz tails could play a relevant role [Moretti and Muñoz, , 2013].

Will the phase transition remain only at one point in presence of disorder or will the critical point split? If so, will the critical behavior change quantitatively, giving new universality classes with new critical exponents, or even qualitatively with new non-power law scalings at criticality? In this section we will review some of these questions in systems with quenched disorder i.e. depending on the spatial variables but that does not evolve in time [Ódor, , 2004; Vojta, , 2006]. One of the most common ways of introducing quenched disorder in a system is the dilution, that is, the absence of spins in some fixed places of the lattice. The dilution reduces the tendency towards magnetic long-range order in the system. Therefore, the critical value of the control parameter (typically the inverse of the temperature) for the pure model (without noise), $b_{c,pure}$, moves into the ordered phase, $b_{c,q}$ (Fig. 3.1). That is, the transition to an ordered system is at a lower temperature.

On the other hand, in the case of infinite systems, as happens in the thermodynamic limit, it is possible to find regions without vacancies of an arbitrary size, regardless of the impurity concentration, that is the number of vacancies in the whole lattice. When the value of the control parameter is between $b_{c,pure}$ and $b_{c,q}$, although the whole system is in the disordered phase these pure regions can exhibit some local order, which means a local non-vanishing value of the magnetization. These pure spatial regions are called rare regions and the probability of finding them decreases exponentially with its size $V_{RR}$ and the impurity concentration, $p$. In addition, the dynamics in these regions is very slow since a coherent change (fluctuation) is needed in order to flip all the spins therein.

The interval between $b_{c,pure}$ and $b_{c,q}$, in the disordered phase, (Fig. 3.1) is the so-called Griffiths Phase, because it was Robert B. Griffiths the first who showed the possible existence of a singularity in the free energy in this region [Griffiths, , 1969]. Its main characteristics are the generic divergences of the susceptibility, as a consequence of the singularity in the free energy, and an anomalously slow relaxation to zero of the order parameter. Other time-dependent quantities also relax specially slow, mainly as a power law or a stretched exponential, in contrast with the fast decay typical of pure systems, usually exponential.

In a given system with impurity concentration $p$, the probability of finding a rare region decreases exponentially with its d-dimensional volume, $V_{RR}$. Calling $P_{RR}$ the probability of finding a rare region, up to constant factors, it is

[TABLE]

so it is very unlikely to find large rare region. In contrast, the long-time dynamics inside the Griffiths phase is dominated by these regions. Consider for instance the temporal evolution of the order parameter, typically magnetization in magnetic systems, $m$. As its long time behavior is governed by rare regions, it is

[TABLE]

where $\xi_{t}$ is the flipping time of a rare region, and increases exponentially with its size, so it is $\xi_{t}(V_{RR})\propto\exp(bV_{RR})$. The integral in Eq. (3.2) can be solved using a saddle-node integration to obtain the slow relaxation of the magnetization typical of Griffiths Phases

[TABLE]

where $\phi$ is a non-universal exponent.

3.3 Mean-field theory of $Z_{2}$-symmetric models with temporal disorder

Interacting particle models, such as ecosystems, evolve stochastically over time. A useful technique to study such systems is the mean-field (MF) approach, which implicitly assumes a well-mixed situation, where each particle can interact with any other, providing a sound approximation in high dimensional systems. One way in which the mean-field limit can be seen at work is by analyzing a fully connected network (FCN), where each node (particle) is directly connected to any other else, mimicking an infinite dimensional system.

In the models that we study here, states can be labeled with occupation-number variables $\rho_{i}$ taking a value $1$ if node $i$ is occupied or [math] if it is empty, or alternatively by spin variables $S_{i}=2\rho_{i}-1$, with $S_{i}=\pm 1$. Using these latter, the natural order parameter is the magnetization per spin, defined as

[TABLE]

where $N$ is the total number of particles in the system. Its value varies between $-1$ an $1$, both of them corresponding to a fully ordered state, with the intermediate value $m=0$ reflecting a disordered system with the same amount of spins up and down. The rate equation for the probability $P(m,t)$ of having magnetization $m$ at a given time $t$, is

[TABLE]

where $\omega_{\pm}(m,b)$ are the transition probabilities from a state with magnetization $m$ to a state with magnetization $m\pm 2/N$.

In the limit of $N\rightarrow\infty$ ($dt\rightarrow 0$) Eq. (3.3) becomes a Master equation (Sec. 1.1.1),

[TABLE]

This describes a process in which a “spin” is randomly selected at every time-step (of length $dt=1/N$), and inverted with a probability that depends on $m$ and the control parameter $b$. The allowed magnetization changes in an individual update, $\Delta m=\pm 2/N$, are infinitesimally small in the $N\rightarrow\infty$ limit. In this limit one can perform a Taylor expansion of the rates around $m$, leading to the Fokker-Planck equation

[TABLE]

with drift and diffusion terms given, respectively, by

[TABLE]

From Eq. (3.7), and working in the Itô scheme (as justified by the fact that it comes from a discrete in time equation [Horsthemke and Lefever, , 1984]), its equivalent Langevin equation is [Gardiner, , 1985]

[TABLE]

where the dot stands for time derivative, and $\eta(t)$ is a Gaussian white noise of zero-mean and correlations $\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$. The diffusion term is proportional to $1/\sqrt{N}$, and therefore, it vanishes in the thermodynamic limit ($N\to\infty$), leading to a deterministic equation for $m$ (Sec. 1.1.1).

The drift and diffusion coefficients in Eq. (3.10) depend not only on the magnetization, but also on the parameter $b$. To analyze the behavior of the system when $b$ changes randomly over time, and following previous works [Vazquez et al., , 2010, 2011], we allow $b$ to take a new random value, extracted from a uniform distribution, in the interval $(b_{0}-\sigma,b_{0}+\sigma)$ at each MC step, i.e. every time interval $\tau=1$. Thus, we assume that the dynamics of $b(t)$ obeys an Ornstein-Uhlenbeck process

[TABLE]

where $\xi(t)$ is a step-like function that randomly fluctuates between $-1$ and $1$, as depicted in Fig. 3.3a. Its average correlation is

[TABLE]

where the bar stands for time averaging. The parameters $b_{0}$ and $\sigma$ are chosen with the requirement that $b$ takes values at both sides of the transition point of the pure model (see Fig. 3.3b), that is, the model with constant $b$. Thus, the system randomly shifts between the tendencies to be in one phase or the other (see Fig. 3.4).

The model presents both intrinsic and extrinsic fluctuations, as represented by the white noise $\eta(t)$ and the colored noise $\xi(t)$, respectively. Plugging Eq. (3.11) for $b(t)$ into Eq. (3.10), and retaining only linear terms in the noise one readily obtains

[TABLE]

where $f_{0}(m)\equiv f(m,b_{0})$, $g_{0}(m)\equiv g(m,b_{0})$ and $j_{0}(m)$ is a function determined by the functional form of $f(m,b)$, that might also depend on $b_{0}$. To simplify the analysis, we assume that relaxation times are much longer than the autocorrelation time $\tau$, and thus take the limit $\tau\to 0$ in the correlation function Eq. (3.15), and transform the external colored noise $\xi$ into a Gaussian white noise with effective amplitude $K\equiv\int_{-\infty}^{+\infty}\overline{<\xi(t)\xi(t+\Delta t)>}\,d\Delta t=\tau/3$. Then, we combine the two white noises into an effective Gaussian white noise, whose square amplitude is the sum of the squared amplitudes of both noises [Gardiner, , 1985], and finally arrive at

[TABLE]

where $\langle\gamma(t)\rangle=0$ and $\langle\gamma(t)\gamma(t^{\prime})\rangle=\delta(t-t^{\prime})$.

In the next two sections we analyze the dynamics of the kinetic Ising model with Glauber dynamics and a variation of the voter model (the, so-called, q-voter model) –which are representative of the Ising and GV transitions respectively– in the presence of external noise. For that we follow the strategy developed in this section to derive mean-field Langevin equations and present also results of numerical simulations (for both finite and infinite dimensional systems), as well as analytical calculations.

3.4 Ising transition with temporal disorder

We consider the kinetic Ising model with Glauber dynamics [Glauber, , 1963], as defined by the following transition rates

[TABLE]

The sum extends over the $2d$ nearest neighbors of a given spin $i$ on a $d$-dimensional hypercubic lattice, and $b=J\beta$ is the control parameter. $J$ is the coupling constant between spins, which we set to $1$ from now on, and $\beta=(k_{B}T)^{-1}$. Note that $b$ in this case is proportional to the inverse temperature.

3.4.1 The Langevin equation

In the mean-field case, the cubic lattice is replaced by a fully-connected network in which the number of neighbors $2d$ of a given site is simply $N-1$. Then, the transition rates of Eq. (3.18) can be expressed as

[TABLE]

which implies $\omega_{\pm}(m,b)=\frac{1\mp m}{2}\,\Omega_{\pm}(m,b)$ for jumps in the magnetization. Following the steps in the previous section, and expanding $\Omega_{\pm}$ to third order in $m$, we obtain

[TABLE]

where $b_{0}$ is the mean value of the stochastic control parameter, $a_{0}\equiv b_{0}-1$, and $c_{0}\equiv b_{0}^{3}/3$.

The potential $V(m)=-\frac{a_{0}}{2}m^{2}+\frac{c_{0}}{4}m^{4}$ associated with the deterministic term of Eq. (3.20) has the standard shape of the Ising class, that is, of systems exhibiting a spontaneous breaking of the $Z_{2}$ symmetry (Fig. 3.5). A single minimum at $m=0$ exists in the disordered phase, while two symmetric ones, at $\pm\sqrt{a/c}$ exist below the critical point.

3.4.2 Numerical results

In this section we study the behavior of the mean crossing time, that has been shown to be relevant in the presence of temporal disorder [Vazquez et al., , 2011]. The crossing time is the time employed by the system to reach the disordered zero-magnetization state for the first time, starting from a fully ordered state with $|m|=1$ (see Fig. 3.6). Crossing times were calculated by numerically integrating Eq. (3.20) for different realizations of the noise $\gamma$ and averaging over many independent realizations. These integrations were performed using a standard stochastic Runge-Kutta scheme (note that, the noise term does not have any pathological behavior at $m=0$ as occurs in systems with absorbing states, for which more refined integration techniques are required [Dornic et al., , 2005]) . Results are shown in Fig. 3.7.

To estimate the critical point, we calculated the time evolution of the average magnetization $\langle m\rangle(t)$ by integrating the Langevin equation (3.20), and also by performing Monte Carlo simulations of the particle system on a fully connected network. At the critical point $b_{0,c}$ the magnetization decays to zero as $\langle m\rangle\sim t^{-\beta}$. We have estimated $b_{0,c}=1$, which coincides with the pure case critical point $b_{c,pure}=1$: the critical point in the presence of disorder in mean-field is not shifted with respect to the pure system, in agreement with the analytical calculation in appendix B. At this critical point, as it is characteristic of TGPs [Vazquez et al., , 2011], a scaling of the form $T\sim[\ln N]^{\alpha}$ is expected. The numerically determined exponent value $\alpha\simeq 2.81$ for $\sigma=0.4$ is higher than the exponent $\alpha=2$ of the asymptotic analytical prediction Eq. (B.53), probably because of the asymptotic regime in $\ln N$ has not been reached. Instead, the behavior for arbitrary values of $N$ appears to be a second order polynomial in $\ln N$, as we can see in Eq. (B.44). Indeed, the numerical data is well fitted by the quadratic function $a\,(\ln N)^{2}+b\,\ln N+c$ (see inset of Fig. 3.7).This is to be compared with the standard power-law scaling $T\sim N^{\beta}$ characteristic of pure systems, i.e. for $\sigma=0$. Moreover, a broad region showing algebraic scaling $T\sim N^{\delta}$ with a continuously varying exponent $\delta(b_{0})$ ($\delta\rightarrow 0$ as $b_{0}\rightarrow b^{+}_{0,c}$) appears in the ordered phase $b_{0}>b_{0,c}$. Both $\alpha$ and $\delta$ are not universal and depend on the noise strength $\sigma$. Finally, in the disordered phase the scaling of $T$ is observed to be logarithmic, $T\sim\ln N$.

We have also performed Monte Carlo simulations of the time-disordered Glauber model on two- and three-dimensional cubic lattices with nearest neighbor interactions. The critical point was computed following standard methods, that is, by looking for a power law scaling of $\langle m\rangle$ versus time, as we mentioned above. In $d=2$, a shift in the critical point was found: from $b_{c,pure}=0.441(1)$ in the pure model to $b_{0,c}=0.605(1)$ for $\sigma=0.4$. However, the scaling behavior of $T$ with $N$ resembles that of the pure model, with $T\sim N^{\beta}$ at criticality (with an exponent numerically close to that of the pure model [Marro and Dickman, , 2008]), and an exponential growth $T\sim\exp(cN)$, where $c$ is a positive constant, in the ordered phase (Arrhenius law) ) (see Fig. 3.8 (Left)). Thus, no region of generic algebraic scaling appears in this low-dimensional system. On the contrary, in $d=3$, results qualitatively similar to mean-field ones are recovered (see Fig. 3.8 (Right)). The critical point is shifted from $b_{c,pure}=0.222(1)$ (calculated in Heuer, [1993]) to $b_{0,c}=0.413(2)$, with a critical exponent $\alpha(d=3)=5.29$ for $\sigma=0.4$, and generic algebraic scaling in the ordered phase. In conclusion, our numerical studies suggest that the lower critical dimension for the TGPs in the Ising transition is $d_{c}=3$. This is in agreement with the analytical finding in Alonso and Muñoz, [2001], establishing that temporal disorder is irrelevant in Ising-like systems below three dimensions. This result is to be compared with $d_{c}=2$ numerically reported for the existence of TGPs in DP-like transitions [Vazquez et al., , 2011] (observe, however, that temporal disorder, in this case, affects the value of critical exponents at criticality in all spatial dimensions). Further studies are needed to clarify the relation between disorder-relevance at criticality and the existence or not of TGPs.

3.4.3 Analytical results

Let us consider the Langevin equation (3.17) in the thermodynamic limit ($g_{0}(m)=0$). Given that the remaining intrinsic noise comes from a transformation of a colored noise into a white noise, the Stratonovich interpretation is to be used to obtain its associated Fokker-Planck equation (see e.g. [Horsthemke and Lefever, , 1984])

[TABLE]

Starting from this N-independent Fokker-Planck equation (3.4.3), we will next provide analytical results on the mean crossing times. An effective dependence on $N$ is implemented by calculating the first-passage time to the state $m=|2/N|$ rather than $m=0$. This is equivalent to the assumption that the system reaches the zero magnetization state with an equal number $N_{+}=N_{-}=N/2$ of up and down spins when $|m|<2/N$, that is, when $N/2-1<N_{+}<N/2+1$. The mean-first passage time $T$ associated with the Fokker-Planck equation (3.4.3) obeys the differential equation (Chapter 1) [Redner, , 2001]

[TABLE]

with absorbing and reflecting boundaries at $|m|=2/N$ and $|m|=1$, respectively. The solution, starting at time $t=0$ from $m=1$ is given by

[TABLE]

where

[TABLE]

Computing these integrals (see Appendix B) we obtain

[TABLE]

These expressions qualitatively agree with the numerical results of Fig. 3.7, showing that $T$ grows logarithmically with $N$ in the absorbing phase $b_{0}<1$, as a power law in the active phase $b_{0}>1$, and as a power of $\ln N$ (i.e. poly-logarithmically) at the transition point $b_{0,c}=1$. The exponents $\delta=6(b_{0}-1)/\sigma^{2}$ do not agree well with the numerically determined exponents. This is due to the fact that we have neglected the $1/\sqrt{N}$ term by taking $g_{0}=0$, which becomes of the same magnitude as the $j_{0}$ term when $|m|$ approaches $2/N$. This was confirmed by testing that analytical expressions Eq. (3.28) agree very well with numerical integrations of Eq. (3.20) performed for $g_{0}=0$, and setting the crossing point at $m=2/N$ (See Fig. 3.9). In summary, this analytical approach reproduces qualitatively –and in some cases quantitatively– the above reported non-trivial phenomenology.

3.5 Generalized Voter transition with temporal disorder

We study in this section the GV transition [Dornic et al., , 2001], which appears when a $Z_{2}$-symmetry system simultaneously breaks the symmetry and reaches one of the two absorbing states. A model presenting this type of transition is the nonlinear $q$-voter model, introduced in Castellano et al., [2009]. The microscopic dynamics of this nonlinear version of the voter model consists in randomly picking a spin $S_{i}$ and flipping it with a probability that depends on the state of $q$ randomly chosen neighbors of $S_{i}$ (with possible repetitions). If all neighbors are at the same state, then $S_{i}$ adopts it with probability $1$ (which implies, in particular, that the two completely ordered configurations are absorbing). Otherwise, $S_{i}$ flips with a state-dependent probability

[TABLE]

where $x$ is the fraction of disagreeing (antiparallel) neighbors and $b$ is a control parameter. Three types of transitions, Ising, DP and GV can be observed in this model depending on the value of $q$ [Castellano et al., , 2009]. Here, we focus on the $q=3$ case, for which a unique GV transition at $b_{c}=1/3$ has been reported [Castellano et al., , 2009].

3.5.1 The Langevin equation

In the MF limit (FCN) 111For a fully connected network the number of neighbors has no meaning. However the MF limit of the model refers to the use of the probability given by Eq.(3.29), the fractions of antiparallel neighbors of the two types of spins $S_{i}=1$ and $S_{i}=-1$ are $x=(1-m)/2$ and $x=(1+m)/2$, respectively. Thus, the transition probabilities are

[TABLE]

Following the same steps as in the previous section, we obtain the Langevin equation

[TABLE]

Let us remark that the potential in the nonlinear voter model (Fig. 3.10) differs from that for the Ising model. Owing to the fact that the coefficients of the linear and cubic term in the deterministic part of Eq. (3.5.1) coincide (except for their sign), the system exhibits a discontinuous jump at the transition point, where the potential minimum changes directly from $m=0$ in the disordered phase to $m=\pm 1$ in the ordered one. Furthermore, the potential vanishes at the critical point [Al Hammal et al., , 2005].

3.5.2 Numerical Results

The ordering time, defined as the averaged time required to reach a completely ordered configuration (absorbing state) starting from a disordered configuration, is the equivalent of the crossing time above. We have measured the mean ordering time $T$ by both, integrating the Langevin equation (3.5.1) and running Monte Carlo simulations of the microscopic dynamics on FCNs and finite dimensions. In Fig. 3.11 we show the MF results. We observe that $T$ has a similar behavior to the one found for the mean crossing time in the Ising model, and for the mean extinction time for the contact process [Vazquez et al., , 2011]. That is, a critical scaling $T\sim[\ln N]^{\alpha}$ at the transition point $b_{0,c}=1/3$, with a critical exponent $\alpha=3.68$ for $\sigma=0.3$, a logarithmic scaling $T\sim\ln N$ in the absorbing phase $b_{0}<b_{0,c}$, and a power law scaling $T\sim N^{\delta}$ with continuously varying exponent $\delta(b_{0})$ in the active phase $b_{0}>b_{0,c}$.

Monte Carlo simulations on regular lattices of dimensions $d=2$ and $d=3$ revealed that there is no significant change in the scaling behavior respect to the pure model (not shown). The critical point shifts in $d=2$ and remains very close to its mean-field value in $d=3$, but results are compatible with the usual critical (pure) voter scaling $T_{2d}\sim N\ln N$ and $T_{3d}\sim N$. In the absorbing phase $T$ grows logarithmically with $N$, while in the active phase $T$ grows exponentially fast with $N$, as in the pure-model case. Therefore, in these finite dimensional systems we do not find any TGP nor other anomalous effects induced by temporal disorder, although we cannot numerically exclude their existence in $d=3$. Such effects should be observable only in higher dimensional systems (closer to the mean-field limit).

3.5.3 Analytical results

The ordering time $T$ can be estimated by assuming that the dynamics is described by the Langevin equation (3.5.1), and calculating the mean first-passage time from $m=0$ to any of the two barriers located at $|m|=1$. It turns out useful to consider the density of up spins rather than the magnetization

[TABLE]

$T$ is the mean first-passage time to $\rho=0$ starting from $\rho=1/2$. The Langevin equation for $\rho$ is obtained from Eq. (3.5.1), by neglecting the $1/\sqrt{N}$ term and applying the ordinary transformation of variables (which is done employing standard algebra, given that Eq. (3.5.1) is interpreted in the Stratonovich sense) is

[TABLE]

with

[TABLE]

where $a_{0}=1-3b_{0}$.

Now, we can follow the same steps as in Sec. 3.4.3 for the Ising model, and find the equation for the mean first-passage time $T(\rho)$ by means of the Fokker-Planck equation. The solution is given by (see Appendix B)

[TABLE]

These scalings, which qualitatively agree with the numerical results of Fig. 3.11 for the $q$-voter, show that the behavior of $T$ is analogous to the one observed in the Ising transition of section 3.4 and in the DP transition found in Vazquez et al., [2011]. Therefore, we conclude that TGPs appear around GV transitions in the presence of external varying parameters in high dimensional systems.

For the GV universality class the renormalization group fixed point is a non-perturbative one [Canet et al., , 2005], becoming relevant in a dimension between one and two. A field theoretical implementation of temporal disorder in this theory is still missing, hence, theoretical predictions and sound criteria for disorder relevance are not available.

3.6 Summary and conclusions

We have investigated the effect of temporal disorder on phase transitions exhibited by $Z_{2}$ symmetric systems: the (continuous) Ising and (discontinuous) GV transitions which appear in many different scenarios. We have explored whether temporal disorder induces Temporal Griffiths Phases as it was previously found in standard (DP) systems with one absorbing state. By performing mean-field analyses as well as extensive computer simulations (in both fully connected networks and in finite dimensional lattices) we found that TGPs can exist around equilibrium (Ising) transitions (above $d=2$) and around discontinuous (GV) non-equilibrium transitions (only in high-dimensional systems).

Therefore, we confirm that TGPs may also appear in systems with two symmetric absorbing states, illustrating the generality of the underlying mechanism: the appearance of a region, induced by temporal stochasticity of the control parameter, where first-passage times scale as power laws of the system size. The algebraic scaling of the crossing time, compared with the exponential one observed in pure system, shows that temporal disorder makes the ordered/active phase less stable. This implies that the system becomes highly susceptible to perturbations. This appears to be a rather general and robust phenomenon and an relevant result with applications in ecology. At a given system size, considering the effect of the environmental variablity considerably reduces the crossing times, that in some cases can be identified with population lifetimes. This is a consequence of the change in the growth of this magnitude with the system size, from exponential to algebraic, when external fluctuations are introduced. From a conservational point of view, the power-law scaling of the crossing times implies that population lifetimes would decrease less, compared with the exponential behavior of the pure model, when going from larger to smaller systems, which is the main consequence of habitat fragmentation.

Additionally, although it is of secondary interest to our focus on ecological systems and it has been not shown in this chapter, we have also confirmed that the response function of the system (i.e. susceptibility) diverges in a finite region close to the critical point. This complementary result has been also obtained in Vazquez et al., [2011] and confirms the phenomelogical similarities between Temporal Griffiths Phases and Griffiths Phases in systems with quenched disorder.

It also seems to be a general property that TGPs do not appear in low dimensional systems, where standard fluctuations dominate over temporal disorder. In all the cases studied so far, a critical dimension $d_{c}$ –at and below which TGPs do not appear– exist ($d_{c}=1$ for DP transitions, $d_{c}=2$ for Ising like systems, and $d_{c}\simeq 3$ for GV ones). Calculating analytically such a critical dimension and comparing it with the standard critical dimension for the relevance/irrelevance of temporal disorder at the critical point (i.e. at the renormalization group non-trivial fixed point of the corresponding field theory) remains an open and challenging task.

Future research might be oriented to the effect of temporal disorder on the formation and dynamics of spatial structures.

Appendix A Itô-Stratonovich discussion.

The integration of stochastic differential equations with multiplicative white noise presents some problems because the integral of the noise is not well defined. These problems are solved choosing either the Itô or the Stratonovich definition of the integral.111Any definition can be chosen or even made, but these two are the most often used. We have chosen one or the other depending on the origin of the noise term. This Appendix explains how Itô calculus works, and the connection between Itô and Stratonovich schemes. We will finish discussing what of them is more suitable in each situation.

A.1 Stochastic integration.

Let us start providing a precise definition of the second integral in

[TABLE]

that is

[TABLE]

The integration interval $[0,t]$ is divided into n subintervals,

[TABLE]

and intermediate points in each interval $\tau_{i}$ defined

[TABLE]

The stochastic integral in Eq. (A.2) is defined as the limit of the partial sums,

[TABLE]

where the Itô vs Stratonovich dilemma resides in the fact that the limit of $S_{n}$ depends on the particular set of points $\tau_{i}$ that are used. Itô stochastic integral is defined taking $\alpha=0$, so Eq. (A.5) becomes

[TABLE]

that is, the known function $g(x(t))$ is evaluted on the beginning point of the interval while Stratonovich is obtained if $\alpha=1/2$ and

[TABLE]

A.2 Itô’s formula.

In spite of being much more elegant from a mathematical point of view, Ito’s prescription is not always the most suitable choice for physical interpretation. Calculus we are used to does not work in this scheme, and a different change of variables must be considered. Let us take an arbitrary function $a[x(t)]$ with $x(t)$ obeying the SDE

[TABLE]

where $\xi(x,t)$ is a white Gaussian noise. Consider

[TABLE]

where higher terms in $dt$ have been neglected. Now, replacing $dW^{2}(t)=dt$ (see [Gardiner, , 1985] for a proof),

[TABLE]

which is known as the Itô’s formula and shows that change of variables is not given by ordinary calculus unless $a[x(t)]$ is linear in $x(t)$.

A.3 From Stratonovich to Itô.

As may be expected, both interpretations of the stochastic integral are somehow related. To show it, consider an stochastic differential equation

[TABLE]

where $\eta(t)$ is a white, zero mean, Gaussian noise. Integrating, it is,

[TABLE]

where $S$ denotes that a Stratonovich integration is used. We will derive the equivalent Itô stochastic differential equation.

Assuming that $x(t)$ is a solution of

[TABLE]

the corresponding $\alpha[x(t),t]$ and $\beta[x(t),t]$ will be deduced. The first step is to compute the connection between $S\int_{0}^{t}\beta[x(s),s]dW(s)$ and $\int_{0}^{t}\beta[x(s),s]dW(s)$, where the lack of notation in the second integral means an Itô interpretation. Then,

[TABLE]

Taking into an account

[TABLE]

in the Stratonovich integral, then

[TABLE]

Now, the Itô SDE (A.13) is used in order to write

[TABLE]

Using Itô’s formula given by Eq. (A.10) as well as simplifying the notation writing $\beta(t_{i-1})$ instead of $\beta[x(t_{i-1}),t_{i-1}]$, Eq. (A.16) becomes,

[TABLE]

Finally, substituing into Eq. (A.14), neglecting terms in $dt^{2}$ and $dWdt$ and setting $dW^{2}=dt$,

[TABLE]

or going back to integrals,

[TABLE]

which means that the stochastic integral in Stratonovich representation is equivalent to a stochastic integral in Itô’s and a drift term. It is also important to remark that this formula gives a connection between both integrals of function $\beta[x(s),s]$, in which $x(s)$ is the solution of the Itô SDE (A.13). It does not give a general connection between the Itô and Stratonovich integrals of arbitrary functions.

[TABLE]

Or

[TABLE]

There are many consequences of this transformation formula, but the more important are

• •

It is always possible to change from the Stratonovich to the Itô interpretation of a SDE by adding $\frac{1}{2}\beta(x,t)\partial_{x}\beta(x,t)$ or in the inverse direction subtracting a similar term.

• •

In the case of additive noise, i.e., $g(x,t)=const.$ in Eq. (A.8) there is no difference between the Itô and Stratonovich integral.

• •

In the case of multiplicative noise, i.e., $g(x,t)\neq const.$ in Eq. (A.8), where the influence of the random force depends on the state of the process, the correlation between both the random force and the state of the process is implicit in the Stratonovich integral. It gives raise to the noise induced drift when moving to Itô appearing in the deterministic part of the equation.

• •

The Stratonovich calculus obeys the classical chain rule, Itô’s formula derived in Section A.2 plays a similar role on Itô’s calculus.

A.4 Stratonovich / Itô dilemma.

The long controversy in the physical literature about what is the right definition of the stochastic integral has created some confusion on this topic. That’s why, although a much more mathematically rigorous and longer discussion can be found in the references, [van Kampen, , 2007; Horsthemke and Lefever, , 1984; Jazwinski, , 2007] some hand waving arguments will be given in this section.

First of all it is important to say that this kind of ambiguity when working with SDE only yields for the particular, but most common, case of differential equations with multiplicative white noise222Cases where the rapidly fluctuating external force depends on the state of the system.. As a first approach, it is natural to tend to believe that due to invariance of the equations under “coordinate transformation” $y=u(x)$ when working on Stratonovich scheme it is the proper choice. However, it means nothing but it obeys the classical calculus rules we are familiar with. The only quantities that have to be invariant under a transformation $u=y(x)$, where $u$ is one to one, are the probabilities,

[TABLE]

and this is of course guaranteed in both calculi. They lead to a consistent calculus.

It looks sensible, then, to change the question. The matter is not what is the right definition of the stochastic integral, but how do we model real systems by stochastic processes. That is, in which situation either Itô’s or Stratonovich’s choice is the most suitable.

On the one hand, if the starting point is a phenomenological equation in which some fluctuating parameters represented through colored noise terms are approximated by Gaussian white noise, then the most appropiate process is the one that is defined by the Stratonovich interpretation of the equation.

On the other hand, in many systems the appropiate starting point is a discret time equation, as it happens, for instance, in biology when working with populations of insects. In these cases the equation reads

[TABLE]

where $t_{i}=t_{i-1}+\Delta t$ in every time step and $Q_{i}$ are Gaussian independent random variables with expected values $<Q(t_{i})>=0$ and $<Q^{2}(t_{i})>=\Delta t$.

If times considered are longer compared to $\Delta t$, the continuous time limit can be taken. Then the system is described by

[TABLE]

which is also a SDE where $W(t)$ is the Wiener process. However, due to the asymmetric form of Eq. (A.24) with respect to time it is much more appropiate the stochastic process defined according to the Itô interpretation in this case.

To sum up, as a take to home message from this section, two different cases can be considered when working with SDE. When the white Gaussian noise limit is considered as the limit of a colored noise when the correlation time tends to zero, the Stratonovich interpretation is more sensible, when Itô’s is more suitable when it represents the continuous limit of a discrete time problem. In any case, there are no universally valid theoretical reasons why one or the other interpretation of an SDE should be preferred and the ultimate test must be the confrontation of the analytical (or numerical) results with the experimental facts.

Appendix B Analytical calculations on the escape time for the Ising Model

We will show here all the analytical calculations done to obtain the result of Eq. (3.28). To make the integrals analytically solvable we take $g_{1,0}(m)=0$, so the Langevin equation is

[TABLE]

with

[TABLE]

where $a_{0}=b_{0}-1$, $c_{0}=b_{0}^{3}/3$ and $\gamma(t)$ is, again, a white Gaussian noise defined by its autocorrelation function $<\gamma(t)\gamma(t^{\prime})>=\delta(t-t^{\prime})$, and its mean value $<\gamma(t)>=0$.

The Langevin equation (B.1) presents one absorbing state in $m=0$ induced by the simplification done when neglecting thermal fluctuations.

Working in the Stratonovich scheme 111Because the noise term comes from taking the white noise limit in a colored one, the associated Fokker-Planck equation is

[TABLE]

where

[TABLE]

According to Gardiner, [1985] and Redner, [2001], the escape time from an starting point $m$ obeys,

[TABLE]

As the size of the system does not appear naturally in the problem because of the simplification done when taking $g_{1,0}=0$, the mean escape time will be defined as that needed to pass through $m=2/N$, which is the length of the jumps of the Brownian particle to whose movement the problem has been mapped. Then, taking into an account that there is an absorbing barrier in $m=0$ and a reflecting one in $m=1$ and the initial condition, the solution is [Gardiner, , 1985]

[TABLE]

with

[TABLE]

which involves $6^{th}$ and $4^{th}$ order polynomial functions.

To make the integral simpler, we expand the functions up to $3^{rd}$ order, and takre the low integration limit in Eq. (B.7) at $1$ instead of $2/N$. This change can be done because $\psi(z)$ appears both in the numerator and the denominator of $T(m)$, so the contribution of the lower limit vanishes, allowing to take it in our interest. The first assumption leads to

[TABLE]

where it has been defined $\omega\equiv\tau\sigma^{2}/3$; $r\equiv a_{0}+\omega/2$; $s\equiv(c_{0}+2\omega b_{0}^{2})$. The size of the system will be rescaled too, so the lower limit in the expression of the escape time Eq. (B.6) is $1/N$. This simplifies the notation and does not affect the qualitative behaviour of the results in the asymptotic limit (only a constant factor appears).

Now, it can be written,

[TABLE]

where $\alpha\equiv 2r/\omega$ and $\beta\equiv s/\omega$.

Lets now define the function

[TABLE]

which presents a singularity when $\alpha=1$ as can be seen integrating by parts. With the definition made of the parameters, it can be shown that it corresponds to $b_{0}=1\equiv b_{0,c}$.

Considering the definitions of Eq. (B.9) and Eq. (B.10), the mean escape time is given by

[TABLE]

Each case will be studied separately.

B.1 Case $\alpha\neq 1$

Integrating by parts Eq. (B.10)

[TABLE]

where the new integral can be solved again integrating by parts. Working recursively this way,

[TABLE]

or

[TABLE]

The mean escape time is given now by

[TABLE]

where

[TABLE]

Integrating $I_{1}(N)$ by parts (taking again the exponential part as $u$ and the rest as $dv$) and following the same procedure as in Eq. (B.10) it is obtained

[TABLE]

while $I_{2}(k,N)$ is easily solved

[TABLE]

At the end, an expresion for the mean escape time is achieved

[TABLE]

In the asymptotic limit $N\rightarrow\infty$ two different cases must be considered.

B.1.1 $\alpha<1$

Under this prescription, $\alpha-1-2l<0$ when $l\geq 0$ so in $I_{1}(N)$

[TABLE]

which leads to

[TABLE]

Finally, for the mean escape time,

[TABLE]

what means,

[TABLE]

B.1.2 $\alpha>1$

It is taken as a starting point

[TABLE]

where considering that $N^{\alpha-1}\gg N^{\alpha-1-2l},\forall l>0$, only the first term in Eq. (B.27) is relevant. It implies

[TABLE]

and in the mean escape time

[TABLE]

B.2 Case $\alpha=1$ Critical point.

It has to be solved now

[TABLE]

using the expansion of the exponential function and integrating it is

[TABLE]

It makes the mean escape time to obey, taking the form of $I(y)$ Eq. (B.31) into Eq. (B.11)

[TABLE]

where

[TABLE]

First of all, lets consider the solution of $I_{3}(N)$ integrating by parts ($u={\rm e}^{\beta y^{2}}$ and $dv={\rm ln}(y)y^{-1}dy$) so,

[TABLE]

where the new integral is solved again integrating by parts taking

[TABLE]

It leads to a solution behaving like

[TABLE]

so finally,

[TABLE]

which scales in the asymptotic limit as

[TABLE]

Secondly, lets focus on $I_{4}(N)$, where, again, an expansion of the exponential function has to be done

[TABLE]

which can be easily solved

[TABLE]

The leading behavior when the size of the system is big enough $(N\gg 1)$ is

[TABLE]

The last integral to be solved, also using the expansion of the exponential function, is

[TABLE]

It finally leads to an expression for the mean escape time in the critical point

[TABLE]

In the limit of very big systems ($N\gg 1$) the mean escape time scales as

[TABLE]

with asymptotic behaviour

[TABLE]

To sum up, it has been obtained analytically the finite size scaling of the mean escape time, defined as the time taken by the system for reaching $m=0$ from an initial condition $m_{i}=1$. It is

[TABLE]

or in terms of the original parameters

[TABLE]

Part V  **Conclusions and outlook

**

CHAPTER   3    Conclusions and outlook

This thesis has addressed a series of ecological problems from the point of view of statistical physics, that has provided the theoretical framework to develop different mathematical models.

The origin of the regular structures of vegetation that are observed in many regions around the world has been studied. They appear in landscapes where there is a limited amount of rainfall during the year, regardless of the type of soil and vegetation. This scarcity of water is an important constraint for the establishment of new plants. Traditionally, it has been thought that the emergence of the patterns comes from the presence of facilitative and competitive interactions among plants acting simultaneously but at different spatial scales. This phenomenon has been referred to as scale-dependent feedback in the literature. The findings presented in this thesis, using mathematical models that contain only competitive interactions, suggest that facilitative interactions could be superfluous if the finite length of the roots is considered in the equations. As an alternative to the scale-dependent feedback, we have introduced the concept of exclusion areas. They are regions, typically between two maxima of vegetation density, where the competition is so strong that it cannot be overcome by new plants. The extension and the location of these areas are given by the interaction kernel of the different species of plants, which is intimately related to the length of their roots. This concept allows to know in which regions the vegetation will disappear and in which it will remain given an initial ditribution of plants. Determining the existence of exclusion areas could also have important implications on the design of farming strategies that minimize the competition in the crops. This would allow an optimal exploitation of the water resources, mainly in arid regions.

In addition, the proposed models follow previous results that allow the use of these patterns as early warning signals of desertification in arid regions, allowing the development of conservation strategies by anticipating the death of vegetation. As the amount of rainfall decreases, the shapes of the distributions show a universal sequence of gaps, stripes and, finally, spots of plants when the water is very limited. This sequence is independent on the species in that particular landscape. Unveiling the basic mechanisms that drive the formation of these structures becomes essential to change the natural tendency that would lead arid regions to dessert states.

The collective searching behavior of some animal species with communication skills has been also studied. Although the main focus of the work is on foraging strategies, our results could be extended to many other situations, such as mating or predation. The influence of different classes of random movements on the results is also analysed. This work constitutes one of the first theoretical approaches tackling the effect that animal interactions have on the duration of their daily tasks. The main result is that the effect of communication on searching times is maximum when they send information at intermediate length scales. Longer communication ranges, that suppose interacting with more individuals, overwhelm the searchers with too much information coming from all the directions. On the other hand, short ranges do not provide all the population with information enough to expedite the search. In both extreme situations the displacements of the foragers lose directionality to the targets. This result is robust against changes in the type of movement, either Brownian jumps or Lévy flights. As a general result, Lévy strategies give faster searches, but Brownian jumps are more influenced by a communication mechanism.

An application of this model to the foraging behavior of the gazelles inhabiting the Eastern Steppe of Mongolia is also presented. The steppe is one of the largest remaining grasslands in the world, where gazelles have to find each other and small areas of good resources. In addition, because of the orography of this landscape, sound can travel long distances therein. This, together with the strength of gazelle’s vocal tract, allows them to communicate accoustically over long distances. The model predicts an optimal search for resources when the communication is on a frequency of $1.25$ kHz, a value that lies in the range measured for gazelles in the wild ($0.4$ kHz to $2.4$ kHz). This result not only confirms the robustness of the model against changes in the communication channel, but also gives realistic values for the measured quantities. This is our central finding, and suggests that, during its evolution, the species could have optimized its vocal tract to efficiently communicate in the steppe. This work aims to open new research lines in the interrelation between communication, optimal search and mobility patterns. From a theoretical point of view, we propose a new collective searching strategy that offers a wide range of potential fields of applicability, even far away from an ecological context. Similar algorithms, based on collective animal behavior, have been recently implemented in collectivities of robots to tackle different problems [Penders et al., , 2011; Werfel et al., , 2014]. Due to its simplicity, our model could be applied to several searching processes, optimizing the first hitting times if the individuals are enforced to communicate over intermediate lengths. Furthermore, the comparison between Brownian jumps and Lévy flights makes possible to choose the mobility strategy that better works in a given scenario.

Finally, it has been studied the effect of external variability on the diversity, robustness, and evolution of many interacting particles systems. The behavior of the crossing times changes substantially when driven by fluctuating environments. It appears a finite region around the critical point where this time scales with the system size as a power law with continuously varying exponents. These results have clear implications in the mean lifetimes of species in an ecological context (species coexistence, competition…) and also allow to extend the concept of Temporal Griffiths Phases, originally found in epidemic spreading models, to a larger variety of systems.

There are several open questions to be tackled in the future. Many of them are related with the influence of patterns and the role that different scales play on its formation. Most of the future challenges in ecology involve many spatial and temporal scales. In fact, most of the natural systems do not have a characteristic scale and the observed spatial strucutres most of the times have their origin in phenomena that take place at smaller scale. The key for understanding and predicting lies in unveiling the mechanisms underlying these structures [Levin, , 1992].

In this dissertation, vegetation pattern formation has been addressed developing mathematical models with a single spatial scale. However, trees also present facilitative interactions, related to the size of the tree canopy, that act at a shorter scale than competition which are mediated by the roots. Although one of our main results is that positive interactions could be superfluous in the formation of patterns, they could have further implications on its shape and stability. In addition, nature is full of examples where many interaction scales are involved, as it is the case of the regular distributions observed in many mussel beds. Beyond spatial degrees of freedom, ecological systems also show different organizational scales. In the particular case of plants, they are not isolated in the landscape but in interaction with many other species that influence its evolution. This is the case of termites or some microbials, that are known as ecosystems engineers. Investigating its influence on the evolution and formation of the vegetation distributions constitutes a promising challenge.

Establishing relationships between vegetation distributions and animal mobility presents also many challenging questions. Most of them should focus on merging both research lines, addressing the influence that grazing could have on the patterns, and how the formation of groups of animals could modify their shapes or destroy them. While larger groups have clear benefits in terms of group defense and predator swamping, they also lead to a faster degradation of the vegetation. This is the problem of foraging influencing the vegetation patterns which should be treated in the future.

The study of how an information flow can modify collective searching processes is attracting more attention the last years. In this work, we have studied how the ranges at which the information is shared modify the duration of foraging. However, many other questions, such as how informed individuals in a population may adapt its mobility pattern in order to increase the success of the group remain still open.

In summary, coming years promise an intense activity trying to answer these and more open questions. Statistical physics is now much more than a discipline devoted to the study of the macroscopic properties of thermal systems, and theoretical ecology is a well established quantitative field. Their development during these years has brought them to a common point, from where natural environment can be better described and understood.