Multiple regimes of robust patterns between network structure and biodiversity

Ecological networks such as plant-pollinator and host-parasite networks have structured interactions that define who interacts with whom. The structure of interactions also shapes ecological and evolutionary dynamics. Yet, there is significant ongoing debate as to whether certain structures, e.g., nestedness, contribute positively, negatively or not at all to biodiversity. We contend that examining variation in life history traits is key to disentangling the potential relationship between network structure and biodiversity. Here, we do so by analyzing a dynamic model of virus-bacteria interactions across a spectrum of network structures. Consistent with prior studies, we find plausible parameter domains exhibiting strong, positive relationships between nestedness and biodiversity. Yet, the same model can exhibit negative relationships between nestedness and biodiversity when examined in a distinct, plausible region of parameter space. We discuss steps towards identifying when network structure could, on its own, drive the resilience, sustainability, and even conservation of ecological communities.


Results
A rule-based framework to identify distinct domains of biodiversity-nestedness relationships. We are interested in studying coexistence of multiple strains in virus-bacteria systems and its relationship to model parametrization. We define biodiversity as the fraction of surviving strains in the system and simulate the interactions between strains using the model described in the Methods. In general, coexistence of all the strains in a system, i.e., a biodiversity value of 1, requires the existence of a steady state with positive densities for all strains. For a given interaction matrix, we refer to a steady state with positive densities for all strains as a feasible steady state and to the associated parameter set as a feasible parameter set.
The mathematical conditions for feasibility can be formulated in terms of the parameters of a model: r, bacterial growth rate; K, carrying capacity; a, bacterial competition; as well as φ, m, and β; denoting virus adsorption rate, decay rate, and burst size respectively. As an example, consider a low-dimensional system with two bacteria, two viruses, and a nested infection network (virus 1 infects bacteria 1 and 2 and virus 2 infects only bacteria 1). To simplify the analysis we assume that virus adsorption rates and burst sizes are independent of bacterial strains (i.e. φ ij = φ j , β ij = β j ) and that intra-specific and inter-specific competition between bacterial strains is equal (i.e. a ii′ = 1). Under these assumptions the conditions for feasibility are: . These conditions represent a feasible "volume" in the 9-dimensional parameter space, i.e., a region of parameter space where coexistence of all 4 strains is possible. Further, we see that the condition for feasibility involving growth rates divides the corresponding two-dimensional parameter subspace into two areas: a feasible region where coexistence of all strains is possible(r 1 > r 2 ), and an infeasible region where coexistence of all strains is not possible (r 2 < r 1 ). As should be evident, randomly sampling parameter space would affect conclusions regarding the total potential biodiversity in the system. We can extend the analysis of this low-dimensional example to the case of a perfectly nested infection network with 10 bacteria and 10 viruses (Fig. 2a). The rules for feasibility are an extension of the rules from the previous 2 virus, 2 bacteria example, namely: > > … > r r r . These rules can be generalized for a perfectly nested network of arbitrary size with equal number of bacterial and virus strains 17 . In the case of 10 bacteria and 10 viruses, the rules define a feasible region in the 41-dimensional parameter space of the model. We sample parameters from two different regions of parameter space to illustrate the effect of model parameterization on biodiversity-nestedness relationships. First, we sample parameters from the feasible region of the perfectly-nested network (Fig. 2a) where all strains coexist when the interaction matrix is nested. Second, we sample parameters from the feasible region of a low-nestedness network (Fig. 2b); details of the network and its feasible region are presented in the Methods for the rule-based framework. Given the sampling regions, numerical simulations of model dynamics are averaged over an ensemble of 100 infection networks representing a gradient in nestedness while conserving the connectance of each network in the ensemble. When parameters are chosen from the feasible region of the perfectly nested network (Table 1), we observe a positive biodiversity-nestedness relationship (Fig. 3a). In contrast, when parameters are chosen from the feasible region of a low-nestedness network (Table 2), we see a negative biodiversity-nestedness relationship (Fig. 3b). The perfectly nested network was included in the ensemble of 100 networks. Sampling parameters from its feasible region resulted in survival of all the strains in the system (average biodiversity of 1, Fig. 3a). The low-nestedness network used to obtain the second sampling region was not included in the ensemble because it has a different connectance than the perfectly nested network and the rest of the networks. We found that the trends in biodiversity were robust to changes in initial conditions when sampling values from biologically plausible regions. The key point of this analysis is that the relationship between biodiversity and nestedness differs qualitatively when examining two distinct feasible sets of parameters.
A feasibility-based framework to identify distinct domains of biodiversity-nestedness relationships. For a general interaction network the conditions for feasibility are more complicated than the conditions presented for the perfectly-nested network derived in the previous section. As a consequence, we developed an alternative framework for selecting parameter sets that maximize biodiversity in Eqs. (1)-(2) for any non-trivial infection network. This alternative framework does not require finding the rules for feasibility explicitly. Instead, we choose a subset of the parameters randomly from a biologically plausible region (Table 3), specify the target steady state densities of the bacteria and viruses, and then solve for the rest of the parameters using the steady state equations. In this way it is possible to obtain a particular feasible parameter set for any infection network for which all nodes have at least b a Perfectly nested network (NODF = 1). White cells denote that a virus in that column is able to infect the bacteria in that row, whereas blue cells denote the absence of infection. one link (see Methods for more details). We selected three particular parameter sets for which biodiversity is maximized for three different infection networks of 10 bacteria and 10 viruses. We calculated average biodiversity for an ensemble of 100 matrices which included the three focal networks (Fig. 4). Figure 4 shows that not only is it possible to maximize biodiversity for different networks, but it is also possible that the resulting trends are qualitatively different. We show that, by maximizing biodiversity for a low-nestedness matrix, we obtain a negative trend of biodiversity vs. nestedness (Fig. 4 left). In contrast we obtain a positive trend of biodiversity vs. nestedness by selecting parameters that maximize biodiversity for a perfectly nested networks (Fig. 4 right). It is also possible that the quantitative strength of the relationship can differ, approaching the case where the is no significant relationship between biodiversity and nestedness ( Fig. 4 middle). This analysis supports the prior conclusions, i.e., that biodiversity-nestedness relationships depend on model parameterization. Table 1. Parameter ranges used to obtain feasible parameter sets for the perfectly nested network. Parameters were sampled from uniform distributions and sorted according to the rules of feasibility for the perfectly nested network presented in rule-based framework of the Results. The limits of the distributions are specified in the column titled Range. A fixed value was used for K.  (Table 1). (b) Parameter sets sampled from the feasible region of a low-nestedness network (  Table 2. Parameter ranges used to obtain feasible parameter sets for the low-nestedness network. Range denotes the limits of the uniform distribution used to generate parameters. The limits were calculated such that the parameters satisfy the rules of feasibility for the low-nestedness network presented in the feasibilitybased framework as shown in the Methods. Fixed values were used for φ, β, and K. Table 3. Parameter and target steady state density ranges used in the feasibility-based framework. Bacteria growth rates, r i , and virus decay rates, m j , were derived using the steady state equations and the parameters presented in this table (see Methods, given feasibility-based framework). The range denotes the limits of the uniform distributions used to generate parameters.  the original parameter set which maximized biodiversity for a particular interaction network (as explained in the previous section). The size of the region was set by the parameter δ such that for a given value of a parameter, θ, the random parameter values were sampled from the interval θ δ θ δ ( − ), ( + ) [ 1 1 ]. Here θ is a dummy parameter denoting any of the life history traits, r, K, m, β and φ that impact Eqs. (1)-(2). Figures 5(a-f) show evidence for a robust, negative biodiversity-nestedness relationship that persists even as the interval width increases. We see that the negative trend is robust to selecting parameters from (Panels (a-f) red) Average biodiversity for 100 parameter sets sampled from a uniform distribution centered around a parameter set that maximizes biodiversity for a low-nestedness network. (Panels (g-l) blue) Average biodiversity for 100 parameter sets sampled from a uniform distribution centered around a parameter set that maximizes biodiversity for a perfectly nested network. Each plot corresponds to a different value of δ which determines the size of the region used for sampling. The slope of a linear fit α and coefficient of determination R 2 are presented (p < 10 −5 for all the fitted lines). a region, but the trend becomes weaker as the size of the intervals are increased. Similarly, Figures 5(g-l) show evidence for a robust, positive biodiversity-nestedness relationship observed in a different area of parameter space that also persists as the sampling interval is increased. In summary, qualitatively distinct biodiversity-nestedness relationships persist, i.e., are robust, when parameter values are chosen at random from a region of parameter space.

Discussion
We analyzed a nonlinear model of phage-bacteria dynamics (Eqs. (1)-(2)) as a means to investigate the entanglement between interaction network structure, life history traits, and biodiversity. Given the complexity in the space of possible networks, we considered ensembles of networks that varied in a particular structural feature -nestedness -such that interaction ranges differ in the extent to which they form partially ordered subsets of one another. We found that there is not one globally applicable, positive relationship between biodiversity and nestedness in this model (Fig. 3). Instead, we identified distinct regions in parameter space where there are contrasting relationships, both positive and negative (Figs 4 and 5).
Elevated nestedness is a common feature of interaction networks spanning both plant-pollinator and phage-bacteria systems 12,13 . Moreover, prior theoretical work has suggested that ecological communities whose interaction networks are nested are more likely to have higher relative biodiversity 12 . Our results highlight the need to understand the life history context underlying a given biodiversity-nestedness relationship. The totality of parameter space includes a subspace of biologically plausible values. Such subspaces often have relatively uninformative prior distributions. Therefore, using biologically plausible regions to restrict the parameters of a model is not a strong enough restriction to uniquely define the effect of network structure on community dynamics.
Recently, it has been pointed out that different model parameterizations can lead to different biodiversity levels and consequently to contradictory results 20 . We expand on this point to show that this is also the case for whole regions of parameter space. This is important because studies using numerical simulations often average over different parameterizations. Indeed, Rohr et al. 20 examine to what extent parameters can vary for a given network structure before the community make-up changes. In contrast, our work examines how changes in network structure affects community make-up for a fixed set (or region) of parameters. These approaches are related, but they are not the same. Our results show how averaging over parameterizations is not sufficient to account for the effects of life history traits and that a more systematic study of parameter space is necessary. Additionally, we make a stronger case for the effect of parametrization by showing that it is possible to not only maximize biodiversity for specific networks but to obtain completely different trends of biodiversity vs. nestedness (Fig. 6).
In our view, statistical inference from numerical simulations can be informative and even advantageous, so long as certain precautions are kept in mind. The key point is that systematic analysis of parameter space is necessary whenever a numerical approach is used to characterize network structure-biodiversity relationships in nonlinear ecological systems. Studies that rely on analytical methods to estimate biodiversity or related features usually focus on fixed-point equilibrium states that are stable either locally or globally. This could be problematic for two reasons. First, general analytical solutions could be hard to find or interpret. Second, coexistence in high-dimensional ecological models could be characterized by non-equilibrium steady states. In such circumstances, fixed-point analyses would overlook configurations that are ecologically relevant.
In the case of phage-bacteria dynamics, our study highlights the value of additional measurements of life history traits, complementary to the recent focus on methods for characterizing who infects whom 21 . We used a particular phage-bacteria model to highlight the importance of systematically studying model parametrization in distinct regions to better understand the relationship between network structure and biodiversity -yet the findings are relevant to a wider debate. The current findings point to the need to revisit the relationship between network structure, life history traits and biodiversity in other systems and given other kinds of network patterns. Optimistically, the systematic study of model parametrization could be of service in resolving ongoing debates concerning the relationship, or relationships, between biodiversity and network structure in plant-pollinator systems 14,15,20,22 .

Methods
Model. The dynamics of virus-bacteria systems can be modeled using systems of coupled, nonlinear differential equations 23,24 . Here we use a system of equations that extend the basic Lotka-Volterra equations 25,26 to incorporate multiple types of virus and bacteria and competition between bacteria 17 : where there are n h bacteria types, each with density H i and n v virus types, each with density V j . In this system: r i is the growth rate of bacteria i in the absence of phage and other bacteria, a ii′ is the competitive effect of bacteria i′ on bacteria i (set to 1 for the analysis), K is the system-wide carrying capacity, φ ij is the adsorption rate of phage j when attaching to bacteria i, β ij is the effective burst size of phage j when infecting bacteria i, and m j is the decay rate of virus j. Finally the element M ij denotes which virus infects which bacteria such that M ij = 1 if bacteria i is infected by virus j and is zero otherwise; altogether these interactions can be represented as a network (Fig. 7) where each type is equivalent to a strain 27 .  We are interested in identifying fixed points where all strains have positive densities, such that: We refer to these points as feasible fixed points. In the present model the feasible fixed point is the solution ⁎ H and ⁎ V to Eqs. (3) and (4). Here,  r and m are vectors whose elements are the growth and decay rates and β, φ, and M are matrices whose elements are the β ij , φ ij , and M ij respectively; "⚬" denotes element-wise matrix multiplication.
Network ensemble. We generated ensembles of networks with fixed connectance and fixed size, n h × n v = n 2 , given the further assumption that n = n v = h h . The ensemble generation procedure builds upon that introduced by James and colleages 15 . Starting with a perfectly nested network, randomly chosen interactions are removed from the existing interactions among bacteria and viruses and an equal number of new interactions are added among bacteria-virus pairs where interactions did not exist originally. In the conventional matrix representation this is akin to random removal of interactions in the top left corner of the matrix followed by an equivalent addition of new interactions in the bottom right corner of the matrix. Figure 7 shows an example of three matrices with different values of nestedness generated using this procedure. To generate a large ensemble of matrices with varying degrees of nestedness the number of randomly selected links that where moved was varied from 1 to (n/2) 2 . We selected 100 invertible matrices for our ensemble and used the Non-overlapping and Decreasing Fill (NODF) metric for nestedness 8 .

Estimating biodiversity from numerical simulation of the dynamics. Biodiversity is defined in
this study as the fraction of surviving strains in Eqs. (1)-(2) after a sufficiently long transient. We estimated biodiversity by numerically integrating the system using MATLAB's ODE45 and enumerating the number of strains with densities greater than 10 −10 ml at the end of the simulation. We used this criteria for survival instead of relying on stability analysis of equilibria (or other attractors) because strain coexistence could occur via a stable fixed point, periodic cycles, or chaotic dynamics.
The stopping time of the simulation was determined via a heuristic that evaluates the convergence of time-averaged densities and mitigates inconsistencies introduced by arbitrarily choosing a stopping time.
In developing the heuristic, we leveraged the fact that the average density of each strains is equal to its equilibrium density in Lotka-Volterra systems. The algorithm to determine stopping time is as follows: • Every 500 hours calculate the infection matrix of only those strains with densities greater than 10 −10 .
• If the subsystem is solvable (invertible community matrix) calculate the theoretical interior fixed point. • Check if the average density over the last half of the simulation is within 10 percent of the theoretical prediction of the subsystem. • The simulation stops when the last condition is satisfied or the simulation has run for 40000 hours.
Note that the algorithm only calculates the stopping time and does not alter the densities of the strains. To evaluate the robustness of the stopping criteria, we calculated persistence using the stopping time heuristics and compared the results with ones obtained after doubling the time for all simulations. Figure 8 shows the relationship between persistence and nestedness and the difference in mean persistence. Note that the relationship between biodiversity and nestedness are the same whether using the heuristic stopping time or twice this value (Fig. 8a,b). We find that there are less than 1% differences in the average biodiversity between the two stopping times (Fig. 8c). Altogether, the heuristic captures the long term coexistence of different strains even when the dynamics do not tend to a steady state.
Parameter range selection. We use two different frameworks to select parameter ranges. In doing so, we use the term "plausible" to denote those parameter sets and steady-state densities whose values are consistent with virus-host biology and ecology. In addition, we use the term "feasible" to denote those parameter sets whose associated steady state densities are all positive for a given infection matrix. The steps to generate parameter sets in each of the frameworks are: Ruled-based framework: • Choose a focal matrix.
• Solve for the feasible steady state of the focal matrix implicitly in terms of model parameters.
• Choose plausible regions in parameter space that satisfy the feasible steady state conditions.
• Generate a parameter set by sampling uniformly from feasible and biologically plausible regions. Feasibility-based framework: • Choose a focal matrix.
• Sample one set of values each of carrying capacity, adsorption rate, burst size, virus and host steady state densities (K, φ j , β j , ⁎ H i , ⁎ V j ) for all species from biologically plausible regions. • Solve the steady state equations to find specific values of host growth rate and virus mortality rate, r i and m i , respectively for all species. • If the values of r i and m j are plausible, then these together with the previously sampled K, φ j , and β j form a feasible and plausible parameter set. • Generate a parameter set by sampling uniformly from a region centered around the specific feasible parameter set.
We used three focal networks in the feasibility-based framework (Fig. 7) and two focal networks in the ruled-based framework (Fig. 2). The following sections explain each framework in detail.
Parameter range selection in the ruled-based framework. In the rule-based framework, we start with a particular focal matrix and analytically solve for the steady states. From the steady states we obtain constraints on the life history traits (rules in the form of inequalities) which guarantee feasibility. For example, feasibility in the perfectly nested network requires an ordering of host growth rates ( > > … > ) r r r 1 2 10 amongst other rules. We then choose biologically plausible regions of parameter space that satisfy these rules and look at how biodiversity varies with network nestedness. For this, we draw uniformly distributed parameters from these fixed regions and calculate average persistence for an ensemble of networks that span a large range of nestedness values. We followed this framework for two different focal matrices: one with high nestedness (Fig. 2b) and one with low nestedness (Fig. 2a). This framework has the advantage that the rules for feasibility characterize the entirety of parameter space by dividing it into feasible and non-feasible regions for each focal network. Nonetheless, it is not trivial to generalize this framework to an arbitrary focal network. As a consequence we were not able to find a low-nestedness network with the same connectance as the rest of the matrix in the ensemble and with simple feasibility rules. Instead, the focal low-nestedness network used to show a different trend of persistence vs. nestedness has a slightly larger connectance (three more interactions) than the rest of the matrices in the ensemble.
1. Feasibility conditions for the perfectly nested network. We make the following assumptions regarding Eqs. (1)-(2): there is equal intra-specific and inter-specific competition across bacterial strains (a ii′ = 1), and virus adsorption rates and burst sizes are independent of bacterial strains (φ ij = φ j , β ij = β j ). Using these assumptions we find general expressions for the steady states, H * and V * , in terms of the life history traits for a 10 by 10 perfectly nested network (Fig. 2b)  networks M (Fig. 7), and values K, φ j , β j , H * , V * from random uniform distributions in biologically plausible regions. Then, using Eqs. (3) and (4), we obtained values for m i and r i . With this procedure we obtain a set of life history traits values that yield a feasible fixed point. For each parameter value, θ, local perturbations are made by sampling from an interval of length 2δθ centered around θ. The particular values of delta used are included in the caption of Fig. 5.