Symmetric and Asymmetric Tendencies in Stable Complex Systems.

A commonly used approach to study stability in a complex system is by analyzing the Jacobian matrix at an equilibrium point of a dynamical system. The equilibrium point is stable if all eigenvalues have negative real parts. Here, by obtaining eigenvalue bounds of the Jacobian, we show that stable complex systems will favor mutualistic and competitive relationships that are asymmetrical (non-reciprocative) and trophic relationships that are symmetrical (reciprocative). Additionally, we define a measure called the interdependence diversity that quantifies how distributed the dependencies are between the dynamical variables in the system. We find that increasing interdependence diversity has a destabilizing effect on the equilibrium point, and the effect is greater for trophic relationships than for mutualistic and competitive relationships. These predictions are consistent with empirical observations in ecology. More importantly, our findings suggest stabilization algorithms that can apply very generally to a variety of complex systems.

Complex systems may undergo transitions between alternate stable states of contrasting behavior. Such a transition is called a critical transition or a regime shift in the literature1. Critical transitions are highly non-linear phenomena in that a small change in a controlling parameter such that a critical point is crossed can unexpectedly provoke a huge response (critical transition). Further away from the critical point, such a small change in the controlling parameter would only result in a comparable response without any critical transition. This non-linear response, along with the fact that critical transitions are common in nature123456, makes the study of critical transitions an important one. Critical transitions can happen as a result of instability in the stable state that the system was residing in.

In order to determine the stability of an equilibrium point, the simplest kind of stable state, a commonly used approach in non-linear dynamics is to linearize the dynamical equations describing the system about the equilibrium point. One obtains from this linearization the n × n Jacobian matrix B evaluated at the equilibrium point, with real matrix elements {b: i, j = 1,…, n} for a system with n dynamical variables x = (x1, x2,…, x). The matrix B is also known as the community matrix. The equilibrium point is stable if all real parts of the eigenvalues of B are negative and unstable otherwise. Henceforth in this paper, we may refer to B being stable or unstable when we actually mean the equilibrium point associated with B being stable or unstable respectively.

The matrix element b describes the dependence of dynamical variable x on dynamical variable x, where i ≠ j. Conversely, b describes the dependence of x on x. We may also refer to b as an interaction and its magnitude as its interaction strength. Here, we define the product bb to be the relationship between x and x. The relationship between x and x is mutualistic if b > 0 and b > 0, competitive if b < 0 and b < 0, and trophic if bb < 0. A relationship is symmetrical when b and b are of comparable magnitudes and is asymmetrical otherwise. For example, a measure of asymmetry for mutualistic relationships is |b − b|/max(b,b) from Bascompte et al.7. The main result of this paper involves using eigenvalue bounds to show that stability in B favors mutualistic and competitive relationships that are asymmetrical and trophic relationships that are symmetrical. The analysis presented here stems from a rather old research question: how do the eigenvalues of B depend on its matrix elements?

Unfortunately, there is no exact answer to this question. An approach has been to use Random Matrix Theory (RMT), originally introduced by Wigner to study spectral properties of atomic nuclei8. RMT has since found applications in a wide variety of disciplines including number theory9 and neuroscience10. In ecology, RMT was used by Robert May to study the stability of a large ecological network at an equilibrium point11. In May’s seminal work, B is a random matrix, with off-diagonal matrix elements being independent and identical random variables of mean zero and variance σ2. The diagonal elements, set at −1, represent characteristic return rates for the populations of species when disturbed from equilibrium. May claimed that for large n, B is unstable when . The main criticism with May’s work is that real-world ecosystems are structured unlike the random matrix studied by May121314. Allesina and Tang, relying on recent advances in RMT from the mathematics literature15, recently confirmed May’s claim and further analyzed random matrices with various structures16, alleviating some of the criticisms associated with May’s work. Research in RMT has hinted that high correlation between random variables b and b in mutualistic and competitive relationships has a destabilizing effect whereas high negative correlation in trophic relationships has a stabilizing effect on B171819. Conjectures in RMT typically assume at the least that n is large and that matrix elements or pairs of matrix elements are independently and identically distributed. The significance of the work presented here is the generality of our results: in fact we make no assumptions about B (besides its matrix elements being real). At the same time, we cannot obtain precise conditions for stability or instability beyond the observation that B will eventually become unstable if certain quantities become large enough.

In the next few sections, we will first present the eigenvalue bounds in terms of the matrix elements and the complex parts of the eigenvalues. Then we will show that the system will become unstable when the off-diagonal sum becomes large enough. Next, we will demonstrate a stabilization algorithm on random matrices using a random strategy, a variance-minimizing strategy, and a χoff -minimizing strategy. This will be followed by a description of a model of B with ecologically motivated constraints on the interaction strengths. In that section, we also analyze the effect of dispersion in the interaction strengths on χoff. Finally, a discussion of the results concludes the paper.

Mathematical Formulation and Eigenvalue Bounds

We start with a dynamical system with n dynamical variables described by n arbitrary non-linear differential equations i.e. x = f(x), where i = 1,…, n and x=(x1…x) is a vector of dynamical variables. x* is an equilibrium point if for every i, f(x*) = 0. The local stability of x* may be studied by linearizing the dynamical equations about x*20. The linearization furnishes B, the Jacobian matrix evaluated at x*. The matrix element b, which we described as the dependence of x on x, is the gradient of f (x) along x at x* i.e. . The equilibrium point is stable when any perturbation of x from x* decays with time. Conversely, the equilibrium point is unstable when any perturbation of x from x* grows with time. Stability is determined by the eigenvalues of B. The equilibrium point is stable if the real parts of all eigenvalues are negative and is unstable otherwise. Equivalently, the equilibrium point is stable if the largest real part of all eigenvalues, which we shall refer to as the maximum real eigenvalue, is negative and unstable otherwise. Eigenvalues are the exponential decay rates of small perturbations from the equilibrium point. Thus, eigenvalues that are more negative indicate greater stability along their respective eigenvectors. Solving for the eigenvalues is equivalent to finding the roots of the characteristic polynomial det(B − λI), where λ is an eigenvalue of B and I is the identity matrix. The eigenvalues depend on the matrix elements in a nontrivial fashion in part because there is no general algebraic expression for the roots of polynomials of the 5th degree or higher. This is the Abel-Ruffini theorem and is a well-known result from Galois theory. While others have resorted to RMT for this problem, we use an alternate approach with eigenvalue bounds to glean information about the eigenvalues’ dependence on the matrix elements.

Given the multiset of eigenvalues of B, , an upper and lower bound for the maximum real eigenvalue of B are respectively,

Here, is the mean while sλ is the standard deviation of the real parts of all eigenvalues. The upper bound is more well-known and was probably first discovered by Laguerre but is more commonly known as Samuelson’s inequality2122. The lower bound is due to Brunk23. If the matrix elements are real, the bounds may be given by an expression in terms of {b: i, j = 1,…, n} and the complex parts of the eigenvalues of B (Methods and Supplementary Information),

Here, is the diagonal sum, is the off-diagonal sum, is a function of the diagonal elements, and is a non-negative number that is positive when the imaginary components are non-zero and zero otherwise. n is kept constant throughout our analysis. The mean of the eigenvalues is controlled by the diagonal elements i.e. . Hence, the eigenvalues’ dependence on the diagonal elements is more straightforward and is generally of less interest than the off-diagonal elements. It follows from the bounds that λ+ < 0 is a sufficient condition for stability while λ− > 0 is a sufficient condition for instability.

From Equation (3), we may then draw some conclusions on two cases: (i) the eigenvalues are real numbers, i.e. h = 0 (ii) the eigenvalues are complex numbers, i.e. h ≥ 0. For the first case, it is always possible to achieve stability or instability by decreasing or increasing χoff respectively. The first case shall not be analyzed further because of the strong assumption that all eigenvalues are real. Instead, we consider the second case, which is more general. For the second case, the upper bound becomes less useful but not the lower bound because h ≥ 0; we can still always increase χoff enough such that B becomes unstable. Therefore, while it is still necessary to keep χoff small enough for stability, keeping χoff small alone does not guarantee stability because of h. For this reason, other contributing factors would still have to be considered in order to form a complete picture of how stability arises in B. For example, ecologists have been obsessed with the nestedness, a persistent structural property observed in mutualistic networks12242526. In a nested architecture of a bipartite network, a more specialist species (defined as having fewer mutualistic interactions) would only interact with a proper subset of mutualistic partners of the more generalist species (defined as having more mutualistic interactions). However, there remains some controversy over how important nestedness is to the stability of mutualistic ecological networks24. Instead of delving into the details of specific structural properties, we will focus our efforts on minimizing χoff instead. χoff is the sum of all relationships in B. A natural way to minimize χoff is to make both interaction strengths in mutualistic and competitive relationships very weak. However, mutualistic and competitive relationships are pervasive in nature. Therefore, we need to constrain the interaction strengths in B so that minimization of χoff will not render both interaction strengths in B negligible. Then, minimization of χoff will require mutualistic and competitive relationships to be asymmetric in order to minimize each summand, 2bb i.e. one large and one small interaction strength. Conversely, trophic relationships will be symmetric. In the section on interdependence diversity and symmetric correlation, we will constrain the interaction strengths in B from an ecological standpoint and demonstrate that minimization of χoff will require mutualistic and competitive relationships to be asymmetrical and trophic relationships to be symmetrical.

Pruning Random Matrices for Stability

The results presented thus far suggests that minimization of χoff might provide an efficient route to stabilize an equilibrium point. We employ a simple algorithm on a well-known example, the random matrix studied by May11. For this example, consider the situation where M is a random matrix and its diagonal elements are set at −d while the off-diagonal elements are independently and identically distributed random variables of mean zero and variance σ2. According to RMT, for large n, the eigenvalues of M are contained in a circle of radius centered at (−d, 0) on the complex plane16. For this example, we use n = 20, d = 2, a standard normal distribution for the off-diagonal elements and a modification factor g = 3/2 that we shall introduce in the description of the algorithm.

The algorithm consists following steps: (1) initialize a random matrix M, (2), calculate the eigenvalues of M, (3) choose an off-diagonal matrix element b randomly, (4) if bb < 0, multiply b by a factor g, else if bb > 0, divide b by a factor g, (5) calculate the new eigenvalues of M after the modification, and (6) if the maximum real eigenvalue of M after the modification is larger than before the modification, revert to step (2) using M before the modification; if the maximum real eigenvalue of M after the modification is smaller than before the modification instead, revert to step (2) using M after the modification. This counts as one iteration.

This algorithm employs a χoff -minimizing strategy due to step (4). We compare this algorithm using the χoff -minimizing strategy against the same algorithm using a random strategy and a variance-minimizing strategy. In the random strategy, step (4) is replaced by the following step instead: (4) b is randomly chosen to be multiplied or divided by g with probability ½. In the variance-minimizing strategy, step (4) is replaced by the following step instead: (4) b is divided by g. We compare the three strategies using 50,000 random matrices over 2,000 iterations. The results are shown in Fig. 1. The χoff-minimizing strategy clearly outperforms the other two strategies. Of course, if we were to accept every modification without checking if it reduces the maximum real eigenvalue at every iteration, then the variance-minimizing strategy will eventually reduce all eigenvalues to −d. However, there are two reasons why such an algorithm might be undesirable: (i) the maximum real eigenvalue may at times increase with iteration number, and (ii) the eventual interaction strengths are small unlike the original algorithm with the χoff-minimizing strategy which allows for larger eventual interaction strengths. Sensitivity analysis of the parameter g reveals that the χoff-minimizing strategy still outperforms the other two strategies for the various values of g tested (Supplementary Information). Statistics of the matrices at the end of the iterations and figures displaying matrices after implementations of the stabilization algorithm can be found in the Supplementary Information (Table S1 and Figure S2).

Figure 1

Pruning unstable random matrices.

Results of the stabilization algorithm employed on 50,000 unstable 20 × 20 random matrices for the three different stabilization strategies (random, variance-minimizing and-minimizing and χoff-minimizing) described in the main text. Data points at each iteration indicate the sample average over the 50,000 simulations. Standard error of the mean estimates are on the order of 10−3 for both figures. (a) The maximum real eigenvalue at the end of each iteration. (b) The proportion of all matrices with a decreased maximum real eigenvalue from the previous iteration number.

Interdependence Diversity and the Symmetric Correlation

Clearly, it is not realistic to stabilize B by rendering the interaction strengths in B negligible since interactions are ubiquitous in nature. Therefore, we now consider ecologically motivated constraints of the interaction strengths in B. To do this, we first consider two sets of similar variables, and , where x is the set of all variables and . These two sets of variables are so defined to delineate interactions of a particular type between variables from the two sets. For example, if the interaction types are consumption and pollination, then the variables in y could represent populations of pollinators while the variables in z will represent populations of plants. We now formulate equations of constraints that allow variables in y to depend on various weighted combinations of the variables in z and vice versa. For notational convenience, let us denote Y(x) to be and Z(x) to be for all k and l. Then, we may find in B the matrix elements , which is the dependence of species y on species z, and , which is the dependence of species z on species y, for all k and l. Here, d and e are weights such that , , , , 0 < d < 1, and 0 < e < 1. α and β are real numbers and because of the bounded weights, their absolute values are the maximum interaction strengths possible for the respective interactions (e.g. consumption and pollination) and species (y and z) they pertain to. These constraints on the interaction strengths can be motivated by the constant interacting effort hypothesis which states that interaction strengths should be stronger, on average, for species interacting with a smaller number of resource species27. This hypothesis was postulated based on the fact that there is a limited amount of time a species has to interact with other species. Hence, if the interaction strengths are to be proportional to the time spent on interacting with other species, and that there is a fixed amount of time for a type of interaction, then we arrive at the constraints and for all k and l. While a species can spend a different amount of time for a type of interaction at the expense of other types of interactions, it is sufficient to fix the total amount of time spent on a type of interaction for the purpose of analyzing the effect of dispersion in the weight distributions on the off-diagonal sum. Additionally, we might also require that the total dependencies on any particular species (y or z) be constrained ( and respectively). Finally, differences in α for all k and differences in β for all l lead to biases in weight distribution amongst the variables in y and z when minimizing χoff. Since our goal is to evaluate the effect of the dispersion in the weights on the off-diagonal sum, we may simplify the problem by assuming that α = α for all k and β = β for all l.

The relationship between any variable y and any variable in z is either mutualistic or competitive if αβ > 0 and trophic if αβ < 0. If we arrange the weights d into an m × m matrix D such that d is a matrix element of D, then αD is a submatrix of B. Similarly, e is a matrix element of E and βE is a submatrix of B. This convenience allows us to define a quantity C that we shall call the symmetric correlation,

The symmetric correlation is bounded 0 < C < 1 (Theorem S2) and contains information about both the variance of the weight distributions and the correlation between the matrix elements of D and the corresponding matrix elements of ET. For example when C = 1, all weights are either zeros or ones and they fulfill d = e for all l and k whereas when C = 0, all weights fulfill de = 0 for all l and k. When the variance of the weight elements in D and the variance of the weight elements in E are fixed, we may use C as a relative measure of symmetry and correlation between the matrix elements in D and the corresponding matrix elements in ET (Methods). We find that χoff contains the summand mαβC i.e. χoff = 2mαβC+…, and that the weight elements d and e for all k and l are contained exclusively in the summand 2mαβC of χoff. Thus, minimizing χoff for the relationship αβ requires minimizing C for mutualistic and competitive relationships, and maximizing C for trophic relationships i.e. mutualistic and competitive relationships will be asymmetric whereas trophic relationships will be symmetric (specifically, the weights associated with trophic relationships will be symmetric).

Next, we define the interdependence diversity,

a measure of the diversity of dependencies among the y and z variables. S is simply the sum of all the squared weight elements. Furthermore, due to the weight constraints, S is bounded −2m < S ≤ −2. When S = −2m at minimum interdependence diversity, then all weights are either zeros or ones. When S = −2 at maximum interdependence diversity, then all weights are equal to 1/m. The interdependence diversity defined here is similar to the Herfindahl index in economics28 or the Simpson index in ecology29. We denote max(C)| and min(C)| to be respectively the maximum and minimum C under any variation of weights and under a fixed S (without violating the weight constraints). The relationships of max(C)| and min(C)| with S are shown using numerical calculations in Fig. 2. Analytical calculations can be found in the Supplementary Information. Minimizing χoff means that mutualistic and competitive relationships will reside on the min(C)| curve while trophic relationships will reside on the max(C)| curve. Both max(C)| and min(C)| are monotonically decreasing and increasing functions of S respectively (Fig. 2). Hence, the capacity of B to minimize χoff decreases with increasing interdependence diversity for mutualistic, competitive and trophic relationships. Additionally, because max(C)| is more adversely affected than min(C)| with increasing interdependence diversity (for m > 2, Fig. 2), this effect is more pronounced in trophic relationships than mutualistic and competitive relationships. Essentially, trophic relationships are more affected than mutualistic and competitive relationships with increasing interdependence diversity because there exists many more possibilities in the network to minimize C. Only when the network is fully connected with equally weighted one-directional links does min(C)| start increasing with S (Proposition S2). For the more general case where x and y contain a different number of dynamical variables, we still expect trophic relationships to be more adversely affected by an increasing interdependence diversity than mutualistic and competitive relationships (Supplementary Information Section S3).

Figure 2

Symmetric correlation and interdependence diversity.

This graph shows the boundaries of values possible for C and S. The blue line is the maximum C attainable under fixed S. The red line is the minimum C attainable under fixed S. Both line plots are calculated by numerical optimization techniques with m = 5 (Methods). An analytical calculation for m ≥ 1 is provided in the Supplementary Information. As a further note, the set of C for a fixed S is not necessarily continuous within the boundaries (e.g. at S = −2m).

Discussion

In this work, we have derived eigenvalue bounds for the maximum real eigenvalue of B in terms of the matrix elements and the complex eigenvalues. From these bounds it follows that a necessary condition for stability is that χoff is small, for χoff can always be increased enough such that B will become unstable. The generality of this result and subsequent calculations allows us to consider different types of interactions in concert, something that was limited in previous studies with RMT due to assumptions on matrix elements being independently and identically distributed. Additionally, we show that two observations, increasing interdependence diversity causing decreasing χoff and this decrease in χoff being more pronounced for trophic than mutualistic and competitive relationships, can both be explained as a result of B losing its capacity to accommodate symmetric and asymmetric relationships.

In the course of implementing the stabilization algorithm, the maximum real eigenvalue will be a monotonically decreasing function of iteration number. From the eigenvalue bounds, we generally expect h + χoff to also decrease with iteration number. Indeed, statistics of the matrices after 2,000 iterations of the stabilization algorithm reveal that the average change in h + χoff is negative for all three strategies, with the effect of decreasing from the initial random matrices (Table S1). In particular, h and the standard deviation of the off-diagonal elements increase for the χoff-minimizing strategy, with the mean of the off-diagonal elements remaining constant. This suggests that while there is a certain risk in increasing h when increasing the interaction strengths, it is possible that the increase in h can be mitigated and overcome by a larger decrease in χoff such that the system can be stabilized with increasing interaction strengths.

The interplay between h and χoff is an important factor to consider when minimizing χoff to stabilize B. We have shown, under a general framework of ecologically motivated constraints, that minimization of χoff will result in trophic relationships being more adversely affected than mutualistic and competitive relationships with increasing interdependence diversity. The validity of this result for lowering the eigenvalue bounds of an ecological community will depend on whether a not minimization of χoff will necessarily give rise to an increase of h larger than the decrease in χoff in every minimization scenario under the general framework of constraints employed. Empirical observations in the ecology literature suggest that this may not the case for most communities. Our result is consistent with empirical observations if we allow the interdependence diversity defined here to be used as a proxy for connectance, a measure that is well known in the ecology community. The connectance is the proportion of non-zero dependencies in B. Hence, we generally expect an increasing connectance to also result in an increasing interdependence diversity as the number of interactions increases and as the interaction strengths become more distributed among the dynamical variables. Thébault and Fontaine found trophic networks to have a lower connectance than mutualistic networks in a meta-analysis of real-world pollination (mutualistic) and herbivory (trophic) networks while controlling for n14. Therefore, our derived result that increasing interdependence diversity having a destabilizing effect being more pronounced in trophic relationships than mutualistic and competitive relationships could provide a plausible theoretical explanation for this empirical observation.

The prediction that mutualistic and competitive relationships are symmetric whereas trophic relationships are asymmetric also agrees with empirical observations. It has been known for some time that mutualistic ecological networks like plant-pollinator networks consist of highly asymmetric interactions between plant and pollinator71230. For example, the manduvi tree relies almost exclusively on the toco toucan for seed dispersal, but the toco toucan is not limited to the manduvi tree’s fruits in its diet31. Overall, consistency of our calculations with empirical observations demonstrates our approach to be promising for further investigations of stability in B.

Our results highlight the importance of asymmetry in mutualistic and competitive relationships, and of symmetry in trophic relationships to the stability of a complex system. Identifying and understanding the contributing factors to stability can be used to help design algorithms to stabilize real-world systems on the verge of critical transitions. For example, the stabilization algorithm described in this paper could be a starting point for future investigations into the stabilization of real-world systems. In a successful realization of such an algorithm, critical slowing down signals could be used to measure the change in stability at every iteration (step (5) of the algorithm). Critical slowing down signals are statistical signals that can be used to detect if a stable state is becoming more unstable. These signals have been detected in a wide variety of real-world systems1. They are based on the premise of a slower return rate to the stable state after a perturbation as the stable state becomes more unstable32. While there have been ample studies on the detection of critical slowing down signals, more research needs to be conducted on the stabilization of potentially unstable stable states.

Stabilization is one way to deal with critical transitions. A recent attempt at this problem involves smoothening the non-linearity of a critical transition33. Network properties not covered in this work can also be very important in dealing with instability. For a formerly stable equilibrium point, initial instability occurs when the maximum real eigenvalue goes above zero. The eigenvector(s) of the maximum real eigenvalue determine the initial directions of instability and which variables will be initially affected by this instability. As the system transitions away from the previously stable equilibrium point, more and more variables might be affected depending on their dependence on the initially and subsequently affected variables in what is known as a cascade of failures. Whether a not such an initial instability will eventually lead to system-wide instability depends on a multitude of factors including the structure of the network connecting these variables and how the system responds to this initial instability. For example, in a load bearing network with a heterogeneous degree distribution, the failure of a single node with a large number of dependencies can cause a large cascade of failures34. In ecological mutualistic networks, the right and left leading eigenvectors not only determine the species affected by perturbations to the system and the size of these perturbations, they also positively correlate with a few network properties like the degree centrality and the page-rank centrality35. The effect of initial instability or failure on the whole system is a topic of great interest in network science343637. The interplay between the factors that determine stability is still an important research topic to be explored in greater detail. While there remains a host of factors that ultimately determine stability in a complex system, the generality of our results suggests that asymmetry in mutualistic and competitive relationships and symmetry in trophic relationships should be universally observed and not restricted to ecology.

Methods

Derivation of eigenvalue bounds

The polynomial equation is . We may express both bounds in terms of c1, c2 and h using Viète’s formulas and the complex conjugate root theorem. The relation between the matrix elements and the coefficients c1 and c2 can be found by expanding the Leibniz formula for matrix determinants. This gives us the bounds in terms of h and the matrix elements of B. A more detailed derivation may be found in the Supplementary Information.

Obtaining the numerical results of Fig. 2

To obtain max(C)|, we (1) construct 5 × 5 matrices Dmax and Emax at max(C)| when S is at the minimum of −2m; Dmax and Emax are initial starting points for a nonlinear constrained optimization (maximization) algorithm implemented in MATLAB (fmincon function with sqp algorithm where the constraints for the optimization problem are the weight constraints 0 < d,e < 1, , , , and , and the interdependence diversity constraint , while the objective function is C), (2) carry out the optimization for the starting point and constraints, and (3) use the solution as the new weight matrices for the starting point of the next optimization where the interdependence diversity is fixed at a positive increment ε = 0.001 from the previous optimization. Steps (2) and (3) are repeated until the maximum interdependence diversity is reached at −2. To obtain min(C)|, we use the same steps, replacing the initial starting point with Dmin and Emin at min(C)| when S = −2m and using the same optimization algorithm but with minimization instead.

The symmetric correlation as a relative measure of symmetry and correlation

Let represent a sequence of the matrix elements of D and represent the corresponding sequence of the matrix elements of ET. A measure of correlation between D and E is the Pearson’s correlation coefficient estimate

where , and are sample means, S is the standard deviation of D and is the standard deviation of E. Because of the weight constraints, are constants independent of the weight distribution of D and E. Also, . Hence, the symmetric correlation may be used as a relative measure of symmetry or correlation between matrix elements of D and the corresponding matrix elements of E when S and S are fixed.

Additional Information

How to cite this article: Tan, J. P. L. Symmetric and Asymmetric Tendencies in Stable Complex Systems. Sci. Rep. 6, 31762; doi: 10.1038/srep31762 (2016).