Evolução – Fisherian Optimality Models – Parte 5

Descendo a quantidade de um post por mês quase, mas ainda tenho esperança de reverter esse quadro.

Agora, voltando a evolução básica de Fisher, lembrando que é básica para ele, pro Fisher, não que eu ache isso básico, de forma nenhuma, gostaria, mas não dou conta de achar isso básico ainda. Mas ok, até agora nos assumimos que a medida apropriada para o fitness é o sucesso reprodutivo ao longo da vida, R0. Apesar que essa medida pode ser apropriada para uma população estável, uma medida de fitness mais geral é o parâmetro maltusiano, r, que é igual a taxa de crescimento da população na distribuição de idades estável.

 \int_0^\infty e^{-rt} l(t)m(t)dt =1

 \sum\limits_{t=1}^\infty  e^{-rt}l_t m_t = 1

onde l(t), l_t são as probabilidades de sobrevivência na idade t e m(t), m_t são os nascimentos específicos da idade de fêmeas (= fecundidade/2, assumindo uma razão sexual de 1:1). A diferença da notação usada nas duas equações são geralmente de pouco ou nenhuma consequência ( a equação da diferença pode ser igualmente escrita como a equação da integral) e usada simplesmente para ilustrar como a diferença na notação não deve implicar em diferença na interpretação. Note que a equação da diferença (do somatório) começa no fim do primeiro período de tempo, desde que nos assumimos que a fecundidade é zero no nascimento (claro que poderíamos iniciar do zero se nos simplesmente colocarmos m_0=0). Apesar que essas equações não englobam diretamente a contribuição dos machos, nos poderíamos escrever de forma similar equações para os machos relacionando seu sucesso em encontrar fêmeas a taxa de crescimento. Esse pressuposto referente ao r é que qualquer mutação que aumenta o r devera aumentar em frequência na população. A algum tempo atras, escrevemos sobre isso aqui, esse post foi muito legal, para ver como o r é determinante ao longo do tempo, para ver como um fenótipo ou característica muda de proporção nas populações, mesmo quando por exemplo ta todo mundo decaindo. Outra coisa para ter claro na cabeça olhando esse modelo, que sempre estamos levando em conta o tamanho, é lembrar da seleção sexual. Da para dar uma olhada no artigo do Ronald Aylmer Fisher aqui para ter uma ideia de como é um cara que escreveu seu nome na historia rediz um artigo, alias esses artigos meio que escritos a maquina de datilografar, sem figuras, porque acho que o povo desenhava a mão, mas uma leitura uma vez na vida de uns paper desse não mata ninguém.

Mas ok, vamos voltar ao que interessa, nos vamos considerar os mesmos pressupostos, com exceção que agora a medida de fitness é o r.

Pressupostos gerais.

1. O organismo é iteroparo.
2. Fecundidade, F, aumenta com o tamanho do corpo, x, o qual não muda depois da maturidade
3. Sobrevivência, S, diminui com o tamanho do corpo, x.
4. 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. Fecundidade aumenta linearmente com o tamanho na maturidade, resultando na fecundidade sendo uma função uniforme da idade

 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. Sobre esse pressuposto, a sobrevivência na idade t é dada por

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

Note que para fazer a sobrevivência declinar em função do tamanho do corpo, dado que temos uma função exponencial, nos trocamos o a_S-b_Sx por a_S+b_Sx, outra coisa que pode ser confuso aqui, é falar que é linear e a equação ser uma exponencial, veja que mesmo que tenhamos uma chance p de morrer no tempo 1, no tempo 2, vai ser o p da primeira e o p da segunda, no tempo 3, vai ser o p do tempo 1, o p do tempo 2 e o p do tempo 3, por isso aquele e elevado a equação ali.

4. O Fitness, W, é o parâmetro maltusiano r, pegando o r para ser a medida de fitness, a função de fitness é dada pela solução da equação característica:

 \int_1^\infty e^{-rt} (a_F+b_Fx)e^{-(a_S+b_Sx)t} dt =1

onde o valor inicial da integral é 1, ja que essa é a primeira idade reprodutiva.
Os dois expoentes podem ser absorvidos em um único termo, ja que os dois são “e” elevado a algo, dando

 \int_1^\infty (a_F+b_Fx)e^{-(a_S+b_Sx+r)t}dt=1

Agora, a equação acima tem a mesma forma geral da equação da da parte 4, e a integração dela vai ser mais ou menos a mesma coisa, a gente so esta adicionando a parte do crescimento populacional exponencial.

 1= (a_F+b_Fx) [\frac{-e^{-(a_S+b_Sx+r)t}}{a_S+b_Sx+r}]_1^\infty = 0 + \frac{a_F+b_Fx}{a_S+b_Sx+r} e^{-(a_S+b_Sx+r)}

Plotando a função do fitness

Para plotar r, nossa medida de fitness, como uma função de x, nos precisamos resolver a integral anterior ou a função acima para cada valor de x, ao gosto do cliente. Podemos primeiro examinar os métodos de estimar o r sem nos mesmos resolvermos a integral, e depois usando a integral resolvida.

Usando o R para fazer a integração.

Nos primeiro consideramos a integração numérica para estimar r. A estratégia é de iterar pelo espaço de x e para cada valor e calcular r.

1. integral produz a função a ser integrada, o  (a_F+b_Fx)e^{-(a_s+b_Sx+r)t}

1
2
3
4
5
6
7
8
9
10
11
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}

2. cal_integral chama a função do R integrate para obter a integral, passando ela para o integral que acabamos de definir

1
2
3
4
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}

3. A função encontra_raiz usa a função do R uniroot para achar o valor de r que satisfaz a equação característica para um dado x.

1
2
3
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}

4. Agora é só criar uma matriz de com uma única coluna valores de x e usa a função apply do R para usar encontra_raiz em cada valor de x (linha da matrix x) para pegar o valor requisito de r, guardando no vetor r.

1
2
3
4
5
6
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)

5.O vetor x e r são então usados para plotar a relação entre r(fitness) e x(tamanho do corpo). Rm sequencia, para um dado x, RCALC chama uniroot que chama INTEGRAL que então chama INTEGRAND

1
2
#Figura
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)

Usando nossa solução para a integral

Nos podemos usar a nossa solução de integral também, é mais simples no quesito de código, entretanto, como muitas funções não podem ser integradas analiticamente, ela é menos geral, além do que estamos assumindo que você consegue integrar bem as funções, eu sou péssimo com integrais, saindo daquilo que são os exemplos “triviais”.

#Fazendo a figura com nossa equação integrada

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#Nossa função de integral
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)

E então temos nossa figura:

01

E temos ai, como o fitness varia em relação ao tamanho do corpo, podemos ver que temos aquele bom e velho ponto máximo, o melhor tamanho de corpo para se ter, agora é so iniciar a busca para saber qual valor é esse.

Encontrando o máximo usando cálculo.

Nossa medida de fitness não está mais em um lado da equação, e ela não pode ser simplificada para ser assim. Mas ainda podemos usar a diferenciação implícita ja que o r depende do x. Por “conveniência” (porque senão não anda), nos usamos logaritmos para converter a função em uma que é aditiva.

 0 = ln(a_F+b_Fx)-(a_S+b_Sx +r) - ln(a_S+b_Sx+r)

Veja que temos 3 termos separados por um sinal de menos

 0 = T1 - T2 - T3

Agora, trabalhando com cada termo separadamente temos:

T1 :  \frac{dr}{dx} = \frac{b_F}{a_F+b_Fx}
Esse é fácil, não tem r na parada

T2:  \frac{dr}{dx} = b_S + \frac{dr}{dx}
Aqui veja que temos o r e x

T3:  \frac{dr}{dx} = ( \frac{1}{a_S+b_Sx+r} )(\frac{dr}{dx}) + \frac{b_S}{a_S+b_S+r}

Juntando todos os termos, ficamos com

 \frac{dr}{dx} (1+ \frac{1}{a_S+B_Sx +r}) = \frac{b_F}{a_F+b_Fx} -b_S - \frac{b_S}{a_S+b_Sx+r}

Agora quando  \frac{dr}{dx}=0 , toda a parte da direita se iguala a zero, dado que o termo no parenteses do lado esquerdo não se iguala a zero também, o que geralmente não sera o caso. Assim podemos rearranjar a parte da direita para fazer o r uma função de x:

 r^*= \frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*} -b_Sx^* - a_S

onde r^* é o máximo valor de r, obtido no x^*. Para achar x^* nos substituímos essa equação la na primeira com os logaritmos.

 0=ln(a_F+b_Fx^*)-(a_S+b_Sx^*+r^*)-ln(a_S+b_Sx^*+r^*)

 0=ln(a_F+b_Fx^*)-[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b^Sa_F-b_Sb_Fx^*} -b_Sx^*-a_S]  -ln[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*}-b_Sx^*-a_S]

que pode ser finalmente resolvido numericamente, mas deus me livre ter que fazer isso na ponta do lápis, acho que desde a parte 2, já não rola mais heheh. Mas o R ta ai pra isso.

Usando uniroot

A primeira abordagem é o uniroot. Note que os limites são setados ligeiramente perto dos valores requeridos. Se os limites estiverem muito longes um do outro () a função pode falhar, retornando uma mensagem de erro.

> uniroot(f=funcao,interval=c(1.2,3))$root Error in uniroot(f = funcao, interval = c(1.2, 3)) : f.upper = f(upper) é NA Além disso: Warning message: In log(As + Bs * x + r) : NaNs produzidos

Esse erro mostra a importância do plot original, no “olhômetro”, da para ver quais seriam valores razoáveis e mais próximos do “ótimo”.

Então definimos a função que chegamos usando o calculo

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As # r from eqn (2.40)
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}

e usamos o uniroot

> uniroot(f=funcao,interval=c(1.2,1.8))$root [1] 1.389974

Usando o nlm

Uma forma alternativa é usar o nlm, usando o valor absoluto na função, que no caso o mínimo tem de ser zero.

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}

Veja que o que mudou para a função anterior é somente o abs de valor absoluto ali na saída.
agora usamos a função nlm

> nlm(funcao, p=1.2)$estimate [1] 1.389943 Warning messages: 1: In log(As + Bs * x + r) : NaNs produzidos 2: In nlm(funcao, p = 1.2) : NA/Inf substituido pelo máximo valor positivo

Como x não é restrito a uma região, esse método não é satisfatório e apesar de termos a resposta correta, uma mensagem de aviso é gerada.

Usando optimize

Usar optimize é bem simples no entanto, basta usarmos a mesma função do nlm. Sendo que aqui definimos o mínimo e máximo, que é fácil estimar a partir da figura do fitness.

> optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum [1] 1.389952

O intervalo de nlm de busca é muito grande, veja também que apesar de muito próximos, encontramos valores diferentes, diferentes a partir da quinta casa depois do zero, mas valores diferentes.

Usando métodos numéricos

Nos vamos ver duas abordagens aqui, a primeira, usando a equação integral resolvida e a segunda usando métodos númericos para a integração do modelo original. Esse é o caso mais geral usando métodos númericos mas também o que consome mais tempo. Mesmo que a integral possa ser resolvida, essa é uma boa prática para checar se a solução da integral está realmente correta.

Usando a função integrada.

Aqui nos usamos a função optimize, dando a ele a função que calcula a raiz (raiz_funcao) usando uniroot para achar a raiz da função especificada em “funcao”.

1
2
3
4
5
6
7
8
9
10
11
12
13
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}

E assim nos temos

> optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum [1] 1.389946

Usando integração númerica para a função.

Lembrando que a gente ja teve uma noção do que é integração númerica aqui e aqui, nos podemos pular a parte da integral feita a mão (caso lidemos com um modelo difícil) e integrar e depois achar o valor de r em função de x e depois maximizar o fitness. Como dito anteriormente, mesmo com uma função integrável, essa pode ser uma garantia que tudo está correto.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}

E agora usamos o optimize para achar o máximo da função de fitness.

> optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum [1] 1.389934

Veja que aqui, estamos começando a esbarrar com problemas de precisão númerica, que não tem muito efeito por enquanto, mas pode vir a ser um problema.

Bem é isso ai, 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. Esse livro do Derek Roff é realmente muito legal, estamos ainda no segundo capítulo para quem quiser acompanhar, e queria eu ter a oportunidade de desde a graduação, estar exposto a um conteúdo desses, mas nunca é tarde para estudar e aprender coisas novas. Nas próxima partes, os modelos vão começar a ficar realmente legais com a entrada de variação e mais variáveis a maximizar. Até a próxima espero eu que a próxima chegue mais rápido.

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

R. A. Fisher 1915 The evolution of sexual preference Eugenics Reviews 7(3): 184–192.

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
#Plotando a função de fitness
 
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função para integrar a equação e retornar os valores
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}
 
 
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}
 
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)
#Figura
jpeg("01.jpg")
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)
dev.off()
 
 
#Fazendo a figura com nossa equação integrada
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)
 
#####################################
## Achando o maximo com calculo
#####################################
 
##uniroot
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}
#
uniroot(f=funcao,interval=c(1.2,1.8))$root
 
## Usando nlm
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}
#
nlm(funcao, p=1.2)$estimate
 
## Optimize
optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum
 
#####################################
## Usando métodos númericos
#####################################
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum
 
#####################################
## Usando integração númerica
#####################################
 
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum

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

Evolução – Fisherian Optimality Models – Parte 3

E vamos continuar olhando os modelos ótimos baseados na ideia de Fisher. A gente já viu dois exemplos aqui, na parte 1 e parte 2, agora vamos ver um caso um pouco mais complicado, do ponto de vista do cálculo pelo menos, ao mudar somente um pressuposto matemática em relação a parte 2, mais especificamente o pressuposto 3.

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 = \sum\limits_{t=1}^\infty F_{t} S_t = (a_F + b_Fx)e^{-(a_S+b_Sx)t}

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.

Gráfico da função de fitness

Antes de tentar encontrar um valor ótimo, nos primeiro plotamos a função do fitness para verificar que há um ponto máximo, que a função não tem uma forma esquisita de forma que faça as rotinas para encontrar o ponto máximo falharem, independente do ponto de inicio. Nos vamos assumir os mesmo parâmetros do exemplo da parte 2 (a_F=0, b_F=4,a_S=1, b_S=0.5), assim a função

W = \sum\limits_{t=1}^\infty 4xe^{-(1+0.5x)t}

O somatório acima vai até o infinito, o que geralmente não é uma opção para a análise numérica. Assim, primeiro, nos temos que decidir para quantas gerações, n, nos precisamos considerar esse somatório. Para isso nos setamos o tamanho do corpo, x, em algum valor arbitrário mas razoável, vamos dizer x=1, e fazemos o seguinte:

Geramos uma sequencia de inteiros de 1 a n e atribuímos isso a um vetor chamado idade.

A partir dai nos criamos outro vetor chamado
Wt, que é o componente especifico da idade da equação, 4xe^{-(1+0.5x)t}

Por último nos calculamos a soma do vetor Wt usando a função sum.

1
2
3
4
5
6
Somatorio <- function(n,x=1) {
    #print(x)
    Idade <- seq(from=1, to=n)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}

Então para o peso = 1, temos a seguinte figura usando este código:

1
2
3
4
nmax <- 20
n <- matrix(seq(from=1, to=nmax))
W <- apply(n,1,Somatorio)
plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3)

01

O somatório rapidamente se aproxima de sua assintota, seu valor assintótico, assim se fizer uma somatório até a idade 20, esse valor deve ser adequado para todos os valores de x.

Apenas para ilustrar um pouco mais, o que aconteceria se escolhêssemos outros valores de x, vamos fazer mais algumas tentativas.

Usamos o seguinte código:

1
2
3
4
5
6
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
for(i in 1:4) {
    W <- apply(n,1,function(n) {Somatorio(n,x=i)})
    plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3,
         main=paste("Tamanho do corpo",i),frame=F)
}

Um pequeno detalhe, se alguém ficou curioso o porque do print na função somatório, quando eu fiz esse for, eu fiquei preocupado com o escopo do código, se ele ia pegar o i certo, dai imprimi o i para garantir, veja que dentro do for o R pegou o valor de i que estava neste escopo e o n do apply.

02

Apesar de diferente o valor da assintota, a gente sempre tende pra ela. Por curiosidade, a gente pode se perguntar, porque diabos isso acontece? Bem basta olhar que a função que está no pressuposto 3, veja que a mortalidade aumenta com o tempo, ou seja a sobrevivência diminui ao longo do tempo, então a cada tempo, a cada tempo que se passa, a contribuição da idade para a sobrevivência é menor, um ano a mais para alguém “velho” não faz mas muita diferença, porque ele já vai morrer e não vai mais contribuir para mais descendentes. Podemos olhar como a função que estamos somando, como ela é para cada tempo.

03

Então conseguimos nos convencer que não é preciso somar a função de fitness até o infinito, até o 20 ta bom, ja que a cada novo valor do somatório, a gente soma um numero menorzinho, menorzinho, cada vez mais insignificante, então procedemos para tanto.

04

E olhando assim, da para ver que temos um ponto máximo para o tamanho do corpo, agora é só encontrá-lo.

Encontrando o máximo usando cálculo.

Agora que nos temos certeza que há um fitness máximo, nos sabemos seu valor aproximado pelo olhômetro da figura anterior e nos podemos procurar esse valor. O somatório nesse caso é resolvível, e para tanto nos vamos tentar chegar a solução exata usando cálculo.
Uma série que ocorre frequentemente em modelos de historia de vida, como o que temos aqui, são as séries geométricas.

\sum\limits_{t=1}^\infty = 1+a+a^2+a^3+ \dots = \frac{1}{1-a}

onde \left| a \right| < 1[/latex] (valor absoluto é menor que 1). No nosso caso, a função de mortalidade é uma série geométrica. Para simplificar, vamos estabelecer que [latex]a_F + b_F x = A[/latex] e [latex]a_S+b_S x = B[/latex], o que nos da [latex]W= \sum\limits_{t=1}^\infty F_tS_t = \sum\limits_{t=1}^\infty (a_F+b_Fx)e^{(a_S+b_Sx)t} = \sum\limits_{t=1}^\infty (A)e^{(B)t}[/latex] E vendo dessa forma, como é uma multiplicação ali dentro do somatório, podemos arrancar ele fora, como fizemos na parte 2 [latex]W = A \sum\limits_{t=1}^\infty e^{(B)t}[/latex] Agora para mudar essa equação, para a formula fechada do série geométrica, note que [latex]e^{-Bt}[/latex] pode ser escrito como [latex]e^Be^{B(t-1)}[/latex] e assim podemos arrancar esse cara também para fora do somatório [latex]W = A \sum\limits_{t=1}^\infty e^{(B)t} = Ae^{-B} \sum\limits_{t=1}^\infty e^{-B(t-1)}[/latex] E agora mudamos o somatório para a formula fechada [latex]W = Ae^{-B} (\sum\limits_{t=1}^\infty e^{-B(t-1)}) = Ae^{-B} (\frac{Ae^{-B}}{(1-e^{-B})})[/latex] Agora nos temos uma função "simples", simples no sentido que ela é diferenciável, porque ta complicado esse monte de equação. Mas é possível seguir em frente e diferenciar para achar o 0 da derivada. Podemos usar a regra da cadeia para encontrar a derivada, que fica. [latex]\frac{dW}{dx} = \frac{(b_F)e^{-(a_S+b_Sx)}}{1-e^{-(a_S+b_S x)}} + \frac{(a_F+b_F x)(-b_S e^{-(a_S+b_S x)})}{1-e^{-(a_S+b_S x)}} + \frac{(a_F+b_Fx)e^{-(a_s+b_Sx)}(-1)(b_S e^{-(a_S+b_Sx)})}{(1-e^{-(a_S+b_S x)})^2} [/latex] Da para destrinchar passo a passo, mas vamos deixar isso para outro post. Agora nos precisamos achar o valor de x no qual [latex]\frac{dW}{dx}=0[/latex], de forma que podemos simplificar essa equação final da derivada um pouco já que [latex]\frac{e^{-(a_S+b_Sx)} }{ 1-e^{-(a_s+b_Sx)}}[/latex] é um elemento comum aos três termos que temos, assim a equação que temos que resolver é [latex](b_F) + (a_F+b_Fx)(-b_S) + \frac{(a_F+b_Fx)(-1)(b_Se^{a_S+b_Sx})}{1-e^{a_S+b_Sx}} = 0[/latex] E para resolver isso no R podemos usar a função uniroot. Agora é simples, é só escrevermos essa função no R, usando os valores que definimos para os parâmetros e chamar a função uniroot. Então definimos a função

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

Como nos estamos apenas interessado nos valores positivos, pesos negativos não fariam nenhum sentido biológico, o limite para a busca da raiz é colocado como 0, excluindo possíveis raízes negativas, Se houver duas raízes positivas, a raiz mais próxima do menor limite sera encontrada e devolvida pela função. Nos sabemos que para esse caso, a partir da figura do fitness, que existe apenas uma raiz positiva, e assim não precisamos procurar por mais valores.

Aproximação numérica

Agora, vamos dar uma olhada como seria para encontrar o valor máximo a partir de uma aproximação numérica.

Nos vamos usar o mesmo somatório anterior, exceto pelo fato que vamos tomar o resultado como um valor negativo, porque a função nlm do R acha o mínimo de uma função, e ao colocar o resultado negativo, nos deixamos a função de ponta cabeça, mas ainda acharemos o valor que nos interessa. Precisamos de uma valor inicial também para iniciar a busca, e setamos esse valor como 1, já que sabemos, no olhômetro do gráfico, que esse é um valor razoável.

Então alterando a função para retornar menos o somatório

1
2
3
4
5
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(-sum(Wt))
}

E encontrando a estimativa de valor ótimo para o peso

> nlm(Somatorio, p=1) $minimum [1] -1.268755 $estimate [1] 1.68281 $gradient [1] -2.111179e-09 $code [1] 1 $iterations [1] 6

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, principalmente nas equações, que é fácil passar error, ou mande um e-mail e até mais. Olhando a solução por métodos numéricos, realmente da uma tentação em ignorar o uso do cálculo, mas não custa nada olhar como faz. O livro do Roff de evolução é realmente muito legal e não terminamos nem o inicio, mas é devagar e sempre.

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
68
69
##################################
## Figuras
##################################
Somatorio <- function(n,x=1) {
    #print(x)
    Idade <- seq(from=1, to=n)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}
 
nmax <- 20
n <- matrix(seq(from=1, to=nmax))
W <- apply(n,1,Somatorio)
 
#Figura 1
plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3)
 
 
#Figura 2
layout(matrix(1:4,nrow=2,ncol=2,byrow=T))
for(i in 1:4) {
    W <- apply(n,1,function(n) {Somatorio(n,x=i)})
    plot(n,W,type="l", xlab="Idade, n", ylab="Somatório do peso, Wt", las=1,lwd=3,
         main=paste("Tamanho do corpo",i),frame=F)
}
 
#Figura 3
curve(4*1*exp(-(1+0.5*1)*x),0,2,frame=F)
 
#
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(sum(Wt))
}
 
 
x<- matrix(seq(from=0,to=5,length=100))
W <- apply(x,1,Somatorio)
 
#Figura 4
plot(x,W,type="l", xlab="Tamanho do corpo, x", ylab="Fitness, W",las=1,
     lwd=2,lty=2,col="red",frame=F)
 
 
##################################
## Usando Calculo
##################################
funcao <- function(x){
    (4)+
        (0+4*x)*(-0.5)+
            (0+4*x)*(-1)*(0.5*exp(-(1+0.5*x)))/ (1-exp(-(1+0.5*x)))
}
 
#Raiz
solucao <- uniroot(funcao, interval=c(0,4))
solucao$root
 
##################################
## Usando aproximação numérica
##################################
 
Somatorio <- function(x) {
    Idade <- seq(from=1, to=20)
    Wt <- 4*x*exp(-(1+0.5*x)*Idade)
    return(-sum(Wt))
}
 
nlm(Somatorio, p=1)

Evolução – Fisherian Optimality Models – Parte 2

Bem a gente começou a ler sobre evolução aqui.

Agora vamos adicionar um pouco mais de complexidade ao primeiro cenário que avaliamos, mas será que maior complexidade implica em mudança? Nos vamos aumentar a complexidade introduzindo uma estrutura de idade, que é algo que podemos esperar que mude a solução anterior. Mas de fato tudo vai continuar na mesma. A introdução da complexidade não vai mudar tanto as conclusões tanto qualitativa quanto quantitativamente.

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. Sobrevivência até a idade da primeira reprodução, sendo 1 aqui, diminui linearmente com o tamanho do corpo e é constante por unidade de tempo e independente do tamanho do corpo.
A sobrevivência na idade t é dada por:

S_t=(a_S-b_{S}x)e^{-M(t-1)}

onde M é a taxa de mortalidade instantânea que é independente da idade.

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

W = R0= \sum\limits_{t=1}^\infty F_{t} S_t

Logo:

W= \sum\limits_{t=1}^\infty (a_F+b_{F}x)(a_S+b_{S}x)e^{-M(t-1)}

Apesar da equação acima parecer muito mais complicada que a equação do cenário 1, ela na verdade não contribui com nenhuma nova informação que muda o tamanho ótimo, já que os componente que dependem do tamanho não dependem da idade e assim podem ser movidos para fora do somatório.

Observe o seguinte, suponha a seguinte somatória:

 \sum\limits_{x=1}^{3} 3x
Isso é o mesmo que
 (3 \cdot 1) + (3 \cdot 2) + (3 \cdot 3)  = 18
Podemos por o 3 em evidência, tirando ele de dentro de cada soma certo?
 3 \cdot ((1) + (2) + (3))  = 18
Agora esse 1+2+3 é o somatório certo?
 3 \cdot (\sum\limits_{x=1}^{3} x)  = 18

Logo para nossa equação:

W= \sum\limits_{t=1}^\infty (a_F+b_{F}x)(a_S+b_{S}x)e^{-M(t-1)}

Fazemos o mesmo e tiramos a parte que não depende do tempo pra fora do somatório

W=  (a_F+b_{F}x)(a_S+b_{S}x) \cdot \sum\limits_{t=1}^\infty e^{-M(t-1)}

E podemos desenvolver essa parte como no cenário 1

W=  (a_F a_S -b_F b_S x^2 + (a_S b_F - a_F b_S)x) \cdot \sum\limits_{t=1}^\infty e^{-M(t-1)}

Dai a gente sabe que aquela somatória vai dar um valor, sei la qual, um valor que vamos chamar de constante, ja que não muda em função de x, e ele vai multiplicar toda nossa equação

W=  (a_F a_S -b_F b_S x^2 + (a_S b_F - a_F b_S)x) \cdot constante

E o que acontece com o valor ótimo quando a equação de fitness é multiplicada?

Bem a gente tem o mesmo valor ótimo, so teremos um intervalo mais estreito de valores viáveis.

Veja uma figura de um exemplo mais simples, imagine que a gente tivesse chegado na equação
W =  -Tamanho^2

O valor ótimo seria 0

01

Agora com a equação

W =  (-Tamanho^2) \cdot constante

A gente continua com o mesmo valor ótimo

02

O que nos interessa é apenas o valor ótimo aqui nessa análise. Consequentemente nos temos que a equação de fitness especificada é a mesma do cenário 1 multiplicado por uma constante, assim o tamanho ótimo continua o mesmo.

Esse é um exemplo interessante de pensar, porque muitas vezes conseguimos avaliar algumas características na evolução independente de outras, nem sempre é assim, mas aqui é um exemplo de que é plausível.

Bem é isso ai, apos um grande lag entre o ultimo post e esse, mais alguma coisa para ler, 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.

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
x<-seq(-2,2,0.1)
jpeg("01.jpg")
plot(x,-x^2,type="l",frame=F,ylab="Fitness",xlab="Tamanho do corpo",main="Sem constante multiplicando",ylim=c(-13,0))
abline(v=0,lty=3)
text(-1,-10,expression("Fitness=-Tamanho^2"))
dev.off()
 
jpeg("02.jpg")
plot(x,3*-x^2,type="l",frame=F,ylab="Fitness",xlab="Tamanho do corpo",main="Com constante multiplicando",ylim=c(-13,0))
abline(v=0,lty=3)
text(-1,-10,expression("Fitness=-Tamanho^2"))
dev.off()

Evolução – Fisherian Optimality Models – Parte 1

Com certeza, uma disciplina muito legal dentro da biologia é a evolução.

Como diz o celebre trabalho do Theodosius Dobzhansky de 1973:

“Nothing in Biology Makes Sense Except in the Light of Evolution”

“Nada na biologia faz sentido a não ser sobre a luz da evolução”

Mas nem tudo é alegria, o que essa frase anima a gente, as vezes outras são meio perturbadoras, como a do John Maynard Smith:

“If you can’t stand algebra, keep out of evolutionary biology”

“Se você não suporta álgebra, fique fora da biologia evolutiva”

Essa frase estava no inicio do livro Algebraic Statistics for Computational Biology do Lior Pachter e Bernd Sturmfels. Eu nunca dei conta de ler esse livro inteiro, porque é meio tedioso as vezes. Mas recentemente eu comecei a ler um livro muito legal chamado Modeling Evolution – an introduction to numerical methods do Derik a. Roff.

É o tipo de livro que eu acho legal, porque no próprio prefácio, o autor fala que o livro veio de uma disciplina que ele dava em Harvard, mas é o livro que ele queria ter lido na graduação. Essa são palavras chave pra gente ver um livro que o cara escreveu com gosto, e que não pega nem tão pesado na biologia como na matemática.

O livro coloca bem como começar do zero modelos em evolução, e como ir procedendo, desde os mais simples ao mais complexos. Eu pelo menos nunca tinha estudado muito a fundo como é a análise de modelos em evolução, o que é ruim porque a gente acaba simplesmente acreditando em tudo que le, ja que a análise passa a ser uma barreira.

Mas já que estamos aqui, vamos dar uma olhada em como são os modelos mais básicos, da forma como o Ronald Fisher começou a analisar as coisas em evolução, depois das explicações do livro, deu para dar uma melhorada. O próprio Darwin deu uma baqueada no fim da sua vida, pensando que sua teoria da evolução tinha talvez muitas coisas que não conseguia explicar, mas depois veio o Fisher e companhia com a teoria neodarwinista e começaram a matar a charada.

Os Fisherian Optimality Models são baseados na abordagem do Fisher para estudar questões evolutivas nas quais ele pegou o parâmetro malthusiano, r, como a medida de fitness. Essa medida pode ser ruim quando uma distribuição etária estável não pode ser assumida, existe uma dependência da densidade ou a frequência é um fator importante na evolução das características em questão, a população é caótica ou o fitness depende de interações sociais. Mas mesmo com tantas restrições, esse tipo de modelo é muito interessante para a análise de questões evolutivas e na geração de predições testáveis, alias essa é uma critica comum quando a pesquisa no Brasil, trabalhos que não tem predições a priore, não testam hipóteses, é coletar dados e ver o que da pra fazer. Agora da para tentar começar a usar esse tipo de abordagem para formular hipóteses, predições a priore.

Mas continuando, três dos parâmetros populacionais mais comumente usados nesse tipo de modelo são a taxa de crescimento populacional r, sucesso reprodutivo ao longo da vida (ou sucesso reprodutivo liquido) R0 e um componente do fitness o qual na sua maximização também maximiza o fitness. Dos três primeiros, o primeiro é o mais comumente usado. O parâmetro malthusiano é a taxa de crescimento da população que está com uma distribuição de idades estável (vou ficar devendo um post sobre isso), e é dado pela equação:

\int_0^{\inf} e^{-rt}l(x)m(x)dx=1

Onde o t é a tempo, l(x) é a taxa de sobrevivência no tempo x, e m(x) é o número de nascimentos por fêmea na idade x (normalmente a razão sexual é 1:1, mas isso não é obrigatório). Na ausência da dependência de densidade, é meio obvio que mutações que aumentem r devem aumentar sua frequência como a gente ja viu aqui. Se existe uma dependência da densidade, então a medida apropriada de fitness é o número relativo de indivíduos que passam do período no qual a dependência da densidade age. Por exemplo, suponha que a dependência da densidade ocorra na fase imatura e que nos queremos calcular o fitness de dois tipos de indivíduos, A e B. Se o número de adultos de A excede o número de adultos de B então A é necessariamente o com melhor fitness. A análise de tal situação pode se tornar mais complicada se o sucesso nesse período não depende simplesmente da densidade, mas também na frequência dos tipos na população.
Se a população não aumenta em tamanho, uma medida plausível do fitness é a taxa reprodutiva liquida.

R_0 = \int_0^{\inf} l(x)m(x)dx

Essa medida e algumas vezes é chamado de expectativa de sucesso reprodutivo ao longo da vida (“expected lifetime reproductive success”). A premissa aqui é que se A tem uma R_0 maior que B então A vai ter um fitness maior. Dependência da densidade é assumida implicitamente como não afetando as características de interesse. Quando existem flutuações estocásticas que movem a população para fora de uma estrutura estável de idade, nem r nem o R_0 podem ser assumidos como medidas apropriadas para medir o fitness.

Em alguns casos, particularmente estudos de comportamento, nosso interesse é em uma característica em particular e para simplificar a análise nos assumimos que a maximização de alguns componentes da historia de vida é equivalente a maximização do fitness. Por exemplo, nos podemos querer determinar a alocação ótima de tempo de forrageamento entre manchas que variam quanto a qualidade de seus recursos, a gente comentou sobre isso aqui, aqui e aqui. Mas uma premissa comum é que a maximização da obtenção de recursos por unidade de tempo é equivalente a maximização do fitness, mas precisamos de cuidado, porque ignorar outros componentes da historia de vida pode produzir conclusões erronias. Dois exemplos famosos são o paradoxo de Cole e a hipótese de Lack.

No seu clássico artigo de 1954, Cole sugeriu que havia um aparente paradoxo porque

“For an annual species, the absolute gain in intrinsic population growth which could
be achieved by changing to the perennial reproductive habit would be exactly equivalent to adding one individual to the average litter size”

“Para uma espécie anual, o ganho absoluto intrínseco no crescimento populacional que poderia ser obtido alterando o habito reprodutivo perenial seria exatamente equivalente a adicionar um individuo ao tamanho médio da ninhada”

Se a conclusão de Cole fosse verdadeira, nos esperaríamos que plantas perenes deveriam ser raras, o que certamente não é o que observamos. O problema foi que cole não incluiu diferenças na sobrevivência entre adultos e jovens, sua análise assumia implicitamente taxas de sobrevivência de 1 para ambos os estágios.

Ja Lack(1947) teve a ideia que em aves…

“the average clutch-size is ultimately determined by the average maximum number of young which the parents can successfully raise in the region and season in question”

“o tamanho médio da ninhada é ultimamente determinado pelo numero médio de jovens que os pais podem criar com sucesso na região em questão”

A hipótese de Lack assumia que apenas interações negativas dependentes de densidade entre irmãos dentro de um ninho eram importante para predizer o tamanho ótimo de filhos a se criar em um ninho para os pais, logo a quantidade de ovos que vemos em um ninho é o máximo que os pais podem criar. Lack não observou a possibilidade que a sobrevivência na fase a adulta dos filhotes da ninhada pode ser afetada pelo numero de irmãos com o qual cresceram. Experimentos manipulando o tamanho de ninhadas em aves, mamíferos, repteis, peixes e plantas demostrou efeitos negativos de uma ninhada maior na sobrevivência futura dos adultos e no seu sucesso reprodutivo. A incorporação de tais efeitos prediz que o tamanho de ninhada ótima vai ser menor que o valor proposto por Lack.

Métodos de análises:

O foco da análise é em condições de equilíbrio e não na trajetória evolucionaria até ele (mas a trajetória pode ser analisada também, espero abordar isso depois quando ler todo o livro de Evolução). Aqui a dependência da densidade não é explicitamente considerada. Dependência da frequência é também assumida como ausente. A formulação do modelo que gostaríamos de chegar é

W = f( \theta_1, \theta_2, \dots, \theta_k, x_1, x_2, \dots , x_n)

onde W é o fitness, \theta_1,\theta_2,\dots,\theta_k são os parâmetros e x_1 ,x_2 , \dots , x_{n} são as características.
Essa função é para ser lida como Fitness é uma função de k parâmetros e n características (“Fitness is a function of k parameters and n traits.”).

Uma regra geral é manter o número de parâmetros e características no mínimo. Suponha que tenhamos 5 parâmetros e queremos avaliar a performance de todas as combinações, dividindo cada parâmetro em 10 partes, que é uma número razoável de divisões, isso nos da 10^5, 100000 combinações para examinar. Apesar de ser possível (usando o computador), essa pode não ser o caminho preferível, se puder ser evitado.

A função de fitness vai invariavelmente ser feita de um número de funções componentes, como veremos no exemplo daqui a pouco, no qual o fitness é o produto da fecundidade e da sobrevivência e ambas são função do tamanho do corpo. Essas funções podem produzir curvas suaves de fitness que podem ser analisadas usando calculo (integral e diferencial) ou a função pode ter descontinuidades ou ser difícil de ser analisadas por algum outro motivo (uma superfície irregular) usando o calculo.

É preciso cuidado também examinando modelos com descontinuidades e lugares no qual os componentes do modelo podem tomar valores fisicamente impossíveis (irrelevantes do ponto de vista biológico). No exemplo que a gente vai olhar aqui é dado por uma função linear negativa do tamanho do corpo, o que significa que abaixo de um tamanho particular a sobrevivência vai ficar negativa, o que não é possível no mundo real. No exemplo isso não é um problema porque o fitness será negativo nesse caso. Entretanto, suponha que o modelo contenha o produto de duas funções de sobrevivência que podem ser matematicamente menores que zero. Agora é possível que os dois valores negativos produzam um valor positivo e o valor de fitness não pode ser visto como incorreto.

No fim, existem três casos gerais em que acabamos sempre caindo.

Primeiro: fitness pode ser escrito como uma função de características de interesse e a diferenciação não é um problema.

Segundo: fitness pode ser escrito como uma função de características de interesse mas a diferenciação é problemática.

Terceiro: fitness não pode ser escrito como uma função na qual o fitness pode ser isolado.

Métodos de análises:
W=f(\theta_1,\theta_2,\dots,\theta_k,x_1,x_2,\dots,x_n) é bem comportado.

Assumindo uma função de fitness diferenciável (uma superfície suave de fitness), nos achamos o valores ótimos para as características diferenciando a função com respeito as características e igualando o tudo a zero. Por exemplo, suponha que a função de fitness é uma função quadrada.

W=f(\theta_1,\theta_2,\theta_3,x) = \theta_1 + \theta_2x + \theta_3 x^2

Então:

\frac{dW}{dx} = \theta_2 + 2 \theta_3 x

\frac{dW}{dx} =  0 quando  x = \frac{-\theta_2}{2 \theta_3}

Note que um fitness intermediário requer que  \theta_3 \textless 0 e \theta_2 > 0 . Podemos ficar tentados a esperar apenas que \frac{-\theta_2}{2 \theta_3} \textless 0, entretanto, a condição \theta_3 > 0 e \theta_2 \textless 0 define uma função convexa de forma que o ponto de flexão seja um mínimo e não uma máxima. Sempre que possível, a função deve ser plotada numa figura para garantir que um ponto apropriado de flexão exista. Se a função é complexa e contem muitos termos, a diferenciação pode se tornar muito tediosa e cuidado é necessário para que termos não sejam acidentalmente omitidos ou sinais trocados.

Código de R para derivadas:

A gente já viu sobre derivadas aqui, mas vamos olhar o código do livro para seguir a linha de raciocínio.

Então suponha que a gente quer derivar uma equação de segundo grau

a+bx+cx^2

No R a gente declara ela como uma formula usando o ~(tio) e usa a função deriv

> y <- deriv(~a+b*x+c*x^2,"x") > y expression({ .value <- a + b * x + c * x^2 .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- b + c * (2 * x) attr(.value, "gradient") <- .grad .value })

A derivada é dada na linha .grad[, “x”] <- b + c * (2 * x). Para calcular o valor em algum gradiente, a gente precisa de um pouco mais de código. Em alguns casos a saída é dividida em expressões separadas. Por exemplo ao derivar a expressão [latex]e^{ax}+bx+cx^2[/latex] A sintaxe no R é a mesma usando a função deriv

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

Mas aqui a derivada é composta, veja que a derivada

.grad[, "x"] <- .expr2 * a + b + c * (2 * x)

precisa da expr2 que está aqui

.expr2 <- exp(a * x)

O primeiro cenário:

Esse cenário ilustra a análise no caso em que não existe uma estrutura de idade e o fitness é simplesmente uma função quadrada da caraterística sendo estudada. Por causa da sua simplicidade, a análise e visualização são facilmente completadas. Apesar da simplicidade, é valido começar com algo simples e ir aprimorando, onde podemos adicionar complexidade aos poucos. Essa abordagem permite visualizar a importância de algumas premissas em particular, por exemplo aqui consideramos o caso de organismos semelparos, que a característica sendo analisada varia positivamente com um componente do fitness e negativamente com o outro componente. A priori, podemos ser tentados a assumir que isso vai necessariamente resultar em uma função que tem um máximo em algum valor intermediário da característica. Mas isso não é necessariamente verdade, suponha por exemplo duas funções de fitness

y_1=e^{ax}
y_2=e^{-bx}

Onde y_1 e y_2 são características como fecundidade e sobrevivência, a e b são constantes e x é a característica sendo estudada, nesse caso o tamanho do corpo. Vamos supor que o fitness seja o produto de y_1 e y_2, então o fitness W é dado por

W = y_1y_2=e^{(a-b)x}

Que não tem um valor ótimo intermediário. No primeiro cenário, um valor ótimo não é garantido, mas ocorre para um particular set de valores de parâmetros.

Pressupostos gerais:

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

Os pressupostos acima descrevem uma situação bem geral. Para fazer uma previsão, nos precisamos converter esses pressupostos em expressões matemáticas. Mas mesmo antes nos podemos fazer uma afirmação geral. Como um componente do fitness, fecundidade, aumenta com o tamanho do corpo, enquanto o outro componente, sobrevivência, diminui com o tamanho do corpo, em muitas circunstancias, mas não todas, nos esperamos que um tamanho ótimo do corpo maximize o fitness. Para investigar isso, nos vamos converter esses pressupostos em expressões matemáticas.

Pressupostos matemáticos

1. Fecundação aumenta linearmente com o tamanho do corpo:

F=a_F + b_F x

onde a_F e b_F são constantes.

2. Sobrevivência diminui linearmente com o tamanho do corpo:

S= a_S - b_S x

3. O fitness, W, é a expectativa reprodutiva ao longo da vida, R_0, dado como o produto da fecundidade e da sobrevivência.

W = R_0 = FS
W = (a_F + b_F x)(a_S - b_S x)
W = a_F a_S -b_F b_S x^2 + (a_S b_F - a_F b_S)x

A equação acima descreve uma parábola que é concava para baixo, isso é, tem um valor máximo em algum tamanho de corpo, digamos x^*. Apesar de sabermos de antemão que existe um valor ótimo de tamanho de corpo, nos não sabemos se ele ocorre em um valor plausível de tamanho de corpo. O primeiro passo na análise é plotar W com x para confirmar visualmente que o máximo existe e mostrar que ele ocorre para a espécie em estudo.

Plotando a função de fitness

Normalmente é uma boa ideia plotar a função para examinar visualmente o seu comportamento. Enquanto nos sabemos que neste caso a função é quadrática e concava para baixo, nos não sabemos se o valor ótimo para a característica avaliada é plausível biologicammente dados os parâmetros usados. Assim nos definimos valores dos parâmetros, vamos supor que valores razoáveis são a_F = 0, b_F=4, a_S = 1 e b_S = 0.5.
A equação do fitness pode então ser escrita como:

W = 0 \cdot 1 - 4 \cdot 0.5 x^2 + (1 \cdot 4 - 0 \cdot 0.5)x
W = -2x^2 + 4x

E no gráfico:

01

A partir da figura nos podemos ver que é um máximo em torno de 1.


Achando o máximo usando calculo

Para obter x^* usando calculo nos diferenciamos a equação, igualamos a zero e achamos o valor de x que satisfaz a condição.
Nesse caso a diferenciação é muito simples

\frac{dW}{dx} = -4x+4

\frac{dW}{dx} = 0 quando -4x+4 =0

logo

x=1

Diferenciação simbólica no R

Diferenciação simbólica, ou também chamada de diferenciação algorítmica (porque tem algorítimo que pensa nisso pra você) no R pode ser feita usando deriv

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

O output é meio “bagunçado”, mas isso porque ele não é pensado como um sumário, mas como uma entrada funcional para outras funções.
Com essa saída, a gente pode determinar em que valor a derivativa é zero usando a função uniroot. Nos setamos a derivativa como uma função separada e chamada pela função uniroot.

1
2
3
4
5
6
7
8
# 
FUNC <- function(w) {
    y <- deriv(~-2*x^2+4*x,"x")
    x <- w
    z <- eval(y)
    d <- attr(z,"gradient")
    return(d)
}
> B <- uniroot(FUNC, interval= c(-2,4)) > B$root [1] 1

Achando o máximo usando uma aproximação por métodos númericos.

Em muitos casos, a função pode não ser assim tão facilmente diferenciada, por exemplo, a função pode consistir de uma modelo de simulação ou ela pode ter discontinuidades, nos vamos encontrar comumente, mas esse cenário pode servir de exemplo para ilustrar como é o procedimento.
As rotinas disponíveis tipicamente localizão o valor mínimo da função. No nosso caso queremos o máximo. Para usar rotinas de minimização nos podemos simplismente pegar o valor negativo da equação. Logo, ao invés de procurar o máximo de 4x-2x^2, nos procuramos o mínimo de -4x+2x^2. Existe muitas rotinas prontas no R para achar o mínimo, aqui a gente vai usar o nlm (uma altenativa ao optimize).

1
2
3
FITNESS <- function(x) {
    return(2*x^2-4*x)
}
> nlm(FITNESS, p=-2) $minimum [1] -2 $estimate [1] 0.9999995 $gradient [1] 2.220446e-10 $code [1] 1 $iterations [1] 1 > nlm(FITNESS, p=-2) $minimum [1] -2 $estimate [1] 0.9999995 $gradient [1] 2.220446e-10 $code [1] 1 $iterations [1] 1

Agora que temos o resultado, é só sair medindo indivíduos dessa espécie na natureza e ver se eles medem algo em torno de 1 para saber que a teoria está correta, senão alterar os pressupostos e começar denovo. Mas o legal é que temos agora um resultado esperado a priori. Tudo bem que faltou definir unidades direito, precisamos pensar também como são os parâmetros da espécie em questão, mas da pra começar a pensar na metodologia para tentar aplicar. Gostaria de ter tido uma cabeça melhor para ter estudado esse tipo de coisa na graduação. Mas Modeling Evolution – an introduction to numerical methods é um livro muito legal, ele nem existia antes mesmo, então não tinha como ler, mas conforme for lendo mais, vou ir postando aqui porque se eu não escrevo, eu não aprendo.

Bem é isso ai, 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.

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
y <- deriv(~a+b*x+c*x^2,"x")
y
 
#
y <- deriv(~exp(a*x)+b*x+c*x^2,"x")
y
 
#
x <- seq(0,2,length=1000)
W <- (-2*x^2 + 4*x)
jpeg("01.jpg")
plot(x,W,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, W",las=1,lwd=3,
     frame=F)
dev.off()
 
 
y <- deriv(~-2*x^2+4*x,"x")
y
 
# 
FUNC <- function(w) {
    y <- deriv(~-2*x^2+4*x,"x")
    x <- w
    z <- eval(y)
    d <- attr(z,"gradient")
    return(d)
}
 
#
B <- uniroot(FUNC, interval= c(-2,4))
B$root
 
 
FITNESS <- function(x) {
    return(2*x^2-4*x)
}
 
nlm(FITNESS, p=-2)