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

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 , between sites in neighboring columns, and the vertical couplings , between sites in neighboring rows, may be different. The “spin” at site in the anisotropic 2D Ising model is a binary random variable, , and the probability of any configuration of spins, denoted , is given by the Gibbs distributionwhere is a so-called “reduced energy function” that can be written aswith and denoting the sets of all horizontal neighboring spins and all vertical neighboring spins, respectively. The third term in the expression describes the effect of an external magnetic field acting on all spins. The overall negative sign means that the reduced energy is lower (higher) when neighboring spins are more (less) aligned. Because the probability of any given configuration of local spin states over the lattice is given by the Gibbs distribution, configurations with a high degree of local alignment, or ordered configurations, are typical if both and are sufficiently large. Alternatively, if the external field is large, the spins will be aligned with the field. Low levels of local alignment, or disordered configurations, are typical if , , and are sufficiently small. Critical transitions between order and disorder occur for and intermediate values of and . In the current context of masting trees, each spin in the Ising model represents the phase of production of a tree. If , the tree at site produces a large crop in even years and vice versa for . Thus, the spin variables are related to first differences in yearly production (). The couplings and represent endogenous interactions that favor phase locking and the field represents the synchronizing effect of global exogenous forces such as a weather event (Moran effect). Thus, the external field is typically a random function of time.

In the context of the Gibbs distribution, ordered (synchronized) states can arise either from strong couplings () or strong fields (); 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 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).

Fig. 1.

Obtaining high-resolution spatial data on individual pistachio tree production. The Lost Hills orchard in Kings County, CA is planted on a precise grid with 5.2 m between columns and 6.4 m between rows. (A) Annual harvesting is fully automated and proceeds down the rows (Movie S1). The “shaker” shakes each female tree near the resonance frequency of the individual fruits on their small stems. The “receiver” collects the fallen pistachios. (B) The receiver was equipped with an electronic scale and on-board computer (25), allowing for accurate measurement of the production of individual female trees occupying a 9.3-hectare study area (Movie S2).

Fig. 2.

Individual pistachio tree production within study area. Color bar indicates production of individual female trees in each year. Black squares indicate the location of male trees, which do not bear fruit. White squares indicate missing data. The homogeneity of the spatial pattern in 2005–2007 indicates that the exogenous forcing of the Moran effect dominates endogenous forces in regulating production. Indeed, this is the expected pattern of production given that pistachio tree genetics, soil conditions, and management practices are highly uniform across our 9.3-hectare study area. This makes the multiscale heterogeneities observed in 2003 and 2004 surprising. Rather than being entrained to large-scale regional growing conditions, the emergent spatial patterns in 2003 and 2004 indicate that endogenous forces dominate the Moran effect in producing spatial correlations.

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 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.

Fig. 3.

Regional patterns of synchrony. Normalized first differences of total pistachio production time series are plotted for our study area, the surrounding county, and the surrounding four-county region. The alternate-bearing noisy two-cycles of pistachio production are clearly synchronized at all three spatial scales in the first differences , , and , providing evidence that the exogenous forcing, also known as the Moran effect, dominates localized endogenous forces in those years. By contrast, the first differences for our study area and the surrounding county are out of sync with regional production, suggesting a much weaker Moran effect. In more detail, let denote the time series of average pistachio production per tree within the study area. Let denote the SD of X. Then the normalized first difference (NFD) in year is defined as . We use the same formula to calculate normalized first differences for the time series of total pistachio production in Kings County, , and across the region spanning Kings, Fresno, Kern, and Tulare Counties, . is the sum of total production time series for each of the four counties, which is calculated from the publicly available data by multiplying average pistachio production per harvested acre by the total number of harvested acres.

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 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) 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 first difference field (Fig. 4, Fig. S1, and ). 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.

Fig. 4.

Scale-invariant spatial correlations and inverse Ising fit. (A) The connected pairwise spatial correlations of the first difference field exhibit a power-law scaling consistent with the critical point of the 2D Ising model. These data also demonstrate a strong asymmetry between the horizontal correlations of trees within rows and the vertical correlations of trees within columns. A, Inset shows first difference data. The 2D connected pairwise correlations on the nonbinarized data are given by Eq. after replacing with . If m is the vertical distance between rows and m is the horizontal distance between columns, then the distance between trees spaced by and is . We distinguish horizontal, vertical, and diagonal correlations by different colors and plot them as a function of , after normalizing by the value of . (B) Phase diagram of the anisotropic 2D Ising model with typical snapshots in the disordered, ordered, and critical phases. Best Ising model fit to the binarized first difference data is shown as a red circle. Note the similarity between the critical Ising snapshot and the pistachio data displayed in A, Inset. The pistachio best-fit point and error ellipse are obtained by the inverse Ising inference methods discussed above. In C, Insets are representative samples of the anisotropic 2D Ising model. The critical line of the anisotropic Ising model in 2D is known analytically in terms of an implicit function of the horizontal coupling, , and vertical coupling, , .

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 and such that spatial correlations in the pistachio data are well matched by spatial correlations in the 2D anisotropic Ising model. Since global synchrony is absent during these years, we set . On the data side, the 2D lattice is the study area containing female trees. Each of the 66 rows is indexed by an integer and each of the 106 columns is indexed by an integer . Let and denote the pistachio production of individual female trees in 2003 and 2004, respectively. Likewise, let and denote the average production per tree within the study area. We define the first difference (46) field asIf site or site is a male tree, or if data are missing for production at site or site in 2003 or 2004, we say that is “undefined.” The binarized version of is denoted , where if and otherwise. Given that any two trees are separated vertically by rows and horizontally by columns, the 2D connected pairwise correlations are defined by a double sumskipping any terms in which the value of or is undefined and denoting by the total number of terms in the sum.

On the modeling side, for any given values of and , we use the Wolff algorithm (47) to draw independent samples of the anisotropic 2D Ising model on a lattice with free boundary conditions. The size of this lattice is equal to the size of the isolated orchard subplot that contains the study area. We calculate the 2D connected pairwise correlations within the study area just as in Eq. . Repeating this procedure for all 2D Ising samples gives us the mean value for each pairwise correlation, , and the SD, .

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 study area generate large systematic deviations from the observed power-law correlations in Fig. 4. A good fit is obtained only when we simulate the 2D Ising model on the full grid of the surrounding orchard subplot, which is bounded by wide access roads, on all four sides, that inhibit root grafting with trees at the boundaries of neighboring subplots. A cutoff to the power-law correlations in small systems is a known property of the 2D Ising model that, in this context, provides additional evidence that the observed scale invariance of pistachio production in 2003 and 2004 is indeed an emergent phenomenon regulated by endogenous interactions of each individual tree with its four nearest neighbors.

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.