Crescimento populacional com atraso para o efeito de dependência da densidade

Voltando ao crescimento populacional, tópico que tem alguns posts aqui já, so dar uma olhada nessa busca aqui do blog.

Mas quando a gente olha o crescimento populacional em modelos como o do crescimento logístico.

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t} )

Lembrando que aqui, o alpha (\alpha) é o efeito de um indíviduo sobre a população, que é a mesma coisa que um sobre a capacidade suporte, ( k= {1\over \alpha} ).

Bem podemos escrever essa formula como normalmente é vista nos livros, com capacidade suporte ao invés de efeito da população sobre ela mesma.

  {dN\over dt} = r\cdot N\cdot (1-{N\over K})

Uma coisa legal que eu ouvi, é que a capacidade suporte não precisa ser pensada somente no efeito da própria população nela mesma, como falta de espaço de uma população de alga para crescer numa placa de petri, mas podemos pensar nela também como um proxy de outras condições que mudam no ambiente, como a “deterioração do ambiente”, que tem seu valor próprio, mas pense que esse valor tem uma relação direta com o tamanho populacional, ele ferra a população, logo ao invés de medir ele, podemos usar o próprio tamanho populacional, que representa esse valor, e simplifica em muito a modelagem, diminuindo uma variável por exemplo.

Mas ok, o ponto que queremos chegar é que aqui

 (1-{N\over K})

é a hora que acontece a diminuição do crescimento de acordo com o tamanho da população, mas isso é baseado no tamanho populacional aqui e agora.
Mas e se, pensando numa população de ratinhos na mata, quando o ratinho tem sua ninhadinha, logo que ele teve os filhos, não temos ainda o efeito da deterioração do ambiente por exemplo, esses filhotes tem que crescer, e viver e se reproduzir para deteriorar o ambiente, mas na vida deles eles não sentem essa deterioração, quem sente isso são seus descendentes. Quem vai lidar com os problemas causados são os filhos, netos e bisnetos. Quando agora os ratinhos estão la desfrutando da fartura e comendo todas as sementes, quem vai se ferrar são as próximas gerações, que se todas as sementes são comidas, não serão germinadas novas arvores e arbustos que deveriam fornecer comida para outras gerações.
Agora como será que isso impacta a dinâmica populacional?

Bem vamos fazer uma simulação aqui para tentar ver o efeito de dependência da densidade atrasado (Delayed density dependence).

Então a gente precisa dos parâmetros populacionais que usamos no modelo logístico:

1
2
3
4
5
N0 <- 1 #População inicial
r  <- 0.1 #Taxa de crescimento Populacional
K  <- 100 #Capacidade Suporte
t  <- 200 #Número de gerações
tau<- 20  #Atraso de resposta para o efeito dependente de densidade

E mais um parâmetro especial, que é o tau, ou o lag, o quanto vai demorar para a população sentir o efeito do tamanho populacional de agora.

E vamos implementar uma simulação discreta da seguinte forma:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
for(i in 1:length(generation)) {
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        } else {
            if(i <= tau) {
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }
        }
    }

Bem temos dois if nessa simulação discreta (para noção da diferença entre discreto e continuo, olhe aqui e aqui).

O primeiro, aqui

1
2
3
4
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        }

A gente so está tratando da geração 0.

1
2
3
4
5
6
        } else {
            if(i <= tau) {
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } 
          }

Depois aqui a gente ta tratando do inicio da dinâmica, lembre-se que falar que temos um lag de 10 gerações, ou 10 tempos, significa que o tamanho populacional de agora afeta a taxa de crescimento de daqui 10 gerações, so que não da para isso acontecer na primeira geração, nem na segunda, nem na decima, porque não tem 10 gerações para trás ainda, então temos que dar um jeito nisso, que aqui é usar o crescimento logístico sem lag nesse inicio.

1
2
3
4
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }

E por ultimo, a gente trata do caso quando começamos a poder implementar o lag, no caso do lag 10, a partir do momento que ja temos 10 gerações para trás.

Veja que a diferença está somente aqui na formula

 (1-{N\over K})

O efeito da diminuição da taxa de crescimento. A gente usa a população do tamanho lag para traz

 (1-{N_{t-lag}\over K})

ou seja

  N_{t+1} = r\cdot N_{t} \cdot (1- \alpha * N_{t-lag} )

A taxa de crescimento, vezes o tamanho populacional agora, mas o efeito na diminuição da taxa de crescimento é da população t-lag gerações atrás.

Bem eu vou colocar tudo isso numa função para facilitar a simulação, que é o que vai estar no script la em baixo, aqui a gente so estava dando uma olhada.

Mas o resultado disso, para um lag de 10 gerações fica o seguinte, sobre a dinâmica populacional.

01

A população consegue passar a capacidade suporte, no inicio, quando a responsabilidade bate, e vem o efeito dessa população alta, a taxa de crescimento fica negativa e ela cai, depois com a população abaixo da capacidade suporte, ela sobe denovo, mas esses altos e baixos vai diminuindo, até haver uma estabilização na capacidade suporte.

A gente pode olhar a taxa de crescimento de acordo com o tamanho populacional

02

Veja que legal, eu achei essa figura muito legal, a gente fica rodando perto taxa de crescimento zero, que é o que acontece quando a gente chega na capacidade suporte, veja que quando mais próximo da estabilidade, mais perto de parar em zero a gente fica rodando. Algo similar ao que acontece com o crescimento per capita, a contribuição que cada individuo da população pode fornecer.

03

Certo então temos o efeito dependente de densidade atrasado na dinamica populacional.

04

Que ela vai ter um momento inicial diferente do crescimento logístico, mas vai tender a uma estabilidade também, na capacidade suporte do ambiente.

Mas será que isso vale para todos os valores de atraso? Bem se a gente diminuir o lag para 3 por exemplo.

05

A gente vai se aproximar do crescimento logístico, o comportamento do crescimento logístico, onde a gente só diminui a taxa de crescimento per capita e temos um crescimento grande no começo e depois que chegamos a metade do tamanho populacional, começamos a ter uma queda na taxa de crescimento.

Mas e se a gente aumentar bastante o lag? Algo como lag = 20.

06

Que legal, a gente entra num ciclo, onde a população sobe muito acima da capacidade suporte, e depois cai bruscamente quando esse efeito vem a tona, mas depois como fica com um tamanho pequeno, ela vai e cresce bastante denovo, e fica girando nesse ciclo, diferente de um lag médio quando ela oscilava no começo mas depois estabilizava.

Bem é isso ai, eu estou vendo umas coisas legais em um curso do Edx chamado Quantitative Biology Workshop, esse modelo estava perdido no meio de um exercício, mas eu achei muito legal. No começo o curso estava meio chato, ficar fazendo exercícios em matlab (eu que nunca tenho grana sobrando, quiça para comprar o matlab) não me atrai muito, mas eles tem exercícios com simulações bem legais, e é apresentado aula dos professores que dão esse curso la no mit no final da semana, que é bem legal de ver também, mas então, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#Post
#http://recologia.com.br/2015/01/crescimento-populacional-com-atraso-para-o-efeito-de-dependencia-da-densidade
 
#Parâmetros de entrada
N0 <- 1 #População inicial
r  <- 0.1 #Taxa de crescimento Populacional
K  <- 100 #Capacidade Suporte
t  <- 200 #Número de gerações
tau<- 20  #Atraso de resposta para o efeito dependente de densidade
 
simu_loglag<-function(N0=1,r=0.1,K=100,t=200,tau=10) {
    generation = seq(0,t,by=1)
    dNdt = rep(0,0,t)
    population = rep(0,0,t)
 
    for(i in 1:length(generation)) {
        #Aqui é o inicio da simulação, como não temos nada de população
        #ainda, usamos a população inicial
        if(i == 1){
            population[i] <- N0
            dNdt[i] <- r*population[i]*(1-population[i]/K)
        } else {
            if(i <= tau) {
                #Aqui, é o caso onde o lag ainda é maior que o tamanho da
                #populacao, lembre-se que se a gente tem um lag de 10, até o
                #tempo 10, não da para olhar 10 tempos atrás, diminuir 10 do
                #i, então a gente tem que usar o crescimento sem lag no
                #começo
                population[i] <- population[i-1]+dNdt[i-1]
                dNdt[i] = r*population[i]*(1-population[i]/K)
             } else {
                 population[i] = population[i-1]+dNdt[i-1]
                 dNdt[i] = r*population[i]*(1-population[i-tau]/K)
             }
        }
    }
    saida<-data.frame(generation,population,dNdt)
    return(saida)
}
 
 
simu_lag10<-simu_loglag(tau=10)
 
#
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag10,type="l",lwd=2,frame=F)
 
#
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag10,type="l",lwd=2,frame=F)
 
#
plot((dNdt/population)~population,data=simu_lag10,type="l",lwd=2,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita")
 
#Lag = 10
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag10,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 10")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag10,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag10,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")
 
#Lag = 3
simu_lag3<-simu_loglag(tau=3)
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag3,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 3")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag3,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag3,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")
 
#Lag = 20
simu_lag20<-simu_loglag(tau=20)
layout(matrix(1:3,ncol=1,nrow=3))
plot(population~generation,xlab="Tempo",ylab="Tamanho da população",
     data=simu_lag20,type="l",lwd=1,frame=F,
     main="Tamanho populacional ao longo do tempo, Lag = 20")
plot(dNdt~population,xlab="Tamanho da população",
     ylab="Crescimento populacional",data=simu_lag20,type="l",lwd=1,frame=F,
     main="Crescimento vs Tamanho populacional")
plot((dNdt/population)~population,data=simu_lag20,type="l",lwd=1,frame=F,
     xlab="Tamanho da população",ylab="Crescimento populacional per capita"
     ,main="Crescimento per capita vs Tamanho populacional")

Métodos comparativos filogenéticos – Simulando o movimento brownian

Muito tempo atrás, eu comecei a olhar sobre os Phylogenetic comparative methods (PCMs), neste post aqui mais especificamente. Bem, quando realmente tentei usar o tal dos métodos comparativos filogenéticos, vi que tudo era muito mais complicado. Dai o jeito foi começar a tentar simular dados (tentando ir na onda do livro do Marc Kery, que é muito legal), para entender melhor como tudo isso ta funcionando de verdade.

Se a gente olhar como é o código da função compar.gee, a gente vai ver que ela pega a árvore que a gente manda para a análise e constrói uma matriz de correlação entre as espécies, usando a função vcv do pacote ape, e essa função faz a matriz de correlação entre as espécies baseada no modelo browniano (único disponível na função até agora, tenho medo de pensar nos outros modelos, se esse ja ta difícil entender).

Então vamos tentar entender melhor esse primeiro passo. A ideia é a seguinte. O movimento browniano é o movimento aleatória de partículas suspensas em um fluido, resultado da colisão com átomos ou moléculas que formam o fluido, seja um gás ou líquido. Na série Cosmos, que eu vi recentemente a versão com o Neil deGrasse Tyson, diz foi assim que veio a sacada que os átomos estavam la, mesmo sem a gente ver, mas só assistindo, porque eu não sei explicar. Mas o termo “Movimento Browniano” também se refere ao modelo matemático usado para descrever movimentos aleatórios (tem esse post aqui pra dar uma olhadinha também), que aqui vem aplicado a evolução.

Mas como é aplicado a evolução?

Primeiro vamos fazer uma simulação de como ele está acontecendo, num contexto evolutivo. Vamos pensar numa característica (fenótipo), que por exemplo não está sofrendo nenhuma seleção, não está sob pressão seletiva, sei la tamanho de sapinhos (é um exemplo). Essa característica pode variar ao longo do tempo, tempo medido em gerações, e como não existe pressão seletiva, pode variar tanto pra mais como para menos. Isso vai seguir a ideia do movimento browniano em uma dimensão apenas, então podemos fazer isso assim:

1
2
3
4
5
6
7
8
#Tempo
t <- 0:100
#Variação
sig2 <- 0.01
## Primeiro simulamos um conjunto de desvios aléatorios
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
## Depois nos computamos a soma comulativa
x <- c(0, cumsum(x))

A gente seleciona um número de gerações (100 nesse caso, veja que o zero é o ponto de início) que queremos simular, estabelecemos o quanto de variação é possível por geração, depois sorteamos aleatoriamente a variação por geração (usando rnorm, veja que o n tem o -1 porque a gente vai começar do zero) e depois vemos como essa variação acumula, usando o cumsum, que é uma função que faz a soma acumulada de um vetor, ou seja soma o primeiro com o segundo, esse resultado com o terceiro elemento, o resultado da soma dos 3 primeiros com o quarto elemento, e devolve um vetor com essa soma, vetor esse que vai ter o mesmo tamanho do anterior.
Agora podemos dar uma olhada melhor no que estamos fazendo.

01

E é isso ai, o pontinho vermelho é onde a gente começa a acumular a variação, e o valor que a gente observa depende sempre da geração anterior, lógico não? Seu descendente vai ser um pouquinho maior ou um pouquinho menor. Assim você vai poder chegar a qualquer altura no final da simulação, salvo o limite que você pode variar (é pouco provável que depois do ponto de saída, você rapidamente fique gigante, ou diminuto, mas pode acontecer, mas é pouco provável).

Agente pode pensar melhor sobre isso tentando fazer muitas simulações como essa e vendo o que acontece.

02

Veja que ao avaliar várias possibilidades, o que a gente vai ver é que todo tipo de coisa está acontecendo (estamos vendo todo tipo de tamanho), mas se a gente olhar como é a distribuição de tamanhos na última geração, veja quem reconhece essa forma!

03

Isso mesmo, teremos uma distribuição normal de tamanhos, lembrando que estamos pensando aqui da seguinte forma, pensa que essa é a variação, ou seja pense numa média, e some esse valor, e terá os tamanhos, ou que os valores que observamos aqui são valores diminuídos do tamanho médio, bem a ideia é só representar a variação.
Mas se a gente olhar essa variação ao longo de gerações…

04

A gente pode ver que a variação está aumentando ao longo das gerações (o tempo), isso porque quanto mais tempo a gente tem, mais é possível acumular variação, mas veja que a média, o centro da distribuição, tende a permanecer no mesmo lugar.

E se a gente usar uma variação menor, teremos um crescimento mais lento da variação.

06

Mas de qualquer forma, sempre esperamos uma relação onde a variação aumenta com o tempo.

07

Veja que até ai não entramos muito bem em o que isso tem haver com a correção feita nos métodos comparativos filogenéticos, ou a correção que tentamos fazer.
Mas agora vamos la, vamos pensar em uma determina árvore filogenética, que nos diz como as espécies foram separando. Com o pacote phytools a gente gera uma árvore aleatória.

1
2
3
4
5
6
library(phytools)
t <- 100  # Tempo Total
n <- 30  # Número de espécies
b <- (log(n) - log(2))/t
tree <- pbtree(b = b, n = n, t = t, type = "discrete")
tree$tip.label <- paste("Sp",sub("t","",tree$tip.label))

Aqui então pensamos no tempo total, em gerações discretas, mas isso pode ser contínuo também. Bem uma especiação discreta é na verdade pouco realista, mas é mais fácil para fazer as contas. Falamos o número de espécies final da árvore e b é a chance de especiação por geração, para o modelo discreto. Esse valor é bem complicado de estabelecer, porque pense se a gente coloca um número muito grande (próximo de 1), a gente vai “gerar” todas as espécies, fazer todas as ramificações logo no início da árvore, e ficar o resto das gerações sem simular nada, por outro lado, se a gente pega um valor muito baixo, próximo de zero, a gente chega a 100 gerações sem o número de espécies que queremos, assim a simulação tem que suar a camisa pra dar conta de resolver isso dependendo do número b escolhido, mas vamos la, então geramos nossa árvore com 30 espécies.

08

Bem agora podemos usar ela para fazer uma simulação do movimento browniano para essas 30 espécies, mas agora vamos levar a árvore em consideração, como assim? Bem vamos fazer uma figura e olhar pra ela que vai dar pra entender melhor o que está acontecendo.

09

Aqui que está a jogada, pense que no inicio, só existia uma espécie, e no tempo zero, que é quando começamos a observar, ela especiou e então começamos a acompanhar duas espécies diferentes. E somente essas duas espécies que estão acumulando variação, e não as 30 espécies na verdade, e elas vão para algum lugar na variação, sei la, vão evoluir para ficar maior ou menor, ou se manter na mesma, qualquer coisa pode acontecer. Mas pouco tempo depois, temos 2 novas bifurcações, ou seja, novos eventos de especiação fazem essas 2 espécies virarem 4, só que cada espécie herda da sua espécie ancestral o tamanho que tinha(veja que aqui eu to falando espécie ancestral, não estamos observando a espécie que temos viva ainda, isso ta la no final da árvore, no final do eixo do tempo), assim, se uma espécie estava ficando grande, positiva, seus dois descendentes vão partir desse ponto, e não do zero, como na segunda figura, da segunda simulação que fizemos. Assim, novamente, denovo e denovo vamos tendo eventos de especiação, e criando os ramos da árvore, só que sempre o acumulo de variação parte do ponto que o ancestral comum parou. Ou seja, quando vamos pensar nas nossas 30 espécies, primeiro elas não tem 100 gerações de acumulo de variação livre, mas dependente dos ancestrais, agora se a gente olhar a distribuição de como estão as espécies no final.

10

Não temos mais aquela cara de distribuição normal, mas algo mais uniforme, isso porque a espécie tem uma dependência de como estava o tamanho do ancestral, ou seja as ramificações da árvore geram uma dependência entre espécies irmãs, que estão no mesmo clado, no mesmo ramo, uma correlação. Nossa senhora, tudo isso para finalmente chegar a essa palavra, correlação entre espécies.

Aqui fizemos uma simulação discreta, mas o pacote phytools tem uma função chamada fastBM que faz a simulação de variação acumulada devido ao movimento browniano, e alias faz a simulação de forma contínua também, o que é mais difícil mas mais real.

11

Mas o código da simulação pode ajudar a entender melhor o que está acontecendo, veja que a gente tem que ir fazendo a simulação para as 2 primeiras espécies, e ir pegando quando elas se dividem e dai pegar o valor da espécie ancestral e continuar para as filhas. Para o caso contínuo isso fica bem mais difícil.

Mas é isso ai, a primeira dificuldade que eu tive para começar a entender melhor esse tipo de método, como fazer essa correção usando uma árvore, era como a árvore filogenética entrava na parada, na análise. E é assim que ela vai entrar, a gente tentando corrigir essa variação dependente da ramificação da árvore. Talvez até não corrigir, porque da a impressão que a evolução é um problema, é melhor falar que vamos levar em conta na análise essa correlação entre espécies, vamos levar em conta que as espécies, quando mais aparentadas, devem ter valores mais próximos.

Bem, aqui foi uma primeiro olhada pra tentar entender melhor o que está acontecendo, mas vamos ver um pouco mais ainda continuando esse post, grande parte do código eu peguei aqui, mas alterei um pouco os gráficos pra ficar mais fácil de visualizar. Na minha opinião, mas esse exercício ajuda bastante a entender melhor como deve começar a entrar o tal da correção da filogenia nas análise, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#http://www.phytools.org/eqg/Exercise_4.1/
set.seed(123)
 
#Tempo
t <- 0:100
#Variação
sig2 <- 0.01
## Primeiro simulamos um conjunto de desvios aléatorios
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
## Depois nos computamos a soma comulativa
x <- c(0, cumsum(x))
 
#Figura1
jpeg("01.jpg")
plot(t, x, type = "l", ylim = c(-2, 2),frame=F,xlab="Tempo (Gerações)",ylab="Variação fenotípica",lty=3)
points(0,0,pch=19,col="red")
dev.off()
 
#Agora vamos simular varias espécies
nsim <- 100
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2)), nsim, length(t)-1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
 
#Figura 2
jpeg("02.jpg")
plot(t, X[1, ], xlab = "Tempo", ylab = "Variação fenotípica", ylim = c(-2.5, 2.5), type = "l",frame=F,lty=3)
apply(X[2:nsim, ], 1, function(x, t) lines(t, x,lty=3), t = t)
points(0,0,pch=19,col="red")
dev.off()
 
#Figura 3
jpeg("03.jpg")
hist(X[,101],xlim=c(-3,3),breaks =seq(-3,3,0.2),main="Após 100 gerações",freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
dev.off()
 
#Figura 4
jpeg("04.jpg")
par(mfrow=c(2,2))
hist(X[,11],xlim=c(-3,3),breaks =seq(-3,3,0.2),main=paste("Gerações = 10, Desvio Padrão = ",round(sd(X[,11]),2)),freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
hist(X[,26],xlim=c(-3,3),breaks =seq(-3,3,0.2),main=paste("Gerações = 25, Desvio Padrão = ",round(sd(X[,26]),2)),freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
hist(X[,51],xlim=c(-3,3),breaks =seq(-3,3,0.2),main=paste("Gerações = 51, Desvio Padrão = ",round(sd(X[,51]),2)),freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
hist(X[,101],xlim=c(-3,3),breaks =seq(-3,3,0.2),main=paste("Gerações = 101, Desvio Padrão = ",round(sd(X[,101]),2)),freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
dev.off()
 
#Simulação usando for
X <- matrix(0, nsim, length(t))
for (i in 1:nsim) {
    X[i, ] <- c(0, cumsum(rnorm(n = length(t) - 1, sd = sqrt(sig2))))
}
 
#Figura 5
jpeg("05.jpg")
par(mfrow=c(1,1))
plot(t, X[1, ], xlab = "Tempo", ylab = "Variação fenotípica", ylim = c(-2.5, 2.5), type = "l",frame=F,lty=3)
for (i in 1:nsim) {
    lines(t, X[i, ],lty=3)
}
points(0,0,pch=19,col="red")
dev.off()
 
#Menor Variação
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2/10)), nsim, length(t) -   1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
 
#Figura 6
jpeg("06.jpg")
plot(t, X[1, ], xlab = "Tempo (Gerações)", ylab = "Variação fenotípica", ylim = c(-2, 2), type = "l",lty=3,frame=F)
apply(X[2:nsim, ], 1, function(x, t) lines(t, x,lty=3), t = t)
points(0,0,pch=19,col="red")
dev.off()
 
#Variação vs Geração
nsim <- 10000
X <- matrix(rnorm(n = nsim * (length(t) - 1), sd = sqrt(sig2)), nsim, length(t)-1)
X <- cbind(rep(0, nsim), t(apply(X, 1, cumsum)))
v <- apply(X, 2, var)
 
#Figura 7
jpeg("07.jpg")
plot(t, v, type = "l", xlab = "Tempo (Gerações)", ylab = "Variancia entre gerações",frame=F)
dev.off()
 
#
library(phytools)
t <- 100  # Tempo Total
n <- 30  # Número de espécies
b <- (log(n) - log(2))/t
tree <- pbtree(b = b, n = n, t = t, type = "discrete")
tree$tip.label <- paste("Sp",sub("t","",tree$tip.label))
 
#Figura 8
jpeg("08.jpg")
plot(tree)
dev.off()
 
 
## Simulando evolução ao longo de cada vertice (espécie ancestral)
X <- lapply(tree$edge.length, function(x) c(0, cumsum(rnorm(n = x, sd = sqrt(sig2)))))
## Reordenando os vertices da arvore em pre-ordem transversal
cw <- reorder(tree)
## Agora simulando na arvore
ll <- tree$edge.length + 1
for (i in 1:nrow(cw$edge)) {
    pp <- which(cw$edge[, 2] == cw$edge[i, 1])
    if (length(pp) > 0) {
        X[[i]] <- X[[i]] + X[[pp]][ll[pp]]
    } else {
        X[[i]] <- X[[i]] + X[[1]][1]
    }
}
## Pegando o ponto de inicio e fim de cada vertice para a figura
H <- nodeHeights(tree)
 
 
#Figura 9
jpeg("09.jpg")
#Linhas da simulação
plot(H[1, 1], X[[1]][1], ylim = c(-2.5,2.5), xlim=c(0,120), xlab = "Tempo (Gerações)", ylab = "Fenótipo",col="red",pch=19,frame=F)
for (i in 1:length(X)) {
    lines(H[i, 1]:H[i, 2], X[[i]],lty=3)
}
#Pontos internos da arvore onde começa o movimento browniano
for (i in 1:nrow(H)) {
    if(H[i, 2]!=100){
        points(H[i, 2], X[[i]][length(X[[i]])],pch=19,col="red")
    }
}
## Adicionando os nomes das espécies
yy <- sapply(1:length(tree$tip.label), function(x, y) which(x == y), y = tree$edge[,2])
yy <- sapply(yy, function(x, y) y[[x]][length(y[[x]])], y = X)
text(x = max(H)+5, y = yy, tree$tip.label,cex=0.6)
dev.off()
 
#Separando os valores para o histograma dess simulação
saida<-rep(NA,n)
contador<-1
for (i in 1:nrow(H)) {
    if(H[i, 2]==100){
         saida[contador] <- X[[i]][length(X[[i]])]
         contador <- contador + 1
    }
}
 
#Figura 10
jpeg("10.jpg")
hist(saida,xlim=c(-3,3),breaks =seq(-3,3,0.2),main="Após 100 gerações",freq =F
    ,ylab="Densidade",xlab="Classes de variação fenotípica")
dev.off()
 
## Simulando a evolução Browniana em uma arvores usando o fastBM
x <- fastBM(tree, sig2 = sig2, internal = TRUE)
 
#Figura 11
jpeg("11.jpg")
phenogram(tree, x, spread.labels = TRUE, spread.cost = c(1, 0))
dev.off()

A evolução da cooperação e o torneio de Axelrod

A primeira vez que li o livro O gene egoísta, do Richard Dawkins, apesar de muito legal o livro, foi meio triste. Porque com todas aquelas teorias de evolução, do gene egoísta e tal, da a impressão que que evoluímos para todos terem um comportamento egoísta, a luta pela sobrevivência. A gente acaba vendo poucos exemplos onde ser “legal” é eficiente na natureza. Segundo a teoria neodarwinista.

Antes de ler os próximos livros eu fiquei preocupado com isso. No fim do livro o Dawkins, ele ainda fala que uma das características que nos distingue das outras espécies é podermos escapar do “Gene egoísta” porque pensamos, mas ainda sim, isso não era muito reconfortante, isso não explicava muito bem porque vemos tanta gente legal no mundo, essas pessoas não podem ser o “fitness” ruim.

Mas então, lendo os outros livros, como O Relojoeiro Cego, eu li a primeira vez sobre o torneio de Axelrod, e pela primeira vez uma explicação que fazia sentido, em como ser legal pode ser evolutivamente estável.

Muitas perguntas em biologia, como em outras ciências, tem o problema de serem muito difíceis de se fazer hipóteses testáveis. So pensar em evolução, como diabos se testa algo em evolução, separar grupos tratamento e controle e esperar organismos evoluírem nunca é uma proposta viável, dado o tempo que leva para algo acontecer, nossa vida é curta para essa escala de tempo. Talvez de para acompanharmos alguns organismos, como bactérias, para algumas perguntas específicas, mas como testar algo como cooperação, evolução da cooperação.

Um cara chamado John Maynard Smith teve uma sacada. Ele resolveu usar na biologia a teoria dos jogos para testar coisas. Quase todo mundo ja ouviu falar de teoria dos jogos, é o que aquele cara do filme uma mente brilhante estuda, o John Nash.

Mas então, dento desse cenário, o Robert Axelrod, que é um cientista politico e o William D. Hamilton publicaram um artigo na revista science que deu o melhor exemplo, de como a cooperação aparece. Mas não a cooperação como em insetos eu-sociais, tipo formigas, abelhas, etc. Cooperação entre indivíduos independente do seu parentesco.

O que o Axelrod fez, e eu não sei exatamente aonde o W.D. Hamilton entrou na jogado, acho que deve ter entrado para escrever o artigo, analisar os dados, já que o Axelrod fez o Axelrod Tournament. Em teoria dos jogos, as pessoas quantificam jogos, cada jogada tem um custo ou retorno aos seus participantes. Talvez o jogo mais famoso e primeiro exemplo da teoria dos jogos é o dilema do prisioneiro.

O Dilema do prisioneiro é relativamente simples, existem dois caras presos, o cara A e o cara B. Eles são colocados separados em celas diferentes. Dai vem um policial para cada um, e pergunta se ele quer confessar o crime ou não. Então três coisas podem acontecer (na verdade quatro, mas ja vemos isso).

Uma primeira possibilidade é todo mundo é dedo duro, o A entrega o B, e o B entrega o A. Agora o policial tem a confissão dos dois e manda os dois para a cadeia para cumprir dois anos cada.

A segunda possibilidade é quando os dois ficam em silencio, se ambos ficam em silencio, o policial não tem nada contra eles, a não ser as provas que já tinha antes e só consegue colocar cada um por um ano na prisão.

A terceira possibilidade é quanto temos um dedo duro e um cara legal, ai o policial vai ter as provas que tem mais o testemunho do parceiro de crimes, colocando quem foi dedurado por três anos na prisão e liberando o dedo duro, porque ele ajudou a criar provas contra o outro prisioneiro. Mas aqui tanto o A pode ser o dedo duro e o B ficar em silencio como o B ser o dedo duro e o A ficar em silencio, em ambos os casos um se ferra e o outro se da bem.

Agora veja, podemos resumir essas possibilidades em uma matriz onde as linhas são as ações de um prisioneiro e colunas as ações do outro prisioneiro.

 \begin{array}{ccc}  & {Prisioneiro\ B\ fica\ em\ silencio} & {Prisioneiro\ B\ entrega\ o\ parceiro\ A} \\ {Prisioneiro\ A\ fica\ em\ silencio} & {1/1} & {3/0} \\ {Prisioneiro\ A\ entrega\ o\ parceiro\ B} & {0/3} & {2/2} \\ \end{array}

Esses valores podem variar de acordo como a gente quiser, mas a ideia do jogo é essa, temos participantes, aqui o prisioneiros, e cada decisão deles leva a custo, ou ganho, só que esse custo/ganho não depende da ação de um dos prisioneiros apenas, o custo/ganho sempre vai ser dado pela combinação da ação dos dois.

Agora vamos pensar um pouco, se a gente pega somente o prisioneiro A. Vamos supor que o Prisioneiro B ja tenha tomado sua decisão de ficar em silencio, aqui a decisão do prisioneiro A é entre também ficar em silencio e ficar um ano preso ou entregar o B e não ficar nenhum ano preso. Agora vamos supor que diferente disso, o Prisioneiro B tenha escolhido entregar A, agora A enfrenta a decisão de passar três anos presos ao ficar em silencio, ou 2 anos se denunciar B. Veja que para cada um dos participantes, se a gente pensar independente do outro prisioneiro, compensa sempre denunciar o outro prisioneiro, ser o dedo duro, sempre você é o dedo duro, você passa menos tempo na cadeia, é a escolha entre 1 ano ou nenhum ano ou a decisão entre 3 anos ou 2 anos.
Para esses valores, compensa sempre ser dedo duro, nunca cooperar, quem coopera com os outros sempre pode se ferrar bonito, passar 3 anos na cadeia se ficar em silencio, ou se der sorte 1 se o outro prisioneiro também ficar em silencio, mas se cada prisioneiro está em uma cela separada, incomunicável, quem garante que o outro prisioneiro não vai te sacanear? Sendo que ele so tem a “ganhar”.

Seguindo essa analise, nunca esperaríamos que o altruísmo aparece nas espécies. Para testar essa hipótese, o Axelrod teve uma idéia, ele bolou um torneio de computador, onde as pessoas submetiam programas de computador para jogar o dilema do prisioneiro, mas de forma iterada. Você não joga somente uma vez, você joga varias “partidas”, e ao invés de anos de cadeia a gente contabiliza pontos seguindo a seguinte matriz.

 \begin{array}{ccc}                      & {B\ Coopera}   & {B\ nao\ coopera} \\   {A\ Coopera}       & {1/1}          & {2/0} \\   {A\ nao\ Coopera}  & {0/2}          & {-1/-1} \\   \end{array}

Agora a gente tira ponto quanto os dois não cooperam, da um ponto pra cada quando os dois cooperam, e se somente um cooperar e o outro não cooperar, que não cooperou ganha 2 pontos e o outro não ganha nada.

E o jogo vai ser jogado por varias partidas, só que cada jogador sabe o historio do outro jogador com quem já interagiu, então a primeira vez que dois jogadores, vamos chamar de agentes, a primeira vez que eles interagem é novidade, a segunda cada um já lembra o resultado da primeira interação.

Mas além disso, teremos vários agentes, com diferentes estratégias. A gente vai fazer essa simulação do torneio do Axelrod.

Então a gente vai fazer varias funções, que serão os agentes, a usaremos TRUE e FALSE para cooperar ou não cooperar respectivamente.

Mas são essas funções? Vamos pegar um exemplo, a estratégia do ditador agressivo:

# Ditador Agressivo
ag.ditador <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}

O ditador agressivo começa o jogo não cooperando.

> ag.ditador(c()) [1] FALSE

Se ele nunca interagiu com você, ele solta logo um “não coopero” (FALSE). Mas ao longo do jogo, se em média você começar a nunca mais cooperar com ele.

> ag.ditador(c(T,T,F,F,F)) [1] TRUE

Veja que o outro jogador cooperou as duas primeiras rodadas e parou de cooperar, ai o ditador agressivo volta atras, e começa a cooperar com você pra ganhar sua amizade denovo.

> ag.ditador(c(T,T,F,F,F,F,T,T,T)) [1] FALSE

Mas ai a hora que você perdoou o ditador agressivo, e começou a cooperar com ele, ele vai e te ferra denovo e volta a não cooperar, e assim vai.

Essa é uma estratégia, podemos fazer funções para as estratégias mais simples possíveis, como sempre cooperar, ou sempre não cooperar, como estratégias elaboradas, para tentar prever o padrão do outro agente e tomar proveito disso, porque sempre que você não coopera e o outro coopera, você ganha 2 pontos.

E na nossa simulação, como eu fiz poucas estratégias, vamos colocar mais de um agente com a mesma estratégia, vamos usar 20 agentes de 11 estratégias diversificadas, e vamos jogar por 500 partidas.
Basicamente cada partida a gente faz 10 pares de agentes, sorteando os pares ao acaso, e fornece para função o histórico de interação entre esses agentes, sendo somente o historio entre esses agentes, e não todas as interações do agente.

Depois disso transformamos as interações em pontos, e podemos ver a pontuação ao longo do tempo.

01

Veja que na minha simulação, eu tive três agentes usando a estrategia tit4tat, essa estrategia foi a ganhadora do torneio, nas duas edições que aconteceram. Quando o Axelrod fez o torneio, ele pediu para um monte de cientista da área de teoria dos jogos, optimização, para enviarem programas com estratégias pro campeonato, as que estão nessa simulação, não são nem de longe tão complexas como a que as pessoas tentaram. Tinha de tudo, no entanto, a estrategia extremamente simples chamada tit4tat venceu. Além disso, quando acabou a primeira edição, o Axelrod publicou o resultado do torneio e deixou publico todos os programas enviados. Mas na segunda edição, mesmo todo mundo podendo estudar os programas dos outros, estudar o resultado do jogo, podendo bolar agora maneiras de como explorar as estratégias das outras pessoas, o tit4tat venceu denovo na segunda edição.

Mas como funciona o tit4tat?
O Tit for Tat foi feito por um russo chamado Anatol Rapoport. Basicamente ele começa cooperando, e sempre coopera enquanto o outro agente também cooperar, se em alguma rodada o outro agente não cooperar, então ele começa a não cooperar também, mas ele só vai não cooperar até que o outro agente volte a cooperar, a partir dai ele volta a cooperar também. Simples assim.

Vamos lembrar o esquema de pontuação.

 \begin{array}{ccc}                      & {B\ Coopera}   & {B\ nao\ coopera} \\   {A\ Coopera}       & {1/1}          & {2/0} \\   {A\ nao\ Coopera}  & {0/2}          & {-1/-1} \\   \end{array}

Agora pense nessa estratégia contra as mais simples, que é de nunca cooperar e sempre cooperar.
Quando o tit4tat encara um agente que nunca coopera, como na primeira rodada ele vai cooperar e o outro agente não vai cooperar, o outro agente vai ganhar 2 pontos e o tit4tat não ganha nada, dai pra frente, os dois nunca mais vão cooperar entre si, sempre ganhando -1 ponto cada um e a pontuação vai somente despencar, mas ainda sim o outro agente sempre estará 2 pontos a frente do tit4tat, ou seja, sozinho o tit4tat nunca ganha desse agente.
Por outro lado se o tit4tat for jogar com alguém que sempre coopera, os dois vão sempre cooperar e cada um vai ganhar sempre um ponto, por quantas rodadas forem o jogo, isso faz com que o tit4tat consiga empatar com esse agente, mas nunca ganhar.

Mais ou menos da para ver que tit4tat nunca ganha, no melhor caso ele empata, no entanto os três agentes deles estão na frente na nossa simulação.

Acontece que temos que pensar no resultados em grupos, o tit4tat não esta ganhando na verdade, ele começa a se dar bem, quando interage com outros agentes que também cooperam, ganhar um ponto é sempre melhor que perder 1 ponto. Veja que algumas estratégias, pouco depois do inicio começam a despencar na pontuação, por exemplo o laranja, que é o sempre trapaceia cai fácil, porque como os programas tem a memoria de como foram as interações anteriores, alguns agentes são fácil de decidir que nunca vão cooperar, sempre vão tentar passar a perna, então quando todo mundo começa a não cooperar com ele, ele sempre perde pontos.
Por outro lado algumas estratégias são mais confusas. Pegue o ditador agressivo por exemplo, como ele muda de estratégia de acordo com a sua, sempre tentando te convencer a voltar a ser cooperativo e então trapaceando denovo, ele acaba se mantendo bem, mas veja de temos 3 agentes dele também, mas somente 2 continuam bem no final, um começa a descambar.

No artigo do Axelrod e do Hamilton, apesar do tit4tat ter sido a estratégia vitoriosa, eles tenta ver o que fez algumas estrategias serem melhores que muitas outras. E o que observaram é que cooperar pode não ser bom, quando se está no meio de um monte de trapaceiros, mas quando você encontra outros agentes que também cooperam, esses dois começam a melhorar a pontuação juntos.
É claro que se você ficar cooperando com todo mundo, você vai se ferrar, como a gente o agente que sempre coopera, ele não vai bem, mas se você cooperar com quem coopera, e trapacear quem trapaceia, logos os trapaceiros começam a não conseguir mais somar pontos, e somente os agentes cooperando entre si conseguem pontos.
So que as vezes alguns agente deixam de cooperar, como o agente ditador agressivo, mas quando ele volta a cooperar, se você perdoa-lo, vocês podem voltar a a fazer pontos juntos.

O que o Axelrod e Hamilton viram foi isso, características como sempre cooperar com quem coopera com você, trapacear os trapaceiros e perdoar quem quer mudar de estrategia gerava a melhor pontuação. Logo ser legal com pessoas legais pode ser uma estratégia evolutivamente estável.

Essa conclusão é muito legal né. Claro que comportamento em espécies é muito mais complexo que essas regrinhas simples, mas por exemplo em primatas, que vivem em grupos e a gente ve um tirando os piolhos dos outros é um exemplo. Talvez seria melhor deixar alguém tirar seu parasitas mas não perder tempo fazendo isso pra ninguém. Mas logo ninguém mais vai querer interagir com você no grupo e você vai encher de piolhos, logo todo mundo tem que ser legal um com os outros, se for querer viver sem piolhos.

Outro exemplo é você pensar em dois insetos, onde cada um poe um ovo, se os dois fizerem um ninho, eles podem revesar, onde um cuida dos ovos e o outro come, e depois eles invertem os papeis, mas se na sua vez de revesar, você abandonar o ovo do seu parceiro, ele não deixara descendentes mas você, nessa geração você se deu bem, mas logo na próxima não haverá ninguém que cooperara. Assim apesar de em uma geração você ter uma grande vantagem, na próxima você se ferra. Como no jogo, no inicio todos os agentes estão embolados, mas no final muitos deles despencam. A estratégia que sempre trapaceia é funciona assim. Veja que a linha laranja, até a rodada 50 por ai, estava entre as melhores pontuações, somente depois que ele começa a cair, ai sem ninguém para trapacear, so perde pontos.

Isso tudo é muito legal, vermos como uma simulação pode trazer tantas conclusões assim. O W.D. Hamilton tem modelos muito legais, além de ser uma personalidade muito interessante. No obituário que o Dawkins escreveu quando o Hamilton morreu, ele comenta que o Hamilton, professor titular de Cambridge ia de bicicleta trabalhar todo dia e não ligava muito para muitas coisas do nosso mundo comum, como dinheiro ou luxo. Muito legal.

Existe muita coisa para escrever sobre esse tópico, mas uma hora o post tem que acabar, e melhor acabar por aqui. O script para a simulação está embaixo, é legal ver ela varias vezes, variando outras caracteristicas como números de agentes, ou mesmo escrever outras estrátegias, e ver como as vezes o tit4tat chega a perder, mas quase sempre ele está bem colocado, e aumentando os pontos. Tanto o número de agentes como a constituição deles influencia muito o resultado, tente colocar muitos ditadores ou ainda muitos agentes imprevisiveis, que da resultado bem legais para se pensar. Mas uma evidencia para a cooperação, para um porque as pessoas serem boas é muito legal para dormir tranquilo, pensando sempre que o mundo não deve ter somente pessoas ruins, ja vivemos com tantos problemas na vida, ler esse artigo traz algum conforto na vida, pelo menos pra mim. Os scripts sempre no repositório recologia e eu ja tinha pego um script parecido de onde copiei essa simulação a muito tempo atrás, mas nem lembro mais o link, mas foi outra pessoa que bolou esse esquema para a simulação, so não lembro o nome para dar o devido crédito, peço desculpa, mas tenho certeza que todo mundo fica feliz em ver seus script e funções rodando por ai. E é isso, desculpem os erros de portugues, ja que escrevi muito rapido esse texto.
Uma ultima coisa, eu ouço um podcast sobre ciência chamado fronteiras da ciência, o episodio 19 da quinta temporada é sobre a teoria dos jogos, são fisicos falando disso, que devem entender teoria dos jogos melhor que eu, então vale muito a pena ouvir, caso tenha mais interesse no assunto, e tem esse programa que eu ja vi aqui, que também reproduz o torneio e é bem legal, para quem quer só brincar e não mexer com script, e até o próximo post.

# Simulação do torneio de Axelroad
set.seed(1415)
 
# Retorno para:
# Cooperação mutua
coopera.coopera <- 1
# Você trapaceia mas seu oponente coopera
trapaceia.coopera <- 2
# Você coopera mais seu oponente trapaceia
coopera.trapaceia <- -2
# Todo mundo trapaceia
trapaceia.trapaceia <- -1
 
# Função que retorna o valor de resultado da interação para o participante
retorno <- function(voce,oponente) {
    saida <- 0*voce
    saida[(voce==1)&(oponente==1)] = coopera.coopera
    saida[(voce==0)&(oponente==1)] = trapaceia.coopera
    saida[(voce==1)&(oponente==0)] = coopera.trapaceia
    saida[(voce==0)&(oponente==0)] = trapaceia.trapaceia
    return(saida)
}
 
# Exemplo da função para tres jogos
retorno(c(0,0,0),c(1,1,0))
 
 
#Os tipos de agentes.
 
#Os agentes simples:
sempre.coopera   <- function(h) { return(TRUE) }
sempre.trapaceia <- function(h) { return(FALSE) }
 
 
# Imprevisivel
imprevisivel <- function(h) { return(1==rbinom(1,1,.5)) }
 
#tit4tat
 
tit4tat <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(rev(h)[1]==1)
    }
}
 
 
# tit4tat inicialmente agressivo
ag.tit4tat <- function(h) {
  if (length(h)==0) {
      return(FALSE)
  } else {
      return(rev(h)[1]==1)
  }
}
 
# tit4tat oportunista
op.tit4tat <- function(h) {
  if (length(h)==0) {
      return(TRUE)
  } else {
      return(((n<n.iteracoes-n.agentes)&(rev(h)[1]==1)))
  }
}
 
 
# fairbot
fairbot <- function(h) {
  if (length(h)==0) {
      return(TRUE)
  } else {
      return(mean(h,na.rm=T)>=.5)
  }
}
 
# fairbot inicialmente agressivo
ag.fairbot <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)>=.5)
    }
}
 
# Ditador covarde
ditador <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}
 
# Ditador Agressiveo
ag.ditador <- function(h) {
    if (length(h)==0) {
        return(FALSE)
    } else {
        return(mean(h,na.rm=T)<=.5)
    }
}
 
# Frenemy
frenemy <- function(h) {
    if (length(h)==0) {
        return(TRUE)
    } else {
        return(mean(rev(h)[1:2],na.rm=T)==1)
    }
}
 
ag.ditador(c())
 
ag.ditador(c(T,T,T,T,T,T))
ag.ditador(c(T,T,F,F,F))
ag.ditador(c(T,T,F,F,F,F,T,T,T))
 
 
#Parametros da simulação
n.iteracoes <- 500
n.agentes   <- 20
 
lista.agentes <- c("sempre.coopera", "sempre.trapaceia", "tit4tat" , "ag.tit4tat",
               "fairbot"    , "ag.fairbot"  , "ditador", "ag.ditador",
               "op.tit4tat" , "imprevisivel"  , "frenemy")
 
# Selecionando agentes para a simulação
agentes <- sort(sample(lista.agentes, n.agentes, replace=T))
 
# Matriz da historia de cada jogador
agentes.numero <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
# Matriz de ações dos agentes
agentes.acoes.recebidas <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
agentes.acoes.aplicadas <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
# Pontuação historicas
agentes.pontos <- matrix(NA, ncol=n.agentes, nrow=n.iteracoes)
 
for (n in 1:n.iteracoes) {
  while (sum(is.na(agentes.numero[n,]))>0) {
  sempar <- 1:n.agentes
    while (length(sempar)>1) {
      i <- sempar[1]
      sempar <- sempar[sempar!=i]
      i.adversario <- sample(rep(sempar,2), 1)
      sempar <- sempar[sempar!=i.adversario]
      agentes.numero[n,i] <- i.adversario
      agentes.numero[n,i.adversario] <- i
    }
  }
  for (i in 1:n.agentes) {
    a.encarado <- agentes.numero[n,i]
    a.numeros <- agentes.numero[,i]
    a.acao <- agentes.acoes.recebidas[,i]
    h <- a.acao[(a.numeros==a.encarado)&(n>1:n.iteracoes)]
    response <- get(agentes[i])(h)
    agentes.acoes.aplicadas[n,i] <- response
    agentes.acoes.recebidas[n,a.encarado] <- response
  }
}
 
# Calculando os scores por rodada
placar.interacao <- retorno(agentes.acoes.aplicadas, agentes.acoes.recebidas)
placar.acumulado <- apply(placar.interacao, 2, cumsum)
 
 
#Figura 1
cores <- rainbow(length(lista.agentes))
plot(x=c(1,n.iteracoes+25),y=c(min(placar.acumulado),max(placar.acumulado)),type="n",frame=F,
     main="Performance ao longo do tempo", xlab="Encontro", ylab="Pontuação")
for (i in 1:n.agentes) {
    lines(1:n.iteracoes, placar.acumulado[,i], col=cores[agentes[i]==lista.agentes], lwd=2)
}
legend("topleft",lty=1,legend=lista.agentes,col=cores,bty="n",lwd=2)

Forrageamento ótimo: Polinizadores!

Ao falar de polinizadores, agora eu lembro dessa imagem.

av0dbmd_700b

Isso porque polinização é talvez o melhor exemplo de serviço ambiental prestado.
Serviço ambiental é quando você ta recebendo um benefício do ecossistema em volta de você, e podemos até quantificar isso em unidade de dinheiro, reais ou dólar se preferir.

É isso ai, é só pensar de forma simples, sem abelhas, muitas culturas produzem menos sementes, ou até nenhuma semente. E ninguém quer pagar alguém pra ficar pegando pólen em uma plantinha e levando para outra. Nem temos muitas maneiras artificiais eficientes de fazer isso, pensando numa escala de produção de sementes em massa. Então polinização é de maior interesse para todos, mesmo que você não acha de primeira.

Bem todos os organismos forrageiam de forma a manter o seu sucesso maximizado, a gente vem vendo isso aqui e aqui, porque se você não faz isso, você tem menos energia para reprodução ou vai morrer no meio da sua vida e como vimos aqui, quem segue a mesma estratégia falha que você vai rodar. Agora podemos pensar em forrageamento ótimo não somente da perspectiva de comida, por exemplo flores forrageiam polinizadores para conseguir produzir sementes, e otimizam estratégias para isso, já que elas querem atrair polinizadores, mas não predadores.

Do outro lado, a gente pode pensar também numa colonia de abelhas, que de modo geral consiste de uma rainha e trabalhadoras em diferente castas com serviços diversos na colônia. Uma abelha forrageadora é uma trabalhadora com uma tarefa simples, coletar pólen e néctar. A eficiência dessa nossa amiguinha e a quantidade de comida que ela coleta é decisiva para o sucesso reprodutivo da colônia. O pólen e néctar é trazido de volta para a colônia, e guardado como reserva energética para a colônia. Alias, mel é um dos alimentos mais fáceis de estocar, não precisa de geladeira, não funga e até em piramides, encontraram mel intacto depois de mais de mil anos estocado, se a discovery não mentiu pra mim. Mas todos os processos metabólicos, incluindo crescimento, sobrevivência e taxa de reprodução são dependentes da quantidade de néctar e pólen estocados. Assim a seleção natural deve favorecer estratégias que maximizam da taxa líquida de ganho de energia enquanto forrageando atrás de néctar.

Vamos pensar num grupo de abelhas que tem pequenas colonias, pra facilitar o entendimento aqui.
As mamangaba, ou bumblebee em inglês. O nome do robozinho do transformes que era um fusquinha no desenho mas virou camaro nos filmes, já que a Alemanha não mandou dindin pro filme, mas whatever é uma abelha que estabelece colonias em para algumas de suas espécies. Vamos ver um modelinho tentando ver como funciona a interação dessa abelhas mamangaba forrageando no tremoceiro.

No tremoceiro, a flor la de baixo da inflorescência é mais velha, abre primeiro e contém normalmente mais néctar. Conforme o tempo passa, a gente corre na linha do tempo da fenologia da espécie, a flor fica velha demais e “morre”, vai virar fruto e a flor em cima dela passa a ser a última, até acabar as flores. Mas em um dia durante a fase da floração, há uma relação linear entre a posição da flor e o conteúdo calórico dessa flor, para o nosso tremoceiro aqui, a gente pode descrever essa relação da seguinte forma.

{C_i}={Calorias\ por\ flor} - {Inclinacao_i}

Aqui o índice i é referente a posição da flor na inflorescência, 1 é a flor mais de baixo, e se por exemplo a gente tem 20 flores, 20 é a la de cima, veja que a inclinação é negativa, isso acaba por impor um limite no número de flores possíveis, isso é, num teremos uma flor com quantidade negativa de néctar, uma flor não toma néctar da abelhinha, mas é só ser sensato ao usar essa relação, sempre lembrar dos limites dela.

Podemos representar essa relação assim.

01

A nossa mamangaba tem o comportamento de primeiro ir na flor de baixo, que contem mais néctar ou mair retorno em calorias, e vai indo na de cima, depois outra, até o topo da haste da inflorescência. Mas normalmente ela deixa a planta antes de visitar todas as flores daquela inflorescência.

Agora suponha que podemos coletar informação sobre a o consumo de energia de uma espécie de mamamgaba em um dia em particular, em um determinado local.

Matemagicamente, nos podemos modelar a eficiência energética baseado nas informações coletadas. Assim podemos estudar o movimento das mamangabas de uma planta para a próxima no ato de coletar néctar e dai analisar se ela está forrageando de uma maneira ótima ou não. O modelo matemático para calcular a taxa liquida de ganho energético para cada flor vai ser dada pela seguinte formula.

{Net_i} = {{P \sum \limits_i^n (\alpha+\beta \cdot i) - eb - P \cdot ew \cdot ( i - 1) - P \cdot i \cdot ef - ee \cdot (1 - P)} \over{tb + P \cdot tw \cdot (i - 1) + P \cdot i \cdot tf + te \cdot (1-P)}}

Onde

Net_i é o total de energia ganho visitando até a flor i

P é a probabilidade da inflorescência não ter sido drenada de néctar naquele momento, estar vazia.

eb é a energia gasta voando de uma inflorescência para outra

tb é o tempo gasto voando de uma inflorescência para outra

ew é a energia gasta voando de uma flor para outra na mesma inflorescência

tw é o tempo gasto voando de uma flor para outra na mesma inflorescência

ef é a energia gasta coletando todo o néctar da flor.

tf é o tempo gasto coletando todo o néctar da flor

ee é a energia gasta experimentando se a flor está vazia

te é o tempo gasto experimentando se a flor está vazia

n é a posição da ultima flor visitada pela abelha

O numerador é o total de energia ganha pela visita até a flor i e o denominador da formula descreve o tempo gasto forrageando na planta e movendo para a próxima planta.

Quando a razão é maximizada, nos podemos dizer que a abelha no seu forrageamento ótimo.

Vamos implementar essa formula como uma função no R

bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,alpha,beta) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+(alpha+beta*i)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}

Nessa função, entramos com os parâmetros das espécie vegetal e a saída é o ganho energético por visitar até a i^th flor, ou seja, da for 1 até a i.

Ai para ver o forrageamento ótimo, é só ver quando a abelha deveria “parar”.

02

Temos ali o ponto critico, se a nossa amiguinha mamangava abandonar a inflorescência antes, ela poderia ter saído daquela inflorescência com um melhor retorno energético, se ela sair depois daquele ponto, ela poderia ter tido um ganho maior, mais ao ficar, o gasto com deslocamento para a próxima flor, para inspecionar a próxima flor foi maior que o retorno, então começou a dar prejuízo energético continuar naquela flor. O ponto indicado é exatamente a hora de parar para ter o maior lucro com o menor investimento, energético.

Mas estamos olhando o forrageamento para um set fechado de parâmetros, e se esses parâmetros fossem um pouco diferentes?

Vamos pensar primeiro na probabilidade da inflorescência não ter sido drenada de néctar naquele momento, estar vazia.
Se essa probabilidade diminui, ou seja, a chance é maior de a flor estar vazia, isso vai diminuir nosso retorno energético como um todo não? Só imaginar o quanto isso vai resultar em encontrar plantas que não tem néctar, mas o gasto energético para ir olhar já foi, essa energia não volta mais.

03

Note ai que o forrageamento ótimo se move para a direita, ou seja, começa a valer mais a pena pegar mais flores, mesmo elas tendo cada vez menos néctar, porque é melhor garantir do que remediar, já que temos grande incerteza quando a próxima inflorescência.

Agora vamos se a gente manter a probabilidade da inflorescência não ter sido drenada de néctar naquele momento constante e alterar a tempo gasto voando de uma inflorescência para outra. Imagine que a vegetação é mais densa, existem mais obstáculos o que gera a necessidade de mais manobras no voo.

04

O ponto de forrageamento ótimo não se altera, mas o ganho energético diminui conforme as plantas estão mais espalhadas no ambiente. Então aglomeração de planta rendem mais energia, o que vai atrair mais polinizadores, ja que esses buscam o melhor ganho. E se a gente pensar, mesmo em muitas espécies, mesmo quando elas não estão exatamente num sistema como esse, tende a estar assim, aglomeradas numa mancha. É manchas podem ser bem mais rentáveis energeticamente do que espalhar uniformemente no campo as flores.

Agora vamos voltar o tempo gasto voando entre inflorescências para o normal e vamos varia o quanto uma flor a mais de néctar do que a próxima, mais especificamente, vamos homogeneizar um pouco mais as flores, diminuir aquela relação linear que falamos la em cima, alterando o coeficiente de inclinação beta.

Bem se a próxima flor tem o mesmo tanto de néctar que a anterior, vamos continuar um pouquinho mais nessa planta, já que da para encher os bolsos, ou perninhas nessa caso, com mais néctar.

05

Veja que o beta é um valor negativo, mas quanto mais reto, mais próximo de zero, onde zero seria todas as flores terem a mesma quantidade de néctar, mais vale a pena passar mais tempo numa inflorescência.

Ja o contrario.

06

É melhor sair antes, e quanto mais flores eu insisto em visitar, mais energia eu desperdiço trabalhando e menos eu levo para casa.

Agora a gente ainda pode fazer uma ultima perguntinha. E se a quantidade de néctar é imprevisível? Se é tudo na loucura e cada hora cada flor tem uma quantidade de néctar. Mesmo a planta fazendo o mesmo investimento energético que faria se fosse previsível a quantidade de néctar inflorescência.
A gente pode simular isso, pegando a média e desvio de conteúdo de néctar para a inflorescência com alta previsibilidade e fazer flores com conteúdo de néctar ao acaso, e jogar na nossa equaçãozinha e ver o que ocorre.

07

Cada linha cinza é uma simulação, uma determinada quantidade de néctar, e a linha vermelha é a média de 250 simulações, o caso médio para todas as simulações.

Olha ai que bolado, primeiro, vemos que o forrageamento ótimo quando é previsível a quantidade de néctar pela abelha, é bem menor para a abelha. Mas agora veja que no finalzinho, na décima flor, a curva de forrageamento ótimo estava subindo para flores com quantidade de néctar ao acaso, e passando as flores com quantidade de néctar previsível. La na frente nesse gráfico, o forrageamento em flores com quantidade de néctar ao acaso superaria até o retorno energético ótimo do forrageamento em flores previsíveis, pense que a próxima flor, se tiver próxima, sempre pode estar cheia até a boca do delicioso néctar, existe essa possibilidade, mas e ai, você teve que passar muito mais tempo forrageando, o que pode ser perigoso, você, uma abelhinha mamangaba pode ser comida no caminho, pode sei la ter um infarto de tanto trabalhar. Prevendo a quantidade de néctar de uma flor, a abelha pode trabalhar bem menos e já conseguir um bom retorno, sendo esse retorno seguro. Agora eu pergunto, como pode a abelha prever a quantidade de néctar?
Ai entra a comunicação entre espécies. Não é a toa que abelhas muitas vezes evoluem para serem especialista, se pensarmos em características das flores como sinais, coloração das pétalas, posição dessas, além de outras características podem funcionar como sinais de que ali tem néctar ou não, pensando por esse ponto de vista claro, existem outros. Além do que se a abelha tem “garantias”, ela pode preferir sempre a mesma espécie de planta, o que garante que o pólen sera transferido de forma mais eficiente, o que é interessante para a planta.

Veja que esse modelinho de polinização é bem legal, e podemos fazer ainda muitas outras modificações pequenas para se adequar a outras situações, por exemplo, trocar inflorescências por flores únicas, adicionar uma distância média entre flores, sorteando valores ao invés de fixo, se as plantas estão ao acaso no campo, ou modificar um pouco a ideia de inflorescência para manchas, e assim vai.

Bem ficamos por aqui, o script está no repositório recologia, polinização é bem legal, mas temos mais um post ainda sobre forrageamento ótimo, os predadores senta-espera maximizando seu forrageamento.

Referência:
Bernstein, R. 2003 – Population Ecology – An Introduction to Computer Simulations. John Wiley & Sons, Ltd 159pp.

Heinrich, B. 1983 – Do bumblebees forage optimally, and does it matter? American Zoologist 23: 273-281.

bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,alpha,beta) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+(alpha+beta*i)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}
 
#Calorias
i<-seq(0,10,1)
calories<-20-1.7*i
 
#figura 1
plot(i,calories,type="l",frame=F,ylab="Número de calorias da flor",xlab="Sequência de flores em um ramo da planta",
     main="Relação linear entre posição da flor e conteúdo calorico da flor",xlim=c(0,10),ylim=c(0,20))
 
#figura 2
saida<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, tb=5, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7, P=0.4",font.main=1,cex.main=0.9)
abline(h=max(saida),v=which(saida==max(saida)),lty=3,cex=2)
text(which(saida==max(saida))+0.5,max(saida)+0.025,"Forrageamento ótimo")
 
#Aumentando a chance da flor ter sido exaurida antes da visita
saida<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, tb=5, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.2,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.1,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("P=0.6","P=0.4","P=0.2","P=0.1"),lty=2,col=1:4,lwd=2,bty="n")
 
#Aumentando distancia entre flores
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, beta=-1.7",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=7,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=9,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=11,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("tb=5","tb=7","tb=9","tb=11"),lty=2,col=1:4,lwd=2,bty="n")
 
#Flores maturam num intervalo de tempo menor
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, tb=5",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.5)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.3)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.0)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("beta=-1.7","beta=-1.5","beta=-1.3","beta=-1.0"),lty=2,col=1:4,lwd=2,bty="n")
 
#Flores maturam num intervalo de tempo menor
jpeg("06.jpg")
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
plot(1:10,saida,type="l",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, alpha=20, tb=5",
      font.main=1,cex.main=0.9)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.9)
points(1:10,saida,type="l",col=2,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-2.1)
points(1:10,saida,type="l",col=3,lty=2,lwd=2)
saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-2.3)
points(1:10,saida,type="l",col=4,lty=2,lwd=2)
legend("topright",legend=c("beta=-1.7","beta=-1.9","beta=-2.1","beta=-2.3"),lty=2,col=1:4,lwd=2,bty="n")
 
#Experimento
saidaoriginal<-bumb(P=0.6,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,alpha=20,beta=-1.7)
 
 
bumb<-function(P,eb,ew,ef,ee,tb,tw,tf,te,n,media,desvio) {
    netE<-vector()
    sumi=0
    for(i in 1:n) {
        sumi=sumi+rnorm(1,media,desvio)
        net=(P*sumi-eb-P*ew*(i-1)-P*i*ef-ee*(1-P))/(tb+P*tw*(i-1)+P*i*tf+te*(1-P))
        netE[i]<-net
    }
    return(netE)
}
 
#figura 7
plot(1:10,saida,type="n",frame=F,xlim=c(1,10),ylab="Taxa de ingestão de energia líquida em cal/seg",
     xlab="Última flor visitou antes a mosca deixa a planta",ylim=c(0,1),col=1,lty=2,lwd=2,xaxt="n")
axis(1,at=c(1:10))
title("eb=0.1, ew=0.03, ef=0.02, ee=0.01, P=0.4, tw=0.05, tf=15, te=9, n=10, tb=5",font.main=1,cex.main=0.9)
 
matriz<-matrix(NA,ncol=10,nrow=250)
for(i in 1:250) {
    saida<-bumb(P=0.4,eb=0.1,ew=0.03,ef=0.02,ee=0.01,tb=5,tw=0.05,tf=15,te=9,n=10,media=mean(calories),desvio=sd(calories))
    matriz[i,]<-saida
    points(1:10,saida,type="l",col="gray",lty=3,lwd=0.7)
}
 
points(1:10,colSums(matriz)/250,type="l",col="red",lty=3,lwd=3)
points(1:10,saidaoriginal,type="l",col="black",lty=3,lwd=3)
legend("bottomright",legend=c("Conteúdo de néctar conhecido","Conteúdo de néctar ao acaso"),lwd=3,lty=3,
       col=c("black","red"))

Crescimento Populacional e a Seleção Natural

A teoria da evolução é talvez a teoria mais importante para a biologia. Ela unificou muitos ramos científicos atribuídos a biologia que antes dela eram praticamente áreas distintas, como fisiologia, taxonomia, comportamento, etc.

Temos a frase celebre do nosso colega Theodosius Dobzhansky.

Nada em biologia faz sentido se não for à luz da evolução

Mas será que a gente realmente entende como a evolução funciona?

Quero dizer, provavelmente, todo mundo tem uma noção básica. O que a gente aprende na escola é tranquilo, mas o problema é que muitas vezes nos lembramos do exemplo básico e fazemos generalizações muito amplas, ou pior ainda, generalizações curtas e não conseguimos ver o processo como ele realmente ocorre.

Na verdade a teoria da evolução é muito mais do que vamos ver aqui, mas vamos olhar como funciona o negócio de variação no número de descendentes produzidos. Quem não lembra do exemplo da mariposa e da revolução industrial. A Biston betularia. Essa ai é aquela mariposa que existia de dois tipo, ou dois fenótipos, como preferir.

Tinhamos a versão, ou fenótipo mais claro.
claro

E a versão, ou fenótipo, mais escuro.
escuro

Ai antes da poluição, a mais clarinha era comum, e a escura era muito rara, com o acumulo de fuligem nos troncos das arvores, a clarinha começou a não ser muito boa para se camuflar, a escura, da cor da fuligem, começou a se dar bem e o cenário se inverteu, com a clarinhas passando a ser a mais rara e a escura sendo a mais comum.

Daqui pra frente vamos chamar de fenótipos sempre. Se a gente fosse la no campo e conta-se todas as mariposas existentes, poderíamos dizer que existe uma quantidade N de mariposas, essa quantidade sendo independente do fenótipo, que é a cor para as mariposas.

E elas vão vivendo, se reproduzindo a uma certa taxa, e bem isso nos faz lembrar da equação de crescimento populacional.

N_t = N_0 \cdot r^t

Lembre-se que estamos falando do crescimento discreto ai em cima, não confunda com a formula N_t = N_0 \cdot e^{r \cdot t} , a gente já discutiu sobre crescimento discreto versus crescimento contínuo aqui, de uma olhada para relembrar e ver a diferença.

Mas continuando, então temos nosso crescimento populacional, vamos ficar com ele por enquanto para facilitar o raciocinio, se começar a meter aquele número e de Euler vai ficar tudo mais difícil.

N_t = N_0 \cdot r^t

So que essa equação é valida para toda a população dessa determinada espécie, todos os indivíduos de uma região que compõem uma população de uma certa espécie, nossa mariposinha por exemplo.

Vamos pensar nos fenótipos claro e escuro como fenótipo A e fenótipo a, cada um deles então vai ter uma frequência na população, que chamaremos de f_A e f_a e cada um deles vai ter sua taxa de crescimento R, r_A e r_a, que também é conhecido como o fitness absoluto (Absolute fitness), que a princípio vamos considerar como igual a taxa de crescimento da população como um todo, escrevendo isso de forma mais clara, vamos assumir que r_A = r_a = R.

Primeiro, vamos pensar que o número inicial de indivíduos A seria igual a sua frequência no total inicial da população.

N_{A,0} = N_0 \cdot f_A

Certo, se a população inicial é N_0, sabemos que uma fração dela é de mariposas claras, a fração é f_A, por isso se multiplicar o total pela fração N_0 \cdot f_A da N_{A,0}

O mesmo vai ser valido para a população de a,

N_{a,0} = N_0 \cdot f_a

Lembrando que f_A + f_a = 1, tem que ser o total da população.

Vamos em frente, A e a são a população total dessa espécie.

N_0 = N_0 \cdot f_A + N_0 \cdot f_a

Sem muita mágica até aqui.
Certo, agora como o crescimento vai funcionar aqui?
Muito simples, para a próxima geração (vamos pensar apenas na diferença de uma geração para a próxima) temos a quantidade inicial vezes a taxa de crescimento, so que por cada um dos fenótipos.

Para A temos

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a

Mas como r_A = r_a = r, o fitness dos dois são iguais, então poderíamos simplesmente dizer.

N_{A,1} = N_0 \cdot f_A \cdot r

e

N_{a,1} = N_0 \cdot f_a \cdot r

Então podemos ver que a população total no tempo 1 N_1 vai ser

N_1 = N_0 \cdot f_A \cdot r + N_0 \cdot f_a \cdot r

Veja que o termo N_0, que é a população total inicial e o termo rocorrem nas duas multiplicações, então podemos simplificar isso.

N_1 = N_0 \cdot r \cdot (f_A + f_a)

Lembre-se que f_A + f_a = 1 já que são frequência, e esses dois A e a tem que ser o todo.

Então apos uma geração a frequência de A será a quantidade de indivíduos A pelo total de indivíduos da população, a frequência f_A será o seguinte.

{f_{A,1}} = {{N_0 \cdot f_{A,0} \cdot r}\over{N_0 \cdot r}}

Como é tudo uma multiplicação podemos cortar os N_0 e r e assim ficamos somente com.

f_{A,1}= f_{A,0}

A frequência na próxima geração não mudou nada.

Mas é isso que esperamos não? Porque nesse exemplo r_A = r_a = r, ou seja os dois fenótipos eram igualmente bons, mas e se r_A \ne r_a, um deles pode ser melhor que o outro, sei la r_A > r_a ou r_A < r_a, mas ambos os casos r_A \ne r_a, como fica esse mesmo esquema?

A principio, a próxima geração vai ter a mesma formula do crescimento.

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a,

Ou ainda, generalizando para o tempo t, temos

N_{A,t} = N_0 \cdot f_A \cdot {r_A}^t

E para a temos

N_{a,t} = N_0 \cdot f_a \cdot {r_a}^t

Ou seja a população N_t será

N_{t} = N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t

Agora vamos fazer aquela conta denovo, após t gerações, a frequência f_A vai ser a quantidade de individuos A sobre o total N

{f_{A,t}} = {N_0 \cdot f_A \cdot {r_A}^t \over{N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t}}

Agora nos podemos avaliar a frequência de a depois de t gerações assim como fizemos antes, mas podemos simplificar um pouco essa equação. Podemos começar tirando esse N_0 em cima e em baixo na divisão.

A frequência f_A vai ser o total dela sobre o tamanho total da população.
{f_{A,t}} = {f_A \cdot {r_A}^t \over{f_A \cdot {r_A}^t + f_a \cdot {r_a}^t}}

Agora podemos dividir em cima e em baixo por r_A e ficamos com o seguinte

{f_{A,t}} = {{f_A}\over{{f_A} + {f_a \cdot {{r_a}^t\over{{r_A}^t}}}}}

E finalmente temos a nossa equação para acompanhar como a frequência gênica muda ao longo das gerações.
Se a gente pensar um pouco, quando a razão {r_a}^t\over{{r_A}^t} é 1, ambos os fenótipos são igualmente bons, essa razão da um, e teremos  {f_A}\over{f_A + f_a \cdot 1}, como  f_A + f_a = 1, teremos  f_A dividido por 1 e nada nunca muda. Agora qualquer desvio na razão {r_a}^t\over{{r_A}^t} para mais ou para menos que 1, teremos mudanças.

Mas de qualquer forma, vamos jogar alguns números e fazer uma simulação no R.

Precisamos estabelecer algumas condições iniciais, e pensar por quantas gerações queremos simular.

# ************************************************** # Condições iniciais # ************************************************** N0 <- 1000 # População inicial rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a" fA <- 0.3 # Frequência do alelo "A" max_gen <- 20 # Número de gerações a simular

Veja que baseado naquelas quantidades, o resto dos ingredientes são derivados

# ************************************************** # Calculando variáveis derivadas # ************************************************** fa <- 1.0 - fA # Frequência do alelo "a" NA_0 <- N0 * fA # População inicial do alelo "A" Na_0 <- N0 * fa # População inicial do alelo "a"

Mas vamos la, criamos uma função para com aquelas condições iniciais e usando as equações que estabelecemos acima para ver o que acontece.

evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}

E testamos nossa funçãozinha.

> evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10) Nt NA_t Na_t fA_t fa_t [1,] 1000.000 300.0000 700.000 0.3 0.7 [2,] 1440.000 432.0000 1008.000 0.3 0.7 [3,] 1728.000 518.4000 1209.600 0.3 0.7 [4,] 2073.600 622.0800 1451.520 0.3 0.7 [5,] 2488.320 746.4960 1741.824 0.3 0.7 [6,] 2985.984 895.7952 2090.189 0.3 0.7 [7,] 3583.181 1074.9542 2508.227 0.3 0.7 [8,] 4299.817 1289.9451 3009.872 0.3 0.7 [9,] 5159.780 1547.9341 3611.846 0.3 0.7 [10,] 6191.736 1857.5209 4334.215 0.3 0.7 [11,] 7430.084 2229.0251 5201.059 0.3 0.7

Temos como entradas as quantidades iniciais, que repetimos na primeira linha, a geração zero, e temos na saída o tamanho da população total como primeira coluna, como ela vai mudando sendo cada geração uma linha, como cada fenótipo muda, bem como suas frequências ao longo das gerações, mas como olhar para uma matriz de números é meio chato, vamos representar esses dados em dois gráficos, um de como está a população e outro de como estão as frequências.

01

Certo, bem simples não, a população está crescendo exponencialmente, e as frequências genicas estão constantes, como previmos, porque as taxas de crescimento iniciais que colocamos para essa simulação são iguais.

rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a"

Agora lembre-se que nem tudo na vida são flores, se tivermos taxas de crescimento abaixo de 1, a população estará diminuindo.

Vamos pegar esse caso aqui.

#N0 = 1000 rate_A = 0.7 rate_a = 0.7 fA = 0.3

02

Mas ainda temos as frequências dos fenótipos iguais, aqui ta todo mundo se ferrando igual. Agora quando a gente diz que está havendo evolução é quando as frequência começam a mudar, quando aparece um novo mutante que derrepente se da melhor, se alguém tem a taxa de crescimento maior que sua contraparte, vemos o seguinte.

#N0 = 1000 rate_A = 2.0 rate_a = 1.5 fA = 0.02

03

Veja que aqui, nenhum dos fenótipos está necessáriamente se ferrando, só temos um fenótipo que é melhor que o outro, mas ninguém está diminuindo em quantidade ao longo do tempo, a população está crescendo como um todo, todos os fenótipos estão crescendo, apenas um mais rápido que o outro, mas isso faz com que a frequência dos fenótipos mude. Quando eu digo que as vezes a evolução pode não ser tão simples, porque como no exemplo das mariposas, a gente não pensa nesse caso assim, a gente sempre pensa que se alguém é pior, ele tem que estar se ferrando, como no próximo exemplo

#N0 = 1000 rate_A = 1.2 rate_a = 0.9 fA = 0.02

04

Essa é a situação que a gente tem na cabeça para as mariposas, pelo menos a primeira que vinha na minha, se a gente pensar na linha azul como as mariposas brancas e vermelho as mariposas escuras, a gente quer que alguém seja muito comum e alguém raro, ai esse cara mais comum, as mariposas brancas, começam a se dar mal por causa da revolução industrial, as mariposas escuras, que antes eram raras e ruim, começam a se dar bem, e temos uma inversão de quem é mais raro e quem é mais comum, mas vemos que a população decaiu um pouco e depois voltou a crescer. Mas como vimos, a evolução pode acontecer mesmo quando todo mundo está se dando bem, com a população crescendo, e além disso ainda, podemos pensar em outra situação.

#N0 = 10000 rate_A = 0.8 rate_a = 0.6 fA = 0.02

05

E é isso ai, ambos os fenótipos estão se ferrando, ambos estão diminuindo sua frequência. Cada vez que você olha a mata, você ve menos e menos daquela espécia que você observa, mas isso não quer dizer que a evolução não possa estar agindo ali. Veja que aqui não é quem se ferra e quem se da bem, aqui a situação é quem se ferra menos, quem pode perdurar aos tempos difíceis, talvez essa espécie seja extinta, mas talvez não, e quando os tempos mudarem, ela pode voltar a crescer, mas as frequências dos fenótipos estão alteradas já na próxima largada para o crescimento populacional.

Mas, não é a toa que nos humanos somos péssimos geradores de números aleatórios, pelo menos eu, sempre tendia a pensar em evolução, no caso de uma espécie que tem um fenótipo se tornando mais raro, porque a população dele está diminuindo, e outro se tornando mais comum porque a população dele está aumentando, mas esteja a população estacionada, aumentando ou diminuindo, as frequências de alelos, fenótipos, snp, como queria pensar, podem estar mudando, e todos esses cenários tem que estar na nossas cabeças.

So que a primeira coisa que vem na minha cabeça é um fenótipo se dando bem e outro se dando mal. Mas existem muitos mais possibilidades. Então é preciso se policiar para não tomar conclusões precipitadas, e conduzir uma análise de dados que de possibilidade de todos os casos serem avaliados corretamente.

Alias, para teste em sistemas mais complexos, sempre esses teste envolvem simulações com parâmetros aleatórios, porque muitas vezes o que podemos pensar não contempla todas as possibilidades, pensa em quantas árvores podemos desenhar com n espécies, eu nunca pensei que existissem tantas possibilidades assim, e é fácil se perder nas contas. Mas ao tentar algumas simulações, passar tudo para um modelo e ir tentando várias possibilidades de parâmetros, podemos talvez ver aquilo que nunca tinha passado pela nossa cabeça.

Bem é isso ai, as equações la do começo do post eu peguei de um curso do coursera, o Computational Molecular Evolution, ótimo curso alias, uma excelente seleção de tópicos para a ementa e conceitos bem complicados são explicados de maneira bem simples, pena que a maior parte pratica é usando o paup, mas agora é ir adaptando tudo pro R com as ferramentas que tivermos a mãos. Scripts no repositório recologia e aqui em baixo, e até o próximo post.

# **************************************************
# Condições iniciais
# **************************************************
 
N0      <- 1000 # População inicial
rate_A  <- 1.2  # Taxa de crescimento do alelo "A"
rate_a  <- 1.2  # Taxa de crescimento do alelo "a"
fA      <- 0.3  # Frequência do alelo "A"
max_gen <- 20   # Número de gerações a simular
 
# **************************************************
# Calculando variáveis derivadas
# **************************************************
 
fa   <- 1.0 - fA  # Frequência do alelo "a"
NA_0 <- N0 * fA   # População inicial do alelo "A"
Na_0 <- N0 * fa   # Polulação inicial do alelo "a"
 
# **************************************************
# Simulação
# **************************************************
 
evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}
 
evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10)
 
saida<-evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 50     rate_A = 1.2  rate_a = 1.2  fA = 0.3
saida<-evosimu(N0=50,rate_A=1.2,rate_a=1.2,fA=0.3,max_gen=20)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 0.7  rate_a = 0.7  fA = 0.3
saida<-evosimu(N0=100,rate_A=0.7,rate_a=0.7,fA=0.3,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topright",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 2.0  rate_a = 1.5  fA = 0.02
saida<-evosimu(N0=1000,rate_A=2.0,rate_a=1.5,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 1.2  rate_a = 0.9  fA = 0.02
saida<-evosimu(N0=1000,rate_A=1.2,rate_a=0.9,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 10000  rate_A = 0.8  rate_a = 0.6  fA = 0.02
saida<-evosimu(N0=10000,rate_A=0.8,rate_a=0.6,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")