Propriedades em rede de interações, como descrever as comunidades de plantas e seus polinizadores?

Eu ja falei uma vez sobre como fazer gráficos de rede aqui, mas podemos explorar também várias propriedades de redes interessantes para a biologia.

Primeiro vamos simular aqui quatro redes para ajudar a pensar, como eu não sei muito bem como começar, vamos usar o modelo de Erdos para a simulação da rede de interações. Alias Paul Erdos era um cara muito interessante, diria até fascinante quanto ao seus estilo de vida, vale a pena uma olhadinha na história dele.

Mas por enquanto vamos olhar as redes que criamos onde bolinhas verdades são espécies de plantas e quadradinhos azuis são espécies de polinizadores.

01

Bem, eu fiz 4 redes, todas tem 15 polinizadores (linhas) e 15 espécies de plantas (colunas). As plantas dependem dos polinizadores para reprodução, e os polinizadores também dependem das flores para viver (de forma geral, existem excessões). Isso nos leva a conclusão óbvia que todo mundo que a gente ta vendo ai tem pelo menos 1 interação, cada polinizador visita pelo menos um espécie de flor.

Mas além disso podemos pensar em alguma possibilidades pensando nos polinizadores. Uma possibilidade é o polinizador ser generalista, visitar várias flores. Isso é bom já que ele deve encontrar flores mais facilmente, porém ele tem que estar preparado para ser um “pau para toda obra” já que cada flor tem uma forma (morfologia), as vezes disponibiliza seus recursos de uma forma diferente umas das outras, o que pode tornar esse polinizador meia sola em muitas espécies de plantas.

Outra estratégia que pode ser interessante se pensarmos bem é ser especializado. Ja que melhor que pegar só um pouquinho de comida de porta em porta que se bate, é ganhar um montão de comida em uma porta. Mas diferente do generalista, o especialista só consegue o que quer de uma única espécie de planta.

Em termo de energia ou tanto de comida que você consegue coletar, as duas estratégias podem ser equivalente, já que passar bastante tempo voando para pegar bastante pólen em uma única flor ou voar menos e ir pegando um pouquinho de pólen em cada flor pode dar na mesma no final das contas, mas existe outra coisa a se levar em conta, um polinizador generalista pode chegar em novos ambientes facilmente, já que ele consegue comida em varias espécies, enquanto o especialista só vai chegar no ambiente depois que sua “fonte de comida” chegou.

Até agora só consideramos comida mas podem ter flores outras utilidades para os polinizadores e existem outros tipos de interações. Existem ainda os interesses das plantas em todo esse jogo relações.

Mas talvez antes de ficar pensando em possibilidades e mais possibilidades, casos e mais casos, vamos olhar como é essa rede de interações como um todo.

Vamos olhar primeiro o exemplo 1 da perspectiva dos polinizadores.
Nele a gente vê bastante espécies de polinizadores especialista, que só visitam uma espécie de planta, como 2 casos de exclusividade polinizador e planta no meio do gráfico, 3 casos de grupinho envolvendo 2 polinizadores e 2 plantas e um grupão maior no canto inferior direito mais ou menos.

Mas a gente pode contar quantas plantas cada polinizador visita e fazer o seguinte gráfico;

02

Então vamos la, no eixo x temos o número de interações e no eixo y temos a contagem de polinizadores. Ou seja exitem 15 polinizadores, e como todo mundo interage com pelo menos uma planta (o obvio, senão nem na matriz de interação eles estariam), nos temos que com 1 planta os 15 polinizadores interagem, desses 15 polinizadores, apenas 7 interagem com 2 espécies de plantas e um apenas interagiu com 3 plantas.

Certo mas esse gráfico vai ser dependente do número de polinizadores no sistema, mas podemos resolver esse problema transformando tudo em porcentagem, ou seja dividir pelo número de polinizadores total. Então ficamos com…

03

Exatamente o mesmo gráfico, mas notem que agora o eixo X vai de 0 a 1, ou seja é independente do número de polinizadores, isso facilita se quisermos comparar características que vemos nesse gráfico com outras comunidades, independente do número de espécies de polinizadores. Mas ao invéz de barrinhas de contagem, podemos simplesmente usar bolinhas e ficamos com isso.

04

Agora sim, bolinha é bom para ajustar modelos, mas na verdade ainda estamos olhando a mesma coisa do primeiro gráfico de barras, apenas alguns ajustes aqui e ali.

Nesse caso um dos modelos que ajuda a descrever essa comunidade segundo o Jordano (2003) é o exponencial, é um modelo com dois parâmetros que ele chama de b e gamma. Gamma é a letra grega  gamma

E o modelo é assim:

 p(k) = b cdot exp(-gamma cdot k)

Mas notem que no artigo o Jordano (2003) não cita o b, que não fala muita coisa, mas ele não usa o símbolo de igual, ele usa o símbolo de aproximadamente, acho que para não citar o b, não entendo muito o porque disso, mas vamos em frente, é vivendo e aprendendo.

Bem podemos ajustar esse modelo não linear no R usando a função nls(). E ajustando temos…

05

Agora que porra é essa? Bem esse modelo (a linha vermelha) ta nos dizendo como muda a distribuição acumulada de interações, ou seja, numa comunidade onde todo mundo é especialista como essa do exemplo 1, as barrinhas dos gráficos la em cima descem muito rápido, formam uma escada muito íngreme, de degraus altos, enquanto tivessem mais generalistas, a escadinha seria bem menos íngreme e subir essa escadinha seria fácil. E é isso que o gamma do modelo ta falando.
Podemos fazer o mesmo procedimento para os 4 exemplos de comunidade e ver o gamma neles.

coef(m1.exe1) gamma b 0.9372338 2.5877177 coef(m1.exe2) gamma b 0.5020859 1.6452600 coef(m1.exe3) gamma b 0.2602011 1.3744055 coef(m1.exe4) gamma b 0.1417438 1.2410299

E sem olhar para o b por enquanto, mas somente o gamma que nos interessa agora, sabemos que a comunidade 1 é a que tem mais especialistas e menos generalistas, enquanto a comunidade do exemplo 4 tinha muitos mais generalistas e poucos especialistas, olhem o emaranhado que são as interações no exemplo 4.

Isso se reflete nos valores de gamma, olhem para o exemplo 4, ele é bem próximo de zero, bem baixinho 0.14, mas o exemplo 1, cheio de especialistas e quase nenhum generalista forte temos um valor de gamma de 0.93, quase 1.

Ou seja valor alto de gamma é uma escada íngreme naquele gráfico de barras, sinal de especialistas na comunidade, muitos especialistas.
Valor baixo de gamma é sinal de mais generalistas.
Agora temos um pouco de intuição sobre o que o modelo esta ajustando não?

Mas o Jordano citou três modelos no artigo dele, ajustar o modelo desse jeito vai falhar em um caso que imagino. No caso de todo mundo ser generalista. Isso eu imagino pois a gente pode tentar olhar o que acontece com as curvas de acordo com a mudança de gamma.
Como olhando a formula não me traz muitos insights pois sou péssimo de matemática, eu resolvo ela para tudo que é valor de gamma e fico olhando o que acontece.

06

E olha ai, vemos como quanto maior o valor de gamma, mais íngreme é a subida, e biologicamente isso é mais gente com 1 interação apenas, ou seja, especialistas dominam na comunidade. Mas o gamma vai ser no máximo 1, quando ele for 1 o modelo vai ser uma reta inclinada, e no caso de somente haver generalistas, precisaríamos de uma reta não inclinada, ou seja esse modelo ia falhar. Mas o Jordano foi legal deixou outras opções. Mas investigaremos elas depois 🙂

Tudo muito legal mas foram um milhão de comandos, eu não sei fazer tudo isso denovo. Tire a preocupação do seu coração que o senhor Carsten foi muito gentil em resumir todo esse processo numa função do pacote bipartite. La é só digitar degreedistr(sua comunidade) que toda essa magia se realiza para os polinizadores e para as plantas de uma vez so, apesar que não abordamos plantas aqui mas a idéia continua a mesma.

O Resultado vai sair assim e com um gráfico que podemos comparar com o que fizemos até aqui, mas lembrando que aqui estamos olhando para o exemplo 4, a função do bipartite não vai te deixar ajustar o modelo se vc tiver somente 3 barrinhas como nos temos no exemplo 1, o que é fato ja que são poucos graus de liberdade, mas que dificultaria nossa didática aqui.

07

degreedistr(exe4, plot.it=T) $`lower trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1417439 0.01486731 1.210787e-05 0.9717753 -18.632825 power law 0.4722450 0.09654721 1.206915e-03 0.8745991 -4.366322 truncated power law -0.3072618 0.12245915 4.045036e-02 0.9856348 -23.467338 $`higher trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1377984 0.01249272 1.573018e-06 0.9746622 -23.366621 power law 0.4814095 0.08664691 3.536921e-04 0.8791545 -7.092832 truncated power law -0.2837947 0.11314399 3.646958e-02 0.9863076 -28.219390

E lembrando que polinizadores que eram nossas linhas são o “lower trophic level” nesse caso, e olhem que o modelo exponential tem um estimate de 0.14…

coef(m1.exe4) gamma b 0.1417438 1.2410299

E nosso estimate é o mesmo valor, 0.14, ohhhhh, será que porque eu copiei o procedimento dele? 🙂

Bem é isso, fica no script todo o procedimento que fiz para gerar as figuras e a estimativa, depois eu preparo mais posts com os outros modelos, principalmente que essa propriedade tem outro grande interesse em biologia, se um gamma maior significa mais especialistas, um gamma maior também traz algo a respeito da susceptibilidade da comunidade a extinção em cascata.

#Semente para os resultados serem reproduziveis
set.seed(123)
###################################################
#Função para gerar rede de interação baseado no modelo de Erdős–Rényi
#Creditos a Edwin Grappin, vi o post no R-bloggers
#http://probaperception.blogspot.com.br/2012/12/have-you-tried-to-understand-your.html
###################################################
generateER = function(n = 100, p = 0.5){
map = diag(rep(1, n))
link = rbinom(n*(n-1)/2, 1,p)
t = 1
for(j in 2:n){
for(i in 1:(j-1)){
map[i,j] = link[t]
t = t + 1
}
}
dimnames(map)<-list(paste("Pol",1:n),paste("Sp",1:n))
return(map)
}
###########################
#Minha função para transformar matriz em lista para o igraph
##############################
mat2lis<-function(matriz) {
lista<-dimnames(matriz)
dados<-data.frame(Pol=rep(lista[[1]],length(lista[[1]])),
SP=rep(lista[[2]],each=length(lista[[2]])),
int=c(matriz))
return(dados[dados$int==1,1:2])
}
############################
library(igraph)
 
#Gerando exemplos
exe1<-generateER(n=15,p=0.1)
exe2<-generateER(n=15,p=0.2)
exe3<-generateER(n=15,p=0.4)
exe4<-generateER(n=15,p=0.8)
 
#Transformando na forma de dados que o igraph usa
exe1.l<-mat2lis(exe1)
exe1.g <-graph.data.frame(exe1.l, directed=F)
 
exe2.l<-mat2lis(exe2)
exe2.g <-graph.data.frame(exe2.l, directed=F)
 
exe3.l<-mat2lis(exe3)
exe3.g <-graph.data.frame(exe2.l, directed=F)
 
exe4.l<-mat2lis(exe4)
exe4.g <-graph.data.frame(exe4.l, directed=F)
 
#Plots das redes
par(mfrow=c(2,2))
 
plot(exe1.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 1")
 
plot(exe2.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 2")
 
plot(exe3.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 3")
 
plot(exe4.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.5,
vertex.label.color="black"
)
title("Exemplo 4")
 
#Ajustando o modelo exponencial para polinizadores
#Jordano, P., Bascompte, J. and Olesen, J. M. (2003)
#Invariant properties in coevolutionary networks of
#plant-animal interactions. Ecology Letters 6, 69–81
 
#Rede do Exemplo 1
#Número de interações
k<-sum(exe1)
#Número de interações por espécie de polinizador
pol.n<-rowSums(exe1)
#Número de polinizadores por quantidade de interação
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
#Juntando os dados
dados.pol1 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
#Ajustando modelo exponencial aos dados
m1.exe1<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol1)
 
#Grafico 2
barplot(P.pol,names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 3
barplot(dados.pol1[,2],names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 4
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
 
#Grafico 5
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
curve(coef(m1.exe1)[2]*exp(-coef(m1.exe1)[1]*x),1,5,add=T,lwd=2,lty=2,
cex=1.5,col="red")
 
#Exemplo 2
k<-sum(exe2)
S<-sum(dim(exe2))
pol.n<-rowSums(exe2)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol2 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe2<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol2)
 
#Exemplo 3
k<-sum(exe3)
pol.n<-rowSums(exe3)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol3 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe3<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol3)
 
#Exemplo 4
k<-sum(exe4)
pol.n<-rowSums(exe4)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol4 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe4<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol4)
 
#Extraindo coeficientes
coef(m1.exe1)
coef(m1.exe2)
coef(m1.exe3)
coef(m1.exe4)
 
#Niveis de Gamma
graus<-seq(0.1,0.9,0.1)
 
#GRafico 6
jpeg("06.jpg")
plot(0,0,type="n",xlim=c(0,10),ylim=c(0,1),frame.plot=F,
xlab="Número de interações",
ylab="Distribuição Cumulativa do número de interações")
for(i in 1:length(graus)) {
curve(1*exp(-graus[i]*x),0,10,add=T,col=i,lty=2,lwd=2)
}
legend("topright",lty=2,lwd=2,col=c(1:length(graus)),legend=paste("Gamma é ",graus))
dev.off()
 
#Grafico e ajuste usando o bipartite para o exemplo 4
library(bipartite)
degreedistr(exe4, plot.it=T)

Leave a Reply

Your email address will not be published. Required fields are marked *