# Beyond the Aggregated Paradigm:
Phenology and Structure in Mutualistic Networks

Clàudia Payrató-Borràsa,b, Carlos Gracia-Lázarob,
Laura Hernándeza, and Yamir Morenob,c,d

###### Abstract

Mutualistic interactions, where species interact to obtain mutual benefits,
constitute an essential component of natural ecosystems. The use of ecological
networks to represent the species and their ecological interactions allows the
study of structural and dynamic patterns common to different ecosystems.
However, by neglecting the temporal dimension of mutualistic communities,
relevant insights into the organization and functioning of natural ecosystems
can be lost. Therefore, it is crucial to incorporate empirical phenology -the
cycles of species’ activity within a season- to fully understand the effects
of temporal variability on network architecture. In this paper, by using two
empirical datasets together with a set of synthetic models, we propose a
framework to characterize phenology on ecological networks and assess the
effect of temporal variability. Analyses reveal that non-trivial information
is missed when portraying the network of interactions as static, which leads
to overestimating the value of fundamental structural features. We discuss the
implications of our findings for mutualistic relationships and intra-guild
competition for common resources. We show that recorded interactions and
species’ activity duration are pivotal factors in accurately replicating
observed patterns within mutualistic communities. Furthermore, our exploration
of synthetic models underscores the system-specific character of the
mechanisms driving phenology, increasing our understanding of the complexities
of natural ecosystems.

a Laboratoire de Physique Théorique et Modélisation, UMR CNRS, Université de
Cergy-Pontoise, 2 Avenue Adolphe Chauvin, F-95302, Cergy-Pontoise Cedex,
France
b Institute for Biocomputation and Physics of Complex Systems (BIFI),
University of Zaragoza, Spain
c Department of Theoretical Physics, Faculty of Sciences, University of
Zaragoza, Spain
d Centai Institute, Turin, Italy

##  1 Introduction

Mutualistic relationships, which for centuries had been mainly overlooked as
fascinating but marginally relevant, are today known to play a crucial role in
shaping natural ecosystems[1]. In mutualistic communities, individuals from
different species interact through activities that provide them mutual
benefits such as food, shelter or an increase in reproductive success –being
pollination of plants by animals the most characteristic example. The
explosion of network science in the 90s fostered the development of a holistic
understanding of these communities using the language of ecological
networks[2], which allowed investigating the existence of widespread
structural or dynamical patterns despite ecosystem’s differences in species
composition, climate or location[3, 4, 5].

In the context of community ecology, the network’s nodes typically represent
different species, while the links (which might be weighted or binary) depict
the ecological interactions among them. In the particular case of mutualistic
networks, the nodes may be moreover separated into two different sets or
guilds –e.g. plants and pollinators– such that the links representing the
mutualistic interaction only exist between nodes of different guilds,
resulting in a so-called bipartite network. This representation can then be
enhanced, for instance, by adding other kinds of ecological interactions like
competition for common resources, which then leads to a bilayered depiction of
the community[6].

Here, we will attempt to understand how by neglecting the temporal dimension
of real mutualistic communities, that is, by working with the aggregated
version of the interaction networks, we can lose relevant insight into the
organization and functioning of natural ecosystems. In this sense, we will
incorporate into the network formalism the information about empirical
phenology, namely the biological activity cycles of species that, to a certain
extent, constrain and articulate how ecological relationships occur among
individuals. This endeavor, although largely overlooked, is not entirely
new[7, 8, 9, 10, 11] and, in fact, during recent years the need of moving
towards a more realistic depiction of ecological communities has been stressed
more than once[12, 2]. All in all, aggregated networks are still commonly used
in the characterization of the structure and dynamics of mutualistic systems,
probably due to the scarcity of data covering different years with information
both about the observed interactions and the phenology. Along this paper we
will address this problem by examining the effects on the perceived network
architecture of taking into account -at least partly- the temporal dimension
of mutualism, using information extracted from two empirical datasets provided
by[13] and [14], which do include simultaneously the dynamical information
mentioned above. Moreover, we will not only consider their consequences for
mutualistic relationships but also for intra-guild competition for common
resources.

Accordingly, the present paper is organized as follows: first, we complete
this introduction by reviewing some basic notions and recent advances in the
study of phenology; next, we propose a framework to characterize empirical
phenology and we propose a set of synthetic models in order to study the
effect of temporal variability beyond the scarce number of open datasets;
third, we present a structural characterization of two real datasets taking
into account the variation of their interaction networks during the season and
compare these observations with the explorations produced by synthetic models
that include stylized characteristics of the phenology; and finally, we
conclude with a discussion of these results and some unresolved, future
challenges.

###  1.1 From cherry blossom to digital cameras

The bloom of flowers in spring, the arrival of the first migratory birds, and
the shedding of leaves in fall are all simple yet beautiful examples of
phenological events that undeniably shape, since ancient ages, our perception
and narration of time, especially the passage of the seasons. Our language and
culture are full of references to this sort of phenomenon and their timing,
from words like ‘late bloomer’ to Aristotle’s famous phrase ‘one swallow does
not make a summer’. Naturally, hence, we can find as well ancient instances of
documented phenological patterns, not only related to agricultural needs but
also in cases in which their timing marked the date of religious or
traditional festivities. A fascinating illustration is the viewing of cherry
blossoms in Japan, known as hanami, a celebration that started in the eighth
century as an elitist ceremony and gradually became a general festivity.
Nowadays, it is so popular that the cherry blossom is nightly forecasted in
Japan, and the advance of the front of blossom across the country -the so-
called sakura zensen\- is keenly followed (see Fig. 1). As a result, it is
possible to find extensive datasets on the phenology of various species of the
cherry tree over the years –mainly the Prunus serrulata-, comprising more than
seven centuries of information on blossoming dates[15]. In fact, this kind of
sequential data has been used to reconstruct springtime temperatures of
undocumented old periods. Similarly, analogous studies have been carried out
using alternative species’ phenology around the globe, such as the grape plant
to produce wine in France[16].

Figure 1: In a), picture called Asukayama hanami no zu (Cherry-blossom viewing
at Asuka hill) by Utagawa Hiroshige, dated about 1831. In b), forecasting of
the dates of blossoming across Japan in 2007. The numbers represent the month
and day of the blossoming, being the darker areas the ones where blossoming is
expected earlier. Source: Wikimedia Commons.

While these datasets were constructed by relying, basically, on visual and
small-scale observations, during the last decades technological advances have
permitted a multiplication and diversification of the means to measure the
phenology of both plants and animals. Such methodologies range from citizen
science projects where participants collect and share individual observations
of phenology in their local area[17], to the use of sophisticated remote
sensing tools based on satellite data[18]. A particularly ingenious case among
these novel approaches is the work by[19], in which they process the images
recorded by public cameras connected to the Internet -primarily related to
traffic surveillance or national parks- to gather information on vegetation
phenology across North America, exploiting a free source of glances of trees
and plants unintendedly caught by the cameras. Nurtured by these
methodological advances as well as spurred by the investigations on how
climate change influences phenology[20, 21], nowadays we face an undeniable
increase in the quantity and the quality of documented phenology.

Paradoxically, though, this copious amount of phenological data comes at
relative help when trying to better understand mutualism over time. This is
due to the fact that public empirical information on both the network of
interactions and the timing of their mutualistic activity is rather scarce.
Consequently, previous studies examining the effect of phenology on
mutualistic communities either focus on a very reduced number of highly-
resolved networks[7, 22, 10], either they succeed in keeping the big numbers
of phenological data at the expense of roughly approximating the patterns of
mutualistic relationships –that is, by neglecting the real complex networks of
interactions[23]. A third approach, yet, to address the lack of data, is to
explore instead synthetic models either for phenology or for the network of
interactions[24, 8, 25]. In this work we chose to use the latter in order to
overcome the problem of data scarcity.

###  1.2 Ecological and evolutionary determinants of phenology

Although phenology certainly plays a part in different types of ecological
mutualism, in this work we will focus on its implications for plant-pollinator
communities, both because it is a paradigmatic example in which seasonality is
strongly marked and also because, despite the aforementioned data limitations,
the phenology of plants and pollinators is, in general, better documented than
that of other mutualistic species[26].

The study of the phenology involves all temporal aspects of species’ life
cycle, regarding on the one hand the timing of their various stages of
development, i.e. egg/seed, larvae, adult, etc.; and on the other hand the
onset and duration of different biological processes such as flowering,
germination, pollination, leaf shedding, etc. Such timing is determined by a
myriad of factors, similar to what occurs with mutualistic interactions in
general[27]. Indeed, both biotic and abiotic forces shape the phenology of
species, which, in addition, is thought to be subject to evolutionary
change[28, 29]. To complicate things further, several sources of intra-species
phenological variability occur simultaneously: inter-annual[7, 30] –that is,
among different seasons–, geographical[31] –for the same species but on
diverse sites– and last but not least, individual, i.e. among different
individuals of the same population, even those coexisting on the same site and
at the same season[29].

All in all, from a statistical viewpoint a few general patterns have been
identified. In particular, we will focus on two fundamental phenological
quantities of plant-pollinator systems that are well characterized: the
starting dates, that is, the time at which the flowering –in the case of
plants– or the pollination –for animals– begins, and secondly the so-called
periods, that is, the length or duration of this active state, during which
mutualistic interactions are possible. Certainly, this selection of
phenological indicators is far from being exhaustive, and other works have
modeled the different stages of plants and pollinators in more realistic
detail[10]. In our approach we will approximate this complex landscape of
temporal variability by focusing only on the adult phase of the species, and
moreover, reducing the whole possible set of biological processes and stages
to a couple of states: active, namely when the species could potentially hold
a mutualistic relationship –i.e. pollinate or be pollinated–, or inactive.
This sort of simplification is not novel and other works have adopted
analogous approaches[20, 8, 13]. In Fig. 2 we depict this schematic
representation of a hypothetical plant and its pollinator.

Figure 2: Schema of the phenology of a plant and a pollinator, represented by
their active states. Each species is defined by its period and its starting
date, called here $t_{s,f}$ for the flowering plant and $t_{s,p}$ for the
pollinator. The portion of time during which the two species overlap and can
actually interact is highlighted by the gray square.

Regarding the periods of activity, in general, both plants’ and pollinators’
periods tend to follow right-skewed distributions, with a small number of
species exhibiting a long phenophase while a large number of species are
active only during a short time[32, 24]. Particularly, Kallimanis et al.[24]
examined the statistical distribution of period lengths in a Mediterranean
scrub community observed over four years by Petanidou et al.[33] and concluded
that the distribution of pollinator’s periods could be fitted by a decreasing
exponential, while plants’ followed a lognormal distribution. On the other
hand, a lognormal shape for both plants and pollinators’ periods and a normal
distribution for plants in one of the years of their observation has been
identified in a plant-pollinator system in the Arctic[7].

In what concerns the onset of flowering and pollinating activity, a
complementary measure that is often used is the middle date. In this sense,
evidence suggests that middle dates tend to be relatively synchronized[28,
34]. For example, both in Ref.[7] and[32] it was observed that most active
periods temporally coexist at a peak, probably due to a similar reaction to a
common set of physicochemical stimuli such as temperature, photo-period,
humidity, etc[28]. In statistical terms, this would correspond to a scenario
where the middle dates are relatively clustered, as modeled for instance in
Ref.[24] using a normal distribution. At the same time, genetic factors seem
to play an important role in determining the timing of flowering or
pollinating activity, which contributes to explaining the heterogeneity of
starting dates. On the other hand, it has also been hypothesized that species
may spread along the season in order to minimize competition. Although the
evidence supporting this hypothesis is controversial[28], some authors have
proposed that such minimization could still occur within a temporal range
determined by genetic and environmental constraints[35].

This short introduction illustrates the complexity of the temporal dimension
of plant-pollinator communities. Indeed, temporal variability not only appears
at different scales, from days to decades and from individuals to species, but
it is furthermore regulated by a multiplicity of factors. Moreover, the
relevance of phenology goes beyond the individual characterization of species’
traits, since it regulates as well the ecological interactions within the
community, affecting both their occurrence and intensity. Understanding how
this may modify the structure of ecological interactions beyond the limits of
the aggregated paradigm is crucial if we aim to produce a refined depiction of
natural systems. In order to do so, in the next section we define a framework
to model the effects of phenology on mutualistic networks using real data,
whose outcomes we will explore in the results section.

##  2 Definitions and Models

Along this section we introduce the tools and models to explore the effects of
phenology on a plant-pollinator community: first, we define the mutualistic
and competitive networks of interactions; second, we describe how we may
account for empirical phenology using overlaps and a coarse-grained temporal
representation of networks; finally, we present a set of basic toy-models to
produce synthetic configurations of phenology using both mechanistic and
statistical principles.

###  2.1 Mutualistic and competitive network

In order to assess the effects of phenology on the plant-pollinator community,
we will not only consider the mutualistic relationship among different species
but also indirect interactions like competition for shared mutualistic
resources[36], which naturally emerge among species of the same kind. The
negative effects of these antagonistic interactions are known to coexist with
mutualism, leading to a trade-off between costs and benefits[6, 37].

Therefore, we start by defining the biadjacency matrix of mutualistic
interactions $B_{i,k}$ of size $N_{P}\times N_{A}$, where $N_{P}$ is the
number of plants and $N_{A}$ is the number of animals, such that:

$\displaystyle\small\mbox{ if }B_{i,k}=1,\mbox{ then plant }i\mbox{ is
pollinated by animal }k;$ $\displaystyle\mbox{ if }B_{i,k}=0,\mbox{ then
species }i\mbox{ and }k\mbox{ do not share any mutualistic interaction.}$ (1)

In order to be able to calculate the consequence of temporal variability also
among competitors, we need to infer the network of competitive interactions
based on shared mutualistic resources. We do so following Ref.[6], namely,
projecting the empirical biadjacency matrix of mutualistic interactions
$B_{i,k}$ into the subspace of intra-guild interactions such that:

$\displaystyle\small\mbox{ if }B_{i,k}B_{j,k}=1,\mbox{ then plants }i\mbox{
and }j\mbox{ compete for the pollinating services of animal }k;$
$\displaystyle\noindent\mbox{ if }B_{i,k}B_{j,k}=0,\mbox{ then }i\mbox{ and
}j\mbox{ do not compete for the resources of pollinator }k,$ (2)

which leads to a competitive network between all pairs of plant species $i$
and $j$. Analogously, we can define the competitive interactions among two
pollinator species $k$ and $l$ as follows:

$\displaystyle\small\mbox{ if }B_{i,k}B_{i,l}=1,\mbox{ then pollinators
}k\mbox{ and }l\mbox{ compete for the flower resources of plant }i;$
$\displaystyle\mbox{ if }B_{i,k}B_{i,l}=0,\mbox{ then }k\mbox{ and }l\mbox{ do
not compete for the resources of plant }i.$

Altogether, this provides a bilayered depiction of the ecological community,
such that the observed network of pollinating contacts not only mediates the
explicit mutualistic interactions, but also the implicit, competitive
relationships among species of the same guild[6].

###  2.2 Empirical overlaps

The most straightforward consequence of introducing phenology into an
aggregated network is the modification of the amount of temporal coexistence
–the so-called overlap– among species. Indeed, we find a continuum of possible
scenarios ranging from full concurrence to absence of overlap, as shown in
Fig. 3. These overlaps concern either two species in the mutualistic case or
three species –two species of the same kind and the shared resource– when
considering the competition.

Figure 3: Diversity of phenological configurations giving rise to different
types of mutualistic and competitive overlap. The panels on the left represent
the mutualistic overlap between the pollinator species $k$ and the plant
species $i$, from the viewpoint of the pollinator species $k$. In a) there is
full mutualistic overlap, in b) there is only partial overlap, and in c) there
is no overlap at all. Right Panels depict the competitive overlap between
three species, particularly two plants, $j$ and $i$, competing for a shared
pollinator species $k$. In this case, the overlap is referred to the first
species, $j$. In d) the competitive overlap is full, in e) there is only
partial overlap, and in f) there is no overlap at all between the competitors.

Within this framework, we can introduce a set of phenological coefficients
$\\{\Phi\\}$ and $\\{\Omega\\}$ to quantify the effect of the phenological
overlap on, respectively, the mutualistic and the competitive interactions. In
the particular case of a plant species $i$ interacting mutualistically with a
pollinator $k$ and competitively with another plant $j$, we define these
coefficients as follows:

$\Phi^{P}_{ik}=\frac{\tau_{ik}}{\tau_{i}},$ (4)

$\Omega^{P}_{ijk}=\frac{\tau_{ijk}}{\tau_{i}},$ (5)

where $\tau_{i}$ stands for the period of flowering of the plant species $i$,
$\tau_{ik}$ represents the temporal overlap between the period of the plant
species $i$ and that of its pollinator $k$ (see Fig. 3 a-c for a graphical
representation), and finally $\tau_{ijk}$ stands for the overlap between
plants species $i$ and $j$ and their pollinator $k$ (see Fig. 3 d-g for an
example). An analogous set of coefficients can be drawn for the pollinators.

Note that, as can be seen in Eqs. 4-5, the overlaps are pondered by the period
of activity of the species –in the example above, the plant $i$–, which means
that each coefficient is always referred to a certain plant or animal species.
This particular normalization implies, moreover, that the effect of the
temporal overlap is non-symmetric among the interacting partners and hence, in
general, $\Phi^{P}_{ik}\neq\Phi^{A}_{ki}$.

###  2.3 Coarse-grained temporal sequence of networks

In order to better understand how the perceived structure of the network
changes if we take the empirical phenology into account, we construct as well
a discrete temporal sequence of networks. Here, each element corresponds to a
different snapshot of the system, taken at a different moment of the season.
In particular, given that the degree of detail of the empirical phenologies
under study is narrowed to days, we construct a set of daily networks. That
is, each network represents the mutualistic interactions observed among plants
and pollinators on a given day.

To construct this sequence of networks, we use the empirical information on
the starting and ending dates of activity of plants and pollinators, then
remove from the daily network certain mutualistic links among species whenever
they do not coincide in time – i.e., when their phenological overlap, as
defined in Eq. 4, is zero. Additionally, we remove inactive species as well as
active species with no interactions –i.e. zero degree.

Figure 4: Schematic representation of the coarse-grained description of
temporal variation. Each network in the grey scale corresponds to a different
day of the season, where certain species and interactions might be absent.
Summing up these interactions produces the aggregated network, depicted in
color, where the information on the turnover of species and interactions is
lost. Despite this aggregation process is indeed quite simple, a pertinent
question is whether the structure and dynamics of the temporal snapshots are
comparable to those of the aggregated network.

On the whole, this provides a coarse-grained description of the community as
shown in Fig. 4, that, even if it is not a pure temporal network, it does
account to a great extent for the intra-annual temporal variation, moving
beyond the static network formalism.

###  2.4 Synthetic toy-models

In order to evaluate how our results depend on the details of the empirical
dataset, we develop a group of toy models that, under different assumptions,
produce synthetic configurations of phenology compatible with a given
mutualistic network. As aforementioned, in our framework, the phenology of
each species is defined by two elements: the starting date and the period of
activity. Here, we will focus on determining the starting dates while
preserving the empirical periods. Therefore, the starting dates will be
treated as parameters to be set by the synthetic model, under the condition
that mutualistic partners share a non-zero mutualistic overlap. This
requirement ensures that the binary mutualistic network is kept equal to the
observed one, constituting a formalism that focuses explicitly on modeling the
phenology.

Before entering into further details about the different synthetic models, let
us briefly discuss the particularities and consequences of preserving the
network of interactions. As aforementioned, we will construct a model under
the requirement that the phenological overlap among mutualistic partners is
comprised in the range $(0,1]$, excluding the zero in order to ensure that the
two species still interact. This means that, strictly speaking, we will only
keep the binary mutualistic network. In other words, the general idea consists
in constructing a synthetic phenology that is compatible with the given
network treated in binary terms.

In this line of thought, the constraint that mutualistic partners share non-
zero phenological overlap may be written as an inequality for each species’
starting date. For instance, let us take a plant species $i$ with period
$p_{i}^{P}$ and unknown starting date $t_{0,i}^{P}$, which interacts with an
animal species $j$ with period $p_{j}^{A}$ and an unknown starting date
$t_{0,j}^{A}$. In order to ensure that their periods overlap, it must be
fulfilled that:

$t_{0,j}^{A}-p_{i}^{P}+1\,\leq\,t_{0,i}^{P}\,\leq\,t_{0,j}^{A}+p_{j}^{A}-1\,,$
(6)

where we have assumed a daily coarse-graining, such that the minimal
overlapping unit is one day. In order to ensure that all interacting species
share phenological overlap, we may write an inequality analog to Eq. 6 per
each link in the bipartite network. Then, if we take the starting dates as
unknown variables, we can view the problem of constraining the starting dates
as a multivariate system of linear inequalities composed by $N$ unknown
variables and $L$ inequalities (where $N$ and $L$ are respectively the number
of nodes and the number of links of the network). Note that both
$\\{t_{0}^{P}\\}$ and $\\{t_{0}^{A}\\}$ are unknown, which means that the
system of inequalities is coupled.

Furthermore, in order to ensure that the resulting distribution is realistic,
we may introduce an additional set of constraints. In particular, we impose a
lower and upper bound for the starting dates of activity, which ensures that
the community exhibits some seasonality, a characteristic that has been
empirically observed[28, 24]. If we represent by $l$ the lower bound and by
$u$ the upper bound delimiting the season, then the condition on the starting
date of an arbitrary species $i$ can be written as:

$l\,\leq\,t_{0,i}\,\leq\,u-p_{i}+1.$ (7)

Overall, Eqs. 6-7 represent a coupled system of linear inequalities, where the
parameters to be determined are the starting dates. Although a problem of this
kind might not generally have a feasible solution, in our particular case the
existence of some trivial solutions –e.g. the same starting dates for all
species– warrants that the system is consistent.

Within this framework, we present a group of approaches that, under different
assumptions, permit determining the set of synthetic starting dates while
fulfilling the conditions in Eqs. 6-7. In particular, we explore four possible
scenarios, corresponding to different fundamental hypotheses about the main
factor governing the species’ phenology, namely: (i) species’ periods tend to
be synchronized, (ii) species seek to organize their activity as to minimize
their intra-guild competition for mutualistic resources, (iii) species periods
get maximally dispersed along the season, (iv) in the absence of further
information, the starting dates can be determined using a minimal model of
maximum entropy.

As can be observed, each of these assumptions relies on a single dominant
factor. Indeed, hypothesis (i) and (ii) are mechanistic, in the sense that
they place the emphasis on environmental and ecological drivers of the
phenology. Instead, the scenarios (iii) and (iv) function as non-mechanistic
models based on statistical properties of the distribution of starting times,
regarding either their dispersion or their entropy. Of course, since biotic
and abiotic forces do not act separately, one could construct a more complex
model where various of these factors operate simultaneously. In this sense,
these synthetic models constitute a first-order approach, which helps to
understand the consequences of each of these hypotheses concerning the
phenology. In what follows, we define each of these models one by one and
describe how we numerically implemented them.

####  2.4.1 Synchronized phenology

This model is based on the assumption that the level of synchronization in the
community is high, hence most species’ activity temporally coincides with the
peak of the season. As we reviewed in the introduction, empirical evidence
suggests that seasonal systems usually exhibit non-negligible synchronization,
in the sense that the flowering and pollination of species tend to occur
simultaneously. Several hypotheses have been proposed to account for this
pattern, from the presence of abiotic pressures such as the temperature, the
precipitation, or the photoperiod –disregarding genetic differences or biotic
constraints–[29], to the existence of facilitative interaction among
concurrent species, which obtain an ecological benefit of overlapping their
periods of activity. All in all, it is still not clear how the mentioned
factors interplay among themselves and with other possible forces to give rise
to the observed heterogeneity in starting dates[29].

In order to attempt to reproduce this pattern, we construct a model where
species’ periods are approximately centered. To do so, we set the middle dates
around the peak of the season and then slightly perturb them, following the
numerical procedure detailed in the Supplementary Material. This leads to
phenological configurations where species are considerably –but not perfectly–
synchronized. Indeed, introducing this imperfection is crucial for two
reasons: first, to account for the aforementioned heterogeneity in starting
and middle dates of activity; second, because a perfectly synchronized system
is trivial, in the sense that all species fully overlap, recovering hence the
aggregated case.

####  2.4.2 Minimization of the intra-guild competition

This second model is built upon the hypothesis that the periods of the species
are located along the season so as to avoid, as far as possible, overlapping
with other species that exploit a common mutualistic resource –i.e.
pollinating services in the case of plants or flower availability for
pollinators. As discussed above, the evidence for this type of mechanism is
controversial, and, when it does occur, its effect is most probably
subordinated to other genetic and environmental factors. Nevertheless, we will
still explore the consequences of this assumption, given that it is the
natural counterexample of the synchronized scenario.

In order to construct a phenological configuration that minimizes the
competition among species, we implement a global search algorithm that
optimizes an objective function under the constraints in Eqs. 6-7 (see the
Supplementary Material for the numerical details). In this particular model,
the objective function corresponds to the total competitive phenological
overlap among species.

####  2.4.3 Maximization of the variance

We now turn our attention to models that focus directly on the statistical
properties of the system’s phenology. As a first approach, we propose a model
based on maximizing the dispersal of the periods along the season. In
particular, we maximize the variance of the middle dates of activity. Indeed,
if we define the middle time $t_{M}$ of a species $i$ as:

$t_{M,i}\,=\,t_{0,i}\,+\,\frac{p_{i}-1}{2}\,,$ (8)

then its variance Var$(t_{M})$ is an approximate measure of the dispersion of
the activity of both plants and pollinators along the season. Once again, the
maximization of this quantity is done under the condition that mutualistic
interactions present in the aggregated network still occur. The numerical
procedure to do so is explained in the Supplementary Material.

####  2.4.4 Maximization of the entropy

Following a similar line of thought, we propose a second model to produce
phenological configurations that is based on their statistical properties. In
particular, we refine the former approach by characterizing the distribution
with a more sophisticated and robust indicator, namely, entropy. This is a
plausible alternative to maximizing the variance which may especially
outperform the latter when the distribution is multimodal.

Let us define the entropy within the constructed framework for describing
phenology. We start by defining a random variable $x_{i}$ as the number of
species whose middle date of activity $t_{m}$ takes place at a given date
$t_{i}$. Then, its corresponding probability $p(x_{i})$ will be:

$p(x_{i})\,=\,\frac{\mbox{number of species with }t_{m}=t_{i}}{\mbox{total
number of species}}.$ (9)

Given this random variable, it is possible to define its Shannon-Gibbs entropy
as follows:

$S\,=\,-\sum_{i}\,p(x_{i})\,\mbox{ln }p(x_{i}).$ (10)

It is straightforward to see that in a scenario of perfect synchronization,
the middle dates of activity are the same for all species. That is, if we name
such a date by the index $k$, then $p(x_{k})=1$, and the information entropy
is zero.

Maximizing the entropy defined in Eq. 10 under the constraints imposed by the
network of interactions results in a maximum entropy model, that randomizes
the distribution of periods of the species yet using only the minimal amount
of necessary information. The particular details of how this model was
implemented can be found in the Supplementary Material.

##  3 Datasets

In order to assess the consequences of phenology on the perceived network
structure, we analyzed two public datasets [13, 14] which contain observations
of mutualistic interactions in plant-pollinator systems together with their
phenology. We detail the peculiarities of each dataset below.

###  3.1 The Illinois dataset

This dataset [13] corresponds to a set of woodland sites in Carlinville,
Illinois (USA), observed during the springs of 2009 and 2010 from March to
May. The observed network of mutualistic interactions contained originally 26
spring-blooming herbaceous plants and 54 bee species. However, 2 plant species
did not exhibit any interaction in the empirical network, so we obtained
eventually a $24\times 54$ network.

The dataset also included the observed phenology, specifically the starting
and ending dates of the active stages of both plants and pollinators. Using
those, we extracted as well their periods of activity and the middle dates.

###  3.2 The Lesvos Island dataset

This dataset is based on a two-year study conducted in Aglios Stefanos, Lesvos
Island (Greece) during the springs of 2011 and 2012 from April to July[14]. It
includes information on both the network of mutualistic interactions and the
phenology. Aggregated over the two years, the resulting network is composed of
41 plant species and 168 pollinators.

Remarkably, a large portion of the pollinator species is present in only one
of the seasons. In particular, out of the total 168 pollinator species
observed over the two years, only 67 species ($39.9\%$) were persistently
found in both seasons, while 53 species ($31.5\%$) were observed just during
the first year, and 48 species ($28.6\%$) just during the second year of
observation. This implies that not only there is an important turnover of
interactions among the two consecutive years, but also there is a considerable
species turnover, as suggested as well by previous studies[7, 11].

In order to obtain the network of interactions specific to each year, we
removed those interactions present in the aggregated network whose periods of
activity do not overlap according to the information concerning phenology. In
such cases, it is clear that the interaction was not possible due to the lack
of phenological overlap and hence we can remove it from the yearly network.
After doing so for each season, we removed as well those species that no
longer hold mutualistic interactions –i.e. they have degree zero. As a result,
both year-specific network sizes’ are reduced to 34 plant species and 113
pollinators in the first year and 38 plant species and 104 pollinators in the
second year of observations.

##  4 Results

###  4.1 Characterization of the empirical datasets

We start by examining the distribution of the duration of the activity of
plants and pollinators. In particular, we fit a variety of functional forms on
the cumulative distributions using a maximum likelihood estimation. Then, we
test the quality of the fit by performing a Kolmogorov-Smirnov test by
bootstrapping, as explained in detail in the Supplementary Material.

In detail, we test three different functional forms: lognormal, exponential,
and beta function. As aforementioned, the first two distribution types have
been claimed to correctly describe some empirical observations, while we
introduce the latter as a generic fitting. The results of the p-value for the
KS-test are reported in Table 1, where we highlight in bold the good-quality
fits. As can be seen, the Illinois dataset is better fitted by the beta
function, while in the Lesvos Island dataset, we can fit a lognormal to the
distribution of plants’ periods but we found no unique fit for the
pollinators. In general terms, the distribution of periods of the Illinois
dataset is less heterogeneous, while the Lesvos Island distributions are
clearly righter-skewed. This is especially notorious for the periods of the
pollinators, among which we find a large proportion of species with very short
periods.

Table 1: Results of fitting different functional forms, disentangled into plants and pollinators. We show the p-value of the fit, and highlight in bold the good quality fits, that is, those which do not significantly differ from the estimated distribution (p-value $>0.05$). The p-values are obtained by performing a Kolmogorov-Smirnov test between the fitted distribution and the empirical sample, then comparing it to the corresponding K-S distribution sampled by bootstrap as explained in the Supplementary Material. Dataset | Lognormal fit | Exponential fit | Beta function fit
---|---|---|---
| Plants | Pollinators | Plants | Pollinators | Plants | Pollinators
Illinois | 0.021 | 0.018 | 0.001 | 0.001 | 0.63 | 0.18
Lesvos Island 1st year | 0.67 | 0.014 | 0.029 | 0.001 | 0.41 | 0.001
Lesvos Island 2nd year | 0.50 | 0.001 | 0.018 | 0.001 | 0.02 | 0.001

###  4.2 Temporal variation of the structure

Now that we have characterized the empirical phenology of the two datasets, we
turn our attention to the question: how does the perceived structure of the
ecological network change if we take phenological information into account? In
particular, we measure a set of fundamental structural features for the
different coarse-grained networks of the daily sequence: the number of active
nodes, the number of active links, and the maximum degree. We do so both for
the mutualistic and the competitive network. Fig. 5 provides a summary of the
variation of these properties over time, for each of the empirical networks in
our dataset.

Figure 5: Phenological configurations and temporal evolution of three
fundamental structural features. On the top, a schema of the temporal
arrangement of the periods of activity of plant species (green) and
pollinators (purple) along the season. Below, is the temporal variation of
three basic structural parameters, measured on a daily basis: number of active
nodes, active links, and maximum degree. Each column corresponds to a
different dataset as indicated on the top of the figure.

To begin with, our results show that in both cases the temporal measures are
significantly smaller than their aggregated counterparts. However, we observe
a non-negligible quantitative difference among datasets: in the Illinois
dataset, the maximum of the peak represents around the 75-80$\%$ of the
aggregated counterpart, whereas in both instances of the Lesvos Island
dataset, this quantity descends to less than 30-40$\%$. Additionally, the
comparison of the results for the two systems shows that the variation of
their network structural properties exhibits a strikingly different pattern
over the season. Indeed, while in the Illinois dataset, the three quantities
show a marked peak approximately at the middle of the season, in the Lesvos
Island dataset they stay relatively stable, showing an almost flat shape.

###  4.3 Comparison with synthetic models

We now proceed to test the performance of the different generative models of
synthetic phenologies proposed above. To do so, we will adopt as input the
observed network of interactions –which is to be preserved–, and the empirical
distribution of periods. The latter will allow us to focus on the quality of
the methodology to produce the starting dates alone, without introducing
additional sources of noise.

In order to examine the adequacy of each model to reproduce the empirical
phenological patterns, we start by extracting the synthetic starting dates
corresponding to each of our three datasets (see the Supplementary Material
for the explanation of the numerical methods). Then, we characterize the
synthetic phenologies with the methods and definitions described in section 2.
That is, we calculate, on the one hand, the distribution of phenological
coefficients defined in Eqs. 4-5, and, on the other hand, we couple the
phenology to the network structure and compute the change of some basic
structural features along time. Along the rest of this section, we will
explain and discuss the results of comparing these synthetic measures with the
empirical observations.

####  4.3.1 Distribution of overlaps

We first study the distribution of the phenological coefficients defined in
Eqs. 4-5. In order to compare the empirical and the synthetic distributions,
we perform two different statistical measures: (a) a two-sample Kolmogorov-
Smirnov test and (b) the calculation of the Kullback-Leibler divergence (see
the Supplementary Material for the numerical details).

Using (a), we analyzed separately the competitive and mutualistic overlaps on
the one hand, and the distributions for plants and for pollinators on the
other hand, as summarized in Table 2. In particular, for the Illinois dataset,
the null hypothesis that both samples come from the same distribution is
rejected for the models which minimize the competition and maximize the
variance, which is hence incompatible with the empirical data. The model with
centered times and maximum entropy is generally compatible with the observed
phenology, except for the distribution of coefficients associated with the
competition among pollinators. Instead, for the Lesvos Island dataset, the
four models are compatible with the data.

Table 2: Results of the K-S two sample test among the synthetic distribution and the empirical distribution. We show the KS distance between the two samples, and highlight in bold the cases in which the null hypothesis is rejected, that is, those where the two samples do significantly differ (p-value $<0.05$). The smaller the KS distance, the smaller the discrepancy between the two samples. The numerical implementation to calculate these distances is described on the Supplementary Material. Model /Dataset | Mutualistic overlap | Competitive overlap
---|---|---
Synchronized | Plants | Pollinators | Plants | Pollinators
Illinois |  0.052 |  0.051 |  0.008 |  0.021
Lesvos Island 1st year |  0.015 |  0.010 |  0.002 |  0.002
Lesvos Island 2nd year |  0.009 | 0.009 |  0.001 | 0.0017
Minimum competition | Plants | Pollinators | Plants | Pollinators
Illinois |  0.100 | 0.103 |  0.036 |  0.031
Lesvos Island 1st year |  0.012 |  0.013 |  0.002 |  0.0015
Lesvos Island 2nd year |  0.007 | 0.008 |  0.001 | 0.001
Maximum variance | Plants | Pollinators | Plants | Pollinators
Illinois |  0.066 | 0.069 |  0.027 |  0.028
Lesvos Island 1st year |  0.015 |  0.014 |  0.004 |  0.003
Lesvos Island 2nd year |  0.011 | 0.011 |  0.003 | 0.002
Maximum entropy | Plants | Pollinators | Plants | Pollinators
Illinois |  0.040 | 0.042 |  0.010 |  0.015
Lesvos Island 1st year |  0.009 |  0.007 |  0.001 |  0.001
Lesvos Island 2nd year |  0.006 | 0.007 |  0.001 | 0.0016
Table 3: Results of the K-L divergence among the synthetic distributions and the empirical distributions of overlap, either among mutualistic or competitive partners. The smaller the KL-divergence, the smaller the difference between the two samples. Further details on the definition and numerical implementation of this measure can be found in the Supplementary Material. Dataset / Model |  $D_{KL}$ mutualism | $D_{KL}$ competition |
---|---|---|---
Illinois | Plants | Pollinators | Plants | Pollinators | Average
Synchronized | 0.035 | 0.038 | 0.002 | 0.014 | 0.022
Minimum competition | 0.095 | 0.101 | 0.026 | 0.033 | 0.064
Maximum variance | 0.042 | 0.056 | 0.010 | 0.019 | 0.032
Maximum entropy | 0.022 | 0.032 | 0.003 | 0.008 | 0.016
Lesvos Island 1st year | Plants | Pollinators | Plants | Pollinators | Average
Synchronized | 0.0116 | 0.0053 | 0.0014 | 0.0010 | 0.0048
Minimum competition | 0.0050 | 0.0073 | 0.0009 | 0.0004 | 0.0034
Maximum variance | 0.0111 | 0.0047 | 0.0032 | 0.0008 | 0.0049
Maximum entropy | 0.0055 | 0.0051 | 0.0005 | 0.0005 | 0.0029
Lesvos Island 2nd year | Plants | Pollinators | Plants | Pollinators | Average
Synchronized | 0.0048 | 0.0036 | 0.0008 | 0.0008 | 0.0025
Minimum competition | 0.0018 | 0.0035 | 0.0005 | 0.0005 | 0.0015
Maximum variance | 0.0080 | 0.0036 | 0.0019 | 0.0005 | 0.0035
Maximum entropy | 0.0019 | 0.0035 | 0.0005 | 0.0007 | 0.0016

In order to differentiate which model better reproduces the empirical
phenology in the cases in which, according to the KS test, various models are
statistically compatible, we carry out (b), that is, the calculation of the
KL-divergence. Indeed, the KL distance measures the amount of information lost
when we approximate the empirical distribution by the synthetic one, and
therefore it can be used as a statistical criterion for model selection. In
Table 3 we can observe that the model based on maximum entropy generally
provides the best approximation, except for the second year of the Lesvos
Island dataset, in which it is slightly outperformed by the model based on
minimizing the competition. It is also interesting to observe that the second-
best model is not the same among datasets: while the second-best model for the
Illinois dataset is the one with synchronized phenology, in the Lesvos Island
dataset it is the model which minimizes the competition.

####  4.3.2 Structural features

Secondly, in order to complement the information provided by depicting the
phenological coefficients, we calculate the daily networks resulting from
considering the phenology. This again produces a sequence of synthetic
networks, each corresponding to a given date, as represented in Fig. 4.
Repeating the procedure detailed in section 3.2. now for each synthetic model,
we calculated a few basic structural quantities as a function of time: the
number of active nodes -size-, the number of active links, and the maximum
degree.

With this information at hand, we can perform a comparison between the product
of each synthetic model and the empirical observations, specifically
addressing the change in fundamental structural properties over time. To
quantify these correlations, we calculated the corresponding Pearson
coefficient for different values of time lag between the distributions. By
doing so, we are taking into account the possibility that the null model
reproduces well the empirical scenario but in a delayed or advanced time. That
is, we are relaxing the assumption of perfect matching between the empirical
and the synthetic starting dates, focusing instead on the relative position
among the starting dates of different species. In table 4 we show the maximum
value of the Pearson correlation calculated following this procedure
(technical details can be found in the Supplementary Material).

The study of table 4 reveals that the analysis of the network structure leads
to similar results to the aforementioned analysis of the phenological overlap.
Indeed, the better fitting models for the Illinois dataset are the ones based
on synchronized phenology and maximum entropy, while for the Lesvos Island
dataset, the maximal correlation is found for the model which minimizes intra-
guild competition and, again, the one maximizing the entropy.

Table 4: Maximum Pearson correlation among the empirical and the null expectation of the structural properties along time. Dataset / Property | Maximum Pearson coefficient for each null model
---|---
Illinois | Synchronized | Max Variance | Min competition | Max entropy
Maximum degree | 0.98 | 0.88 | 0.93 | 0.98
Number of links | 0.98 | 0.87 | 0.78 | 0.98
Size | 0.99 | 0.79 | 0.92 | 0.98
Average | 0.99 | 0.85 | 0.88 | 0.98
Lesvos Island 1st year | Synchronized | Max Variance | Min competition | Max entropy
Maximum degree | 0.78 | 0.76 | 0.91 | 0.91
Number of links | 0.74 | 0.71 | 0.89 | 0.83
Size | 0.81 | 0.58 | 0.89 | 0.85
Average | 0.78 | 0.69 | 0.90 | 0.86
Lesvos Island 2nd year | Synchronized | Max Variance | Min competition | Max entropy
Maximum degree | 0.90 | 0.67 | 0.87 | 0.93
Number of links | 0.85 | 0.78 | 0.89 | 0.85
Size | 0.89 | 0.63 | 0.90 | 0.88
Average | 0.88 | 0.69 | 0.88 | 0.88

Overall, these multiple statistical tests reveal at least two general
conclusions. On the one hand, we find that each dataset is better described by
a distinct mechanistic assumption, i.e. the synchronization of the species’
phenologies in the Illinois dataset and the minimization of the competition in
the Lesvos Island dataset. This disparity seems to reflect the substantial
divergences between the two communities’ phenologies, which we described in
detail in section 3.1. On the other hand, the temporal networks’ features are
accurately reproduced by the statistical model that maximizes the entropy
associated with the distribution of the middle dates of activity, despite the
considerable differences between the two datasets. What this latter finding
suggests is that the observed network of interactions together with the
periods of each species provide sufficient information to reproduce, fairly
closely, the observed temporal patterns of activity. This is particularly
interesting given the mentioned differences among datasets, which provide, to
a certain extent, a warranty of generality despite the limited data at our
disposal. This finding does not imply that, forcedly, the starting dates of
activity are set at random, but that the information enclosed in the
corresponding network of interactions and the species’ periods may be
sufficient to reproduce well the main characteristics of the community’s
phenology.

##  5 Conclusions

Along this paper, we have addressed the challenge of moving beyond the
aggregated paradigm by incorporating some of the temporal variability that
plant-pollinator communities exhibit throughout the year. In particular, we
have concentrated on seasonal ecosystems, assessing the daily change in both
inter-guild, mutualistic relationships and intra-guild, competitive
interactions.

The analyses carried out by using both synthetic models and empirical data
have revealed, in the first place, that non-trivial information is lost when
portraying the real network of interactions by a static representation, that
is, neglecting its temporal dimension. Indeed, we have observed that the
consequences of introducing the empirical phenology into the network formalism
are system-dependent, and hence no general pattern can be expected a priori.
Moreover, the process of aggregation not only disregards the richness of
temporal variability, but it also tends to overestimate the value of the main
fundamental structural features, as had been remarked as well by[9]. This is
especially relevant given that structural features like the degree or the
connectivity are consistently used to characterize, respectively, the relative
vulnerability of species[38] or the community stability[5].

In light of the limitations of the aggregated paradigm, and driven by the
scarcity of available datasets, we proposed a group of models to produce,
given a fixed mutualistic network, a compatible hypothetical phenology. The
comparison of these results with the empirical datasets revealed that the
soundness of certain mechanisms is, at least in the first-order approach we
have adopted here, specific to the particular system under consideration. We
have found that the purely statistical assumption of maximizing the entropy
associated with the distribution of middle dates performs generally well, as
we have tested in two dissimilar empirical examples. Importantly, the
remarkable performance of the maximum entropy hypothesis is partly explained
by the fact that, actually, the network of interactions is closely dependent
upon the starting dates, in the sense that the mutualistic contacts observed
corresponded, forcedly, to concurrent species –a condition that, indeed, we
impose to our models. Therefore, preserving the network of interactions is a
strong constraint, which could justify the general adequacy of this model.

In perspective, these synthetic models offer a methodological set that might
prove useful in different aspects. On the one hand, they can be exploited
merely as a group of realistic models to construct synthetic ensembles in the
absence of highly-resolved empirical data, in those cases where the main
driving forces of the phenology are known. On the other hand, they can also be
applied as null models that permit testing a variety of null hypotheses, from
the mechanistic forces shaping the phenology to the existence of temporal,
structural, or dynamical patterns. At this point, it is worth reminding as
well that ours is just a first approach to modeling the temporal variability
of networks. In particular, we considered the description level at which
species are active or inactive during a certain fraction of the season.
However, it could be possible to refine the scale of description to include
the weekly or even daily turnover of interactions, a sort of hyper-realistic
depiction of the temporal variability of the network that is gaining attention
during the recent years[22]. This is due, partly, to the technological
advances that permit monitoring phenology in great detail, and therefore it is
probable that in the coming years, we will find a rising number of this type
of studies.

In the bigger picture, the interest in moving towards a more realistic
portrayal of interacting systems is not exclusive to ecology, and indeed the
study of temporal complex networks has received great attention in recent
years[39]. How this change of paradigm will eventually challenge our
understanding of natural systems is something we are just now beginning to
explore.

## Data availability

The datasets analyzed in this paper were released by [13] and [14] and are
publicly available, respectively, at:

  * •

https://doi.org/10.5061/dryad.rp321

  * •

https://figshare.com/articles/dataset/Data_from_Disentangling
_the_role_of_floral_sensory_stimuli_to_pollination_networks_/5663455

## Acknowledgments

The authors would like to acknowledge and thank A. Kantsa for sharing some
valuable clarifications on the construction of the Lesvos Island dataset, as
well as additional information on the starting and final dates of the
pollinator’s activity. We acknowledge partial support from the Government of
Aragon, Spain, and “ERDF A way of making Europe” through grant E36-23R (FENOL)
and from Ministerio de Ciencia e Innovación, Agencia Española de Investigación
(MCIN/ AEI/10.13039/501100011033) Grant No. PID2020-115800GB-I00 to C.G.L. and
Y.M.

## References

  * [1] Bascompte J, Jordano P.  Mutualistic networks. vol. 70.  Princeton University Press; 2013.
  * [2] Heleno R, Garcia C, Jordano P, Traveset A, Gómez JM, Blüthgen N, et al.  Ecological networks: delving into the architecture of biodiversity.  Biology Letters. 2014.
  * [3] Jordano P.  El ciclo anual de los paseriformes frugívoros en el matorral mediterráneo del sur de España: importancia de su invernada y variaciones interanuales.  Ardeola. 1985;32(1):69-94.
  * [4] Bascompte J, Jordano P, Melián CJ, Olesen JM.  The nested assembly of plant–animal mutualistic networks.  Proceedings of the National Academy of Sciences. 2003;100(16):9383-7.
  * [5] Thébault E, Fontaine C.  Stability of ecological communities and the architecture of mutualistic and trophic networks.  Science. 2010;329(5993):853-6.
  * [6] Gracia-Lázaro C, Hernández L, Borge-Holthoefer J, Moreno Y.  The joint influence of competition and mutualism on the biodiversity of mutualistic ecosystems.  Scientific reports. 2018;8(1):9253.
  * [7] Olesen JM, Bascompte J, Elberling H, Jordano P.  Temporal dynamics in a pollination network.  Ecology. 2008;89(6):1573-82.
  * [8] Encinas-Viso F, Revilla TA, Etienne RS.  Phenology drives mutualistic network structure and diversity.  Ecology Letters. 2012;15(3):198-208.
  * [9] Sajjad A, Saeed S, Ali M, Khan FZA, Kwon YJ, Devoto M.  Effect of temporal data aggregation on the perceived structure of a quantitative plant–floral visitor network.  Entomological Research. 2017;47(6):380-7.
  * [10] Ramos-Jiliberto R, de Espanés PM, Franco-Cisterna M, Petanidou T, Vázquez DP.  Phenology determines the robustness of plant–pollinator networks.  Scientific reports. 2018;8(1):1-10.
  * [11] Chacoff NP, Resasco J, Vázquez DP.  Interaction frequency, network position, and the temporal persistence of interactions in a plant–pollinator network.  Ecology. 2018;99(1):21-8.
  * [12] Ings TC, Montoya JM, Bascompte J, Blüthgen N, Brown L, Dormann CF, et al.  Ecological networks–beyond food webs.  Journal of Animal Ecology. 2009;78(1):253-69.
  * [13] Burkle LA, Marlin JC, Knight TM.  Plant-pollinator interactions over 120 years: loss of species, co-occurrence, and function.  Science. 2013;339(6127):1611-5.
  * [14] Kantsa A, Raguso RA, Dyer AG, Olesen JM, Tscheulin T, Petanidou T.  Disentangling the role of floral sensory stimuli in pollination networks.  Nature communications. 2018;9(1):1-13.
  * [15] Aono Y, Kazui K.  Phenological data series of cherry tree flowering in Kyoto, Japan, and its application to reconstruction of springtime temperatures since the 9th century.  International Journal of Climatology: A Journal of the Royal Meteorological Society. 2008;28(7):905-14.
  * [16] Chuine I, Yiou P, Viovy N, Seguin B, Daux V, Ladurie ELR.  Grape ripening as a past climate indicator.  Nature. 2004;432(7015):289-90.
  * [17] Havens K, Henderson S.  Citizen science takes root.  American Scientist. 2013;101(5):378-85.
  * [18] Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JC, Gao F, et al.  Monitoring vegetation phenology using MODIS.  Remote sensing of environment. 2003;84(3):471-5.
  * [19] Graham EA, Riordan EC, Yuen EM, Estrin D, Rundel PW.  Public Internet-connected cameras used as a cross-continental ground-based plant phenology monitoring system.  Global Change Biology. 2010;16(11):3014-23.
  * [20] Memmott J, Craze PG, Waser NM, Price MV.  Global warming and the disruption of plant–pollinator interactions.  Ecology letters. 2007;10(8):710-7.
  * [21] Hegland SJ, Nielsen A, Lázaro A, Bjerknes AL, Totland Ø.  How does climate warming affect plant-pollinator interactions?  Ecology letters. 2009;12(2):184-95.
  * [22] CaraDonna PJ, Petry WK, Brennan RM, Cunningham JL, Bronstein JL, Waser NM, et al.  Interaction rewiring and the rapid turnover of plant–pollinator networks.  Ecology letters. 2017;20(3):385-94.
  * [23] Duchenne F, Thébault E, Michez D, Elias M, Drake M, Persson M, et al.  Phenological shifts alter the seasonal structure of pollinator assemblages in Europe.  Nature Ecology & Evolution. 2020;4(1):115-21.
  * [24] Kallimanis A, Petanidou T, Tzanopoulos J, Pantis J, Sgardelis S.  Do plant–pollinator interaction networks result from stochastic processes?  Ecological modelling. 2009;220(5):684-93.
  * [25] Duchenne F, Fontaine C, Teulière E, Thébault E.  Phenological traits foster persistence of mutualistic networks by promoting facilitation.  Ecology Letters. 2021.
  * [26] Rafferty NE, CaraDonna PJ, Bronstein JL.  Phenological shifts and the fate of mutualisms.  Oikos. 2015;124(1):14-21.
  * [27] Vázquez DP, Chacoff NP, Cagnolo L.  Evaluating multiple determinants of the structure of plant–animal mutualistic networks.  Ecology. 2009;90(8):2039-46.
  * [28] Rathcke B, Lacey EP.  Phenological patterns of terrestrial plants.  Annual review of ecology and systematics. 1985;16(1):179-214.
  * [29] Forrest J, Miller-Rushing AJ.  Toward a synthetic understanding of the role of phenology in ecology and evolution. 2010.
  * [30] Cirtwill AR, Roslin T, Rasmussen C, Olesen JM, Stouffer DB.  Between-year changes in community composition shape species’ roles in an Arctic plant–pollinator network.  Oikos. 2018;127(8):1163-76.
  * [31] Post E, Steinman BA, Mann ME.  Acceleration of phenological advance and warming with latitude over the past century.  Scientific reports. 2018;8(1):1-8.
  * [32] Bawa KS, Kang H, Grayum MH.  Relationships among time, frequency, and duration of flowering in tropical rain forest trees.  American Journal of Botany. 2003;90(6):877-87.
  * [33] Petanidou T, Ellis WN, Margaris NS, Vokou D.  Constraints on flowering phenology in a phryganic (East Mediterranean shrub) community.  American Journal of Botany. 1995;82(5):607-20.
  * [34] Tébar F, Gil L, Llorens L.  Flowering and fruiting phenology of a xerochamaephytic shrub community from the mountain of Mallorca (Balearic islands, Spain).  Plant Ecology. 2004;174(2):295-305.
  * [35] Kochmer JP, Handel SN.  Constraints and competition in the evolution of flowering phenology.  Ecological monographs. 1986;56(4):303-25.
  * [36] Jones EI, Bronstein JL, Ferrière R.  The fundamental role of competition in the ecology and evolution of mutualisms.  Annals of the New York Academy of Sciences. 2012;1256(1):66-88.
  * [37] Bronstein JL.  The costs of mutualism.  American Zoologist. 2001;41(4):825-39.
  * [38] Dakos V, Bascompte J.  Critical slowing down as early warning for the onset of collapse in mutualistic communities.  Proceedings of the National Academy of Sciences. 2014;111(49):17546-51.
  * [39] Holme P, Saramäki J.  Temporal networks.  Physics reports. 2012;519(3):97-125.
  * [40] Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods. 2020.

Supplementary Material

##  SM-1 Numerical implementation of synthetic
models

In this section, we explain the different numerical methods used to determine
the starting dates of species’ activity, on the basis of different hypotheses
and under the constraint that mutualistic partners share a minimum
phenological overlap.

First, we describe two numerical approaches to solve the problem of
maintaining the mutualistic overlap: (i) by shifting some given starting
dates, (ii) by optimizing a certain quantity within the constraints imposed by
the non-zero overlap. Secondly, we explain how we implemented these methods
for each particular synthetic model.

###  SM-1.1 Shift of starting dates

Here we introduce a method to shift a given set of initial times of activity,
with the condition that a minimal mutualistic overlap is preserved. The
periods of activity, that we will denote as $p^{P}_{i}$ for a plant $i$ and
$p^{A}_{k}$ for an animal $k$, are kept unchanged. We solely modify the timing
of the beginning of the activity of each species, where we note the initial
time by $t^{P}_{0,i}$ for each plant $i$ and $t^{A}_{0,k}$ for each animal
$k$.

The program takes as input the aforementioned information together with the
starting dates of activity, which, importantly, needs to already fulfill the
set of inequalities written in Eqs. 2.9 and 2.10 in the main text. This can be
achieved by taking, depending on the case, a trivial solution –e.g. the same
starting dates for all species– or the empirical dataset –if available. Then,
the proposal of the novel starting dates is done under the condition that a
non-zero amount of temporal overlap between interacting mutualistic partners
must be preserved. In detail, we impose, at least, 1 day of overlap, which is
the minimal unit of phenology in the studied datasets. To warrant that this
occurs, we fix the following boundaries to the new initial times. Let us call
$t^{\prime P}_{0,i}$ the novel time proposed for plant species $i$:

  * •

The upper bound is given by the neighbor species whose activity ceases
earlier. That is, by the minimum of the set of times
$\\{t^{A}_{0,k}+p^{A}_{k}\\}$, where $k$ runs over the indexes of the animal
species interacting with plant $i$ (so, if B is the bipartite matrix, those
fulfilling that $B_{i,k}=1$). Let us illustrate this with the example in Fig.
6. The upper boundary for the initial flowering time of plant number 1 is
given by the pollinator that finishes pollinating earlier, which is animal
number 1. Therefore, as can be seen in the figure, the upper limit for
$\\{t^{P}_{0,1}\\}$ coincides with $\\{t^{A}_{f,1}\\}$.

  * •

The lower bound is found by subtracting the period of plant i to the time at
which starts the last species. Thus, by the maximum of the set of times
$\\{t^{A}_{0,k}\\}$, minus $p^{P}_{i}$. In the example of Fig. 6, the
pollinator that begins its activity the latest is animal number 3, at
$\\{t^{A}_{0,3}\\}$. To $\\{t^{A}_{0,3}\\}$, we need to subtract the period of
plant number 1 (in yellow) in order to recover the lower limit.

According to these criteria, the new starting dates of activity are extracted
from the following uniform distribution:

$t^{\prime
P}_{0,i}\>\subset\>\textit{\Large{U}}\,\Big{(}\,\textbf{max}\big{(}\\{t^{A}_{0,k}\\}\big{)}\,-\>p^{P}_{i}+1\,,\,\textbf{min}\big{(}\\{t^{A}_{0,k}+p^{A}_{k}\\}\big{)}-1\,\Big{)},$
(11)

where the 1 is respectively subtracted and added to the boundaries in order to
ensure at least one day of overlap between mutualistic partners.

Figure 6: Example of a plant that interacts mutualistically with three
pollinator species. We calculate the upper and lower boundaries for the new
initial time of the plant $t^{\prime P}_{0,1}$ using, respectively, the
minimum final time of activity of its pollinators (here
$t^{A}_{f,1}=t^{A}_{0,1}+p^{A}_{1}$) and the maximal initial time (here
$t^{A}_{0,3}$) minus the period of the plant ($p^{P}_{1}$). As can be seen in
the figure, to these boundaries we correspondingly added and substracted one
day in order to maintain at least one day of mutualistic overlap.

Note again that this procedure requires an initial configuration of starting
dates, either empirical or hypothetical, that is then randomized.
Subsequently, the resulting set of starting dates shows some dependency on the
initial configuration.

As explained in the main text, we exploit this feature to generate synthetic
phenologies that are approximately centered. In order to produce a
configuration where species’ periods are relatively synchronized, we proceed
in several steps. First, we set the middle dates by sampling a normal
probability distribution. Once the medium dates are determined, the location
of the period of each species along the seasons is randomized using an
algorithmic procedure that conserves the mutualistic interactions provided by
the empirical network.

In more detail, we programmed a software in Python that operates as follows:

  * •

We select one of the guilds (for instance, the pollinators) and determine
their corresponding starting dates by extracting their middle dates from a
normal distribution with standard deviation $\sigma=0.5$ days and mean
$\mu=120$ days.

  * •

We then set the starting dates for the other guild, by using the randomization
process explained above. In detail, we determine the range -maximum and
minimum dates- within which we can extract a random starting date while
ensuring that the mutualistic overlap will not be lost. This novel starting
date is extracted from a uniform distribution.

  * •

To finish, we re-randomize again the first guild, using the same procedure as
in the previous step, in order to relax the assumption of the normality of
middle times of activity.

Indeed, we repeat steps (ii) and (iii) iteratively several times in order to
improve the randomization. This eventually yields a non-perfectly synchronized
configuration of phenology.

###  SM-1.2 Constrained optimizations

Here we describe a second procedure to generate a set of starting dates that
fulfills the inequalities described in Eqs. 2.9 and 2.10 in the main text. In
contrast with the algorithm described in section SM-1.1, this method does not
require any other input than the network of mutualistic interactions and the
periods, and hence there is no initial condition for the starting dates.
Instead, we seek a solution that, on the one hand, verifies the system of
coupled inequalities, and on the other hand, optimizes a given quantity,
described by an objective function.

In particular, given the constraints in Eqs. 2.9-2.10, we may use linear
programming techniques to find a feasible solution. Since $L>>N$, the problem
is overdetermined, yet we know already that it is consistent because there are
some trivial solutions that satisfy all the constraints -e.g., all species
starting or finishing the same date-. However, we are particularly interested
in non-trivial solutions that can lead to realistic or statistically
interesting configurations of phenology. This leads to the introduction of the
objective function, which we will seek to optimize over the possible set of
solutions. In sections SM-1.3-SM-1.5 we specify the numerical implementation
of different forms of this objective function, but the general structure of
the algorithm is maintained.

In detail, we solve the constrained optimization problem defined by Eqs.
2.9-2.10 through a two-step process. First, we perform multiple local
optimizations with different initial seeds, using the function minimize from
the SciPy package in Python[40] –specifically the method ’trust-constr’.
Second, we set these local solutions as the starting population for a
differential evolution algorithm, which performs a global stochastic search.
In particular, we use the differential_evolution function from SciPy. This
combination of an iterative local search plus a global optimization aims at
ensuring a robust finding of the global optimum.

###  SM-1.3 Minimization of the competition

In order to find the configuration that, under the mentioned constraints,
minimizes the competition, we apply the numerical algorithm described in
section SM-1.2 where the objective function is the total competitive overlap.
In detail, at each iteration of the optimization algorithm, we calculate the
total sum of unnormalized phenological overlaps among all competitors,
including both competition among plants and among pollinators.

###  SM-1.4 Maximization of the variance

We perform a constrained search that maximizes the variance of the middle
dates of both guilds by using the algorithm described in section SM-1.2. In
this case, the objective function corresponds to the variance, and it is
calculated at each iteration of the optimization process.

###  SM-1.5 Maximization of the entropy

In this model, we find the phenological configuration that maximizes the
entropy by using the algorithm explained in section SM-1.2. The objective
function corresponds here to the Shannon-Gibbs entropy of the middle dates, as
defined in the main text, which we explicitly computed using the Python
software.

##  SM-2 Statistical Tests

###  SM-2.1 Quality of a fit using the Kolmogorov-Smirnov test

In order to assess the quality of the fit and obtain its corresponding
p-value, we start by performing a Kolmogorov-Smirnov one-sample test between
the fitted function and the empirical distribution, using the stats.kstest
software from the SciPy package[40] in Python. Next, in order to determine its
p-value we explicitly compare the obtained KS-value with the expected
statistics, sampled using Monte Carlo simulations.

In more detail, we proceed as follows:

  * •

We calculate the K-S test statistics –using the SciPy software–, that we may
call $D_{obs}$, between the empirical distribution and the fitted function,
whose parameters are estimated using the maximum likelihood fitting software
from the stats function in SciPy.

  * •

Using Monte Carlo methods, we generate a sample of size $N_{sampl}$ composed
by synthetic distributions, each sampled from the fitted function and having
the same size as the observed empirical distribution.

  * •

For each synthetic sample indexed by $i$, we fit a function using the same
functional form as in the empirical case and using the maximum likelihood
estimator from the stats function in SciPy.

  * •

We then calculate the K-S test statistics between this novel fit and the
corresponding synthetic sample, which produces a new distance that we will
call $D_{syn,i}$.

  * •

Once this procedure has been repeated for all the synthetic distributions
generated, we can compare the observed K-S statistics $D_{obs}$ with the
distribution of $D_{syn,i}$ in the sample.

After carrying out this computation, the p-value is then calculated as:

$\text{p-value}\,=\,\frac{\text{number of samples with
}D_{syn,i}>D_{obs}\,+\,1}{N_{sampl}},$ (12)

which, as normally, quantifies whether the observed KS distance $D_{obs}$
differs significantly from what would be expected if the observed distribution
of periods was really generated by the fitted function. Hence, when the
p-value is sufficiently large the fit is compatible with this assumption.

###  SM-2.2 Kolmogorov-Smirnov two-sample test

We programmed a Kolmogorov-Smirnov test on two samples by using the ks_2samp
function from the SciPy package. In detail, the K-S test for two samples
permits challenging the null hypothesis that two particular samples come from
the same statistical distribution. The smaller the KS distance, the smaller
the discrepancy between the two samples. In the case where the p-value is
relatively small (e.g. p-value $<0.05$) the difference among distributions is
significant and hence the null hypothesis is rejected, which in our case
implies that the distribution produced by the given synthetic model is
incompatible with the empirical observations.

###  SM-2.3 Kullback-Leibler divergence

The Kullback-Leibler divergence quantifies how one probability distributions
$Q(x)$ is different from a second, reference distribution called $P(x)$. In
particular, the KL divergence from Q to P is defined as:

$D_{KL}\,=\,\sum_{x}P(x)\,log\frac{P(x)}{Q(x)},$ (13)

and it measures the amount of information lost when approximating $P(x)$ by
$Q(x)$. In our case, we calculate this quantity by taking as the reference
function $P(x)$ the empirical distribution of phenological overlap, while
$Q(x)$ corresponds to the distribution of phenological overlaps produced by
the null model under consideration. We implement this calculation in a program
in Python that uses the stats.entropy software from the SciPy package.

###  SM-2.4 Pearson correlation with time lag

In order to calculate the Pearson correlation for different values of time
lag, we take the evolution of a given structural property in the synthetic
model, and shift it alternatively backwards or ahead in time. For each
different value of time lag, we calculate the corresponding Pearson
coefficient between the synthetic, shifted array, and the unchanged empirical
vector. In particular, we use the stats.pearsonr software from the SciPy
package to compute the Pearson correlation. We repeat this procedure for a
wide range of different temporal lags (from zero to a full season), and
eventually keep the maximum value of the Pearson coefficient. This provides
the time lag at which the correlation between the empirical and the synthetic
sequences is maximized, as well as the value of this correlation.