High-order species interactions shape ecosystem diversity — TEI-tagged view

Eyal BaireypersName, Eric D. KelsicpersName, Roy KishonypersName (2016date). “High-order species interactions shape ecosystem diversity.” Nature Communicationstitle·j.

Object of study — random multi-species communities with pairwise, three-way and four-way interactionsterm·key

Research question — how high-order interactionsterm·dom affect the critical community size and the classical cap on diversityterm·key

Contribution — high-order interactionsterm·dom set a lower bound on diversityterm·key, complementing the upper bound from pairwise interactionsterm·dom

Full text — tagged

High-order species interactions shape ecosystem diversityterm·key Abstract Classical theory shows that large communities are destabilized by random interactions among species pairs, creating an upper bound on ecosystem diversityterm·key. However, species interactions often occur in high-order combinations, whereby the interaction between two species is modulated by one or more other species. Here, by simulating the dynamics of communities with random interactions, we find that the classical relationship between diversityterm·key and stabilityterm·key is inverted for high-order interactionsterm·dom. More specifically, while a community becomes more sensitive to pairwise interactionsterm·dom as its number of speciesterm·dom increases, its sensitivity to three-way interactionsterm·key remains unchanged, and its sensitivity to four-way interactionsterm·dom actually decreases. Therefore, while pairwise interactionsterm·dom lead to sensitivity to the addition of species, four-way interactionsterm·key lead to sensitivity to species removalterm·key, and their combination creates both a lower and an upper bound on the number of speciesterm·dom. These findings highlight the importance of high-order species interactions in determining the diversityterm·key of natural ecosystems.

A major challenge of theoretical ecology is explaining the stable coexistence of multi-species communities 1 , 2 , 3 , 4 , 5 . While species interactions are important for maintaining natural diversityterm·key 1 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , their random nature destabilizes large communities. Even a community that satisfies the exclusion principle 13 , or is otherwise inherently stabilized 6 , can have its stabilityterm·key jeopardized by random interactions among species 11 . This challenge was embodied in the seminal work by May 14 , showing that the number of speciesterm·dom that can stably coexist in randomly interacting communities is inversely proportional to the strength of their effective pairwise interactionsterm·mth 15 . In ecological communities with strictly random pairwise interactionsterm·dom, this result predicts an upper threshold on diversityterm·key, as observed in random pairwise interactionterm·dom models 16 , 17 . Constrained properties of interactions, such as trophic level 18 , 19 , intervality, broad degree distributions 20 , positive and negative reciprocity 21 , 22 , 23 , distribution skewness 24 , 25 and connectance between differentially self-regulating species 26 have been suggested as ways to explain the coexistence of large communities. Less is known about the feasibilityterm·key, and possibly even the necessity, of diversityterm·key within the strict context of unconstrained random interactions.

Species can exhibit interactions that inherently involve multiple species 27 . While modelling in ecology often assumes species communities with pairwise interactionsterm·dom ( Fig. 1a ), species can also interact in higher-order combinations, that is, the interactions between two species can be modulated by other species 27 , 28 , 29 , 30 . For instance, while a microbial species can produce an antibiotic that inhibits the growth of a competing species, this pairwise inhibitory effect can be attenuated by a third species that produces an enzyme that degrades the antibiotic 31 , 32 , 33 . This third species thus modifies the interaction between the antibiotic-producing and sensitive species, without having a direct effect on any of them in isolation ( Fig. 1b ). The expression or activity of the antibiotic-degrading enzyme may in turn be inhibited by compounds produced by a fourth species 34 , thereby generating a four-way interactionterm·key ( Fig. 1c ). Another general class of high-order interactionsterm·dom may arise when species exhibit adaptive behaviour 30 , 35 , such as a predator switching its prey when a more preferred prey becomes available 36 , or a prey that reacts to the presence of a predator by decreasing its foraging activity directed at a third species 37 . Such behavioural effects can quickly lead to interactions at even higher order; for example, the third species may subsequently increase its activity against a fourth species, generating an effect that is inherently four-way 28 , 38 . While the paradigm of pairwise interactionsterm·dom can capture the density-mediated effect of a predator on a species devoured by its prey that arises from the predator's impact on the prey's density, it misses the often larger trait-mediated effect 28 , 37 , 39 of the predator on this third species through its impact on the prey's behaviour.

Species can exhibit high-order interactionsterm·dom, whereby the interaction between two species is affected by other species. ( a ) Pairwise interactionsterm·dom in a community: species affect each other directly. For instance: species ‘1' could be producing an antibiotic, inhibiting the growth of species ‘2'. ( b ) Three-way interactionsterm·key in a community: one species can modulate the interactions between two others. For instance, species ‘3' might degrade the antibiotics produced by species ‘1', thereby attenuating the inhibitory effect of species ‘1' on species ‘2'. ( c ) Four-way interactionsterm·key in a community. For example, species ‘4' might produce a compound that inhibits the antibiotic-degrading enzyme produced by species ‘3'.

While the importance of high-order interactionsterm·dom has been recognized, their general impact on critical community size has not been studied. One class of three-way interactionsterm·key, whereby one species attenuates the negative interactions between two others, can stabilize well-mixed communities in specific network configurations 31 . Three-way interactionsterm·key have further been shown to promote diversityterm·key in patchy environments by increasing sensitivity to initial conditions 40 . While increasing the order of interactions weakens their destabilizing effect 41 and can increase the variance of species abundances at stationary states 42 , it is unclear how high-order interactionsterm·dom affect the critical thresholdterm·mth on community size 28 , 43 , and whether the classical result that random interactions impose a cap on community diversityterm·key still applies.

Here, we combine simulationsterm·mth and theoretical analysis of random communities to examine the relation between diversityterm·key and stabilityterm·key in ecosystems with both pairwise and higher-order interactions. We find that high-order interactionsterm·dom impose a lower bound on the number of speciesterm·dom, in addition to the upper bound imposed by pairwise interactionsterm·dom. Altogether, interactions of different orders therefore define an optimally stable range of diversities.

Results Modelling communities with pairwise and high-order interactionsterm·dom We use a replicator dynamics modelterm·mth 16 , 23 , 44 , 45 to simulate communities with interactions of different orders 41 , 42 , 46 . Our model assumes N species with abundances x i , and whose per-capita growth rates are determined by the sum of pairwise interactionsterm·dom A ij x j (where A ij determines the effect of species j on species i), of three-way interactionsterm·key B ijk x j x k (where species j and k have a joint effect on species i) and so on. This model therefore represents the Taylor expansion of a general system that includes trait-mediated effects as well as pairwise interactionsterm·dom: The model assumes random interactions perturbing an inherently stable community. The stabilityterm·key is provided by the first term (− x i ), which reflects the assumption that species are self-limiting in high concentrations. The pairwise interactionsterm·dom matrix is given by , where is a random matrix whose elements are drawn from a Gaussian distribution with mean 0 and variance 1, and α determines the strength (variance) of the pairwise interactionsterm·dom. The three-way interactionsterm·key are set by the tensor , where is a random three-dimensional ( N × N × N) tensor whose elements are drawn from a Gaussian distribution with mean 0 and variance 1, and β represents the strength of three-way interactionsterm·key; and similarly with four-way interactionsterm·key whose strength scales with a parameter γ . As the total biomass is often determined by external resources, such as nutrients, energy and space, we keep the sum of all species abundances constant by subtracting from the individual growth rate of each species the average growth rate of the species in the community.

High-order interactionsterm·dom set a lower bound on diversityterm·key As expected from the analysis by May 14 and subsequent results 16 , 17 , without high-order interactionsterm·dom we find that coexistence is lost when the pairwise interactionterm·dom strength α is increased beyond a critical thresholdterm·mth that decreases with the number of speciesterm·dom. We ran simulationsterm·mth with strictly pairwise interactionsterm·dom ( α >0 with β =0, γ =0) for varying interaction strengths α and different numbers of species N . When the pairwise interactionterm·dom strength α is small the simulationsterm·mth reach a steady state where all species coexist ( Supplementary Fig. 1a ), but when α increases simulationsterm·mth start exhibiting species extinctions ( Supplementary Fig. 1b ). Therefore, for each number of speciesterm·dom N , we can define a critical thresholdterm·mth α c as the strength of pairwise interactionsterm·dom at which community feasibilityterm·dom (positive abundance for all species at equilibrium 24 ) is lost for a given fraction of the simulationsterm·mth (5%, Supplementary Fig. 1c,d ). When the number of speciesterm·dom N is increased, the critical thresholdterm·mth α c decreases as 1/ N , in accordance with the result by May 14 ( Fig. 2a ). Thus, as communities become more diverse they are more sensitive to random pairwise interactionsterm·dom ( Fig. 2b ).

The destabilizing effect of diversityterm·key is inverted when interactions are high-order. ( a ) The critical strength of interactions beyond which the community becomes unfeasible (defined as the value at which 5% of random communities exhibit extinctions; error-bars indicate the range of 2–10%). While the critical strength of pairwise interactionsterm·dom decreases with community diversityterm·key ( α c , red, slope=−1.20±0.05), it remains unchanged for three-way interactionsterm·key ( β c , blue, slope=−0.07±0.04) and increases for four-way interactionsterm·key ( γ c , green; slope=0.95±0.04). ( b ) Example simulationsterm·mth of two communities with pairwise interactionsterm·dom of the same strength but with different numbers of species ( N =7, left; N =20, right; red square and circle in a ). While the small community shows convergence to a stable fixed-point with all species coexisting, the large community exhibits species extinction. ( c ) As shown in b , but for communities with four-way interactionsterm·key. Here the trend is reversed, and the small community exhibits extinctions, while the large community exhibits stable coexistence.

Strikingly, this classic relation between diversityterm·key and stabilityterm·key is inverted when high-order interactionsterm·dom are considered. We repeated the analysis above with strictly three-way interactionsterm·key (only β >0) and strictly four-way interactionsterm·key (only γ >0). Similar to the case of pairwise interactionsterm·dom, when the interaction strengths β or γ are increased beyond a certain level, β c or γ c , feasibilityterm·key is lost and extinctions abound. However, in contrast to pairwise interactionsterm·dom, the critical thresholdsterm·mth for high-order interactionsterm·dom do not decrease with the number of speciesterm·dom. The critical thresholdterm·mth of three-way interactionsterm·key β c remains unchanged, while the critical thresholdterm·mth of four-way interactionsterm·key γ c increases in proportion to the number of speciesterm·dom ( Fig. 2a ). These four-way interactionterm·key communities can therefore withstand stronger interactions the more species they have. For a fixed strength of four-way interactionsterm·key, these communities are stabilized, rather than destabilized, as the number species increases ( Fig. 2c ; increased stabilityterm·key with diversityterm·key can also appear in reciprocity models 21 ). Thus, while pairwise interactionsterm·dom put a higher bound on diversityterm·key, it is a lower bound that is created by high-order interactionsterm·dom. Provided that the total abundance is constant (does not change with the number of speciesterm·dom), these results are robust to various changes in the underlying model: they hold for different distributions of the random coefficients ( Supplementary Fig. 2 ), for connectance values smaller than 1 ( Supplementary Fig. 3 ) and when high-order diagonal terms are set to zero ( C ijk ≠0 only when i≠j≠k≠i; Supplementary Fig. 4 ). Similar scalings of the critical strengths of the different interaction orders with the number of speciesterm·dom also occur when the dynamic equations are replaced by a Lotka–Volterra modelterm·mth ( Supplementary Fig. 5 ).

Diversityterm·key scales differently with different interaction orders To better understand the relation between the diversityterm·key of a community and its sensitivity to interactions of different orders, we examine a modified pairwise interactionsterm·dom matrix that captures the high-order interactionsterm·dom around a fixed point. Considering an approximation for how the interactions of different orders contribute to the effective pairwise interactionsterm·mth considered in the analysis by May 14 (namely, the Jacobianterm·mth), we derive a scaling relationship for the dependence of the feasible range of N on the different orders of interactions (see the ‘Methods' section, Supplementary Fig. 6 ): Intuitively, this condition tests whether the s.d. of the sum of all random interactions affecting a species near a fixed point is smaller than the self-stabilizing term (assuming fixed total abundance, see the ‘Methods' section).

This theoretical analysis explains the different scaling of each of the interaction orders with N in accordance with our simulationterm·mth results ( Fig. 2a ). When high-order interactionsterm·dom are turned off (only α >0), the critical strength of pairwise interactionsterm·dom α c decreases as 1/ N . The critical thresholdterm·mth of exclusively three-way interactionsterm·key β c does not depend on the number of speciesterm·dom, while that of strictly four-way interactionsterm·key γ c is proportional to the number of speciesterm·dom. Thus, interactions of any given order destabilize communities, and the effects of different orders are additive, but they scale differently with diversityterm·key. Small differences between these theoretical slopes and the numerical slopes of Fig. 2a might result from the small species numbers of the simulationsterm·mth, for which the analytic analysis holds only approximately.

Mixed interactions define range of allowed diversities Since communities with a small number of speciesterm·dom are sensitive to high-order interactionsterm·dom, while communities with a large number of speciesterm·dom become sensitive to pairwise interactionsterm·dom, we asked whether the combination of pairwise and high-order interactionsterm·dom might dictate an intermediate range of species numbers that optimizes feasibilityterm·key. To answer this question, we ran simulationsterm·mth with both pairwise and four-way interactionsterm·key for three different community sizes N (for simplicity we set β =0; the effect of β >0 is discussed below). For each value of N , we scanned α and γ and identified the domain where communities are feasible ( Fig. 3a ). Since the critical value of α decreases with N , while the critical value of γ increases, there are regions where an intermediate-sized community would be feasible, while either smaller or larger communities would not (area within the N =8 domain but not within the N =5 or the N =18 domain). To understand how the feasibilityterm·key of specific communities depends on the number of speciesterm·dom, we ran simulationsterm·mth where we gradually add species with random interactions to communities with given values of interaction strengths α and γ : pairwise-dominated, four-way dominated and mixed. When pairwise interactionsterm·dom are dominant, feasibilityterm·key decreases when N is increased beyond a threshold ( Fig. 3b ). In contrast, when four-way interactionsterm·key are dominant, the community becomes sensitive to the removal of species ( Fig. 3d ). In the mixed case, when pairwise and high-order interactionsterm·dom are combined, feasibilityterm·key is maximized for an intermediate number of speciesterm·dom; both the removal and the addition of species destabilize the community ( Fig. 3c ).

The combination of high-order and pairwise interactionsterm·dom determines an optimally stable range of diversities. ( a ) Regions of stabilityterm·key for 5, 8 and 18 species communities in the space of pairwise α- and four-way γ-interaction strengths (assuming no three-way interactionsterm·key, β =0). Dots show the critical thresholdsterm·mth at which 5% of the simulationsterm·mth become unstable when the interaction strengths increase along radial lines in the log α -log γ space (see the ‘Methods' section). ( b – d ) The fraction of stable communities over the initial number of speciesterm·dom is shown for three combinations of pairwise and four-way interactionterm·key strengths (at the points labelled by ‘x' in panel a ); pairwise-dominated in b , mixed pairwise and four-way in c , and four-way dominated in d . When pairwise interactionsterm·dom dominate, diversityterm·key is bounded from above; when four-way interactionsterm·key dominate, a lower bound on diversityterm·key appears; with both types of interactions, the level of diversityterm·key is defined to a narrow range, and within this range the community is sensitive not only to the addition but also to the removal of species.

The combination of pairwise and high-order interactionsterm·dom thus defines a range of stable diversities. Assuming that a given natural ecosystem can be characterized by its average strengths of interactions of different orders ( α , β , γ …), which are characteristics of its species and the environment, the upper and lower bounds on its diversityterm·key N can be derived as the roots of equation 2 (see the ‘Methods' section; Supplementary Fig. 7 ). As the total strength of interactions α + β + γ increases, these upper and lower bounds grow closer and the range of allowed diversities narrows around a defined number of speciesterm·dom N* ( Fig. 4 for β =0; using β = γ gives qualitatively similar results, Supplementary Fig. 8 ). This value of optimally stable diversityterm·key is determined by the relative strengths of the four-way and pairwise interactionsterm·dom ( , see the ‘Methods' section). High-order dominated and strong interactions can thus require communities to maintain a large number of speciesterm·dom to remain feasible.

Range of stable diversities over interaction strength for a community with pairwise and four-way interactionsterm·key. As the total interaction strength is increased, the lower bound on the number of speciesterm·dom increases, while the upper bound decreases, until the range of allowed diversities narrows around a defined number of speciesterm·dom N* . Here, the relative strengths of pairwise and four-way interactionsterm·key is fixed so that both will have a significant effect ( γ =10 3 α ). A similar behaviour appears when accounting for three-way interactionsterm·key ( Supplementary Fig. 8 ). Note that an ecologically relevant lower bound ( N >>1) requires the coefficient of the high-order interactionsterm·dom to be much larger than the pairwise interactionsterm·dom ( γ ⩾ α ).

Discussion The combination of pairwise and high-order interactionterm·dom strengths determines both upper and lower bounds on the diversityterm·key of multi-species communities. While very little is known about the abundance and strength of high-order interactionsterm·dom in nature, these results bring a new perspective to understand much empirical evidence showing that communities are sensitive not only to the addition of species but also to the opposite perturbation, namely the removal of species 47 , 48 , 49 , 50 . Considering the perplexing effects of high-order interactionsterm·dom may allow a better understating of why high diversityterm·key is often a necessity rather than an option for natural ecosystems.

Methods Numerical simulationsterm·mth of community dynamics To find the probability that a random community made up of N species with specified strengths of pairwise and high-order interactionsterm·dom is feasible, we numerically integrate equation 1 with different randomizations of the interaction matrices. For a given N , we generate a set of random pairwise, three- and four-way interactionterm·key matrices , and then for a given choice of α , β , γ we substitute into equation 1 : , and . We then follow the species abundances as they change from a uniform initial value ( x i =1/ N ; in the stable regime, the equations should not be sensitive to initial conditions 45 ) using the MATLAB ode45 solver with default integration properties. Simulationsterm·mth were ended when they reach a persistent state (steady-state or bounded oscillations). We repeat this process R times (∼300) with different sets of random pairwise, three- and four-way interactionterm·key matrices.

We concluded that a community had reached a persistent state by demanding that the entropy (−Σ i x i log( x i )) remain constant or bounded in its fluctuations. We evolve the equations for a fixed number of time units (Δ t =10), and after each of those intervals, we calculate the ratio between the minimal value taken by the species distribution entropy in two different periods of time: the last tenth and the last three-tenths of the simulationterm·mth time length. If the ratio between those two values, as well as the similar ratio between the entropy maxima in those two periods, falls within a range of of 1, we concluded the simulationterm·mth. If this condition was not met after a long time (10 4 time units) we concluded the simulationterm·mth as well, and to account for the possibility that a species would have become extinct after this time threshold, we counted those simulationsterm·mth as unfeasible in the calculation of the lower error-bars.

We defined a community as feasible if all the species existed at the end of the simulationterm·mth, with an abundance above a defined threshold (set as 10 −5 / N ). The probability that a community would be feasible was then defined as the fraction of stable communities out of the R random communities simulated. To find the critical strength of interactions for a given number of speciesterm·dom, we used a fixed collection of R sets of pairwise, three- and four-way interactionterm·key matrices and increased the interaction parameters until 5% of the communities exhibited extinctions. In Fig. 2 (as well as Supplementary Figs 2–5 ) this was done by increasing the strength of interactions of a given order, while the rest are kept at 0. In Fig. 3a this was done along radial lines in log ( α ), log ( γ ) space: α =0.001·1.2 r ·cos θ , γ =0.05·1.3 r ·sin θ , such that for each value of θ , we increased r until stabilityterm·key was lost.

Deriving stabilityterm·key criterion by considering effective pairwise interactionsterm·mth The Jacobianterm·mth of the system at a given point is: Where, is the Kronecker delta. At a fixed point, the first term vanishes. We may also neglect the terms in the last line which are derivatives of the normalization, as we will justify below. Assuming stabilityterm·key of all species coexisting, the value of species abundances at this fixed point must scale as x i ≈1/ N (because the total abundance is fixed to 1). Equation 1 can therefore be approximated as pairwise dynamics: with the fixed point will be stable if the eigenvalues of the Jacobianterm·mth all have negative real parts. Within our assumption x i ≈1/ N , we may focus on the eigenvalues of A eff , which differs from the Jacobianterm·mth by a constant factor N that does not affect the signs of the eigenvalues (as previously reported 15 , the actual Jacobianterm·mth in a steady-state deviates somewhat from this simple analysis due to deviation of the fixed-point x i from 1/ N ; Supplementary Fig. 9 ).

The entries of the random component of A eff are sums of independent identically distributed Gaussians with 0 mean, and the variance of such a sum is additive. For simplicity, we assume that only elements with j>k>l are non-zero (such as in Supplementary Fig. 4 ), so the three-way interactionsterm·key term averages over N-1 such elements and its variance is therefore ∼ N times smaller, and similarly with the higher-order terms. The variance of the random component of A eff is thus given by: by Girko's circular law, as N →∞, the eigenvalues of a random matrix M of order N × N will form a disk of radius around the origin on the complex plane. In a large enough community, and taking into account the stabilizing term, the eigenvalues of A eff will therefore be approximately distributed in a disc around (−1, 0) ( Supplementary Fig. 6 ), so that the community will be stable with high probability if the radius of the disc is smaller than 1, or: If we remove the assumption that only elements with j>k>l are non-zero, the high-order interactionterm·dom strengths in the last equation are multiplied by constant factors, but the scaling is unaffected.

This analysis does not depend on the specific distribution of the matrix elements 51 , as long as their variance is fixed, and our simulationsterm·mth indeed show it holds for different distributions ( Supplementary Fig. 2 ). Similar scaling results can be obtained relying on Gershgorin's circle theorem (see ref. 52 for the case of pairwise interactionsterm·dom; generalization to high-order case follows from the same reasoning applied above). The normalization term should not have a significant impact on the instantaneous stabilityterm·key analysis, since by the same reasoning, its contribution to the variance of the Jacobianterm·mth would be N times smaller than the interaction term, as it is itself an average of N interaction terms. However, this analysis does depend on the sum of species abundances being constant (and x i thereby scaling as 1/ N ); a different scaling of critical interaction strengths will appear when the total abundance scales with the number of speciesterm·dom 41 , 42 , 46 .

Note that local stabilityterm·key analysis models exclusively pairwise interactionsterm·dom, since the equation is linearized around a point and the high-order interactionsterm·dom are therefore embedded in the coefficients of the effective pairwise interactionsterm·mth obtained in the linearization. Our analysis therefore extends the result by May 14 in the sense that it suggests how the interaction strength in the Jacobianterm·mth May analyses should scale with the number of speciesterm·dom with respect to the strengths of interactions of different orders.

Deriving stabilityterm·key criterion by analysis of interaction variance The same scaling of diversityterm·key with interactions of different orders can also be obtained by demanding that the expected variance of the random interactions be weaker than the self-stabilizing term. The n -species interaction term is given by: with an n -dimensional interaction tensor whose variance is 1, and α ( n ) the strength of the n -species interactions. Since the value of species abundances near a fixed point with all species coexisting scales as x i ≈1/ N , and from the additivity of the variance of independent random variables, the standard deviation of such a term scales as . Comparing this value with the stabilizing term |− x i |≈1/ N , we get the condition α ( n ) < N n−3 , in agreement with our result ( α <1/ N for pairwise interactionsterm·dom, β <1 for three-species interactions and so forth). The additivity of the variance allows us to extend this reasoning to obtain equation 2 when interactions of different orders are combined.

The optimal community size N * Restricting interactions up to the fourth order, communities are stable if they satisfy the equation: demanding an equality and solving for N gives the lower and upper bounds for number of speciesterm·dom: these two bounds narrow around a defined optimal number of speciesterm·dom when the discriminant vanishes: , which gives .

Code availability MATLAB code for the replicator and Lotka–Volterra modelsterm·mth is available from the corresponding author on request.

Data availability The simulationterm·mth results are available from the corresponding author on request.

Additional information How to cite this article: BaireypersName, E. et al . High-order species interactions shape ecosystem diversityterm·key. Nat. Commun. 7:12285 doi: 10.1038/ncomms12285 (2016).

Supplementary Material Supplementary Figures 1-9 Acknowledgments We thank Ariel Amir, Michael Baym, Guy Bunin, Michael Elowitz, Renan Gross, Simon A. Levin and Joshua Weitz for important feedback and discussions. We would also like to thank Stefano Allesina and additional anonymous reviewers for their constructive input. E.B. acknowledges the support of the Technion Rothschild Scholars Program for Excellence. E.D.K. acknowledges government support by the Department of Defense, Office of Naval Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. R.K. acknowledges the support of European Research Council Seventh Framework Programme ERC Grant 281891, National Institutes of Health Grant R01GM081617, and Israeli Centers of Research Excellence I-CORE Program ISF Grant No. 152/11.

Footnotes References Associated Data Supplementary Materials Supplementary Figures 1-9 Data Availability Statement The simulationterm·mth results are available from the corresponding author on request.