Let’s imagine a competitive network with 2 plant species and 2 pollinator. Plants and pollinators interact mutualistically an comepte within guilds. The system can be represented by the following matrix and graph:
\[A=\begin{bmatrix} -1 & -0.5 & 0.5 & 0.5 \\ -0.5 & -1 & 0.5 & 0 \\ 0.5 & 0.5 & -1 & -0.5 \\ 0.5 & 0 & -0.5 & -1 \end{bmatrix}\]
[,1] [,2] [,3] [,4]
[1,] 2.000000e+00 -1 0 1
[2,] -1.000000e+00 2 1 -1
[3,] 1.665335e-16 1 2 -1
[4,] 1.000000e+00 -1 -1 2 [,1] [,2] [,3] [,4]
[1,] 1.6 -0.4 0.4 0.4
[2,] -0.4 1.6 0.4 0.4
[3,] 0.4 0.4 1.6 -0.4
[4,] 0.4 0.4 -0.4 1.6We want to calculate the structural stability of the system, also called \(\Omega(A)\). For that we can use the R function from (song2018?) that is the classic function to do the job:
One would note that in that function the matrix \(A\) is not used directly, but is used to calculate a variance-covariance matrix of a multivariate normal distribution. This variance-covariance matrix is \((A^{T}A)^{-1}\). The step \(A^{T}A\) produce a symmetric matrix, which is required for the matrix to be variance-covariance matrix, leading to:
\[A^{T}A = \begin{bmatrix} 1.25 & 0.5 \\ 0.5 & 1 \end{bmatrix}\]
tab=expand.grid(seq(-1,1,0.01),seq(-1,1,0.01))
calcN=function(x){sum(-solve(A)%*%unlist(tab[x,1:2])>0)}
obj=sapply(1:nrow(tab),calcN)
sum(obj==2)/nrow(tab)
Omega(A)
A=matrix(c(-1,-0.5,-0.5,-1),2)
obj2=sapply(1:nrow(tab),calcN)
sum(obj2==2)/nrow(tab)
Omega(A)