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)

Todas as estradas terminam na filosofia, mas passam pela matemática

01

As vezes da para se perguntar porque estudar matemática, porque abstair o mundo a formulas. Bem aqui não vamos falar o porque, mas como isso acontece queiramos ou não.

Uma coisa que acontece no Wikipédia legal é que se a gente observar o rede de interações, onde cada artigo é um node e ele tem ligação com outros artigos, para mais de 95% dos artigos, se você começar a clicar no primeiro artigo que ele esta linkado, você sempre vai acabar em filosofia. Mas algo que também é legal, é que você vai passar por matemática. Isso começando de Batman a Batata ou Pólvora.

Isso não tem muita explicação, talvez porque no texto inicial de todo artigo, sempre o foco é justificar o artigo, o que leva a cada vez coisas mais basicas, até fechar em filosofia, e ocorre mais claramente nos artigos de língua inglesa.

Essa descoberta teve uma certa expressão, que deve ter começado como uma piada num site de humor, mas teve grandes desdobramentos, virando noticia em jornal, ganhando o proprio artigo no wikipedia, inclusive um cara fez esse aplicativo legal, que testa o caminho até a filosofia.

Gráficos de Rede no R com o Pacote Igraph

Quando fui olhar pela primeira vez como fazer gráficos de rede no R, a primeira parada foi o pacote bipartite.
O pacote bipartite além de calcular as métricas e estatísticas mais comuns para entender como redes funcionam, tem duas opções para gráficos, as funções plotweb e visweb.

Então quando temos nossa matriz de interações de forma pouco visual:

> table(dados) flor pol Flor 1 Flor 2 Flor 3 Flor 4 Flor 5 Polinizador 1 0 0 1 1 0 Polinizador 2 0 1 0 1 0 Polinizador 3 1 0 0 0 0 Polinizador 4 3 0 0 0 2 Polinizador 5 0 0 2 1 2

Onde o valor é o número de vezes que observamos determinada interação. Assim fica difícil responder facilmente qual polinizador é mais especialista, qual é mais generalista e assim vai. Mas a função plotweb() resolveu bem esse problema.

Mas, apesar de essa figura ficar muito bonita, eu senti dificuldade de conseguir, por exemplo, grudar uma filogenia na rede, já que a posição dos polinizadores não é simétrica, e não é possível controlar de forma simples essa característica. Mas o pacote bipartite tem como dependência o pacote igraph, e fui olhar do que se tratava.

O pacote igraph é exatamente especifico para plotar redes, além de cálculo de métricas e “coisas” de rede. Ele é na verdade uma biblioteca de funções e esta implementado para mais de uma linguagem de computador, como por exemplo para python também.
Mas ele trabalha com uma classe própria de dados, que não é difícil de entender. A matriz acima nele seria representada mais ou menos assim:

> unique(dados) pol flor 1 Polinizador 4 Flor 5 2 Polinizador 1 Flor 3 3 Polinizador 5 Flor 5 4 Polinizador 4 Flor 1 5 Polinizador 2 Flor 4 7 Polinizador 5 Flor 4 8 Polinizador 2 Flor 2 10 Polinizador 5 Flor 3 11 Polinizador 1 Flor 4 13 Polinizador 3 Flor 1

Ou seja, quem interage com quem. Então convertemos esses dados para a classe graph usando a função graph.data.frame() nos dados organizados dessa forma e estamos prontos para fazer os gráficos.

Podemos dar plot nos dados e vemos o que acontece…

Ele faz um gráfico de rede, espalhando todas as plantas e polinizadores e ligando conforme as interações. Na verdade ele tem varias formas de espalhar os nodes, que no nosso caso são polinizadores e flores, e isso é muito interessante pois proporciona varias possibilidades de exaltar varias características da rede, como quem são os especialistas, como são as divisões de flores e polinizadores, além disso é possível modificar formas, cores e muitas outras opções.

Bem mas são funções que fazem esse posicionamento, então podemos ao invés de fornecer uma função, podemos simplesmente dar o posicionamento que desejamos, por exemplo filas emparelhadas como na função plotweb(), claro que a função plotweb ele tem algumas opções de reorganizar a ordem para diminuir a sobreposição entre os riscos que representam as interações. Mas de forma simples podemos fazer no igraph algo similar ao plotweb().

A vantagem que eu achei em usar o igraph é em controlar o tamanho dos nodes, dos símbolos que representam polinizadores e flores, e assim fica mais simples a tarefa de por exemplo colar a filogenia no gráfico de rede, tarefa que não consegui fazer usando o pacote bipartite.

library(bipartite)
library(igraph)
 
#gerando dados de exemplo
set.seed(15)
pol.s<-paste("Polinizador",1:5)
flor.s<-paste("Flor",1:5)
dados<-data.frame(pol=sample(pol.s,15,replace=T),flor=sample(flor.s,15,replace=T))
 
#matriz de interações
table(dados)
 
#Rede com bipartite
plotweb(table(dados))
 
#Convertendo dados
dados.g <-graph.data.frame(unique(dados), directed=F)
 
#Figura 2
plot(dados.g,vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)),vertex.size=25)
 
#Figura 3
par(mfrow=c(2,2))
plot(dados.g, layout=layout.circle,main="Círculo")
plot(dados.g, layout=layout.sphere,vertex.shape="square",main="Esfera")
plot(dados.g, layout=layout.fruchterman.reingold,vertex.color="red",main="Fruchterman Reingold")
plot(dados.g, layout=layout.fruchterman.reingold.grid,
main="Fruchterman reingold Grid",edge.color="red")
 
#Criando posicionamento dos nodes
posicao<-matrix(c(rep(1:5,2),rep(1:2,each=5)),ncol=2,nrow=10)
 
#Rede similar ao bipartite
plot(dados.g, layout=posicao,vertex.color=rep(c("green","yellow"),each=5),
vertex.shape=rep(c("square","circle"),each=5),vertex.size=23,
vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)))
 
#Ligando rede com filogenia
library(ape)
 
par(mfrow=c(2,1),mar=c(1,9,1,9))
plot(rtree(5),show.tip.label=F,direction="downwards",use.edge.length = F,edge.width = 2)
title("Filogenia")
par(mar=c(0,0,0,0))
plot(dados.g, layout=posicao+1,vertex.color=rep(c("green","yellow"),each=5),
vertex.shape=rep(c("square","circle"),each=5),vertex.size=35,
vertex.label=c(paste("Pol",1:5),paste("Flor",1:5)),edge.color="black",
edge.width=3)