Dynamics of a plant–pollinator network: extending the Bianconi–Barabási model

We study the dynamical assembly of weighted bipartite networks to understand the hidden mechanisms of pollination, expanding the Bianconi–Barabási model where nodes have intrinsic properties. Allowing for a non-linear interaction rate, which represents the seasonality of flowers and pollinators, our analysis reveals similarity of this extended Bianconi–Barabási model with field observations. While our current approach may not fully account for the diverse range of interaction accretion slopes observed in the real world, we regard it as an important step towards enriching theoretical models with biological realism.


Introduction
The main patterns of natural ecosystems are determined by the physical environment, with the interactions among the community of species responding to these drivers (Odum 1969).However, human activity impacts the structure, composition, and functioning of ecosystems, leading to direct and indirect effects on environmental health, biodiversity, and overall ecosystem well-being.To understand the influence of human activities on ecosystems, we study the dynamics of species interactions that are the drivers of ecosystem functioning.Currently, theoretical work tries to describe the time-evolution, and indeed the evolutionary dynamics of networks, using differential equations (Bascompte and Jordano 2013;Bastolla et al. 2009;Metz et al. 2023).Alternatively, one can try to extract patterns from observed time-evolution of interaction networks in field data, and describe network dynamics based on generalisations from such observations.This is the topic of the present study.
In ecosystems, species engage in diverse interactions (Diamond and Case 1986).Mutualism describes a relationship in which the organisms involved derive reciprocal benefit from each other.Mutualistic interactions play an important role in the functioning and stability of ecosystems, contributing to overall biodiversity and health of trophic levels along the food web.Animal pollination (here referred to simply as pollination) is a mutualistic relationship where flowering plants provide nectar and pollen as a food source for the pollinators (e.g., insects, birds and bats), while the pollinators transfer pollen (the male reproductive cells) from the male part of a flower (anther) to the female part (stigma) of the same or another flower (Waser and Ollerton 2006).
In the quantitative description of mutualistic interactions, accounting for their interrelated and interdependent components is challenging (San Miguel et al. 2012;Dormann et al. 2017).The most studied factors causing pattern within pollination networks are the morphology of the interacting species (particularly the match of species functional traits, e.g.proboscis size to corolla flower size), and phenology (the temporal overlap in occurrence due to biological life cycles) (Diamond and Case 1986;McCann and Gellner 2020;Valdovinos 2019).In combination, these perspectives have led to the development of theoretical and empirical methods for learning about patterns and stability of pollination network interactions (Bascompte and Jordano 2013;Boccaletti et al. 2006;Metz et al. 2023).However, current models so far fail to provide a clear mechanistic description of the dynamics of real bipartite networks.Instead, they rely on statistical analyses that use machine learning to identify patterns within datasets (Pichler et al. 2019;Terry and Lewis 2020).
Here we extend a theoretical statistical mechanics model to describe a real ecological complex system.In this adaptation, the functional traits and phenology of plant and pollinator species are incorporated in the model to describe the dynamical assembly of a mutualistic interaction network.To do so, we first describe our study system, a plantpollinator interaction network, and its relevant properties.Subsequently, we outline the theoretical framework employed in this study, aimed at developing the Bianconi-Barabási model (Bianconi and Barabási 2001) into a simple yet effective non-random network assembly model that encapsulates the principal mechanisms influencing the temporal evolution of the network.As result, this model is based on simple assumptions, yet produces an improvement on the original model.It allows us to describe the overall dynamical behaviour of the community of interacting species.However, our assumptions also introduce limitations that can be explored and addressed in future research.Burkle et al. (2019Burkle et al. ( , 2022) ) investigated the effects of wildfire severity on ecological drivers of interaction in plant-bee communities, considering the role of bees as pollinators.These studies provides an opportunity to investigate the temporal evolution of real weighted bipartite networks, where nodes possess intrinsic properties such as functional traits.In the experiment, interspecific biotic interactions (hereafter referred to as 'interactions') were quantified by observing floral visitors at open flowers within a 25-m diameter plot (491 m 2 ) during peak activity periods (sunny, calm conditions between 9:00 and 16:00 in 2014).Pollinators, defined as floral visitors that contacted the reproductive parts of flowers and moved among them, were captured using hand nets, collected, and euthanized for species identification.Plant species with flowers that were never visited were excluded from consideration.

Database: pollination network interactions
In this study, we recorded the time each bee visited a flower, allowing the creation of a time series where each step represents a state in the assembly process of a plant-pollination interaction network (Fig. 1).A total of 2589 biotic interactions are distributed over 53 time steps with an average interval of 1.7 days and an average of approximately 49 interactions per day (Fig. 2).We observed a notable increase in the number of interactions over time, peaking towards the end of the season.We attribute this trend to the phenology of both plant and bee species and incorporated it into our model (described in "Method" section).A more comprehensive understanding of this phenomenon is provided in Fig. 3a, where the illustrated plant species initiate their interactions at a later stage than the illustrated bee species.From an ecological perspective, our interest lies in the dynamics of biotic interactions at the species level.In this context, the two disjoint sets of nodes in the bipartite network represent two trophic levels in a food chain.Specifically, there are 191 bee species and 116 plant species connected with links indicating mutualistic relationships between these trophic levels (Fig. 1).While the degree to which a flower visit represents a pollination event remains a topic of research and discussion (Alarcón 2010;King et al. 2013;Young et al. 2021), we assume that mutualism within this network assembly begins with the first recorded interaction.Furthermore, the weight assigned to each link denotes the number of interactions between each pair of nodes (Fig. 1).
Within this dataset, the functional traits of bees (intertegular distance, bee size) and plant species (plant height, number of flowers, and flower head dimensions) was meticulously documented (Burkle et al. 2019).We think that these node-traits are the intrinsic properties crucial for shaping interactions and competition dynamics.Consequently, we hypothesize that the process governing the assembly of weighted bipartite network of ecological interactions can be effectively represented by an extension or modification of the method developed by Bianconi and Barabási (2001)

(described below).
Before proceeding, it is important to note that the process of collecting data on pollination network interactions raises several methodological issues that affect data analysis.Firstly, the number of nodes of these networks is substantially lower than those in the theoretical models upon which our analysis is based.We work in a system where the limits of some assumptions are explored.Secondly, the non-uniform time intervals between data collection rounds deviate from the assumption of constant intervals in our models.Finally, with time resolution of days, the number of interaction between species at each step becomes a stochastic variable.In our simple Fig. 3 a Time evolution of node's strength for plant (filled square and filled circle) and bee species (filled trianlge and filled diamond).The lines, with coefficients of determination of 0.93, 0.94, 0.98 and 0.97 respectively, represent the linear extrapolation of the dynamic evolution of each species.b Cumulative node's degree distribution ( P k = p k ), representing the probability that a node's degree is greater than or equal to k. Lines denote the best fit, with slopes − 1.616 and − 1.886 for plant and bee species, respectively.c Cumulative strength distribution strength distribution, with the best fit lines having slopes of − 1.0281 and − 0.8910 for plant and bee species, respectively model, we represent the number of interactions, m(t 2 ) , as a second-degree polyno- mial dependent on time (Fig. 2).This approach allows us to refine the model while still producing an analytical solution, as described below.

Method
The Barabási-Albert model describes the dynamic assembly of networks by analyzing the number of links connected to nodes, denoted as k i or degree of node i.This model characterizes networks where most nodes have few links, but some nodes have a large number of links (Barabási and Albert 1999).For a better representation of real networks, several authors have proposed modifications and generalisations (Dorogovtsev et al. 2000;Krapivsky et al. 2001;Tadic 2001).Among them, Bianconi and Barabási (2001) proposed a model that associates the dynamics of node degrees with an intrinsic attribute of each node.
In the Bianconi-Barabási model (BB model), the network assembly starts with a set of m 0 connected nodes.At each time step, a new node with links to m (≤ m 0 ) different nodes that are already present in the network is added.The probability, �(k i ) , of a node receiving a new link is described by a fitness factor, η i (where {η i } i∈N has distribution ρ(η) ), which represents the ability of such a node in the competition for new links: BB show that the dynamics of a node with fitness factor η i introduced at time t 0 can be described by multiscaling power law functions, where the value of C is given by: with η max is the maximum value of {η i } i∈N .
To mathematically describe the number of interactions between bee and plant species, we introduce the concept of strength of a node, s i .It is defined as the sum of the weights of the links, w ij , connecting it to its neighbouring nodes, expressed as where N i represents the set of neighbours of node i, or the interacting species of the other trophic level.If, in our model, the fitness factor represents the functional traits of bee and plant species, then even a relatively new species with few links can acquire interactions at a high rate as long as it has a large trait value.
Similar to the Bianconi-Barabási model (Bianconi and Barabási 2001), our network assembly begins with an initial set of species, m 0 , interconnected through interactions.With (1) (2) each subsequent time step a new random species is introduced to the system.However, unlike the BB model, the already connected species will increase their interactions (node's strength).Throughout this process, representing the phenology of species, the dynamics of interactions follow a second-degree polynomial, m(t 2 ) , with parameters derived from the total number of interactions per time step (see Fig. 2), This representation of the phenology, with a coefficient of determination of 0.21, provide a better approximation of the dynamical evolution of the real bipartite network, while still allowing for an analytical solution of the equations.Third-degree polynomials do not yield better results, and piece-wise polynomials (splines) preclude an analytical solution.Additionally, due to the relatively low coefficient of determination, the parameters of Eq. ( 4) can be modified to better describe a specific species of interest without altering the outcomes for the species community, thereby enhancing the model's flexibility and adaptability.
The likelihood that a species establishes new interactions, �(s i ) , has a functional form similar to that described in the BB model (Eq.1).This mathematical relationship is effective, as demonstrated by the fact that 93.3% and 96.9% of the tested plant and bee species, respectively, exhibit interaction dynamics that change over time with a coefficient of determination greater than 0.7 (Fig. 3a).Additionally, the node degree probability distribution ( p k ) of the final stage of the network follows a power law distribution for both bee and plant species, a phenomenon extensively documented in the literature (Boccaletti et al. 2006;Dorogovtsev and Mendes 2003;Pastor-Satorras and Vespignani 2004).This distribution calculates the probability that a randomly chosen node from a network has degree k.For enhanced clarity regarding these findings, Fig. 3b and c present the cumulative degree and strength distribution, respectively.
Due to the definition of links and weights in the network assembly (in "Database: pollination network interactions" section), the relationship between the degree and the strength of the nodes is linear, s i ∝ k i (also called "class I network" in Bianconi 2005).In addition, it is proved analytically that the degree distribution follows a power law (Fig. 3b) if and only if the probability of a new link is linearly proportional to the degree of the node (Krapivsky et al. 2000).Therefore, we expect that the probability of a species to get a new interaction due to the introduction of a new node has a functional form similar to the Eq. ( 1): On the other hand, the likelihood of acquiring a new interaction with a mutualistic partner depends on the individual's ability to compete for interactions.In other words, the probability that the weight of a link between two nodes increases depends on the capacity of the neighbours to increase such weight: (4) m(t 2 ) = m 2 t 2 + m 1 t + m 0 . (5) Here, the ratio w ij /s j defines the importance of this mutualistic relationship over all the links that the species j has already.
In consequence, using a continuum theory in an approach similar to Barabási et al. (2002), Barrat et al. (2004) and Bianconi (2005), the dynamical growth of our bipartite interaction network is described by the contribution made by a new species to the network (Eq.5) and the contributions of the existing species in the system (Eq.6) as follows: To solve this equation, we employ the assumptions of the dynamical assembly described above.Then, we note that � j η j s j � ≈ m(t 2 )tC (for details check Bianconi and Barabási 2001), and j∈N i η j w ij ≈ η max s i .Lastly, with initial conditions s i (t 0 ) = k i , we have that the model that takes into account the influence of functional traits in the competition for interaction between species, as well as the effect of the phenology in the assembly process is approximately:

Results and discussion
The dynamic exponent of Eq. (8), β(η i ) = 2η i +η max C , is influenced by species' trait η i , the maximum trait value among species at the same trophic level η max , and the constant C, which is shaped by the density distribution of the functional traits (Eq.3).Then we expect that the dynamics of a species interacting depends on the community of species at the same tropic level and their ability to interact.Additionally, the representation of plant and bee species phenology (Eq.4) also depends on the community of species that comprise the network (Fig. 2).This second-degree polynomial improves the fit and provides flexibility and adaptability to the model.The initial conditions for the dynamic evolution of a species of interest can always be adjusted with a proper selection of the parameters of such a polynomial (Fig. 4a).
As mentioned before (in "Database: pollination network interactions" section), our interest lies in the dynamics of biotic interactions at the species level.To assess the model, our attention will be directed towards the dynamic exponents β(η i ) of the power law function (Eq.8).For each plant and bee species, we compare the exponent of the power law fitted from the empirical data (the number of interactions over time) with the exponent of the power law predicted by their respective theoretical model.On a logarithmic plot (Fig. 4), we thus compare the slope of the power law curve (represented by black dashed lines) fitted from the empirical data (depicted as black dots) with the slopes of the power laws produced by the theoretical models (depicted by solid lines).
We can extend the BB model from simple links to weighted links because the interactions are integers, and we can consider the weighted bipartite network as a multigraph (network in which nodes can be connected by multiple links).Therefore, with a constant rate of interactions per step ( m = 49 ), we can compare the results of applying the BB (7) . model (Eq.2) with those of our proposed model (Eq.8).In addition, we do not consider the intercept with the axis of the model in the logarithmic plot.This is because, in our model, it can be fitted by an appropriate selection of the parameters of the seconddegree polynomial (Eq.4).
In our model evaluation, we selected three functional traits for bee and plant species, characterized by an approximate density distribution.Conversely, the BB model employed random numbers with density distributions similar to those of the chosen traits.The log-normal density distribution fitted the data best, but exponential and normal also sometimes were suitable (Table 1).Moreover, for each trophic level, the traits of the entire community of bee and plant species were taken into account to compute the dynamic exponents.The model was applied to all species that occurred at least two times in the network assembly.As a result, the model was applied to 75 plant species and 133 bee species, respectively.
For plant species, we chose three traits that describe different mechanisms of the plant to compete for interactions.These traits include the plant height from the soil surface to the top of the canopy, the number of flowers arranged on a stem (number of flowers per inflorescence), and the width of the flower head (capitula's width) (Fig. 5a).For the bee species case, we used intertegular distance (ITD) and body length, two traits used in  entomology as morphological characteristics related to flight ability and hence foraging range, as well as size.Additionally, to test if the functional traits are correlated with the number of interactions observed in the dataset, or if simply the most abundant species are the ones with the highest number of interactions, we use bee species abundance in the final network assembly as traits ( {η i } i∈N ) to evaluate such correlation (Fig. 5b).
In both trophic levels, our model approximately reproduces the dynamics of the species community (Fig. 5).For each trait, the set of dynamic exponents produced by our model has a mean value that closely matches the exponents from the empirical analysis.In contrast, the dynamic exponents produced by the BB model do not reproduce the dynamics of either trophic level for each trait.As depicted in Fig. 4b, the underestimation is due to the fact that the BB model adds a fixed amount of links per step, ignoring the ability of already connected nodes to alter the propensity for links.Overall, our model improves upon the base model at both the species level (Fig. 4) and the community level (Fig. 5).
The range of values of the empirical exponents of the whole community is larger than the range of dynamic exponents produced by the two models within each functional trait (Fig. 5).This discrepancy arises from differences between the assembly process The empirical exponent corresponds to the slope of the best fit of the data in a logarithmic plot.For each functional trait, we compute the dynamic interaction of the species with our model and the BB model, which is our base model of the real bipartite network and the mechanism we assume to generate our model, as described in "Method" section.Particularly, the stochastic number of new species per time step in the network deviates from our assumption of adding just one species to the system per time step.Additionally, the relatively small size of the network produces deviations in the expected behaviour described in the theory, known as edge effect.Here, in particular, the steep decline in interaction rate at the end of season is not matched well by the quadratic function, leading to what appears to be a saturation.
In our model, we note that the values of the dynamic exponents for each trait typically range between 1.61 and 3.55 (Table 1).This range closely aligns with the range commonly reported in the literature, approximately between 2 and 3.Moreover, we observe that the model's fit, for all tested traits, improves as we include more traits in computing the dynamic exponents, underscoring the significance of the final network size.This observation aligns with our earlier discussion where we emphasized that species dynamics are shaped by the interaction within the community at the same trophic level.Additionally, the exponents are sensitive to the density distribution of such traits.
From this analysis alone, we cannot establish a trait that performs significantly better for both trophic levels.The trait that best describes the dynamic interactions of the species depends on each individual species (Fig. 4).The sum of square residual (RSS) and root-mean-square error (RMSE) of the model for each trait are similar between models of the same level in the food chain (Table 1).

Conclusions
Despite the relatively small size of the system we are studying, the model proposed in this paper demonstrates the applicability of complex network theory to ecological networks.Building upon the Bianconi-Barabási model, our approach reveals that a new species in the network with few interactions, such as Melilotus officinalis (in Fig. 4a), can acquire interactions at a high rate if it has a large trait value.This highlights the model's ability to potentially capture the dynamics of species interactions by integrating morphology and phenology.Its simplicity allows for an approximate reproduction of community dynamics across both trophic levels.However, it is worth noting that the constraints imposed by the assumed mechanism for generating the model introduce limitations, preventing a precise fit.Nevertheless, this model approach may be a valuable tool for understanding and simulating complex ecological systems.
Throughout the development of the model, we thoroughly explored various topological features of bipartite networks and alternative assumptions regarding the mechanics of the assembly process.However, none of these models yielded superior results compared to those presented in this work.Moving forward, further development of the model would benefit from datasets representing larger networks with even higher sampling intensities throughout the seasons.Additionally, datasets collected in non-seasonal environments, such as the tropics, would provide valuable insights for enhancing the model's robustness and applicability.

Fig. 1
Fig.1Dynamic assembly of a weighted bipartite network of bee and plant species.This dataset was collected in the Northern Rocky Mountains of Montana, USA.The top/bottom boxes represent bee/plant species.The size of these boxes represents the number of observed interactions for each species in the dataset(Burkle et al. 2022)

Fig. 4
Fig.4Interaction dynamics of Melilotus officinalis.This plant species was introduced to the network at t = 9 .Characterized by a plant height of 56.92 cm, this plant species has a coefficient of determination for the empirical slope of 0.93.We compare the empirical exponent with the dynamic exponent of a the model we propose, and b Bianconi-Barabási model.The inserted figures represent the density distribution ρ(η) of the trait that yields the optimal fit for the data

Fig. 5
Fig. 5 Distribution of dynamic exponents for a flowers and b pollinators with three functional traits each.The empirical exponent corresponds to the slope of the best fit of the data in a logarithmic plot.For each functional trait, we compute the dynamic interaction of the species with our model and the BB model, which is our base model

Table 1
Functional trait information for the bee and plant species included in our analysis Traits of flowers and pollinators used in our analysis.ρ(η i ) refers to the distribution best describing the trait distribution, and hence employed in the Bianconi-Barabási model; minimal and maximal slope estimates as well as error statistics (residual sum of squares: RSS; root mean squared error: RMSE) are provided for each trait