Evolução – Fisherian Optimality Models – Parte 4

Continuando com os modelinho de evolução, seguindo a parte 1, parte 2 e parte 3, na qual usando um somatório para ver o sucesso de uma espécie ao longo do tempo, dado o seu tamanho, nos estamos fazendo uma afirmação sobre quando o censo é feito, no tempo 1, 2, 3, etc. Uma suposição diferente é a que o processo é contínuo (de uma olhadinha aqui, na diferença entre discreto e contínuo para o crescimento populacional) e o uso de uma integral é apropriado. Reiterando as suposições sobre o nosso modelo.

Pressupostos gerais:

1. O organismo é iteroparo
2. Fecundidade, F, aumenta com o tamanho do corpo x, que não muda depois da maturidade atingida.
3. Sobrevivência, S, diminui com o tamanho do corpo x.
4. O fitness, W, é uma função da fecundidade e da sobrevivência.

Pressupostos matemáticos:

1. Maturidade ocorre na idade 1 depois da qual não há mais crescimento.

2. A Fecundidade aumenta linearmente com o tamanho na maturidade, resultando na fecundidade sendo uma função uniforme da idade, t:

F_t=a_F+b_{F}x

3. A taxa instantânea de mortalidade aumenta linearmente com o tamanho do corpo obtido na idade 1 e é constante por unidade de tempo. A partir desse pressuposto, a sobrevivência na idade t é dada por

S_t = e^{-(a_S+b_Sx)t}

Note primeiro que a sobrevivência não é mais independente do tamanho do corpo e que para fazer a sobrevivência uma função que diminui em relação ao tamanho do corpo, dado que a função é exponencial, nos alteramos o a_S-b_Sx com a_S+b_Sx.

4. Fitness, W, é o sucesso reprodutivo esperado ao longo da vida, R_0, dado como o produto acumulado da sobrevivência e da fecundidade

W = R_0 = \int\limits_{1}^\infty F_{t} S_t dt =\int\limits_{1}^\infty (a_F + b_Fx)e^{-(a_S+b_Sx) t} dt

Agora a gente não pode retirar do somatório os efeitos dependentes da idade e assim o tamanho ótimo para o corpo não será como o visto na parte 2.

Plotando a função de fitness

Há duas formas de fazer uma figura de W com x:

Primeiro, nos podemos integrar a equação acima para dar W como uma função de x.
Segundo, nos podemos usar uma rotina de integração numérica.

A primeira seria a melhor escolha, com um resultado exato, mas integração é frequentemente complicada demais ou simplesmente impossível.

Sabendo que \int e^{-Bt} dt = \frac{-e^{-Bt}}{B}, a equação do W pode ser integrada para.

W = (a_F + b_Fx)[\frac{-e^{-(a_S+b_S x)t}}{a_S + b_S x}]_1^\infty = 0 + \frac{a_F + b_F x}{a_S + b_S x} e^{a_S+b_S x}

Nos vamos usar a integração numérica da função integrate, passando pra ela a função que tem que ser integrada, que contém a expressão la de cima (antes de ser integrada). Note também que a integral é sobre a idade, mas o que nos queremos é variar o x, então a gente vai precisar de um loop para variar nele os valores de x. Os valores dos parâmetros que vamos estabelecer serão Af=0, Bf=8, As=1 e Bf=0.5

Assim vemos o seguinte.

01

Bem é fácil ver que existe um valor de tamanho ótimo, o ponto mais alto da curva, e é pra la que a evolução deve levar os indivíduos dessa espécie, mas vamos calcular exatamente qual deve ser esse valor.

Achando o máximo usando calculo.

Achando a derivativa de W em relação a x temos

\frac{dW}{dx} = \frac{b_F}{a_S+b_S x}e^{-(a_S+b_S x)} - \frac{(a_F+b_F x)b_S}{a_S+b_S x}e^{a_S+b_S x} -  \frac{(a_F+b_Fx)b_S}{(a_S+b_S x)^2} e^{-a_S+b_Sx}

e juntando o que todo mundo tem em comum temos

\frac{dW}{dx} = \frac{e^{-(a_S+b_S x)}}{a_S+b_S x} \{  b_F-(a_F+b_F x)b_S - \frac{(a_F + b_F x)b_S}{a_S + b_S x}    \}

Isso pode ser feito com a função deriv do R, mais fácil até.

> y <- deriv(~(0+4*x)*exp(-(1+0.5*x))/(1+0.5*x),"x") > y expression({ .expr2 <- 0 + 4 * x .expr4 <- 1 + 0.5 * x .expr6 <- exp(-.expr4) .expr7 <- .expr2 * .expr6 .value <- .expr7/.expr4 .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- (4 * .expr6 - .expr2 * (.expr6 * 0.5))/.expr4 - .expr7 * 0.5/.expr4^2 attr(.value, "gradient") <- .grad .value })

Agora nos não queremos entender tintin por tintin toda essa equação gigantesca, o que a gente quer é topo do fitness que vimos na figura, e para tanto a gente quer o local onde  \frac{dW}{dx}=0, nos precisamos do valor de x que faz a parte dentro das chaves da derivativa ser zero.

Duas possibilidades surgem aqui, usar diretamente a equação que derivamos, fornecendo ela para a função uniroot e a outra é pegar a saída do deriv e fornecer isso para o uniroot.

No primeiro caso:

1
2
3
funcao <- function(x) {
    return(4+0.5*(0-4*x)-(0+4*x)*0.5/(1+0.5*x))
}
> saida <- uniroot(funcao, interval= c(0,4)) > saida$root [1] 1.236068

Ou podemos usar diretamente a saída do deriv

1
2
3
4
5
6
7
8
9
10
11
funcao2 <- function(w) {
    y <- deriv(~(0+4*x)*exp(-(1+0.5*x))/(1+0.5*x),"x") # calcule a derivada
    #coloque x=w
    x <- w
    #calcule a derivada em w
    z <- eval(y)
    # Atribua o valor para d
    d <- attr(z,"gradient")
    # retorne esse valor
    return(d)
}
> B <- uniroot(funcao2, interval= c(0,4)) > B$root [1] 1.236068

E temos o mesmo resultado, e se conferir com a figura la atrás, ta certinho o local do pico de W em relação a x.

Nesse cenário, alias, nos podemos achar a solução exata. Primeiro nos colocamos todos os termos da derivada de dentro das chaves sobre  a_S+b_S x para termos:

b_F-(a_F+b_F x)b_S - \frac{(a_F + b_F x)b_S}{a_S + b_S x}    = \frac{-a_Fb_Sa_S -x(b_Fb_Sa_S + a_Fb_S^2) - x^2b_Fb_S^2}{a_S+b_Sx}

O numerador é uma equação quadrática (ax^2+bx+c=0), lembre-se que os valores do parâmetros nos temos, logo as raízes podem ser extraídas com a velha formula de bhaskara (x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}), que da 1.236068, o mesmo que ja encontramos. Note que o tamanho do corpo ótimo, x, não é o mesmo encontrado quando usamos o somatório no lugar da integral. Usar uma integral ou um somatório vai depender dos pressupostos biológicos. Nesse caso, a modelo com integral pode ser mais facilmente resolvido e assim não existem motivos fortes para preferir um modelo em detrimento do outro. Mas precisamos pensar muito sobre os pressupostos biológicos, e não sacrificar estes por um modelo mais simples de solucionar numericamente.

Achando o máximo usando uma aproximação numérica.

Nos podemos usar duas funções do R, integrate e nlm. O valor da função para uma dada idade e x, (a_F + b_Fx)e^{-(a_S+b_S x)idade} é determinado de uma função fornecida por nós. Lembre-se que nós temos que retornar um valor negativo de fitness, porque estamos achando o mínimo! Para obter a integral de uma função que fizemos, nos chamamos a função integrate, passando nossa função, o valor ótimo é encontrado passando a saída disso tudo para o nlm.

1
2
3
4
5
6
fitness <- function(idade,x,Af=0,Bf=4,As=1,Bs=0.5)  {
    #Fitness negativo
    return (-(Af+Bf*x)*exp(-(As+Bs*x)*idade))
}
# Função para chamar a integral
funcao <- function(x){integrate(fitness,1,Inf,x)$value}

E minimizando

> nlm(funcao,p=1)$estimate [1] 1.236067

Bem é isso ai, o script vai estar la no repositório recologia, o tempo grande foi duro, mas espero que agora passe e começamos a entrar nos modelos evolutivos mais legais, se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais, e quase tudo aqui vem do texto do livro do Derek A. Roff, se você achou legal esse post, o livro é um “must read”, com código de matlab também.

Referência:
Derek A. Roff 2010 – Modeling Evolution – an introduction to numerical methods. Oxford University Press 352pp.

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
# Função para integração númerica
Fitness <- function(idade,x,Af=0,Bf=8,As=1,Bs=0.5) {
    return ((Af-Bf*x)*exp(-(As+Bs*x)*idade))
}
 
# Figura do fitness
n<- 100
z<- seq(0,3,length=n)
W<- matrix(0,n,1)
 
# Iterando sobre n tamanhos do corpo
for (i in 1:n) {
    #Valor de x
    x<- z[i]
    #Integrando de 1 até o infinito e guardando na saida W
    W[i] <- integrate(Fitness,1,Inf,x)$value
}
 
#Plotando a figura
#jpeg("01.jpg")
plot(z,-W,type="l", xlab="Tamanho do corpo, x", ylab="Fitness, W",las=1,lwd=4)
#dev.off()
 
##################################
## Derivada
##################################
y <- deriv(~(0+4*x)*exp(-(1+0.5*x))/(1+0.5*x),"x")
y
 
##################################
## Achando o maximo
##################################
funcao <- function(x) {
    return(4+0.5*(0-4*x)-(0+4*x)*0.5/(1+0.5*x))
}
saida <- uniroot(funcao, interval= c(0,4))
saida$root
 
#saida do deriv
funcao2 <- function(w) {
    y <- deriv(~(0+4*x)*exp(-(1+0.5*x))/(1+0.5*x),"x") # calcule a derivada
    #coloque x=w
    x <- w
    #calcule a derivada em w
    z <- eval(y)
    # Atribua o valor para d
    d <- attr(z,"gradient")
    # retorne esse valor
    return(d)
}
 
# Encontre a raiz
B <- uniroot(funcao2, interval= c(0,4))
B$root
 
 
##################################
## Usando métodos númericos
##################################
fitness <- function(idade,x,Af=0,Bf=4,As=1,Bs=0.5)  {
    #Fitness negativo
    return (-(Af+Bf*x)*exp(-(As+Bs*x)*idade))
}
# Função para chamar a integral
funcao <- function(x){integrate(fitness,1,Inf,x)$value}
# Minimizando a função
nlm(funcao,p=1)$estimate