In franzikoch/Feedbackloops: Feedback loop analysis for competition networks

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(Feedbackloops)
library(purrr)

#load an example Jacobian matrix
contact_matrix <- read_contact_matrix("Competition_Signy_1.csv")
abundance <- read_abundance("Abundance_Signy_1.csv", contact_matrix)
interaction_table <- interaction_strengths(contact_matrix, abundance, c(-0.1, -0.9, -0.2))
Jacobian <- assemble_jacobian(interaction_table, abundance$species, "a_ij", "a_ji")

Determining stability

The eigenvalues $\lambda$ of the Jacobian matrix can be used to asses the behaviour of the system in the vicinity of the (assumed) equilibrium. Their real parts describe the rate of exponential growth or decay with which a small perturbation would increase (if the rate is positive) or decrease (if it is negative). An equilibrium is only considered stable, when all eigenvalues have negative real parts, so that all perturbations decay over time. The so-called dominant eigenvalue $\lambda_d$, that is the eigenvalue with the largest real part, can be used as an indicator of system stability: If it is negative, the real parts of all other eigenvalues must be negative as well and thus, the system is stable.

We can calculate the dominant eigenvalue $\lambda_d$ of a community matrix using the eigen() function (part of base R):

#calculate eigenvalues
eigs <- eigen(Jacobian)$values
#identify the lambda with the largest real part
lambda_d <- eigs[which.max(Re(eigs))]
print(lambda_d)

Our example matrix has a positive real part, which means that it is unstable.

$Re(\lambda_d)$ tells us whether a given system is stable or not, but it cannot be used to directly compare the stability of different systems, with differing diagonal values. To be able to make this comparison, we use the concept of critical selfregulation $s^$ [@neutel2002StabilityReal]. s represents the factor by which observed intraspecific interaction strength has to be multiplied to obtain stability. The metric s thus describes how far the system is from stability, as a multiplier of the actual amount of self-regulation, in the case of an unstable system ($s > 1$). In the case of a stable system ($s* < 1$), it indicates how much ‘buffering capacity’ a system has, giving the fraction of observed self-regulation that is enough to provide stability.

We determine $s^$ relative to the estimated diagonal values, by mutliplying all diagonal elements with a control parameter $s$ and continously increase $s$ until the matrix becomes stable. The critical amount of self-regulation $s$ is then the smallest s value that causes the real parts of all eigenvalues to be negative.

In order to calculate $s$ we need to replace missing diagonal values. The replace_zeros() function replaces any diagonal element $a_{ii} =0$ with the mean interaction strength of the matrix, multiplied with a factor $f$. Then, we can calculate the critical amount of self-regulation of a community matrix with the find_s()* function:

#replace missing diagonals with mean interaction strength * 0.1
Jacobian <- replace_zeros(Jacobian, f = 0.1)
#calculate s* with a precision of 0.01
s_star <- find_s(Jacobian, step_size = 0.01, max_s = 500)
print(s_star)

As our community matrix is unstable, it needs additional self-regulation to become stable and $s^* > 1$.

Scaling interaction strengths

$Re(\lambda_d)$ can only be used to compare the stability of matrices that with identical diagonal values. However, the strength of self-regulation differs a lot within and between empirical networks. We follow the approach of @neutel2014InteractionStrengths to deal with this problem: The interspecific interaction strengths are normalized (scaled) by dividing each row in the matrix by the absolute value of its corresponding diagonal term. This translates the diagonal structure of a matrix $\mathbf{A}$ onto the off-diagonal structure, which facilitates comparison of stability (at the cost of not knowing whether matrix (in)stability arises from patterns in the diagonal or off-diagonal structures). The resulting scaled matrix $\overline{A}$ has the elements $\frac{a_{ij}}{|a_{ii}|}$ and a uniform diagonal of $-1$.

The critical self-regulation $s^*$ can also be obtained by setting the diagonals of the scaled matrix to 0 ($\overline{A_0}$) and then calculating the real part of the dominant eigenvalue [@neutel2014InteractionStrengths; @thorne2021MatrixScaling].

A given Jacobian matrix can be scaled directly using scale_matrix(). The parameter aii specifies the value of all diagonal elements in the scaled matrix (default is 0):

scale_matrix(Jacobian, aii = 0)

To combine the scaling with randomisation procedures,the values in the interaction tables must be scaled using scale_interaction_table(). The function requires r_factor as an input parameter, that determines how missing diagonal values are replaced.

table_scaled <- scale_interaction_table(interaction_table, r_factor = 0.1)
head(table_scaled)

The function adds two new columns to the dataframe: \$a_ij_scaled and \$a_ji_scaled. These can be used as the input to the randomisation functions.

Analysing feedback structure of the network

Local stability analysis enables us to asses the stability of an equilibrium point by analysing whether a small perturbations that is introduced to the system grows or decays. It is not able to explain why the equilibrium is stable or not. In order to investigate the mechanisms behind the dynamic stability of a network, we have to take a closer look at its structure. How a network reacts a perturbation to one of its elements is tightly linked to the feedback loops that are present in the system.

Feedback loops are circuits of effects, caused by mutual interactions between species. In general, a feedback loop of length $n$ contains $n$ links, each with a different weight that corresponds to the interaction strength. The strength of the whole loop, defined as the product of all links within it, determines its effect on the system's response to a small perturbation. If the product is positive, a disturbance is amplified and the system is driven away from its initial state, while negative loops dampen disturbances and cause the system to move back. Since all species in a network are connected via loops of different lengths, a disturbance to just one species can have direct or indirect effects on all others - but how exactly the system responds depends on the interaction strengths and how they are distributed over positive and negative feedback loops.

Loop weights

The strength of a given feedback loop is defined by the product of all links within it. These values are often very small and hard to compare for loops of different lengths. We use @neutel2002StabilityReal 's loop weight measure to make this comparison easier. The weight of a loop is defined as the geometric mean of all interaction strengths within it.

The function loops(n, A) calculates the strengths and weights of all loops of given length $n$ within the matrix $\mathbf{A}$ and returns a dataframe containing the strength (product of all links) and weights (geometric mean of al links) of all loops:

loops(2, Jacobian)

For long loops (more than 4 links), it takes a while to compute. To save time use the loops_parallel() function in order to run the calculation on more than one core (Does not work on windows at the moment) :

tictoc::tic()
l <- loops_parallel(7, Jacobian, n_cores = 2)
tictoc::toc()

tictoc::tic()
l <- loops(7, Jacobian)
tictoc::toc()

Total feedback values

We use Levin's feedback measure [@levins1974QualitativeAnalysis] to understand how the effects of all feedback loops within a network are combined. His measure allows us to calculate the total feedback for a whole system, as well as for all subsystems within it, that form lower levels of organisation. At any given level $k \leq n$ within a system of $n$ species, the total feedback $F_k$ is defined as the sum of the strengths of all feedback loops of length $k$ and that of all combinations of non-overlapping shorter loops containing $k$ elements.

For larger systems and total feedbacks at higher levels of organisations, the way in which the feedback loops must be combined becomes increasingly complicated. However, it has been shown that the characteristic polynomial of the Jacobian matrix can be used to calculate total feedback values. The roots of the polynomial are the eigenvalues, while the coefficients are formed combinations of matrix elements which correspond to the total feedback values [@brooks2006CoefficientsCharacteristic]. The calculation of the characteristic polynomial is already implemented in the charpoly() function of the R-package pracma.

The function get_Fk(k,A) is a wrapper around this function, that picks the right coefficienct that corresponds to the total feedback at a given level k:

#calculate total feedback at level 3: 
get_Fk(3, Jacobian)

List of references

franzikoch/Feedbackloops documentation built on July 1, 2023, 12:42 p.m.