Duas espécies em metapopulações de Levins

Ja vimos metapopulações para uma espécie aqui, aqui e aqui, mas será que poderíamos adicionar mais uma espécie ao modelo de Levins?

Primeiro, vamos lembrar como é a metapopulação de levins.

\frac{dO}{dt}=cO(1-O)-eO

onde:

  • O – Ocupação
  • c – taxa de colonização
  • e – taxa de extinção

Então basicamente temos dois termos,

cN(1-N)

que faz a população crescer a uma taxa c, mas quanto maior a ocupação, mais propagulos são emitidos, mas veja que quanto maior a ocupação, menor o crescimento, que é a multiplicação por (1-N), conforme o N cresce, o valor tende a zero, porque cada vez teremos 1 menos alguma coisa, que da menos que um, estaremos multiplicando por uma fração menor.

eN

Cada mancha ocupada tem uma chance de extinção e, logo quanto mais manchas ocupadas, mais manchas veremos sendo extintas.

Agora vamos colocar outra espécie, ela vai ocupar as manchas vagas pela primeira espécie, ou seja, vagou, ela tem chance de entrar ocupar, mas a se a espécie um ocupar novamente, ela perde a competição.

Como são duas espécies, precisamos agora acompanhar a mudança na ocupação de cada uma, então precisamos de duas equações

E espécie 1 continua como a equação anterior
\frac{dO_1}{dt}=c_1O_1(1-O_1)-eO_1
mas a segunda fica
\frac{dO_2}{dt}=c_2O_2(1-O_1-O_2)-c_1*O_1*O_2-eO_2

Ou seja, a metapopulação 2 expande sua ocupação, diminuindo a expansão de acordo com o o quanto ja está ocupado pela ocupação da espécie 1 e dela mesmo, e tiramos as manchas que ela perde na competição, além do que ela perde com a extinção.

Agora, quem lembra do pacote desolve, vamos usar ele para ver como isso fica ao longo do tempo, então temos que transformar numa função para usar com a função ode.

1
2
3
4
5
6
7
8
9
10
11
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}

Agora temos que determinar por quanto tempo queremos observar a ocupação e os paramtros definir os parametros de ocupação e extinção

1
2
3
4
5
6
tempo <- seq(0, 200, by =5)
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25

E vemos o seguinte

01

A espécie um rapidamente domina o ambiente, retirando a espécie 2 da jogada.

Mesmo quando mudamos os parametros um pouco, deixando ambas com a mesma taxa de colonização, maior a taxa de extinção, lembrando que para uma espécie, quando a extinção ficava maior que a colonização, ela tendia a extinção, sendo a diferença o quão rápido isso ia acontecer.

1
2
3
4
5
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05

02

Agora se mudarmos o jogo para a espécie 1 extinguir, a espécie 2 consegue dominar, salvo que ela também não tenha a colonização menor que a extinção

1
2
3
4
5
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06

03

E existem também os casos onde podemos ter uma coexistencia

1
2
3
4
5
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05

04

Agora, as possibilidades de combinações de parâmetros são grandes, tornando mais complexo definir quando uma ou outra espécie ganha, ou quando temos a cooexistencia, mas podemos ter um feeling da coisa já. além disso, nossas equações garantem a vitoria da espécie 1, mas isso não precisa ser assim.

Bem é isso ai, lembrando um pouco de equações diferencias, que não falamos a um bom tempo aqui e vendo a parte essa expansão para a metapopulaçã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.

Referência:
A. A. de Oliveira e P. I. Prado – Coexistência em Metapopulações – Roteiro no EcoVirtual
Robert May & Angela McLean 2007 – Theoretical Ecology: Principles and Applications. Oxford University 272 pp

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
##
##install.packages("deSolve")
library(deSolve)
rm(list=ls())
 
##Função
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}
 
##Tempo
tempo <- seq(0, 200, by =5)
 
##Caso 1
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
 
##jpeg("01.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 2
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("02.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 3
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("03.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 4
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05
 
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("04.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()

Biogeografia de ilhas

MacArthur e Wilson propuseram a teoria de biogeografia de ilhas onde o número de espécies de uma ilha oceânica seria uma função da taxa de imigração e extinção de espécies. O número de espécies em qualquer tempo é um equilíbrio dinâmico, resultando de ambos os lentos processos de extinção e continua chegada de novas espécies. Assim, é esperado que a composição de espécies pode mudar ao longo do tempo, gerando um “turnover” de espécies.

Vamos pensar da taxa de imigração, y, como uma função do número de espécies, x, que ja está em uma ilha. Essa relação tem uma inclinação negativa, porque conforme o número de espécies aumenta, a chance de um novo individuo chegar ser de uma nova espécie fica mais baixa, lembre-se de como é a distribuição das abundancias das espécies aqui. Assim, a imigração vai ter seu maior valor quando a ilha está vazia, quando o x é igual a zero, e vai cair para zero quando todas as espécies possíveis tiverem chegado la. Além disso, a inclinação deve desacelerar porque algumas espécies são muito mais prováveis que outras de imigrar. Isso significa que a taxa de imigração decai abruptamente assim que as espécies mais prováveis de imigrar chegam, e depois só aquelas que tem pouca chance de dispersar, ou baixas abundancias estão faltando chegar. Denovo pensando nas distribuições de abundancias de espécies, como Preston observou, algumas espécies são simplesmente mais abundantes que outras, produzindo mais indivíduos para dispersar, sendo “melhor dispersores”
Ou seja, esperamos algo desse jeito.

01

Agora vamos imaginar a taxa de extinção, y, como uma função do número de espécies, x. Essa relação vai ter uma relação positiva com o número de espécies, de tal forma que a probabilidade de extinção aumenta com o número de espécies, fácil pensar nisso, quando mais espécies, mais a gente vai adicionando a chance de extinção delas a taxa geral de extinção.

Essa relação também é predita de ter uma inclinação cada vez maior, mais ou menos pelos menos motivos do declínio cada vez maior da imigração. Algumas espécies são menos comuns que outras, e assim mais fácil de se tornarem extintas devido a estocasticidade demográfica e ambiental, tipo um drift genético, e segundo algumas espécies vão ter um menor fitness por muitas razões. Com o acumulo de espécies, maior sera o número de especies propicias a extinção (em relação a baixa abundancia e fitness) presentes, e que vão se extinguir. Lembrando que estamos falando de extinção localmente. Porque lembro que uma vez que apresentei um seminário sobre meta-populações,a palavra extinção soava agressiva para as pessoas, mas veja que extinção aqui não significa o fim da especie, só o fim local.
Então localmente, a taxa de extinção vai ser algo assim:

02

Assim, a taxa de mudança do número de especies em uma ilha vai ser \Delta R, que vai ser a diferença entre imigração I e extinção E, ou

\Delta R = I-E

Quando \Delta R=0, nos temos um equilíbrio. Se soubermos a forma quantitativa da imigração e extinção, nos podemos resolver para o equilíbrio. Esse equilíbrio vai ser um ponto no eixo x, onde as duas taxas se cruzam.

Na teoria do MacArthur e Wilson de biogeografia de ilhas, essas taxas podem ser determinadas pelo tamanho das ilhas, onde:

  • Ilhas maiores tem menor taxa de extinção por causa do maior tamanho médio das populações
  • Ilhas maiores tem maior taxa de colonização porque são alvos maiores para a dispersão

A distância entre as ilhas e as fontes de propágulos também influência essas taxas da seguinte forma:

  • Ilhas próximas do continente tem maiores taxas de colonização de novas espécies porque mais propágulos podem chegar, eles não morrem na jornada por exemplo
  • O efeito da área seria mais importante em ilhar longe do continente do que ilhas perto do continente, para reverter o continuo processo de extinção

Derivar uma curva de imigração ou extinção é um bocado complicado. Mas para ilustrar a teoria, a gente pode assumir que a taxa de imigração I, pode ser representada como uma função exponencial negativa e^{I_0-iR}, onde I_0 é a taxa de imigração em uma ilha vazia e i é o efeito negativo por espécies, que multiplica R que é quantas espécies tem la.

Então a gente pode criar um exemplo, com o número de espécies na ilha:

1
exemplo<-data.frame(especies=seq(0.1,50,0.2))

E ver a taxa de acordo com o número de espécies:

1
2
3
4
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)

Agora para a extinção, E, está precisa ser zero se não houver espécies presentes. Imaginando que a extinção é uma função da densidade e que a densidade média declina conforme o número de espécies aumenta, então \bar{N}=\frac{1}{R}. então a gente vai ter algo mais ou menos assim

e^{dR}-1

Nos subtraimos esse 1 para ter certeiza que E=0 quando R=0.
O número de espécies, R, vai resultar de \Delta R = 0 = I-E, o ponto onde as linhas se cruzam.

I=e^{I_0-iR}
E=e^{dR}-1
\Delta R = 0 = I-E

03

O exato ponto onde essas linhas se cruzam podem ser encontrada de muitas formas, a gente tem visto isso em evolução aqui, mas vamos criar uma função no R para minimizar, nos vamos minimizar (I-E)^2, porque elevando ao quadrado a diferença nos temos uma quantidade com a conveniente propriedade que o mínimo ira ser aproximado tanto pelo lado positivo como o negativo.

1
2
3
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}

Nos usamos essa função de um parâmetro no optimize, e especificamos que queremos saber o ótimo em um intervalo de 1 a 50.

> optimize(f=deltaR,interval = c(1,50)) $minimum [1] 16.91321 $objective [1] 2.813237e-13

A saída nos diz que o mínimo é algo em torno de 16.9, e podemos colocar isso na figura

04

Agora suponha uma ilha um pouco mais longe, do mesmo tamanho, mas mais longe do continente, veja que a taxa de imigração diminuirá, porque pensa, as espécies vão morrer antes de alcançar a ilha e isso vai resultar um outro ponto mínimo, outro ponto de equilíbrio.

05

A beleza dess teoria é que ela foca nossa atenção em processos a nível de paisagem, geralmente fora dos limites espaciais e temporais das nossas capacidades de amostragem. Ele demostra que qualquer fator que ajude a determinar taxas de imigração ou extinção, incluindo a area de uma ilha ou a sua proximidade a uma fonte de propágulos, vai alterar o equilíbrio do número de espécies. Mas temos que lembrar que a identidade das espécies pode mudar ao longo do tempo, isso é, um turnover de espécies, mas um turnover devagar, porque as espécies que tem uma grande chance de imigrar e menor chance de extinguir vão se manter no ambiente muito provavelmente.

Bem é isso ai, incrível como as teorias mais legais são extremamente simples, mas com um efeito profundo na forma como pensamos sobre as coisas, 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.

Referência:
M.H.H. Stevens 2011, A Primer of Ecology with R

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
library(ggplot2)
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)
 
ggplot(data=exemplo,aes(x=especies,y=imigracao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de imigração")
ggsave("01.jpg")
##curve(exp(IO-b*x),0,50,xlab="n. de espécies",ylab="Taxa Imigração")
 
##Extinção
d<-0.01
exemplo$extincao<-exp(d*exemplo$especies)-1
ggplot(data=exemplo,aes(x=especies,y=extincao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de extinção")
ggsave("02.jpg")
##curve(exp(d*x)-1,0,50,xlab="n. de espécies",ylab="Taxa Extinção")
 
 
##Os dois juntos
exemplo2<-stack(exemplo[,c("imigracao","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,2)
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("03.jpg")
##
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
estavel<-optimize(f=deltaR,interval = c(1,50))
segmento<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,2:3],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("04.jpg")
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
exemplo$isp1<-exp(IO-0.1*exemplo$especies)
exemplo$isp2<-exp(IO-0.14*exemplo$especies)
exemplo$extincao<-exp(d*exemplo$especies)-1
exemplo2<-stack(exemplo[,c("isp1","isp2","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,3)
 
##
deltaR<-function(R,b){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
b<-0.1
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento1<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(2,4)],1,min))))
b<-0.14
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento2<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(3,4)],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração Sp1","Imigração Sp2"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento1,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento1[2,],colour="black")+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento2,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento2[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("05.jpg")

Relação espécie área

A relação entre o número de espécies encontradas em amostras de diferentes áreas é um padrão a muito tempo conhecido, relativamente simples de pensar no processo que gera esse tipo de padrão, já que áreas maiores cabem mais coisas, as vezes meio confuso, principalmente quando a gente realiza que áreas maiores podem ter tanto maior abundancia de espécies assim como mais espécies. Mas areas maiores possem mais de muitos tipo de recursos, o que diminui a chance de extinção, aumentam colonização e muitas possibilidade de manter mais espécies, mas não estamos falando de biogeografia de ilhas, ainda.

Certo, mas como isso funciona? Geralmente é um padrão empírico do número de espécies encontrados em manchas de diferentes tamanhos, plotado contra como uma função dos seus respectivos tamanhos, podendo essas manchas serem ilhas oceânicas ou uma mancha de vegetação numa matriz de campo aberto.

Quantitativamente, essa relação é representada como uma potencia (power law).

R=cA^z

onde o R é o número de espécies em uma mancha de tamanho A, e os c e z são constantes ajustaveis. Mas normalmente a gente ve essa lei plotada como uma relação de log-log, que a torna linear.

log(R)=b+zA

So aplica log dos dois lados.

Podemos gerar um plot aqui, para ter uma ideia melhor.

1
2
3
4
5
6
7
8
A<-10^(1:10)
c<-1.5
z<-0.25
 
fn<-function(x,c=1.5,z=0.25) {
    return(c*x^z)
}
ggplot(data.frame(x=c(10, 10^10)), aes(x)) + stat_function(fun=fn) + xlab(label="Área(ha)") + ylab(label="Número de espécies")

E isso gera a seguinte figura.

01

Essa é a cara da relação, como uma lei de potência, e podemos linearizar, transformando para log os dois eixos.

02

E a parte boa, é que linearizando, tudo fica mais fácil, em quesito de modelagem, porque agora é so passar uma reta. Para tanto, vamos dar uma olhada no artigo do Johnson 1968, na table I do paper.

> exemplo Location Area Species 1 Tiburon Peninsula 5.9 370 2 San Francisco 45.0 640 3 Santa Barbara area 110.0 680 4 Santa Monica Mountains 320.0 640 5 Marin County 529.0 1060 6 Santa Cruz Mountains 1386.0 1200 7 Monterey County 3324.0 1400 8 San Diego County 4260.0 1450 9 California Coast 24520.0 2525

Temos a área e o número de espécies encontradas, em manchas de vegetação. Vamos fazer um plot dos dados e dar uma olhada.

03

Parece que existe ai um padrão de espécie área bem claro, então a gente pode ajustar um power law nesses dados e extrair os valores constantes.

1
modelo<-lm(log(Species,10)~log(Area,10),data=exemplo)

E olhamos o resultado.

> summary(modelo) Call: lm(formula = log(Species, 10) ~ log(Area, 10), data = exemplo) Residuals: Min 1Q Median 3Q Max -0.129947 -0.013386 0.003153 0.041204 0.057274 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.38560 0.05568 42.85 9.84e-10 *** log(Area, 10) 0.21976 0.01917 11.46 8.64e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.06025 on 7 degrees of freedom Multiple R-squared: 0.9494, Adjusted R-squared: 0.9422 F-statistic: 131.4 on 1 and 7 DF, p-value: 8.642e-06

Agora aqui é importante para interpretar, a gente pensar em algumas coisa. A gente usou uma base 10 nos logs, isso porque é mais fácil de pensar, só pensar que log de 10 é 1, de 100 é 2, de 1000 é 3, melhor que aqueles números quebrados do log natural, só facilita visualizar as coisas a base do log aqui. Além disso o que interessa é a inclinação, que representa o nosso valor de potência, o quão rápido essa curva cresce.

Podemos ajustar como um modelo não linear também, para ver como fica o resultado.

1
modelo_nls<-nls(Species~a*Area^z,start=list(a=1,z=0.25),data=exemplo)
> summary(modelo_nls) Formula: Species ~ a * Area^z Parameters: Estimate Std. Error t value Pr(>|t|) a 195.7664 31.4382 6.227 0.000434 *** z 0.2494 0.0186 13.407 3.01e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 122 on 7 degrees of freedom Number of iterations to convergence: 4 Achieved convergence tolerance: 3.028e-06

E antes de ficar olhando os números, vamos ver como essas curvas ficam em relação aos dados.

04

O ajuste fica bem bom, parece predizer bem a quantidade de espécies encontradas, apesar de que os valores não são exatamente iguais, mas olhando os intervalos de confiança.

> confint(modelo) 2.5 % 97.5 % (Intercept) 2.2539475 2.5172501 log(Area, 10) 0.1744247 0.2650924 > confint(modelo_nls) Waiting for profiling to be done... 2.5% 97.5% a 128.9454047 283.572614 z 0.2055227 0.296902

A gente ve que ambos pegam muitos valores semelhantes, apesar do modelo linear ter fica com a média menor.
Bem um nome recorrente quando falamos de espécie área é o Preston, que propos que dado dado a distribuição de espécies log-normal, então áreas maiores devem acumular mais espécies em taxas específicas.

Bem é isso ai, o script vai estar la no repositório recologia, os gráficos com ggplot ainda são meia boca, mas é treinando que se aprende e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Johnson, Mason Raven 1968 – Ecological Parameters and Plant Species DiversityThe American Naturalist Vol. 102(926) pp. 297-306

M.H.H. Stevens 2011, A Primer of Ecology with R

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
library(ggplot2)
 
###Desenhando as curvas da relação espécie área
A<-10^(1:10)
c<-1.5
z<-0.25
 
fn<-function(x,c=1.5,z=0.25) {
    return(c*x^z)
}
ggplot(data.frame(x=c(10, 10^10)), aes(x)) + stat_function(fun=fn) + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("01.jpg")
 
 
fn_log<-function(x,c=1.5,z=0.25) {
    return(log(c,10)+z*x)
}
ggplot(data.frame(x=c(1, 10)), aes(x)) + stat_function(fun=fn_log) + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("02.jpg")
 
###Johnson, Mason, and Raven 1968
exemplo<-data.frame(Location = c("Tiburon Peninsula", "San Francisco","Santa Barbara area", "Santa Monica Mountains", "Marin County", 
                                 "Santa Cruz Mountains", "Monterey County", "San Diego County","California Coast"),
                    Area = c(5.9, 45, 110, 320, 529, 1386, 3324,4260, 24520),
                    Species = c(370L, 640L, 680L, 640L, 1060L, 1200L,1400L, 1450L, 2525L))
exemplo
 
 
ggplot(data=exemplo,aes(x=Area,y=Species))+ geom_point() + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("03.jpg")
 
log(1000,10)
 
modelo<-lm(log(Species,10)~log(Area,10),data=exemplo)
summary(modelo)
 
modelo_nls<-nls(Species~a*Area^z,start=list(a=1,z=0.25),data=exemplo)
summary(modelo_nls)
 
linha<-data.frame(x=log10(seq(1,25000,100)),y=log10(predict(modelo_nls,newdata=data.frame(Area=seq(1,25000,100)))))
 
 
 
ggplot(data=exemplo,aes(x=log10(Area),y=log10(Species),linetype="solid") )  + geom_point() +
    geom_abline(intercept = modelo$coefficients[1], slope = modelo$coefficients[2]) +
    geom_line(data=linha,aes(x=x,y=y,linetype="dashed"))+        
    xlab(label="Área(ha)") + ylab(label="Número de espécies")+
    scale_linetype_manual(name='Modelo',values =c('dashed','solid'), labels = c('Não Linear','Linear'))
ggsave("04.jpg")
 
confint(modelo)
confint(modelo_nls)

Phase Plane

Quando falamos de equilíbrio de duas populações aqui a gente citou a palavra chave State Space, que são as possibilidades que o sistema pode assumir, sejam elas possíveis ou não.

Sistema é um sistema oras, um conjunto de elementos interconectados de modo a formar um todo organizado.

Quando a gente falava de competição entre espécies, o todo é o ecossistema de duas espécies, os elementos são as espécies, a elas são interconectadas pelo efeito que uma tem sobre a outra, a competição.

Lembrando nosso sistema

\begin{cases}  {dN_{1}\over{dt}} = r_{1}\cdot N_{1}\cdot (1-{N_{1}\cdot \alpha_{11}} - {N_{2}\cdot \alpha_{12}}) \\  {dN_{2}\over{dt}} = r_{2}\cdot N_{2}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})  \end{cases}

Que era mais fácil de visualizar com essa figura.

03

Beleza, comumente a gente resolvia esse sistema ao longo do tempo, para ver o que acontecia com as espécies, pra ver o que é possível e o que não é possível.

01

Nesse caso, elas iam crescendo, ai quando a população aumentava, elas começavam competir, as abundâncias iam mudando mas chegava uma hora que tudo estacionava num ponto, um ponto de equilíbrio, e la ficava a não ser que houve-se uma pertubação.

Mas ao invés de olhar dessa forma, podemos fazer um gráfico da abundância da espécie 1 em relação a abundância da espécie 2.

02

E esse que é o tal do Phase Plane, aqui em duas dimensões, porque temos duas espécies, mas poderia ser em três, n dimensões, claro que para mais de três dimensões não da para fazer exatamente um gráfico, mas da para tentar visualizar ele também.

Mais que fluflu para fazer charme, esse tipo de figura nos ajuda a determinar como funciona um sistema, se ele tem algum tipo de equilíbrio ou não, ciclo, essas coisas. A gente viu aqui os casos de equilíbrio para duas populações. Mas sempre tudo ia para algum ponto, meio diretão assim, Mas existem outros equilíbrios bem legais também.

Se a gente olhar o help da função ode a gente conseguir ver que a autora mandou um exemplo pronto do famoso sistema de predador e presa com uma presa com crescimento logístico.

\begin{cases}  {Ingestion} = rIng * Prey * Predator \\  {GrowthPrey} = rGrow * Prey * (1-\frac{Prey}{K}) \\  {MortPredator} = rMort * Predator \\  {dPrey} = GrowthPrey - Ingestion \\  {dPredator} = Ingestion * assEff - MortPredator \\  \end{cases}

O sistema é esse conjunto de equações, que fica mais simples talvez escrita na função do R para usar com o ode

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
 
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
 
    return(list(c(dPrey, dPredator)))
  })
}

Ao resolver esse sistemas a gente tem o seguinte comportamento de predadores e presas ao longo do tempo.

03

Agora se a gente fizer um Phase Plane, para olhar esse sistema (lembrando que definimos um certo set de parâmetros, como vimos para populações, dependendo dos valores dos parâmetros, os alphas la, víamos coisas diferentes).

04

Uma espiral, o início é onde tudo está oscilando mais, mas com o tempo as oscilações entre abundâncias de presas e predadores vai diminuindo até estabilizar em um ponto fixo.

Agora vamos olhar ainda mais um pouquinho o sistema de predadores e presas, mas sem a presa com crescimento logístico, vamos pegar o modelo que aparece no wikipedia nessa página aqui.

O modelo é esse aqui.

\begin{cases}  {\frac{dx}{dt}} = x(\alpha - \beta y) \\  {\frac{dy}{dt}} = - y (\gamma - \delta x) \\  \end{cases}

Que podemos escrever no R assim:

predadorpresa <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        dpresa <- presa*(alpha-beta*predador)
        dpredador <- -predador*(gamma-delta*presa)
        list(c(dpresa, dpredador))
  } )
}

Claro que cada parâmetro tem um significado especial, mas não vamos ficar discutindo isso agora, a gente só quer ver os phase planes legais.
Mas então resolvendo esse sistema para um set de parâmetros e valores iniciais de predador e presa vemos o seguinte phase plane.

05

Uma bolinha, nesse caso aqui, a gente não tem um ponto de estabilidade, mas a gente tem um tipo de ciclo, esse trilho ai e o sistema nunca vai parar, ele vai ficar oscilando entre vários estados, mas os estados vão estar nessa figura ae, nesse circulo. Mas ainda diferente dos sistemas anteriores onde víamos ele “andando certinho”, a gente pode ligar linhas de um tempo para o próximo e ver algo legal.

06

A gente não esta indo suavemente de um ponto para o próximo, a gente ta indo de um canto do phase plane para o outro la longe. Melhor que isso, resolvemos o sistema para 20 tempos, então podemos numerar as bolinhas na ordem, para ver melhor onde tudo esta começando e acabando.

07

Bem mais bagunçado que a competição entre duas espécies não? Mas ainda sim existe um padrão ai, e o equilíbrio é mais algo como um ciclo que um ponto.

Veja só como a população esta oscilando ao longo do tempo, para gerar esse Phase Plane

08

Muito legal, agora a gente sabe que esse tipo de gráfico tem um nome, e pode ajudar, visualmente a entender melhor como as coisas funcionam, e na biologia a gente consegue descrever uma boa quantidade de coisas com sistemas de equações assim. Claro que é mais complicado que isso, entender, a partir formulações matemática de sistemas, como eles funcionam, o que podemos esperar e o que não devemos esperar, mas esse tipo de gráfico é bem legal, para visualizar algumas coisas que são talvez bem complicadas de entender somente olhando as equações.

Bem, acho que é isso ai por enquanto. O script está no repositório do github, além de aqui em baixo, que da pra olhar os parâmetros usados, que não citei la em cima, ainda fica faltando estudar as tal das matrizes jacobianas com mais detalhe, que é uma forma de achar pontos de equilíbrio, eu acho que serve pra isso pelo menos, na verdade falta estudar tanta coisa que nem sei por onde começar as vezes.

Referência:
M Henry H Stevens 2011 A primer of ecology with R. Springer

library(deSolve)
 
#Modelo para crescimento continuo
lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}
 
parms <- c(r1=0.8,r2=0.5,a11=0.010,a21=0.005,a22=0.010,a12=0.005);
initialN<-c(1,1)
out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
 
#Espécies ao longo do tempo
#figura 1
matplot(out[,1],out[,-1],type='l',xlab="Tempo",ylab="Tamanho populacional",
        frame.plot=F,main="Crescimento logístico contínuo para duas espécies")
legend("right",c(expression("Espécie 1 "*(alpha[21]==0.005)),
                expression("Espécie 2 "*(alpha[12]==0.005))),
       lty=1:2,bty='n',col=c("black","red"))
 
#Representação das abundancias das duas espécies num plano
#figura 2
plot(1, 1, type = "n", ylim = c(0, 100), xlim = c(0,100),frame=F,
     ylab = expression("N"[2]),xlab=expression("N"[1]))
points(out[,2],out[,3],type="b",cex=0.5,pch=19)
 
 
##########################################
#  #
##########################################
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
 
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
 
    return(list(c(dPrey, dPredator)))
  })
}
 
pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity
 
yini  <- c(Prey = 1, Predator = 2)
times <- seq(0, 200, by = 1)
out   <- ode(yini, times, LVmod, pars)
 
## User specified plotting
#figura 3
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tempo", ylab = "Abundância",
        main = "Lotka-Volterra", lwd = 2,frame=F)
legend("topright", c("Presa", "Predador"), col = 1:2, lty = 1:2,bty="n")
 
#figura 4
plot(out[,2],out[,3],type="b",pch=19,frame=F,xlab="Presa",ylab="Predador")
 
 
#########################################
 
predadorpresa <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        dpresa <- presa*(alpha-beta*predador)
        dpredador <- -predador*(gamma-delta*presa)
        list(c(dpresa, dpredador))
  } )
}
 
pars <- c(alpha=2,beta=1.1,gamma=1,delta=0.4)
yini<-c(presa=5,predador=1)
times <- seq(0, 20, by = 1)
out   <- ode(yini, times, predadorpresa, pars)
 
#figura 5
plot(out[,2],out[,3],type="p",pch=19,xlab="Presa",ylab="Predador")
 
#figura 6
plot(out[,2],out[,3],type="b",pch=19,xlab="Presa",ylab="Predador")
 
#figura 7
plot(out[,2],out[,3],type="l",pch=19,lty=3,lwd=0.5,xlab="Presa",ylab="Predador")
text(out[,2],out[,3],0:20)
 
#figura 8
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tempo", ylab = "Abundância",
        main = "Predador-Presa", lwd = 2,frame=F)
legend("topright", c("Presa", "Predador"), col = 1:2, lty = 1:2,bty="n")

O elefante de John von Neumann

A gente falou no post anterior, a muito tempo atrás, sobre como pode ser ruim usar muitos parâmetros em um modelo. Ja que muitos parâmetros implicam em ajustar quaisquer dados a um modelo.

Mas aqui no blog a gente sempre fala de vários modelos teóricos em ecologia, modelos que nos fazem ter previsões, que nos fazem achar alguma coisa antes mesmo de começar a coletar dados.

Mas no fim das contas, podemos dizer que ciência é feita de cima pra baixo ou de baixo pra cima, na maioria das vezes.

Como assim?

As vezes a gente tem um monte de dados, e usa estes dados para chutar como o mundo funciona.
Outras vezes a gente chuta como o mundo funciona, depois tenta coletar dados específicos para ver se esse chutes fazem sentido.

Esses são os Top-Down ou Bottom-Up da ciência. Sendo que em ecologia, talvez, a estratégia Top-Down é bem comumente vista, enquanto em, sei la, física, Bottom-Up é uma estratégia mais comum.

E ambas as abordagem podem sofrer desse mal.

Como o Einstein disse:

Everything Should Be Made as Simple as Possible, But Not Simpler

Tudo precisa ser o mais simples possível, mas não mais simples que o possível

Isso deixa uma dúvida eterna, se o modelo que estamos ajustando não tem parâmetros demais, ou parâmetros de menos, e chega uma hora que não tem nada, além do bom senso, intuição do pesquisador, para salvar o dia.

Mas principalmente no Top-Down, que a gente discute bastante aqui, sempre ficamos nessa discussão sobre over-fit, sobre problemas com muitos parâmetros, parcimônia e tal. Talvez no exemplo do post anterior não seja muito elucidativo, um exemplo melhor seria a historinha que o Freeman Tyson contou nesse pequeno artigo na nature aqui.

Freeman Tyson é um físico e matemático ingles, e aqui nesse artigo, ele conta quando foi visitar o Enrico Fermi, que era um pica-grossa na física. Todo mundo quer falar com os pica-grossa, ou bigwig, até porque eles aumentam o fator de impacto dos seus artigos, ou sua chance de publicar. Pelo menos em ecologia isso parece ser um fato, nesse artigo aqui os cara verificaram esse efeito.

Ele diz que uma grande reviravolta na vida dele foi devido a um encontro com o Enrico Fermi em 1953. Em poucos minutos, Fermi de uma maneira educada mas sem piedade destruiu um programa de pesquisa que ele e seus estudantes vinham trabalhando duro a vários anos. Fermi provavelmente salvou todo mundo de muitos mais anos de pesquisa infrutífera na verdade ao destruir o projeto atual deles dizendo apenas a verdade do que achava. Mas a vida é complicada as vezes, como na fabula do passarinho, da vaca e do gato, a gente so entendendo o que está acontecendo depois.

Então, o Tyson era um jovem professor de física teórica na Universidade de Cornell, responsável por dirigir um pequeno exercito de alunos de pos-graduação e pos-docs. Como manda quem pode, obedece quem tem juízo, Tyson, colocou todo mundo para calcular a dispersão de meson-protons, eu nem sei o que é isso, mas o Fermi era um cara que tinha medidas dessa dispersão, e o Tyson queria fazer o modelo teórico que explicasse aqueles dados, por isso o Tyson foi la mostrar ao Fermi os cálculos que estava desenvolvendo.

Em 1948 e 1949, Tyson tinha feito cálculos similares de alguns processos atômicos e tinha encontrado um ajuste muito bom dos seus modelos quanto a medidas retiradas em outros experimentos, o que deixou ele confiante no que estava fazendo.

Mas continuando, Tyson pediu para conversar com o Fermi, na época devia ser via carta, sei la, e um dia foi la na universidade mostrar o trabalho que ele e seus alunos estavam fazendo. Quando ele chegou na sala do Fermi, ele entregou seu trabalho, vários gráficos, e o Fermi mal olhou para os gráficos e convidou o Tyson para sentar. Perguntou da família, como estão as coisas e tal. Depois de trocar uma ideia, ele soltou seu veredito.

There are two ways of doing calculations in theoretical physics, one way, and this is the way I prefer, is to have a clear physical picture of the process that you are calculating. The other way is to have a precise and self-consistent mathematical formalism. You have neither.

Existem duas maneiras de fazer cálculos em física teórica, um deles, e esse é o que eu prefiro, é ter uma visão clara do processo físico que você está fazendo cálculos. O outro modo é ter um formalismo matemático preciso e auto-consistente. Você não tem nenhum dos dois

Qual a primeira reação?
O Tyson ficou puto com certeza, quem gosta de receber críticas, quem gosta de estar errado? Uma resposta dessa é um tapa na cara, mas ele ainda se aventurou a perguntar porque ele não considerava o modelo do Tyson e seus alunos bom, consistente matematicamente. E ele respondeu…

Quantum electrodynamics is a good theory because the forces are weak, and when the formalism is ambiguous we have a clear physical picture to guide us. With the pseudoscalar meson theory there is no physical picture, and the forces are so strong that nothing converges. To reach your calculated results, you had to introduce arbitrary cut-off procedures that are not based either on solid physics or on solid mathematics.

A teoria eletrodinâmica quântica é boa porque as forças são fracas, e quando o formalismo é ambíguo, nos temos uma imagem clara do processo físico para nos guiar. Com a teoria do meson pseudoescalar, não existe imagem física, e as forças são fortes de modo que nada converge. Para alcançar seus resultados, você teve que incluir uma serie de procedimento de corte arbitrários que não são baseados em conceitos sólidos de física nem de matemática

Em desespero, o Tyson ainda perguntou ao Fermi, se ele não tinha ficado impressionado com o excelente ajuste entre o modelo do Tyson e os dados coletados do Fermi. E o Fermi respondeu:

“How many arbitrary parameters did you use for your calculations?” I thought for a moment about our cut-off procedures and said, “Four.” He said, “I remember my friend Johnny von Neumann used to say, with four parameters I can fit an elephant, and with five I can make him wiggle his trunk.”

“Quantos parâmetros arbitrários você usou nos seus cálculos? Eu pensei por um momento sobre nossos procedimentos de corte e disse,”Quatro”. Ele disse, “Eu lembro que meu amigo Jhonny Von Neumann costumava dizer, com quatro parâmetros eu posso ajustar um elefante, e com cinco eu faço ele mexer a tromba”

E ali a conversa terminou, o Tyson agradeceu o Fermi por telo recebido, e pegou o próximo busão de volta a Ithaca para dizer aos seus alunos sobre as mas noticias. Como a gente é obrigado a publicar, eles ainda não abandonaram aquele projeto, e acabaram por publicar na revista Physical Review com o nome de todo mundo, e depois disso todo mundo dispersou para outras linhas de pesquisa, o Tyson foi para Berkley, Califórnia, para começar uma nova carreira em Física de matéria condensada (devo estar traduzindo errado isso).

Mas olhando para traz depois de 50 anos, da pra ver que o Fermi estava certo todo o tempo, mesmo sem saber. A descoberta crucial que colocou sentido nas forças forte foi o quark. Mesons e prótons são pequenos sacos de quarks. Antes de Murray Gell-Mann descobrir os quarks, nenhuma teoria de forças fortes poderia possivelmente descrever o sistema adequadamente. Fermi não sabia nada de quarks, e morreu muito antes dos quarks serem descobertos. Mas de alguma forma ele sabia que alguma coisa essencial estava faltando a teoria dos meson em 1950. Sua intuição física disse a ele que a teoria pseudoescalar do meson não podia estar correta. E foi a intuição do Fermi, e não nenhuma discrepância entre modelo e dados que salvou o Tyson e seus alunos de ficarem presos a uma teoria sem futuro.

Bem, e pra variar, o cara não estava zuando, com quatro parâmetros, realmente da para ajustar um elefante aos dados.

elefantinho

Agora pense, se com 4 parâmetros da para ajustar um elefante, aquelas regressões com 10, 15 parâmetros, ajusta qualquer nuvem de pontos, sempre, agora será que por acaso ou porque os parâmetros realmente tem sentido? Essa é sempre uma boa pergunta!
Bem essa historia é algo para refletir, nem sempre quem encontra defeitos no seu trabalho é seu inimigo. Mas é isso ai, ficamos sem post por um tempo porque eu estava meio ocupado, mas agora espero voltar a frequência normal de posts.

#John von Neumann disse:
#"With four parameters I can fit an elephant,
#and with five I can make him wiggle his trunk".
 
#Quatro parâmetros, eles são número complexos, podemos escreve-los no R assim:
 
p1 <- complex(real = 50, imaginary = -30)
 
#Ou de forma simplificada:
p1<-  50-30i
p2<-  18+8i
p3<-  12-10i
p4<- -14-60i
 
#O modelo vai ser representado por uma série de fourier, que de acordo com o
#wikipedia é:
#Em matemática, uma série de Fourier, nomeada em honra de
#Jean-Baptiste Joseph Fourier (1768-1830), é a representação de uma
#função periódica (muitas vezes, nos casos mais simples, tidas como tendo
#período 2pi) como uma soma de funções periódicas da forma.
 
 
#Então definimos uma função pra calcular a série de Fourier:
#Veja que ela começa definindo um vetor pra descarregar os resultados.
#Em seguida extrai a parte real e a imaginaria do parametro em A e B
#o R faz coerção para numeros reais na maioria das operações o que faz tudo dar
#pau. No loop ali dentro da função ele faz a formula que esta no artigo.
 
fourier<-function(tt,CC){
                         f<-numeric(length(tt));
                         A<-Re(CC);
                         B<-Im(CC);
                         k<-1:length(CC);
                         for(i in 1:length(tt)){
                                                co<-cos(tt[i]*k)
                                                si<-sin(tt[i]*k)
                                                f[i]<-sum(A*co+B*si)
                                                f
                                               }
                        f
                        }
 
 
#Agora que temos uma função para calcular a serie de fourier, fazemos a função
#que é o modelo do nosso elefantinho, nossas observações e 4 parametros
 
 
elefante <- function(tt,p1,p2,p3,p4) {
                                      n<-6
                                      Cx<-complex(n);
                                      Cy<-complex(n);
                                      Cx[1]<-Re(p1)*1i
                                      Cx[2]<-Re(p2)*1i
                                      Cx[3]<-Re(p3)
                                      Cx[5]<-Re(p4)
                                      Cy[1]<-Im(p4)+Im(p1)*1i
                                      Cy[2]<-Im(p2)*1i
                                      Cy[3]<-Im(p3)*1i
 
                                      x<-fourier(tt,Cx)
                                      y<-fourier(tt,Cy)
                                      res<-cbind(y,-x)
                                      res
                                      }
 
 
#Pegamos uma sequencia de observações para introduzir no nosso
#modelo deterministico
observações<-seq(from=0,to=2*pi,length=100)
 
#E vemo o que acontece
elefantinho<-elefante(observações,p1,p2,p3,p4)
 
plot(elefantinho,type="l",col="red",ylab="",xlab="",lwd=2)
 
#Com 5 parametros da pra adicionar um olhinho :)
p5<-40+20i
points(Im(p5),Im(p5),pch=20,col="red")
 
#Artigo original:
#“Drawing an elephant with four complex parameters”
#by Jurgen Mayer, Khaled Khairy, and Jonathon Howard,
#Am. J. Phys. 78, 648 (2010), DOI:10.1119/1.3254017.
#http://java-srv1.mpi-cbg.de/publications/getDocument.html?id=ff8080812daff75c012dc1b7bc10000c
 
 
#Codigo em python que usei de base:
#"""
#Author: Piotr A. Zolnierczuk (zolnierczukp at ornl dot gov)
#
#Based on a paper by:
#Drawing an elephant with four complex parameters
#Jurgen Mayer, Khaled Khairy, and Jonathon Howard,
#Am. J. Phys. 78, 648 (2010), DOI:10.1119/1.3254017
#"""
#import numpy as np
#import pylab
#
## elephant parameters
#p1, p2, p3, p4 = (50 - 30j, 18 +  8j, 12 - 10j, -14 - 60j )
#p5 = 40 + 20j # eyepiece
#
#def fourier(t, C):
#    f = np.zeros(t.shape)
#    A, B = C.real, C.imag
#    for k in range(len(C)):
#        f = f + A[k]*np.cos(k*t) + B[k]*np.sin(k*t)
#    return f
#
#def elephant(t, p1, p2, p3, p4, p5):
#    npar = 6
#    Cx = np.zeros((npar,), dtype='complex')
#    Cy = np.zeros((npar,), dtype='complex')
#
#    Cx[1] = p1.real*1j
#    Cx[2] = p2.real*1j
#    Cx[3] = p3.real
#    Cx[5] = p4.real
#
#    Cy[1] = p4.imag + p1.imag*1j
#    Cy[2] = p2.imag*1j
#    Cy[3] = p3.imag*1j
#
#    x = np.append(fourier(t,Cx), [-p5.imag])
#    y = np.append(fourier(t,Cy), [p5.imag])
#
#    return x,y
#
#x, y = elephant(np.linspace(0,2*np.pi,1000), p1, p2, p3, p4, p5)
#pylab.plot(y,-x,'.')
#pylab.show()