Key gaps to be filled in population and community ecology are predicting the strength of species interactions and linking pattern with process to understand species coexistence and their relative abundances. In the case of mutualistic webs, like plant–pollinator networks, advances in understanding species abundances are currently limited, mainly owing to the lack of methodological tools to deal with the intrinsic complexity of mutualisms. Here, we propose an aggregation method leading to a simple compartmental mutualistic population model that captures both qualitatively and quantitatively the size-segregated populations observed in a Mediterranean community of nectar-producing plant species and nectar-searching animal species. We analyse the issue of optimal aggregation level and its connection with the trade-off between realism and overparametrization. We show that aggregation of both plants and pollinators into five size classes or compartments leads to a robust model with only two tunable parameters. Moreover, if, in each compartment, (i) the interaction coefficients fulfil the condition of weak mutualism and (ii) the mutualism is facultative for at least one party of the compartment, then the interactions between different compartments are sufficient to guarantee global stability of the equilibrium population.
In 2006, the US National Science Foundation convened a panel to discuss the ‘frontiers of ecology’ and to make recommendations for research priority areas in population and community ecology. Several key gaps to be filled were identified, and fundamental questions were raised that, if answered, would substantially advance population and community ecology. One of these questions was understanding the effects of mutualisms in structuring ecological communities in order to assess the degree to which they determine the relative species abundances .
Advances in understanding species abundances in mutualistic webs are currently limited, mainly owing to the lack of methodological tools to deal with the intrinsic complexity of mutualisms in general, and particularly in plant–pollinator communities . The combination of the heterogeneity of their constituent species (size, morphological traits, behaviour, fitness, degree of generalization and phenology) and the large number of interactions between them have made plant–pollinator network systems difficult to analyse. Large sets of parameters for describing the interactions between plants and pollinators are typically introduced, but they are very hard to estimate from the scarce available empirical data. Therefore, theoretical approaches for plant–pollinator networks unavoidably focus on explaining qualitative features rather than making quantitative predictions. Some of these general and qualitative findings provide useful insights into the relationship between community structure and biodiversity (e.g. the fact that nestedness enhances the number of coexisting species  or the strong relationships found among species abundances, nested architecture and community stability [4,5]).
More quantitative issues such as linking species abundances with their interspecific interactions are a pending task. To be able to predict species abundances is of paramount importance in order to assess the impact on a given community of human actions or other external factors (such as climate) that lead, for example, to the loss of existing species or the introduction of new ones.
Models play an important role in understanding ecological systems, because they serve to test the importance and explicative power of different plausible mechanisms or types of interspecific interactions driving plant–animal communities, such as mutualism, competition, predation or parasitism. However, incorporating too much complexity into a modelling effort may be undesirable for various reasons. A large number of tunable parameters permits high-quality fits to data; however, this comes at the expense of complexity that is usually not justified by the empirical evidence . Moreover, such models are prone to measurement errors. Changes in the data will give rise to considerable changes in the fit parameters that are difficult to account for or to predict. Thus, from a modelling perspective, one has to strike a balance between accuracy and a level of complexity that allows for predictability, suggesting that there is an optimal level of complexity for a model [7,8]. In this setting, aggregation is a useful method to simplify a complex system in order to reach a level of description that is appropriate to the amount of data available. Aggregation in ecology has been mostly performed either by common function (functional groups) or by interaction strength into strongly interacting subsystems .
Here, we seek to explain the abundances of plants and pollinators of an ecological community in terms of their mutualistic interactions. It has been observed that the nectar depths of plants and proboscis length of their visitors are the dominant size constraints limiting the mutualistic interactions between different plants and pollinators [10,11]. Aggregating the populations of plants and pollinators by nectar depth and proboscis length into size compartments, we present a modelling framework that leads to a simple population model. We apply our framework to a Mediterranean plant–pollinator community, for which size information is available along with species abundances . We find that with two tunable parameters, our model predicts rather well the 10 size-aggregated populations (five for plants and five for animals) of this ecosystem as its equilibrium solution. As has been stressed recently , a proper discussion of equilibrium abundances predicted by population models should not only address stability, but also specify the conditions for its feasibility. We show that this equilibrium solution is feasible and globally stable.
For our modelling, we consider the dataset Stang and colleagues of [10,11], which contains extensive observations of plant–pollinator interactions, species information, abundances and size data on morphological traits of the interacting species.
We have chosen this dataset because of the simultaneous availability of morphological information, abundances and interactions. To the best of our knowledge, this is the only publicly available dataset containing such information.
The data were collected from a site in the southeast of Spain (38°22′ N, 0°38′ W), consisting of 10 observation plots of 200 m2 each, which were observed over a period of six weeks in March and April 2003 . In each of the plots, all nectar-producing plant species with more than five flowering individuals were selected, resulting in a total of 25 plant species. Samples of 5–10 plants from each species were collected, and the depth of their nectar-holding tube and its width were measured. Flower abundances per species were estimated from counting the (mean) number of flowers per blossom and the number of blossoms per plant individual. Detailed information on the 25 plant species of the study can be found in appendix 2.1 of .
Only plant–pollinator interactions for which it was verified that the pollinator was searching for nectar and/or pollen, and touching the stigma or anthers of the flower [10,11], were recorded. As a result of these observations, a total of 887 pollinator individuals were identified. Out of a total of 111 observed pollinator species, 278 specimens were measured.
Table 1, taken from , shows the data for counts of observed interacting individuals for each pair of plant–pollinators, grouped by size intervals. The rows (columns) labelled 1–5 are the n = 5 size classes for pollinator proboscis length (plant nectar depth). Each size class has a width of 2 mm, and the five size classes span a range from 0 to 10 mm. The leading diagonal from the lower left corner to the top right corner represents the size classes for which nectar depth and proboscis length match.
(a) Choice of model
Our aim is to construct a population dynamics model as simple as possible for plant–pollinators, in terms of their mutualistic interactions, capable of reproducing the measured abundances for the size classes (i.e. the values of a total of 2n different variables). We choose the number of flowers as an estimate of the plant abundance from the pollinators' perspective, because it focuses on the flower as the resource unit for both pollen and nectar . We therefore let Pi denote the abundance of flowers in size class i. Apart from the recorded numbers of plant–pollinator interactions as captured in table 1, Stang et al. did not provide separate estimates for pollinator abundances. We therefore used the sum over numbers of recorded visits to flowers over size classes of flowers as a proxy of abundance (last column of table 1) and denote these populations by Ai. We will refer to a pair (Ai, Pi) of size-compatible pollinator–plant abundances as a compartment.
The simplest choice was to extend the Lotka–Volterra equations for mutualism between two species  to a larger group of species, as considered for example by Goh , taking the general form 3.1and 3.2
In equations (3.1) and (3.2), the index i runs over the n size classes or compartments of our model, superscripts A and P are for animals and plants, respectively, and rA, rP are the growth rates. The expression in parentheses on the right-hand sides of the equations (3.1) and (3.2) are referred to as response functions. and are the carrying capacities in size class i. The terms −Ai and −Pi correspond to intraclass competition, resulting from aggregating both the intraspecific and the interspecific competition for all the species belonging to size class i. The interclass mutualistic parameters stand for the coefficients measuring the effect of plants of class k over pollinators of class i, whereas describe the effect of pollinators of class k on plants of class i.
Note that the simplest linear functional responses were chosen. A rationale for this choice is that there is empirical evidence supporting linear functional responses in plant–pollinator mutualistic interactions [16–19]. Moreover, as we will show below (see also the electronic supplementary material), the non-zero (i.e. feasible) equilibrium populations resulting from our model are globally stable. Therefore, we can view equations (3.1) and (3.2) as the result of a linearization of the response functions around a feasible stable equilibrium point of a more realistic set of responses . Global stability of the former implies that for initial populations sufficiently close to equilibrium the evolution under either set of equations will be qualitatively the same. Because our primary goal here is to predict the feasible equilibrium abundances of plant and pollinators from their interactions, we are not interested in the dynamics of how such equilibria are attained. We therefore chose the growth rates for pollinators and plants, rA, rP > 0, independent of size i.
The feasible stable equilibrium abundances are obtained by setting the response functions on the right-hand sides of equations (3.1) and (3.2) to zero. Thus, in order to predict abundances, we will have to estimate the carrying capacities for pollinators (n parameters) and plants (n parameters), the interaction parameters of (n × n = n2 parameters) and (n2 parameters) (e.g. a total of 60 parameters for the 10 response functions when n = 5). By making neutral assumptions wherever possible and taking into account features of the empirical data, such as interactions ‘forbidden’ by incompatible pollinator proboscis lengths and plant nectar depths , we next successively reduce the number of free parameters, leaving us ultimately with only two.
(b) Parameter reduction
We consider first the reduction of parameters in the response function for the plant populations, equation (3.2). Note that a proboscis length larger than nectar depth is beneficial for the pollinator but not for the flower, because pollination does not generally occur in such an interaction. On the other hand, a proboscis length shorter than the nectar depth is not beneficial to either of them. For successful pollination of a flower, we therefore should require that both the proboscis length and the nectar depth have to be comparable . In terms of the coefficients entering the pollinator response function (3.2), this implies that we should set , whenever i ≠ k, the size classes do not match. This leaves us with the n non-zero coefficients .
Next, we observe an interesting fact that singles out the consolidation of plant and pollinator abundances into n = 5 compartments (i.e. size classes of 2 mm) as an optimal aggregation level. For n = 5, the corresponding compatible pairs of plant and pollinator abundances (Pi, Ai) listed in table 2 exhibit an almost perfect linear relationship, 3.3(figure 1b and electronic supplementary material, figure S1b), which faded away when using either 10 compartments (i.e. size classes of 1 mm; electronic supplementary material, figure S1c), or worsened when grouping into three compartments (i.e. size classes of 3.33 mm; electronic supplementary material, figure S1a). The quality of the linear fit for n = 5 would be nearly perfect, were it not for the excessive abundance of pollinators in size class 3, as can be seen from the second row of table 2. This excess population is due to Apis mellifera (mean proboscis length of 5.95 mm), which had a disproportionate weight in the dataset, representing almost one-third of the total abundance of pollinators across the 111 species (295 of 887 individuals) [10,11]. Apis mellifera is an artificially introduced, domesticated agricultural species, bred and maintained, whose abundance depends mostly on local beekeeping practices [22,23].
The fact that the abundances of the other four plant–pollinator pairs are rather well fitted by equation (3.3), containing only two parameters, KP and αAP, suggests neutrality: the effect of pollinators on plants is generally independent of size class i, that is, and . Seen from this perspective, it is not surprising that the presence of A. mellifera should adversely affect this neutrality. It has been observed that when A. mellifera is extremely abundant it leads to competitive exclusion of native pollinators , implying that when removing this species its competitor species should increase. Because the competitors of A. mellifera are mostly in the same size class, we expect that the removal of A. mellifera will be rather contained, having a small effect on the population of pollinators in the other size classes. For our modelling, this suggests that we should handle model parameters involving size class 3 separately and make neutral choices wherever possible for all other size classes. In the electronic supplementary material, we develop a model along these lines and show that its predictions can be related to those of an equivalent neutral model in which one discounts the pollinator population in size class 3, so that equation (3.3) holds for all size classes. In this setting, A3 will refer to a hypothetical population in which A. mellifera was removed, so that the plant–pollinator community becomes neutral across all size classes. For the sake of simplicity, we proceed with the neutral model, choose n = 5 and assume that a linear relation of the form equation (3.3) holds between all equilibrium plant–pollinator abundances. We obtain the discounted equilibrium population A3 by interpolating the pollinator abundance in compartment 3 to the linear fit of equation (3.3), yielding the value shown in parentheses in the corresponding entry in table 2.
Equation (3.3) serves as a neutral constraint that equilibrium abundances should satisfy. We incorporate this constraint into our model by postulating that the response function for the plant abundances Pi takes the particular form (KP − Pi + αAPAi), and then equation (3.2) becomes 3.4where we choose for KP and αAP the values obtained from the linear fit of the abundances in table 2 to equation (3.3), 3.5Thus, regarding the plant response function we made two assumptions: (i) the neutral assumption that the carrying capacities of plants are independent of size class i (i.e. we set in equation (3.2)), and (ii) that plants can only benefit from pollinators of the same size class and that this benefit is neutral (i.e. independent of the size class i). We turn next to a reduction of parameters in the response function for the pollinator population, equation (3.1). We again make a neutral assumption for the carrying capacities of the pollinators by assuming that they are independent of size class, setting . To estimate the value of mutualistic effects of the plant populations on the pollinators, we resort to the data in table 1. The near absence of non-zero entries below the leading diagonal demonstrates the so-called size threshold rule [10,11]: nectar-searching pollinators will only visit flowers when their mouthparts suffice to reach the nectar. This information implies that the coefficient that describes the effect of plants in size class k on the abundance of pollinators in size class i can be chosen as = 0 if k > i.
(1) Interactions between plants with openly accessible nectar and pollinators with long proboscises were far less frequent than expected, as one can see from table 1 where the number of visits of pollinators of size class 4 or 5 to plants of class 1 or 2 are low (upper left corner of table). Hence, a significant tendency towards size matching was observed; let us call this the size matching rule.
(2) There existed a positive correlation between body mass M and proboscis length λ: λ ∼ M0.61.
Fact (1) suggests the possibility of setting the coefficients to a single parameter αPA weighted by a function wik that decays with the distance to the leading diagonal (the difference of mean traits between size class k and i). A natural and simple choice for wik was the inverse of the mean body mass of class i, because it simultaneously led to (2), and to the general ecological trend, imposed by energetic constraints, of smaller abundance for species with larger body masses  (even though they have more sources to feed from). Hence, we took for all k and q = 1.64. Here, λi = 2i − 1 is the median size (in mm) of size class i. We do not treat q as an additional parameter, as taking q in the range of 1.5–2 yields rather similar results (see electronic supplementary material, figure S2). In this way, the 5 × 5 matrix is given in terms of a single parameter, αPA, and equation (3.1) reduces to 3.6
From the response functions in equations (3.4) and (3.6), we finally obtain the feasible equilibrium abundances as solutions to the system of equations 3.7and 3.8for i = 1, 2, … , n.
Recall that the parameters KP and αAP are obtained by a linear fit (equation (3.3)) to the empirically observed 2n = 10 plant–pollinator abundances listed in table 2, and are therefore fixed and given by equation (3.5). Thus, we are left with only two free tunable parameters, αPA and KA, which we will use to predict the actual values of the abundances.
(a) Feasibility and global stability of equilibrium
Before assessing the fit of the model to the data, let us focus on its structural stability and feasibility. Therefore, we worked out the range of parameter values αPA and KA under which the equilibrium given by equations (3.7) and (3.8) is feasible and globally stable.
A theorem by Goh  asserts for mutualistic Lotka–Volterra equations such as our model that if an equilibrium is feasible and locally stable, then it necessarily is globally stable. As we show in the electronic supplementary material, a sufficient condition for the equilibrium abundances to be feasible requires that 4.1and the resulting equilibrium point is stable if , which, because for all i, becomes 4.2We conclude that the constraints of equations (4.1) and (4.2) are sufficient to ensure that the equilibrium point given by equations (3.7) and (3.8) is feasible and globally stable.
(b) Model–data agreement
We found that our model explains the abundances of plants and pollinators in terms of their mutualistic interactions with great accuracy. This is shown in figure 1b, corresponding to the values of tunable parameters producing the best agreement with empirical data: αPA = 2.816 × 10−4 and KA = 98.05.
(c) Model robustness: sensitivity analysis
We also assessed the robustness of our model by analysing its sensitivity with respect to changes of its two tunable parameters. Figure 1a shows the accuracy of our predictions for the equilibrium abundances, obtained from equations (3.7) and (3.8), as we change αPA (horizontal axis) and KA (vertical axis). The constraints of equations (4.1) and (4.2) for a globally stable equilibrium with non-zero abundances are given by the dashed horizontal and vertical lines. Thus, admissible parameters lie in the top left quadrant established by these lines. The colour of each point is a measure of the root square error, the square root of the sum of the squared errors, between the empirical abundances of table 2 and those predicted, on a colour range from white (indicating low error) to dark (red, large error; see the electronic supplementary material for definition and details). The circle marked (figure 1b) corresponds to the value of parameters for which this error is lowest, whereas the position of the triangles marked (figure 1c–e) indicates parameter choices in its vicinity. Figure 1b–e depicts the empirical abundances (blue boxes) against the predicted ones (red circles) at the corresponding points in the parameter plane indicated in figure 1a.
We noted before that the linear relationship given by equation (3.3) is a good approximation for the empirical abundances (dashed line in figure 1b–e) and is a constraint that the predicted abundances have to satisfy, equation (3.8). Thus, by construction of our model, regardless of the choice of tunable parameters, the empirical and predicted abundances of plant–pollinator pairs will lie on the dashed line. The choice of parameters only affects their positions on the line, as can be seen in figure 1b–e. In each of these panels, the size class i associated with the empirical and predicted abundances increases from upper right (i = 1) to lower left (i = 5). Observe that in all four panels the predicted populations for each compartment follow the spacing of empirical data points. In particular, the abundances of plants and pollinators decrease roughly geometrically with increasing size class i. Moreover, we see that the parameters corresponding to points (b) and (d) in figure 1a yield abundances that are in rather good agreement with the empirically obtained ones. This agreement diminishes for the parameter choices in figure 1c,e.
(a) The optimal level of description: a minimal sufficient model
We have presented a method providing a minimal model that can be extended or refined in several directions, making it more realistic or accurate. We preferred the minimally sufficient description for predicting abundances of pollinators and plants with the available empirical data, because it involves the least amount of parameter fitting (two).
The problem of the estimation of parameters is closely linked with the issue of the optimal level of description, namely the number of compartments or size classes n. By optimal, we mean a model with minimum number of tunable parameters that is at the same time as robust as possible with its predictions. We proposed a simple recipe to find this optimal level for modelling plant–pollinator networks in terms of Lotka–Volterra mutualistic equations: try different levels of aggregation and choose the one that produces the simplest relationship of plant abundances versus animal abundances. This relation is then taken as a constraint on the possible equilibrium abundances, by turning it into a prescription for some of the nullclines (cf. equations (3.3) and (3.8)).
Determining appropriate levels of aggregation or complexity constitutes a major problem in describing ecological systems, because there is no a priori correct level of aggregation. For the particular Mediterranean mutualistic network we analysed, we found that an aggregation into n = 5 compartments by size of proboscis and nectar holder yields a simple relationship between equilibrium plant and pollinator abundances. We thus took this as the optimal level of description. Additionally, this simple relation reinforces the internal consistency of our aggregated model, which incorporates interspecific animal–animal and plant–plant competition through the intraclass competition term, because an interclass competition term would spoil the linear relationship. In other words, it seems that the interspecific competition for either plants or pollinators is more important within size classes than between size classes, so that this latter contribution can be neglected. This is in agreement with the recent finding that competition among pollinator species is important only for species with strong niche overlap .
There are alternative ways in which plant–pollinator communities can be aggregated. Bi-partite networks of plant–pollinator interactions have been extensively compiled and reported . A characteristic feature of such networks is that distinct, yet species-wise similar pollinator species tend to visit the similar plant species and vice versa. In network theory, nodes that are adjacent to the same set of nodes are called structurally equivalent . It is known that statistical inference methods can be successfully applied to networks to aggregate nodes according to their degree of structural equivalence . We have recently found that when this approach is applied to plant–pollinator networks, it produces aggregations that are consistent with a grouping of the species along a one-dimensional niche axis. Moreover, we find that niche position is significantly correlated with trait size .
Thus, using statistical methods to aggregate species of a plant–pollinator interaction web by connection similarity might be an alternative way of reducing the complexity in such ecosystems.
(b) Stability of equilibrium and its robustness as a result of the size threshold rule
In our model, the coupling between different compartments follows the size threshold rule, which constrains the mutualistic effects of plants on pollinators as well as that of pollinators on plants. On the one hand, the population of pollinators of a given size class can only be influenced by plants of the same or smaller size classes. On the other hand, the population of plants of a given size class can be affected only by pollinators of the same size class. This is because if the pollinator has a proboscis length greater than the nectar depth, it fails in general to pollinate the flower (i.e. the pollinator feeds without landing on the plant and thereby does not pollinate it). These two types of constraint give rise to the upper block-triangular structure of the interaction represented schematically in figure 2a.
The conditions for global stability in a mutualistic Lotka–Volterra-type system are rather mild: the presence of an equilibrium point with non-zero abundances along with its local stability is sufficient to ensure global stability . Global stability in our model arises as a simple consequence of the size threshold rule. That is, this rule generates the block-triangular coupling of the compartments, and the block-triangular structure in turn reduces global stability to the requirement that each compartment taken in isolation be locally stable.
Note that the size threshold rule is only one way of establishing a block-triangular structure, meaning that we have quite some freedom in constructing more general mutualistic compartmental models with a globally stable equilibrium. In particular, blocks of the form given in figure 2a can themselves be combined in a way to obtain more complex systems, while retaining the overall block-triangular structure. This is illustrated in figure 2b. If the resulting equilibrium for the full system has non-zero abundances, local stability of the individual constituent compartments will guarantee global stability of the interacting system. Thus, hierarchical nesting of compartmental interactions is one way of propagating stability at the (isolated) compartment level to that of the full (interacting) system.
(c) Type of mutualism predicted by the model
As we have seen, for the global stability of the system with interacting compartments, it is necessary that each of these compartments i taken in isolation is stable. By the latter, we mean setting to zero for each i in turn the abundances Aj, Pj = 0 of all other compartments j ≠ i in the population equations. Because the compartments of our model are pairs of size-compatible plants and pollinators in isolation, they each form a mutualistic two-species Lotka–Volterra system. Quite generally, for such mutualistic two-species systems to have a stable equilibrium with non-zero abundances, it is necessary that at least one of the mutualisms is facultative and that mutualism is weak .
Having set KP < 0, the plants in our model are necessarily all obligate mutualists. Stability then requires that all insects must be facultative mutualists KA > 0. For each compartment, this corresponds to figure 2a.
Obligate mutualism for plants in our model is consistent with the fact that the most abundant plant species are obligate mutualists. Indeed, the flowers of Helichrysum stoechas, Dorycnium pentaphyllum and Phagnalon saxatile, which together represent more than 60% of the total number of flowers for the 25 plant species , have hermaphroditic flowers that need insect pollination , because the plants have mechanisms to prevent self-fertilization. However, one must bear in mind that this is just a result that applies for aggregate variables, reflecting the dominant situations, and that in the ecosystem we are analysing the mutualism between some pairs of plant–animal species can actually be facultative–facultative, facultative–obligate and even obligate–obligate.
Another finding concerns the estimated interaction coefficients for our model; even though they correspond to weak mutualism, their product is very close to one (for the location of the point labelled (b) in figure 1a, 1 − αAPαPA 4.6 × 10−4) and thus at the frontier of strong mutualism, αAPαPA > 1.
A very challenging problem when modelling an ecological community as a dynamical system is that, if S is the number of species, the number of parameters is in general at least equal to S2 (interspecific interaction strengths plus carrying capacities). In the case of mutualistic networks, where S of the order of 100 is not rare, estimating these parameters represents a daunting task. By making reasonable biological assumptions, well supported by the data of a thoroughly studied plant–pollinator network, we aggregated the populations of plant and pollinator species into size classes, thereby ultimately reducing the number of free parameters to two. The resulting compartmental model is robust, explaining the observed abundances of plants and animals in terms of their mutualistic interactions. Moreover, we proved that the estimated parameters lead to a dynamical feasible and globally stable plant–pollinator community.
Our work provides a general methodology that can be used to describe, explain and predict the species abundances for mutualistic plant–pollinator networks encountered in nature, thus contributing to fill an important gap identified in community ecology.
Regarding the usefulness of this method as a predictive tool, increasing problems in many ecosystems are both deforestation  and defaunation , with profound impacts on the diversity and composition of, respectively, pollinator  and plant  communities. We believe that the proposed methodology could be used to assess the impact of species loss or the introduction of new species in a community on the abundances of the remaining species.
A caveat worth mentioning is that our modelling was done using data from a single plant–pollinator network. This limitation is due to the scarcity of datasets of plant–pollinator webs that contain abundances as well as size data about their species. However, datasets containing such morphological traits and phenology are currently being compiled . It would be interesting to apply our methodology to these datasets once they become available.
This work was supported by PEDECIBA-Uruguay, ANII-Uruguay. This work was also supported by grant no. 14B03P6 of Boğaziçi University.
H.F. and M.M. participated in the design and theoretical foundation of the study, and drafted the manuscript. All authors gave final approval for publication.
Conflict of interests
We have no competing interests.
The authors thank Martina Stang for clarifying the presence of the pollinator Apis mellifera. The authors also thank Diego Vazquez for comments about the nature of plant–pollinator mutualisms. M.M. acknowledges the kind hospitality of the Faculty of Sciences of the Univ. de la República.
- Received March 14, 2015.
- Accepted April 14, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.