New forms of structure in ecosystems revealed with the Kuramoto model

Ecological systems, as is often noted, are complex. Equally notable is the generalization that complex systems tend to be oscillatory, whether Huygens' simple patterns of pendulum entrainment or the twisted chaotic orbits of Lorenz’ convection rolls. The analytics of oscillators may thus provide insight into the structure of ecological systems. One of the most popular analytical tools for such study is the Kuramoto model of coupled oscillators. We apply this model as a stylized vision of the dynamics of a well-studied system of pests and their enemies, to ask whether its actual natural history is reflected in the dynamics of the qualitatively instantiated Kuramoto model. Emerging from the model is a series of synchrony groups generally corresponding to subnetworks of the natural system, with an overlying chimeric structure, depending on the strength of the inter-oscillator coupling. We conclude that the Kuramoto model presents a novel window through which interesting questions about the structure of ecological systems may emerge.

JV, 0000-0002-3366-4343; ZH-F, 0000-0002-9092-4500 Ecological systems, as is often noted, are complex. Equally notable is the generalization that complex systems tend to be oscillatory, whether Huygens' simple patterns of pendulum entrainment or the twisted chaotic orbits of Lorenz' convection rolls. The analytics of oscillators may thus provide insight into the structure of ecological systems. One of the most popular analytical tools for such study is the Kuramoto model of coupled oscillators. We apply this model as a stylized vision of the dynamics of a well-studied system of pests and their enemies, to ask whether its actual natural history is reflected in the dynamics of the qualitatively instantiated Kuramoto model. Emerging from the model is a series of synchrony groups generally corresponding to subnetworks of the natural system, with an overlying chimeric structure, depending on the strength of the interoscillator coupling. We conclude that the Kuramoto model presents a novel window through which interesting questions about the structure of ecological systems may emerge.
where ω i is the winding number of oscillator i (the rate of advancement on the circle dictated by the inherent oscillations), K is the intensity of coupling and N is the number of oscillators. This model is commonly used when the phases of the oscillations are taken to be the key dynamical force. In its classical form, the model assumes all oscillators identical and couplings are taken to be global (all to all). This elementary model yields a simple and universal result. With random initiations of Θ, and low coupling, no synchrony occurs. As coupling intensity (K) increases, a critical value exists at which point all oscillators rapidly synchronize and remain synchronized with further increases in coupling intensity. This model has been useful in studying large systems of coupled oscillators [25]. One interesting result is that among groups of strongly synchronous pairs of oscillators, for some arrangements of coupling, there are individual oscillators that fail to synchronize with the rest, creating a so-called chimeric pattern [26], although for finite N, they are now generally thought to be extremely long transients [27].
Based on the reality of coupling types in consumer/resource systems, when two consumers share at least one key resource, they will experience competition from one another. If their sharing is relatively royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 weak, they will converge on a pattern of relative in-phase synchrony with one another [18]. If, contrarily, the coupling is through the competitive interactions of the resources, the oscillators will converge on a pattern of relative anti-phase synchrony with one another [19]. This arrangement has led to some speculations on its meaning for ecological communities in general [19,20] as well as empirical confirmation in the field [23]. Elsewhere it has been shown that the pattern of coupling based on the classical consumer resource model of MacArthur [28] follows precisely the qualitative predictions of coupling patterns from the Kuramoto model [29].
To represent an actual empirical community, we relax the assumption of the universal constant coupling in Kuramoto's model (equation (2.1)), obtaining, where Kuramoto's mean field approach has been disaggregated with the adjacency matrix Γ (with elements γ i,j ) stipulating the coupling of each pair of oscillators. For a non-weighted, non-directed graph, γ i,j is either 0 or 1. The parameter K is retained since it is used as a tuning parameter to study the system as a whole (rather than simply incorporating it into the γs).
The degree of synchrony is frequently measured by Kuramoto's order parameter, which is the absolute value of z where, and N is the number of oscillators.
The behaviour of the Kuramoto model applied to ecological situations under various assumptions about the distribution of the γ i,j has been analysed, albeit infrequently. For example, Banerjee et al. [30] studied the behaviour of the model with γ i,j distributed as a distance related power law in a spatially explicit framework, and Girón and colleagues examined the consequences of allowing some of the γ i,j to be positive and others negative. In an earlier work, Hajian-Forooshani & Vandermeer [29] compared a six-dimensional predator/prey framework modelled in the Lotka-Volterra style to the same framework in the Kuramoto model showing that the synchronization patterns were qualitatively identical. Numerous studies, not necessarily ecological, have focused on a central feature of the model, the chimeric state that inevitably emerges when the γ i,j collectively are too small to engage all the oscillators in complete collective synchrony but too large to permit complete independence of oscillator behaviour [31]. Others have noted the emergence of subnetworks as synchrony groups [32], depending on the distribution of the γ i,j for several theoretical distributions.

A real ecological system
As interesting as the theoretical diversions of the basic model are, thus far there has not been an attempt at applying the Kuramoto framework to a real ecological system. Here, we use the well-documented system of pests and their natural enemies in the coffee agroecosystem [33,34] as a system in which predators/parasitoids/diseases attack four well-known pests of coffee and, therefore, can be viewed as a single system of oscillators. The four pests are (i) the coffee berry borer (Hypothenemus hampei Ferrari), (ii) the coffee leaf miner (Leucoptera coffeella Guérin-Méneville), (iii) the green coffee scale (Coccus viridis), and (iv) the coffee leaf rust fungus (Hamaelia vastatrix Berk. & Broome). We ask, to what extent does the Kuramoto model provide insight into the community structure of this welldefined real community? Will the Kuramoto model reflect what we know about the system? Will it provide any insights regarding the structure of this community? We diagram the basic system in figure 1. It is critical to note that this rendering, although it includes top predators and parasitoids (those that eat the consumers of the prey), is a simplified rendering, emphasizing the direct energy transfer interactions and ignoring the well-known indirect higher order interactions [33]. Anticipated further studies will incorporate those more complex additions of reality.
Applying the model to this well-studied system, we seek to determine whether modules of synchronization appear and whether these modules reflect the reality of what we know about that real system. It is not claimed that the Kuramoto model will provide extra evidence that is not apparent from the fairly obvious structure of this simple community. Rather, we seek to (i) use the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8:

The Kuramoto model and synchrony groups
The dynamics of coupled ecological oscillators is complicated and strongly dependent on the nature of the coupling [8,18,29], with in-phase synchrony resulting from predators sharing prey and anti-phase synchrony resulting from prey competing with one another. Empirical verification of these theoretical propositions is rare [23]. However, the network we seek to investigate (figure 1) emerges from the 'pest' guild associated with the coffee agroecosystem and thus involves only predators consuming prey, and thus are expected, if the coupling is weak enough, to produce nearly in-phase synchrony. More complicated networks that include competition at lower trophic levels would require, minimally, the expansion of equation (2.2) to include both positive and negative K [35]. Here, all of the oscillators are predator/prey (in principle). Based on the simple idea that oscillators are coupled as long as one element is shared (two predators eating the same prey, two prey eaten by the same predator or a trophic triplet [chain]). Thus, for example, oscillators 8 and 3 are coupled because they are part of a trophic triplet whereas oscillators 8 and 5 are not coupled because they neither share a resource or are part of a trophic triplet. Based on these couplings, we established an adjacency matrix, which is reflected visually in the graph of figure 1b. Using this matrix as Γ, we employed the extended Kuramoto model at various coupling strengths (with identical ω i = 0.01), where every coupling of oscillators (figure 1) is at the same strength, K. The general result shows particular patterns of synchrony, occurring in groups of oscillators, what we call 'synchrony groups'. We judge two oscillators, i and j, to be in the same synchrony group if the difference in Θ for the two oscillators is less than C, where the critical value of C is a parameter that may be tuned, and indicates membership  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 in a group. For every pair of oscillators we computed, to compare with C (for all simulations reported herein, C = 0.01 radians). The 'stability' of synchronous groups was determined by randomly initiating the model 20 times and recording the number of times that oscillators i and j fell into the same group. If two oscillators fell into the same group in at least 19 of the 20 runs, they were judged to be in the same 'persistent' synchronous group. Exploratory simulations with 50 runs also run for a few examples to verify that the qualitative results from 20 runs are robust. Calculating the order parameter over a range of values of the coupling coefficient, we obtain the usual rapid increase to full coupling (figure 2). Contrary to the classic Kuramoto model, we do not see the lower plateau at lower values of the coupling coefficient, nor a critical transition to full coupling, as occurs in the classic model. We attribute this to the complexity of the coupling of the system, as well as its small size [36]. With the coupling K = 0, a relatively high value of the order parameter is obtained (approximately 0.4 as an average), so any approach to 1 is obscured by the small range (from less than 0.1 to 1.0), suggesting that the order parameter is not necessarily an efficient measure of synchronization with such a small number of oscillators.
Contrarily, when we extract only the oscillators associated with an ecologically significant subset of the oscillators, such as the coffee berry borer sub-community (oscillators 17-21, in figure 1), we obtain a result that mirrors the classic results of the Kuramoto model ( figure 2b). This suggests that key, wellconnected modules of interacting oscillators within our network behave qualitatively similarly to the classic Kuramoto model (figure 2b), but this dynamic is obscured when considering the full network of oscillators (figure 2a), as is typically done with uniformly global coupling. As with many real networks [37], most nodes in this community have few links, potentially increasing their ability to synchronize and increasing the sum that is the real global order parameter. This result highlights the utility of the criterion |c i,j | < C, from equation (4.1).
The amount of scatter in the order parameter before full synchronization reflects a complicated approach to that synchrony, based on the sequential synchronization of independent synchrony groups. In figure 3, we illustrate the approach to complete synchronization for one exemplary run. The structure of the synchrony graph is based on K = 0.01, and we see a clear division of oscillators in three distinct sets in figure 3a. It is notable that oscillator 3 does not synchronize with any of the three groups. To fully understand the significance of the graph, it is convenient to follow the trajectories of all the oscillators over time as in figure 3b. It is evident that all the oscillators eventually synchronize royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 as they converge to a common oscillatory angle θ. Yet the pattern of synchronization is not a uniform coming together of all the oscillators. It is clear that synchrony groups form rapidly (about t = 5) and subsequently the groups themselves synchronize, perhaps most clearly seen in the three diagrams in the complex plane (figure 3c) placed at approximately the same position as the time series of the angles representing the oscillators. Thus, the organization of synchrony groups is clearly a transient phenomenon, suggesting the time to synchrony as a potentially important characteristic of the synchrony group itself.
In figure 4, we illustrate the basic pattern of development of the stable synchrony groups as the coupling coefficient (K) is increased, harvesting the simulations at t = 20 (an arbitrary designation based on empirical observations- figure 3). Note that the coffee berry borer group, narrowly construed (i.e. oscillators [17][18][19][20][21] begins synchronizing at a very low value of K. Observations of the process reveal that any two of the five oscillators can synchronize first at K ∼ 0.010, while any three of the five at K ∼ 0.014, and any four of the five at K ∼ 0.015. These small synchrony groups are not 'stable,' in the sense that they appear from all initiation points at these low coupling coefficients and  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 probably reflect the random positioning of the oscillators at the beginning of a simulation. Similar observations are evident for the coffee leaf rust group (oscillators 10-13). It is clear, however, from figure 4, that once a synchrony group is formed, it is invariant relative to further increases in K, with the further addition of other oscillators as K increases. For K ∼ 0.6, four distinct synchrony groups are stable. It is also notable that the make-up of the groups corresponds quite well to the general nature of the network with the four groups corresponding to the biological control system of each of the pests, another group associated with the phorid parasitoid, forms at a higher coupling coefficient. As coupling reaches approximately the level of 0.7, the coffee berry borer group merges with the leaf miner group, and the scale insect group merges with the coffee leaf rust fungus group at approximately K = 0.9. Oscillator 3, the predation of a spider on the small parasitoids, synchronizes only at the point that all oscillators are in synchrony.
It is also worth noting that over specific ranges of coupling, there are some oscillators that act something like oscillator 3, for a range of coupling values. For example, in the range K = 0.1-0.4, while there are two synchrony groups in which seven of the 22 oscillators are involved, while the remaining oscillators fail to synchronize either with one of these groups or among themselves. Such isolated oscillators are sometimes referred to as chimeric elements in the overall structural framework. Thus, when K = 0.6, for example, we could describe the 'community structure' as consisting of four specific synchrony groups plus eight chimeric elements, at K = 0.7 we have three synchrony groups plus four chimeric elements, and at K = 0.8, we have four synchrony groups plus one chimeric element.
Thus, a comparison of the model with this real-world ecological network reveals a complex set of patterns within the synchrony group range. Of particular interest is the fact that the larger synchrony groups are not necessarily predictable from a qualitative interpretation of the elementary network structure (figure 1), nor is the timing of synchrony groups obvious. For instance, oscillator 3 is coupled with both oscillator 8 and oscillator 2, and those two oscillators are in different large synchrony groups (i.e. when K is high), explaining why oscillator 3 only joins in synchrony when the large group synchronizes, since even at very high coupling strengths it apparently remains chimeric. Furthermore, the coupling of 8 with 3 connects what would otherwise be completely independent networks ({1-4, 14-22} and {5-13}-see inset in figure 1). Thus, if oscillator 3 were eliminated from the system, the complete synchronization of the entire ensemble would not occur at all, even though that very oscillator is extremely resistant to synchronizing with other oscillators, because of its coupling to elements in the two major independent large synchrony groups.
An alternative look at the basic graph (figure 1b) reveals another obvious division of the graph into independent subgraphs with the elimination of oscillator 22 (leaving the two completely separate groups 1-13 and 14-21). This is an interesting point of departure for examining the consequence of different weightings of the coupling coefficient, since this particular connection is indeed a bit special in the real system. Oscillator 22 involves the keystone ant species Azteca sericeasur and its consumption of the coffee berry borer. There is some debate among experts (see for example, the alternative interpretations of Perfecto & Vandermeer [38] versus Jiménez-Soto et al. [39] as to how much energy transfer occurs in this interaction). Furthermore, there are dramatic trait-mediated indirect effects involved [40][41][42] which ultimately must translate into at least a weak coupling of 22 with all the other oscillators associated with the coffee berry borer. Consequently, we did a series of simulations setting the subset of coupling coefficients involving oscillator 22 to 10% less than all the other couplings. The results of this were that, on the one hand, the basic smaller synchrony groups emerged similarly as to when all coefficients were constant, suggesting that it is the structure of the network that mainly drives the overall community dynamics and structure, in the case with these adjustments to K. On the other hand, subtleties did emerge. For example, in our original analysis, oscillator 3 (an orb-weaving spider, Pocobletus sp., a Linyphiidae catching parasitoids in its web) failed to synchronize with any group unless the coupling became very strong ( figure 4). Yet, after we lowered oscillator 22 couplings, it merged with one of the two major groups [1][2][3][4][5][6][7][8][9][10][11][12][13] at the relatively low general coupling of K = 0.5. Based on the real network structure (figure 1), we infer that when oscillator 22 is tightly coupled with its associated pairs, oscillator 3 is pulled in both major synchrony group directions, first because of its 'indirect' coupling with oscillator 22 (through the couplings with oscillators 2 and 4) pulling it in the direction of group 14-21 simultaneously with its coupling with oscillator 8, also pulling it in the direction of group 1-13. As the system becomes completely coupled, the major transition that previously occurred from four large synchrony groups to two (and eventually to complete synchrony) changes dramatically to include entirely different members when the coupling coefficient of oscillator royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122

Discussion
Reflecting on the prescient ideas of Platt and Denman, we suggest that joining their perspective with the powerfully elegant model of Kuramoto provides a new window through which ecological communities and ecosystems might be examined both theoretically and empirically. Here, we study this framework as applied to a real-world example of four pest species and their biological control elements, all of which are conceptualized as oscillators and coupled together in specific ways according to published studies. We argue that this mode of analysis generates new ways of looking at ecosystems. In particular, rather than a focus on species population densities (or biomasses) and attendant features such as diversity and traditional stability, it might be useful to analyse the system based on synchronization groups and patterns of synchronization. Such a suggestion depends on an initial assumption that the ecosystem is itself relatively permanent, which is to say the traditional questions of stability (Lyapunov or more generally) are not involved. Rather we look, theoretically, at patterns of approach to synchrony and the synchrony groups that emerge, and empirically, at spatial and temporal correlations predicted by the Kuramoto model (i.e. elements co-occurring in a synchrony group would probably be correlated empirically). Somewhat similar to the change in scale of the description of groups of populations into meta-populations, the attraction of our approach is that it permits the prediction of coarser scale synchrony patterns to be expected based on nothing more than a qualitative understanding of consumption patterns.
An issue that has attracted a great deal of attention in the literature surrounding the Kuramoto model is perhaps relevant to ecosystems as well: chimeric elements [25,36,43,44]. The result herein is that increasing the coupling strength of the oscillators in the model may, under at least some circumstances, produce one or more synchrony groups, but also some individual oscillators that 'refuse' to synchronize with any group at all. We refer to those as chimeric elements. In the most extreme case, a large number of identical oscillators synchronize at some critical coupling strength, but other oscillators, identical in every way, do not synchronize and continue wandering in state space, apparently in perpetuity. The idea of chimeric elements is similar in spirit to chaos in that it appears to emerge almost magically under some circumstances, and in physical systems also seems to have some empirical support [45,46].   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 Mathematically there is some question about the reality of chimeric elements, although in practice extremely long transients are quite relevant in practical applications. Given the extensive literature suggesting that chimeric elements, even if just long transients, exist in a wide variety of physical systems, we might also expect them to occur in ecological situations. If this be the case, we can summarize further a vision of community structure as the structure of synchrony groups plus associated chimeric elements. So, for example, the synchrony groups that emerge in the present example are clearly associated with particular pest species, which is not surprising. The parallel question of how they relate to the natural history of the empirical system stems partly from 'competition' from distinct synchrony groups. Thus, oscillator 3 in the present system, by bridging the gap between the two main synchrony groups (at relatively high values of K) is pulled in both directions and, at best, stabilizes in a position between them, but never moves toward one or another.
Yet the idea of a chimeric element perhaps recalls a basic idea in community ecology. It is easy to postulate other, perhaps more complicated, ecosystems, or perhaps with additional complications to the basic Kuramoto framing, as having chimeric elements in addition to synchrony groups. For example, the idea of fugitive species [47,48], a common referent in early literature, fits this idea perfectly, in that we could imagine the normal semi-stable or permanent community, structured along the lines of synchrony groups, but the fugitive species coming and going, apparently without fitting into the basic community structure of the permanent resident species. Other formulations from classical ecology are also possible (e.g. the very nature of the keystone species, if coupled to two relatively equally attractive synchrony groups could also be chimeric, as in the case of oscillator 22).
Our general results may be summarized as follows: first, as the strength of oscillator couplings increases, synchrony groups tend to form, each one of which is associated with a particular pest species. Second, these groups tend to merge, ultimately forming three major groups, (i) the coffee berry borer/miner group, (ii) the coffee leaf rust/scale insect group, and (iii) the phorid group (figures 3 and 4), all mirroring ecological structures within the system. Third, relaxing the assumption of perfectly uniform couplings, allowing the A. sericeasur/coffee berry borer oscillator (number 22) to be only a tenth of the others (based on our knowledge of the system), provokes changes, the main one of which is that the phorid group synchronizes with the coffee leaf rust/scale group rather than the coffee berry borer/miner subgroup as it did previously. We conclude that the general pattern of synchrony group formation is a consequence of the structure of the network, but that particulars of synchrony group formations may vary depending on the heterogeneity of oscillator coupling strengths, and, of course, the inevitable stochasticity of nature is bound to have an effect. All of the synchrony group results that emerge from the Kuramoto model make perfectly consistent sense with what we know of the system in the field, an encouraging result.
It is satisfying that the Kuramoto model does seem to represent the underlying qualitative dynamics of the real system. Since this system is about pests of an important crop, it is obvious of potential practical importance that an approximate toy model such as this one can reflect the basic structure of this system as we have come to know it from over two decades of study [33]. There is a clear relationship between the number of nodes and the timing of the formation of synchrony groups, suggesting that a herbivore guild containing multiple species attacking it is likely to synchronize the whole group quite rapidly. Yet there are other examples in which it seems to be the relative isolation of the group that allows for rapid synchronization as is the case of Pseudacteon and its antagonists, of clear importance to the key spatial patterning of the system as a whole [33,49]. From practical experience on the ground, it is clear that these various synchrony groups do indeed occur in the field, correlated in space and time.
However, it must be noted that this real-world system itself is represented in a simplified fashion, as we know from the large amount of empirical work that has already been published (see references in electronic supplementary material). Most important are the trait-mediated indirect interactions (or higher order interactions) associated, for example, with oscillator 22, which are known to be important empirically [39], and, given the analysis here, are likely to be important in the further detailed study of synchrony groups of the system [41,42]. Incorporating such complexities is a challenge for the future. Yet it is worth emphasizing that the apparent concordance of the qualitative nature of synchrony groups is remarkable especially since the model is very explicitly ignoring some key features.
Nevertheless, even with this simplified version of the real system, certain patterns suggest qualitative predictions. For example, the set of oscillators associated with the scale insect [5][6][7][8] is almost always a synchrony group, suggesting that spatial and temporal correlations among the members of that group should be observable in nature. Such a prediction is, in our practical experience in the field, certainly true [50,51]. Furthermore, our model of the system also provides new predictions about the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210122 spatio-temporal correlations that may be apparent given the modification of the underlying coupling of components of the system. This may be realized in the empirical agroecosystem through several avenues such as coupling strengths being modified due to seasonal factors or under different regimes of management of the agroecosystem.
The present study is in the spirit of many other network studies in framing the question, 'What can I say about the system knowing nothing more than which nodes are connected?' While it is tempting to compare our results with other approaches, for example, foci on connectedness (e.g. [52]) or stability [53], our intent here is to simply offer an example of how the Kuramoto model concords with the basic natural history of a real system, using a non-traditional framing of the network. Since the nodes we propose here are not species or populations but, rather, oscillators, and our focus is not stability or permanence but rather synchrony, perhaps some unique insights emerge. Connecting this theoretical programme with empirical work provides, perhaps, more insight than a direct application of a dynamic model (e.g. a system of ODEs) to the system. The pattern of oscillator couplings (which oscillators are coupled with which), determination of which is a rather easy empirical exercise, suggests questions about qualitative patterns of co-occurrence in the field (which elements are expected to be correlated in their spatial and/or temporal occurrence) which would be subject to empirical verification.