# Spatial patterns of tree yield explained by endogenous forces through a correspondence between the Ising model and ecology

^{a}Department of Environmental Science and Policy, University of California, Davis, CA 95616;^{b}Department of Physics, University of Massachusetts, Amherst, MA 01003;^{c}Department of Plant Sciences, University of California, Davis, CA 95616;^{d}Land Health Decisions, World Agroforestry Centre (ICRAF), Nairobi 00100, Kenya;^{e}Santa Fe Institute, Santa Fe, NM 87501

See allHide authors and affiliations

Edited by William Bialek, Princeton University, Princeton, NJ, and approved January 2, 2018 (received for review November 14, 2016)

## Significance

Explaining correlations across space of cyclic dynamics in ecology is a fundamental challenge. We apply ideas from statistical physics, originally used to explain the behavior of magnets, to a dataset on yield from pistachio trees, obtaining a robust description and potential explanation for the generation of spatial correlations in cyclic dynamics. These results suggest looking for mechanistic underpinnings at the level of interactions between neighboring trees that lead to spatial correlations in dynamics and a surprising correspondence between the descriptions of physical phenomena, magnetization, and ecological dynamics. This work demonstrates with data, and not just models, that correlations in cyclic dynamics can be generated from local interactions and dynamics even in a very noisy ecological system.

## Abstract

Spatial patterning of periodic dynamics is a dramatic and ubiquitous ecological phenomenon arising in systems ranging from diseases to plants to mammals. The degree to which spatial correlations in cyclic dynamics are the result of endogenous factors related to local dynamics vs. exogenous forcing has been one of the central questions in ecology for nearly a century. With the goal of obtaining a robust explanation for correlations over space and time in dynamics that would apply to many systems, we base our analysis on the Ising model of statistical physics, which provides a fundamental mechanism of spatial patterning. We show, using 5 y of data on over 6,500 trees in a pistachio orchard, that annual nut production, in different years, exhibits both large-scale synchrony and self-similar, power-law decaying correlations consistent with the Ising model near criticality. Our approach demonstrates the possibility that short-range interactions can lead to long-range correlations over space and time of cyclic dynamics even in the presence of large environmental variability. We propose that root grafting could be the common mechanism leading to positive short-range interactions that explains the ubiquity of masting, correlated seed production over space through time, by trees.

Spatial dynamics and patterning are central issues in many fields of science (1⇓⇓⇓–5), including ecology, where the interface of theory and data has the potential to provide unique insights into emergent behaviors and underlying causes (6). Among the most dramatic emergent behaviors in ecology is “masting” of perennial plants, in which individual fruit or seed production varies through time and is correlated across space (7⇓–9). In some tree species, individual production is primarily described by a two-cycle alternating of years of high and low production (10, 11). Patterns of spatially synchronized two-cycles are of great interest as they are found not only in trees but also in systems ranging from childhood diseases to small mammals (12⇓⇓⇓⇓–17).

Determining aspects of process from pattern is a fundamental problem because setting up experiments or gathering observations on relevant temporal and spatial scales is logistically challenging, so experimental treatments have been limited to microcosms (17). Here, we take advantage of an agricultural system that is on a much larger spatial scale than typically available experimentally (11). This kind of agricultural dataset then leads to the challenge of developing a modeling approach focused on the relevant details in a very large complex system. Even the first step of developing appropriate measures and expectations for spatial synchrony has proved to be a difficult challenge (17).

There are two approaches to explaining spatiotemporal patterns in fruit production. One could ask, What are the advantages from an evolutionary standpoint? Or somewhat differently one could ask, What are the proximate mechanisms, both from biological and from modeling frameworks, that lead to synchrony in reproduction by plants at different spatial scales (8)? We take the latter approach since the behavior of the trees in the agricultural systems we study is not the subject of direct selection but perhaps reflects evolutionary history. Our motivation is to use modeling and data analysis to infer the plausibility of potential proximate mechanisms that could then be further analyzed.

The basis of our work is that the ubiquity of spatiotemporal patterning in tree fruit production implies that there should be relatively generic, detail-free explanations. It is also plausible that direct interactions are short range, yet correlations are observed over longer spatial scales. Finally, a deterministic approach is not sensible as large environmental noise typically affects both agricultural systems and more general ecological systems. The search for this kind of robust system-independent description leads us to ideas from statistical physics.

Approaches from statistical physics have provided insights into how local behavior and short-range interactions can generate emergent spatial patterns ranging from desertification (18⇓–20) to schooling behavior of fishes and flocking behavior of birds (21, 22), but have not previously been used to understand patterns and causes of spatial synchrony in cyclic dynamics in data from population biology. More generally, statistical physics seeks to explain emergent behavior, such as scale-invariant pattern formation, in large (high-dimensional) complex systems (3⇓–5, 18, 19, 23, 24). A key theme is critical phase transitions between order and disorder (3⇓–5, 23), which lead to identifiable long-range spatial patterns and other characteristic behaviors. Identifying these characteristic signals in data can then be used to infer properties of underlying mechanisms that provide a parsimonious explanation for observed behavior. We previously used this approach to demonstrate, within a modeling framework, how long-range (i.e., power-law decaying) spatial correlations in systems of coupled ecological oscillators are explained by the Ising model, which is a canonical model of statistical physics (24).

## Ising Model Description of Spatial Ecological Dynamics

The Ising model is a standard model in statistical physics for understanding how scale-invariant correlations emerge in a model-independent, or universal, way at critical transitions (3⇓–5). Our previous work on Ising universality of critical transitions in the spatial synchrony of ecological population models uses the isotropic Ising model on a 2D lattice (24). Since the vertical and horizontal spacing of trees differs in our study area, here we use the anisotropic Ising model in which the horizontal couplings *Materials and Methods*). The couplings

In the context of the Gibbs distribution, ordered (synchronized) states can arise either from strong couplings (J) or strong fields (h); however, within a dynamical Ising model it may be possible to distinguish these alternatives. If an ordered state in a large system arises quickly from a disordered state, it must be due to a strong field h rather than strong couplings because the ordering effect of strong couplings requires long times to propagate over large distances.

## Results

Observations of masting based on individual production of over 6,500 trees over a 5-y period in a commercial pistachio (*Pistacia vera*) orchard (11, 25) (Figs. 1 and 2) provide a dataset that is distinguished by its size and underlying homogeneity, presenting a unique opportunity for using ideas from statistical physics to analyze an ecological system. This system is at least an order of magnitude larger than laboratory microcosms (26, 27) and much more homogeneous than epidemiological datasets (16), which have provided much of the basis for understanding spatiotemporal dynamics in ecology. By using ideas from statistical physics to analyze the pistachio dataset, we can begin to assess the relative importance of exogenous forces, such as the number of winter chilling hours for trees, the amount of precipitation falling during the spring flowering season, and irrigation and other management practices, from endogenous local forces, such as chemical signaling, trophic interactions with microorganisms, and root grafting (28, 29).

We do note one potential complication in simply relating nut production to the state of the tree that is important in our analysis. Yield from an individual tree in a year is a function of the state of the tree (e.g., carbon storage) when it awakes from dormancy and environmental factors during flowering and fruit development. Some of the environmental factors, which affect flowering, also affect the future state of the tree, as reduced flowering or pollination means fewer resources are required to meet vegetative and reproductive needs and more are available for storage and use the subsequent year (30). But other factors, such as those producing “June drop” (31) or nut flower bud abscission, reduce yield but do not change the state of the tree as much as if the nuts were to develop and be harvested (32). Thus, in any given year, the state of the tree is a random variable which is a function of the state of the tree the year before, and the yield from a tree can be thought of as a function of the state of the tree and a random component which has large but not complete spatial correlation among all of the trees.

Given the high degree of underlying homogeneity in soil and moisture conditions as well as the use of a single cultivar in the study area and its small spatial extent relative to the typical length scales over which significant variation tends to be found in environmental conditions, the a priori expectation for pistachio production in the study area would be spatially uniform and highly synchronized by long-range, exogenous forces—the Moran effect (33). This is the case for 2005–2007 (Fig. 2). In addition, total study area production in 2005–2007 is strongly correlated with county and regional production (Fig. 3). In the context of the Ising model, a strong field h is present at least in 2006 – 2005 and perhaps subsequent years. Additionally, the June drop phenomenon means that the observation is a function of both the state of the tree at the beginning of the growing season and an additional random factor. Nonetheless, the entrainment of local orchard-level production to regional production and ultimately to climatic conditions is a particularly clear signal of a Moran effect dominating the system. In principle, strong local couplings might also produce the synchrony clearly observed in the 2006 – 2005 orchard data. However, since synchrony is absent in 2003–2004, there is insufficient time for local interactions to synchronize the whole orchard.

This baseline expectation of complete synchrony makes the observed deviations in 2003 and 2004 especially surprising. In our analysis of these 2 y, we find highly heterogenous spatial patterns in the production data (Fig. 2), with clusters of synchronized production on many different length scales in the study area. Simultaneously, we observe a desynchronization of local production from the regional average (Fig. 3). Thus, there is no strong global forcing and, in the Ising setting, the field h is weak or absent. The fractal-like structure of the 2003 and 2004 spatial patterns in pistachio production is strikingly similar to spatiotemporal patterns previously found in theoretical predictions for the behavior of locally coupled noisy two-cycle ecological oscillators (24), in which the emergent fractal patterns were shown to be model independent and described by the 2D Ising universality class of statistical physics. The spatiotemporal use of the pistachio system has all of the ingredients required to observe Ising-like behavior: The pistachios are planted on a 2D grid; the individual production of alternate bearing female trees oscillates in a noisy nonlinear two-cycle (10, 30); and there are potential sources of positive local coupling among neighboring trees, such as root grafting.

This motivates us to estimate the same first difference statistics in the pistachio data that were used to map ecological population dynamics onto the 2D Ising model in our theoretical work (24). Calculating pairwise correlations in the “2004 – 2003” first difference field, we find a result which displays long-range correlations (Fig. 4*A*) significantly different from expectations for an uncorrelated random spatial process or the kind of patchiness described in ref. 18. Indeed, we find strong evidence of long-range correlations in the pistachio data consistent with the 2D Ising critical point where the exponent of the power law decay is 1/4, as seen in the data. Additional evidence of this mechanism is obtained by explicitly demonstrating that a critical point of a 2D anisotropic Ising model, a model with only two free parameters, is consistent with all observed spatial pairwise correlations in the binarized *B*, Fig. S1, and *Materials and Methods*). We also find evidence of a strong asymmetry in horizontal and vertical couplings, which, given the larger spacing between rows than along rows, is consistent with root grafting being the dominant endogenous mechanism driving the emergent fractal pattern in production. The fractal structure of the 2D Ising critical point emerges only for model parameters near the critical line, and a significant asymmetry in horizontal and vertical couplings is required to obtain a good fit. A parsimonious explanation for the patterns observed in these data is the emergence of long-range spatial correlations from short-range interactions in this noisy system.

Despite strong synchronization in 2005 and later years, remnant long-range correlations in the difference fields are observed in these years (Fig. S3). Unlike the 2004 − 2003 correlations, these correlations are visible only in the full production data and disappear if the data are binarized. An explanation of these correlations would thus require going beyond the (binary) Ising model. Fig. S2 shows the time correlations of the first difference fields in later years with the 2004 − 2003 difference fields. This oscillatory pattern reveals both the underlying periodicity and the slow decay of time correlations.

## Discussion

Our results both have implications for ecology and open up additional applications of methods from statistical physics. Although we have previously shown a correspondence between the Ising model and dynamic theoretical spatial models in ecology (24), we did not tie this finding to data. This left open the question of whether these results could describe ecological phenomena, given the simplifying assumptions made in the models. In fitting the Ising model to the admittedly very noisy ecological data we show the potential for a strong correspondence between model and data.

Here, we observe long-range correlations in pistachio data that are consistent with the 2D anisotropic Ising model near criticality. This is a unique demonstration of long-range spatial correlations arising from endogenous forces in population biology. Analogous approaches have been used to study networks of artificial neurons (34) and cell membranes (35), but the discovery of Ising-like long-range correlations in a masting ecological system is unique. Even with the massive, by ecological standards, dataset we use, the limitations of ecological data mean that the analysis cannot be as complete and quantitative as that in these other areas of biology.

Three obvious challenges to finding this kind of phenomenon in ecological systems are that (*i*) ecological systems either are relatively small in spatial extent or have large underlying spatial variability, (*ii*) stochasticity plays a major role, and (*iii*) dynamics of ecological systems are typically emerging over relatively short timescales and the data available will typically be over an even shorter timescale. One further challenge is not a problem with the data we have examined in this paper: Datasets that encompass a large enough spatial extent in a homogeneous environment with sufficiently precise measurements are exceedingly rare. We first discuss our results using this dataset from the standpoint of statistical physics and then turn to ecological implications.

We would not expect a very clean fit to the patterns predicted by the Ising model for the reasons just outlined, so any correspondence is a dramatic finding. We do see that there are both periods when exogenous forces dominate and periods when the spatial pattern of nut production in the pistachio orchard is consistent with critical behavior of an anisotropic Ising model. That the parameters of the inverse Ising fit are so near the Ising critical point suggests the possibility of self-organized or tuned criticality (34, 36). The good fit of the 2004 − 2003 data to the critical point of an anisotropic 2D Ising model suggests that local interactions such as root grafting (29) may generate sufficient coupling among the noisy two-cycle productions of neighboring trees to produce the observed pattern, a possibility previously predicted on theoretical grounds (24). Alternatively, the correlations observed in pistachio production in 2003 and 2004 could conceivably be the manifestation of some hidden process (37), such as underlying correlations in the soil conditions or microbial communities that are locally coupled to the pistachio trees and influence production. Teasing apart the different candidate mechanisms at local scales driving the spatial and temporal dynamics (Fig. S2) requires further study.

Our finding of correspondence between the Ising model and the yield data has several important implications from an ecological standpoint. First, we have found strong evidence that scale-invariant correlations can emerge from endogenous local couplings in producing the masting behavior of plants. For the specific issue of masting in plants, this has important implications as many previous models focusing on this kind of synchrony relied on pollen limitation and limited pollen production in individuals that had limited fruit production in monoecious species (38). But in the dioecious pistachio, there must be another explanation for interactions between trees as male and female production occurs in different trees. At first glance, one might assume that competition for resources could be the explanation. But this would not fit the model here because we need an interaction that promotes local synchrony, which is a positive local interaction. This suggests important directions for uncovering mechanisms that would lead to synchrony. Various mechanisms of plant–plant communication are known to occur, including root- or leaf-emitted volatile organic molecules (39, 40), acoustic signals (41, 42), and mycorrhizal networks between root systems (43, 44) and intraspecific root grafting. A recent study (45) has shown that the coupling due to mycorrhizal networks can be truly substantial in terms of the amount of carbon exchanged among roots.

Intraspecific root grafting (28) is one ubiquitous example of short-range coupling between trees that is present both in our agricultural system and in natural systems. What makes this observation even more intriguing is that root grafting is common in many genera that also exhibit masting (7, 8). In particular, root grafting is common (28) in two genera, *Pinus* and *Quercus*, pines and oaks, respectively, that are also prime examples of masting (7). Our confirmation of the possibility that general short-range positive interactions can lead to correlations in cyclic dynamics at large spatial scales is consistent both with general observations and with the observation of root grafting in the particular orchard from which the data analyzed here come. Of course, mycorrhizal networks are another possible explanation.

More generally, our result may provide keys to understanding complex spatiotemporal patterns in other ecological systems, ranging from the masting in perennial plants to disease dynamics. Obtaining the requisite high-resolution spatiotemporal data as well as extending our ideas to more complex spatial arrangements represents a large challenge, but could greatly increase our understanding of spatiotemporal processes in population dynamics. However, the universality of the behavior of the Ising model holds out hope that the results would be quite robust and would apply even to systems that are quite far from the assumptions underlying our analysis here. Conversely, the success of the approach here would also suggest that applying other ideas from statistical physics to studies of spatiotemporal dynamics in ecology could prove very fruitful both in explaining ecological phenomena and in suggesting directions that could be interesting from the standpoint of physics.

## Materials and Methods

In this section we analyze the data from the years 2003 and 2004, seeking the best-fit Ising model. We seek parameters

On the modeling side, for any given values of **4**. Repeating this procedure for all 2D Ising samples gives us the mean value for each pairwise correlation,

There are other modeling assumptions of ecological interest. First, we cannot ignore the male trees. In fact, the fit presented here assumes that the local interactions of male and female trees are identical. This result suggests that levels of the unidentified endogenous signal in male trees also oscillate in a noisy nonlinear two-cycle and that the local coupling (root grafting is one example) is independent of sex. Second, simulations restricted to the *A*. A good fit is obtained only when we simulate the 2D Ising model on the full

The county data summarized in main text and Datasets S1–S3 were obtained from standardized reports made publicly available by Agricultural Commissions of Fresno, Kern, Kings, and Tulare Counties (Fresno, www.co.fresno.ca.us/departments/agricultural-commissioner/crop-report-history; Kern, www.kernag.com/caap/crop-reports/crop-reports.asp; Kings, www.countyofkings.com/departments/general-services/agriculture-department-measurement-standards/ag-services/crop-reports-1941-2013/test; and Tulare, agcomm.co.tulare.ca.us/default/index.cfm/standards-and-quarantine/crop-reports1/crop-reports-2000-2010/). All measurements of individual tree yield are rounded to the nearest pound. Row 0, column 0 references the female tree in the far southwest corner of the study area. The term “-1” indicates a missing value and “-2” indicates a male tree.

## Acknowledgments

We thank Joe MacIlvaine and Wonderful Orchards for assistance and critical insights. We thank Marcel Holyoak, Sandy Liebhold, Ben Machta, Dan Stein, and Steve Strogatz for productive conversations and two referees for very helpful comments. We thank Uriel Rosa for leading the engineering of the pistachio harvester. Simulations were performed on a high-performance computing cluster housed within the Computer Science and Engineering department at University of California, Davis. We are grateful to Bill Broadley and Terri Knight for their excellent administration of that resource. We thank the Santa Fe Institute for supporting several working groups during which this work was discussed. This work is supported by the US National Science Foundation under Integrated National Science Foundation Support Promoting Interdisciplinary Research and Education Grant 1344187.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: andrew.e.noble{at}gmail.com.

Author contributions: A.E.N., T.S.R., P.H.B., J.M., and A.H. designed research; A.E.N. and T.S.R. performed research; A.E.N. and T.S.R. analyzed data; and A.E.N., T.S.R., P.H.B., J.M., and A.H. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1618887115/-/DCSupplemental.

- Copyright © 2018 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Strogatz SH

- ↵
- ↵
- Goldenfeld N

- ↵
- Sethna JP

- ↵
- Solé RV

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Rosenstock TS,
- Hastings A,
- Koenig WD,
- Lyles DJ,
- Brown PH

- ↵
- ↵
- Rohani P,
- Earn DJD,
- Grenfell BT

- ↵
- ↵
- Björnstad ON

- ↵
- ↵
- ↵
- Rietkerk M,
- Dekker SC,
- De Ruiter PC,
- Van De Koppel J

- ↵
- ↵
- Cavagna A, et al.

- ↵
- ↵
- Bialek W, et al.

- ↵
- ↵
- Noble AE,
- Machta J,
- Hastings A

- ↵
- Rosa UA, et al.

- ↵
- ↵
- Bonsall MB,
- Hassell MP

- ↵
- ↵
- ↵
- Stevenson MT,
- Shackel KA

- ↵
- Kaska N

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Klein T,
- Siegwolf RTW,
- Körner C

- ↵
- Box GEP,
- Jenkins GM

- ↵
- Newman MEJ,
- Barkema GT