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()

Epidemiologia – O Modelo SIR com transmissão dependente da frequência

No post anterior aqui, a gente deu uma olhada no modelo SIR, para avaliar a dinâmica de uma doença parasitaria, sendo que então é mais fácil contabilizar indivíduos infectados ou não, que a população do parasita per se.

O modelo é o seguinte.

\frac{dS}{dt} = -\beta I S
\frac{dI}{dt} = \beta I S - \gamma I
\frac{dR}{dt} = \gamma I

So para facilitar o entendimento, revisando com uma figurinha.

01

São os três compartimentos, S, I e R, e veja que temos as formulas de como é o fluxo de individuos entre esses compartimentos, veja que no grafo ela aparece uma vez, isso porque o \beta I S tira indivíduos do S e põe para o I, por isso que nas equações la em cima, ele é com um sinal de negativo na derivada do S e positivo na derivada do I.
No caso do I, ele perde para individuos para o R pelo \gamma I, ai mesmo esquema, menos no I e mais no R.

Essa conclusão se aplica quando S,I e R são densidades. Isso quer dizer que se a gente aumentar o tamanho da população, mas também a área associada a essa população, então você não aumentou a densidade. Ja se apenas o tamanho da população aumenta, mas a densidade é mantida, então a frequência de interações não aumenta, mas algumas populações podem aumentar em densidade com o aumento da população, se não temos mais área para expandir por exemplo.

Além da transmissão por ação de massa (dependente da densidade), existem outras formas de modelar o problema de forma dependente do tamanho da população, da frequencia desses compartimentos dentro do total da população. Uma das formas mais simples é conhecida como transmissão dependente da frequência (“frequency-dependente transmission”), onde:

\frac{dS}{dt} = -\beta \frac{I S}{N}
\frac{dI}{dt} = \beta  \frac{I S}{N} - \gamma I
\frac{dR}{dt} = \gamma I

A forma correta da função de transmissão depende de como ela ocorre. Imagine duas pessoas em um elevador, uma doente (Infectada) e uma sadia (Suscetível), então o indivíduo infectado espirra, isso resulta em uma uma probabilidade particular \beta que o indivíduo suscetível fique infectado. Agora imagine a mesma situação, mas um indivíduo resistente entre no elevador, será que a adição de indivíduos resistentes modifica a probabilidade que indivíduos suscetíveis fiquem infectados? Note o que mudou e o que não mudou. Primeiro, com a adição de indivíduos resistentes, N aumentou dentro do elevador, foi de 2 para 3, e a prevalência diminuiu, I/N foi de 1/2 para 1/3, no entanto as densidade de I e S continuam a mesma (1 de cada no elevador). O que pode acontecer?

1. Se uma quantidade suficiente de parasitas(imagine a doença que quiser) se espalha igualmente pelo elevador, adicionar um indivíduo resistente não muda a chance do indivíduo suscetível ficar contaminado, e a chance de infeção vai continuar dependente das densidade de I e S, a taxa não vai variar com o declínio da prevalência nesse caso.

2. Se apenas o vizinho imediato do indivíduo infectado é exposto ao patógeno, pense que ele precisa passar pelo toque, pelo contato físico, um esbarrão, então a probabilidade de transmissão diminui com o aumento de indivíduos Resistentes (R) certo? E assim a taxa de transmissão vai diminuir com o declínio da prevalência, porque agora além da chance de \beta, temos também que o indivíduo infectado tem que estar ao lado indivíduo suscetivel, se o resistente fica no meio do caminho, a doença não se espalha, ou seja o N fez o \beta diminuir, e veja la nas formulas acima que o N está dividindo tudo que o \beta multiplica, logo quanto maior o N, menor o \beta.

É simples imaginar alguns cenários diferentes que se enquadram nesses dois resultados, e é muito importante justificar a forma da função.

Transmissão dependente da densidade é independente do número de indivíduos resistentes, ter maior densidade de indivíduos infectados ou mais indivíduos suscetíveis sempre aumenta a taxa de transmissão, assumindo que ambos I, S > 0. Vejamos como isso funciona numa figura.

02

Aqui é relativamente simples, direto, a formula é \beta I S para o fluxo, logo se a gente aumenta a quantidade de I ou de S (não esqueça do I, S > 0 para não zerar a transmissão), a quantidade transmitida aumenta linearmente, ou o que nos interessa, quantos indivíduos ficam doentes.

Porém, em contraste, a transmissão dependente da frequência de indivíduos, como quando os resistentes ficam no meio da transmissão. Isso é, eles reduzem a probabilidade de infectados entrarem em contato com suscetíveis, similarmente, ter mais indivíduos suscetíveis torna mais provável que o vizinho imediato de um indivíduo suscetível seja outro individuo suscetível, o que não gera transmissão também, lembre-se que o N está dividindo, \beta \frac{I S}{N}.

Isso fica assim então

03

Veja que nesse exemplo, temos zero indivíduos Resistentes, e estamos vendo o que acontece no aumento de indivíduos Suscetíveis e Infectados, a transmissão é dependente da prevalência, do tamanho total da população, então por exemplo adicionar muitos indivíduos suscetíveis não aumenta linearmente a transmissão no caso anterior, porque pensamos que eles podem ficar um do lado dos outros e isso não gera transmissão, o que gera transmissão é um infectado junto a um individuo suscetível, mas isso depende do tamanho total da população, então crescemos a transmissão ao adicionar indivíduos infectados e Suscetíveis, veja que a superfície é curvilinear, precisa dos dois para subir.

Mas o que a transmissão dependente de frequência implica sobre a dinâmica? Vamos revolver para dI/dt >0

0 < \beta \frac{S I}{N} -\gamma I[/latex] e [latex]\gamma < \beta \frac{S}{N}[/latex]  Como nos fizemos no <a href="http://recologia.com.br/2015/06/epidemiologia-o-modelo-sir/">post</a> anterior, vamos considerar a situação anterior a explosão da epidemia onde [latex]S \approx N, de forma que \frac{S}{N} = 1. Nesse caso, a taxa básica de reprodução é R_0 = \frac{\beta}{\gamma} que é independente de N. Uma epidemia vai ocorrer se \beta > \gamma, independente da densidade da população, interprete isso como o seguinte biologicamente falando, se todo mundo se cura muito rápido, não existe tempo suficiente para os infectados infectarem muita gente.

Isso está em contraste direto com o modelo de transmissão dependente da densidade, onde a epidemia pode ser prevenida simplesmente reduzindo o tamanho populacional N, para uma densidade baixa o suficiente. Ambos os modelos de transmissão são observados em populações humanas e não humanas, então é importante conhecer como a doença se espalha para conseguir predizer sua dinâmica.

Outro fenômeno interessante com a transmissão dependente de frequência é que a prevalência (I/N) pode declinar com o aumento da tamanho populacional. Dois processos contribuem com esse padrão. Primeiro que a epidemia em uma população totalmente suscetível começa com um único indivíduo infectado, e como consequência, a prevalência inicial é sempre I/N = 1/N. Segundo, e como uma consequência disso, a transmissão é menor em grandes populações já que \beta \frac{S I}{N} é pequeno. Como consequência a prevalência se mantêm baixa por relativamente um longo período.

Veja que quando comparamos a prevalência contra a densidade populacional ao longo do tempo, onde os tipos de linhas são os tamanhos populacionais diferentes, para o modelo dependente da frequência (o tamanho populacional importa):

04

Um maior tamanho populacional maior aqui implica numa menor prevalência geral.

Enquanto para o caso dependente da densidade (aqui o que importa a proporção de infectador para suscetíveis, ou o contrario, o que interessa são essas quantidade apenas).

05

A prevalência aumenta com a densidade.

Mais uma questão interessante ainda, em uma população sazonal, muitos indivíduos em grandes populações se mantêm não infectados depois de um certo tempo, e dependendo do organismo, isso pode ser tempo suficiente para se reproduzir, ai entramos na pressão evolutiva para resistir a doença, ou não pressão, porque depois que você se reproduziu, para uma plantinha sazonal por exemplo, dane-se, é só secar e morrer porque sua parte você já fez, agora é com os descendentes.

07

Em contraste, o modelo dependente da densidade tipicamente mostra um padrão contrário, com uma epidemia mais rápida e extrema explosão de infectados, e maior prevalência em populações densas.

06

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:
Stevens, M. H. – 2011 – A Primer of Ecology with R. Springer

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
##################################
## Grafo
##################################
library(igraph)
 
SIR_grafo<-graph.formula(S-+I,I-+R)
E(SIR_grafo)$label<-c(expression(paste(beta,"IS")),expression(paste(gamma,"I")))
loc<-matrix(c(0,1,1,1,2,1),ncol=2,byrow=T)
 
#jpeg("01.jpg")
plot(SIR_grafo,layout=loc,rescale=F,xlim=c(0,2.2),ylim=c(0,2))
#dev.off()
 
##################################
## Modelos
##################################
library(deSolve)
 
SIR <- function(t, y, p) {
    {
    S <- y[1]
    I <- y[2]
    R <- y[3]
    }
    with(as.list(p), {
                    dS.dt <- -B * I * S
                    dI.dt <- B * I * S - g * I
                    dR.dt <- g * I
                    return(list(c(dS.dt, dI.dt, dR.dt)))
                })
}
 
SIRf <- function(t, y, p) {
    {
    S <- y[1]
    I <- y[2]
    R <- y[3]
    N <- S + I + R
    }
    with(as.list(p), {
                    dS.dt <- -B * (I * S)/N
                    dI.dt <- B * (I * S)/N - g * I
                    dR.dt <- g * I
                    return(list(c(dS.dt, dI.dt, dR.dt)))
                })
}
 
R <- 0
S <- I <- 1000
Ss <- Is <- seq(1,S,length=11)
N <- S + I + R
betaD <- 0.1
betaF <- betaD * N
 
mat1 <- sapply(Is, function(i) {betaD*i*Ss})
mat2 <- sapply(Is, function(i) {betaF*i*Ss/(i+Ss+R)})
 
#jpeg("02.jpg")
persp(mat1,theta=20,phi=15,r=10,zlim=c(0,betaD*S*I),xlab="I",ylab="S",
      main="Dependente da densidade",zlab="Taxa de Transmissão")
#dev.off()
 
#jpeg("03.jpg")
persp(mat2,theta=20,phi=15,r=10,zlim=c(0,betaF*S*I/N),xlab="I",ylab="S",
      main="Dependente da frequência",zlab="Taxa de Transmissão")
#dev.off()
##################################
## Taxa Transmissão
##################################
 
S <- 4^(0:4)
I <- 1
parmsf <- c(B=1, g=0)
parmsd <- c(B=1/16, g=0)
 
Months <- seq(0,8,by=0.1)
outd <- sapply(S, function(s) {
                   out <- ode(c(s,I,R),Months,SIR,parmsd)
                   out[,3]/apply(out[,2:4],1,sum)
               })
 
outf <- sapply(S, function(s) {
                   out <- ode(c(s,I,R), Months, SIRf, parmsf)
                   out[,3]/apply(out[,2:4],1,sum)
               })
 
#jpeg("04.jpg")
matplot(Months,outd,type="l",col=1,ylab = "Prevalência (I/N)",
        main="Dependente da densidade")
legend("bottomright",legend=S,lty=1:length(S),bty="n")
#dev.off()
 
#jpeg("05.jpg")
matplot(Months,outf,type="l",col=1,ylab = "Prevalência (I/N)",
        main="Dependente da frequência")
legend("bottomright",legend=S,lty=1:length(S),bty="n")
#dev.off()
##################################
## Comparação da dinamica
##################################
N <- 1000
I <- 1
R <- 1
S <- N - I - R
 
 
#
#jpeg("06.jpg")
meses <- seq(0, 0.8, by = 0.01)
parms <- c(B = 0.1, g = 4)
SIR.out <- data.frame(ode(c(S, I, R), meses, SIR, parms))
colnames(SIR.out)<-c("tempo","S","I","R")
matplot(meses, SIR.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4,main="Modelo dependente da densidade")
legend("right", c("Resistentes", "Infectados", "Susceptiveis"), lty = 3:1, col = c(4,3,2), bty = "n",lwd=2)
#dev.off()
 
#jpeg("07.jpg")
layout(matrix(1:4,ncol=2,nrow=2,byrow=T))
#
N <- 10000
I <- 1
R <- 1
S <- N - I - R
SIRf.out <- data.frame(ode(c(S, I, R), meses, SIRf, parms))
colnames(SIRf.out)<-c("tempo","S","I","R")
matplot(meses, SIRf.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4,main="N=10000")
#
N <- 1000
I <- 1
R <- 1
S <- N - I - R
SIRf.out <- data.frame(ode(c(S, I, R), meses, SIRf, parms))
colnames(SIRf.out)<-c("tempo","S","I","R")
matplot(meses, SIRf.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4,main="N=1000")
#
N <- 100
I <- 1
R <- 1
S <- N - I - R
SIRf.out <- data.frame(ode(c(S, I, R), meses, SIRf, parms))
colnames(SIRf.out)<-c("tempo","S","I","R")
matplot(meses, SIRf.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4,main="N=100")
#
N <- 10
I <- 1
R <- 1
S <- N - I - R
SIRf.out <- data.frame(ode(c(S, I, R), meses, SIRf, parms))
colnames(SIRf.out)<-c("tempo","S","I","R")
matplot(meses, SIRf.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4,main="N=10")
#dev.off()

Epidemiologia – O modelo SIR

Epidemiologia é a ciência que estuda padrões, causas e efeitos de doenças em populações definidas. Tipicamente, patógenos que causam doenças, como fungos, bactérias e vírus, com alguma especificidade de hospedeiro, e que tem um crescimento populacional dentro do hospedeiro.

O modelo mais simples para doenças, não modela as populações dos patógenos diretamente, mas contabiliza os diferentes estados dos hospedeiros dentro de sua população, primariamente aqueles com e sem sintomas. Desenvolvido por Kermack e McCormick, esse modelo inicial é conhecido como SIR, e basicamente ele registra ao longo do tempo três estados na população de hospedeiros, que são:

Suscetíveis (Susceptible) – Indivíduos que não estão infectados, mas podem vir a serem infectados.
Infectados (Infected) – Indivíduos que estão infectados
Resistentes (Resistant) – Indivíduos que são resistentes a doença, tipicamente aqueles que ja conseguiram adquirir uma imunidade devido a uma infecção passada.

Assim o total da população N é igual a S+I+R. Lembrando que N,S,I e R são densidades. Isto é, nos avaliamos número de indivíduos em uma área fixa. Isso é importante, ja tem um efeito direto sobre a dispersão, ou transmissão da doença.

As doenças se espalham dos indivíduos infectados para os suscetíveis. A taxa depende de certa forma ao número, ou alternativamente, na fração da população que está infectada. Indivíduos resistentes não são tipicamente considerados vetores do patógeno, e dessa forma, conforme vão aumentando, estes diminuem a taxa de transmissão ao diluir os dois outros grupos.

O ponto de partida para entender tudo isso, é uma população com tamanho constante, sem imigração ou emigração.

\frac{dS}{dt} = -\beta I S
\frac{dI}{dt} = \beta I S - \gamma I
\frac{dR}{dt} = \gamma I

No R, podemos criar uma função para descrever esse conjunto de equações.

1
2
3
4
5
6
7
8
9
10
11
12
13
SIR <- function(t, y, p) {
    {
    S <- y[1]
    I <- y[2]
    R <- y[3]
    }
    with(as.list(p), {
                    dS.dt <- -B * I * S
                    dI.dt <- B * I * S - g * I
                    dR.dt <- g * I
                    return(list(c(dS.dt, dI.dt, dR.dt)))
                })
}

Ela está no esquema para ser usada com o pacote DeSolve, para equações ordinárias diferencias (ODE, muito legal esse nome né, equações ordinárias hehe)

Nesse modelo, o coeficiente de transmissão, \beta, descreve a taxa instantânea na qual o número de hospedeiros infectados aumentam por indivíduo já infectado. Lembre-se que ele é baseado na lei de ação de massa, emprestado da física e da química. Ele assume que a taxa na qual um evento ocorre(novas infeções) depende da mistura aleatório dos reagentes (S,I), e que a taxa na qual na qual os reagentes colidem e reagem pode ser descrita por uma simples constante, \beta. Assim como a densidade de qualquer uma das moléculas (S,I) aumenta, aumenta também taxa de transmissão. A taxa de transmissão é então a taxa instantânea de novos infectados po unidade de tempo.

Indivíduos resistentes podem ser resistentes por uma ou duas razões. Eles podem morrer, ou podem desenvolver imunidades. Em ambos os casos, nos assumimos que eles não podem pegar a doença de novo. Como esse modelo assume uma população de tamanho constante, nos continuamos a contar todos os indivíduos R, independente se eles morreram ou se tornaram imunes.

Esses indivíduos se tornam resistentes a doença em uma taxa per capita \gamma. A taxa \gamma é também a média do tempo de residencia, ou duração da doença.

A incidência é o numero de novos afetados ou casos ocorrendo dentro de um intervalo de tempo. Essa definição faz uma versão de tempo discreto da taxa de transmissão. A prevalência é a fração da população que está infectada \frac{I}{N}.

Ao observar essa dinâmica ao longo do tempo, num exemplo onde por exemplo, iniciamos a maior parte dos indivíduos como Susceptíveis, aqui 9998 susceptíveis, um individuo infectado e um resistente.

01

Veja que com apenas um individuo infectado, conseguimos uma epidemia, um “outbreak”.

Uma pergunta comum, talvez a primeira a surgir em ecologia é perguntar sobre quais circunstâncias uma epidemia vai acontecer? Outra forma de perguntar isso é perguntar sobre quais condições esperamos I > 0. Nos podemos colocar \frac{dI}{dt}>0 e resolver para uma condição interessante sobre o que é preciso para uma epidemia.

0<\beta I S - \gamma I[/latex]  e  [latex]\frac{\gamma}{\beta}<S[/latex]  O que isso nos diz? Primeiro, como podemos dividir por I, isso significa que se ninguém está infectado, então uma epidemia não pode acontecer, senão teríamos divido por zero. Mas isso é mais que intuitivo, um equilíbrio instável em 0. Segundo, essa equação nos diz que uma epidemia vai ocorrer se a densidade de indivíduos Susceptíveis é maior que [latex]\frac{\gamma}{\beta}[/latex]. Se nos considerarmos uma pre-epidemia onde [latex]S \approx N[/latex], então diminuindo a abundância da população o suficiente pode conter a expansão da doença. Isso é o porque epidemias tendem a ocorrer em populações com altas densidades, como hospedeiros na pecuária, como o gado, ou historicamente como já vimos em populações humanas ou escolas.  Vacinações são uma forma  de reduzir S sem reduzir N. Se um número suficiente de indivíduos na população são vacinados para reduzir S abaixo de [latex]\frac{\gamma}{\beta}[/latex], isso tende a proteger os indivíduos não vacinados também.   Outra representação comum é chamada de "forca de infeção" ("force of infection") ou taxa de produção básica da doença ("basic reproductive rate of the disease"). Se nos assumirmos que em uma grande população [latex]S \approx N[/latex], então temos que  [latex]R_0 = \frac{\beta N}{\gamma}[/latex]  onde [latex]R_0[/latex] é a taxa básica de reprodução da doença. Se [latex]R_0>1, então uma epidemia é plausível.

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:
Stevens, M. H. – 2011 – A Primer of Ecology with R. Srpinger

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
SIR <- function(t, y, p) {
    {
    S <- y[1]
    I <- y[2]
    R <- y[3]
    }
    with(as.list(p), {
                    dS.dt <- -B * I * S
                    dI.dt <- B * I * S - g * I
                    dR.dt <- g * I
                    return(list(c(dS.dt, dI.dt, dR.dt)))
                })
}
 
N <- 10000
I <- 1
R <- 1
S <- N - I - R
parms <- c(B = 0.01, g = 4)
meses <- seq(0, 1.5, by = 0.01)
 
#install.packages("deSolve")
library(deSolve)
SIR.out <- data.frame(ode(c(S, I, R), meses, SIR, parms))
colnames(SIR.out)<-c("tempo","S","I","R")
 
matplot(meses, SIR.out[, -1], type = "l", lty = 1:3,frame=F,
        xlab="Tempo em Meses",ylab="N",lwd=2,col=2:4)
legend("right", c("Resistentes", "Infectados", "Susceptiveis"), lty = 3:1, col = c(4,3,2), bty = "n",lwd=2)

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")

Outros casos de estabilidade em populações.

Certo, a gente começou a olhar aqui, como modelar duas espécies usando a ideia do crescimento logístico de Lotka-Volterra, depois aqui como avaliar o resultado da competição entre duas desenhando isolinhas, olhamos aqui um pouco melhor como elas funcionam, mas sempre usamos o mesmo exemplo, os mesmos valores para as interações, que resultavam em coexistência. Mas agora vamos olhar algumas das outras possibilidades.

Como vimos nas análises anteriores, o fator determinante parece ser como são as interações entre as espécies, que no caso da coexistência entre duas espécies tinhamos:

 (\alpha_{12} < \alpha_{22}) \land  (\alpha_{21} < \alpha_{11}) [/latex]  Caso alguém não conheça esse simbolo, "[latex]\land[/latex]", esse é o "e" da matemática, vem da <a href="http://pt.wikipedia.org/wiki/%C3%81lgebra_booleana">álgebra booleana</a>, que merece um post aqui cedo ou tarde, mas estamos querendo dizer aqui que isso [latex] (\alpha_{12} < \alpha_{22})[/latex] e [latex] (\alpha_{21} < \alpha_{11}) [/latex] tem que acontecer ao mesmo tempo, não da para tirar uma conclusão olhando uma das comparações sozinhas, temos que olhar as duas e então chegamos a uma conclusão.  Veja que isso pode parecer complexo, mas traduzindo, temos que para as espécies coexistirem estavelmente, os efeitos delas, nelas mesmas precisam ser maiores que os efeitos na outra espécie.  E o resultado é esse gráfico que já vimos várias vezes.  <a href="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg"><img src="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg" alt="01" width="480" height="480" class="aligncenter size-full wp-image-1694" srcset="http://recologia.com.br/wp-content/uploads/2013/10/01.jpg 480w, http://recologia.com.br/wp-content/uploads/2013/10/01-150x150.jpg 150w, http://recologia.com.br/wp-content/uploads/2013/10/01-300x300.jpg 300w" sizes="(max-width: 480px) 100vw, 480px" /></a>   Mas então como são os outros casos?   Quando a gente ver…  [latex] (\alpha_{12} < \alpha_{22}) \land  (\alpha_{21} > \alpha_{11})

A espécie 1 pode invadir quando rara, mas a espécie 2 não pode.
Isso leva a exclusão competitiva pela espécie 1, sempre. A espécie 1 vence sempre. Isso também é citado como equilíbrio de borda, porque uma das espécies ocorre somente na borda do gráfico, do state-space.
A gente não ta falando que a espécie 1 ganha sempre? Mas lembre-se que se ela for zero, só quando algum indivíduo dela invadir que a população vai começar a crescer, ou seja, só na borda de N1=0 que a espécie 2 está na jogada. O equilíbrio onde N1>0 é chamado de equilíbrio interno.

No gráfico vemos:

02

Agora quando a gente ve…

 (\alpha_{12} > \alpha_{22}) \land  (\alpha_{21} < \alpha_{11}) [/latex]  A espécie 1 não consegue invadir o sistema, mas a espécie 2 consegue. O que gera o inverso do caso anterior, aqui a espécie 2 vence, a espécie 2 exclui a 1 por competição. Mas lembre-se que essas equações e suas relações são derivações de simples equações da reta e tem que ser pensado a limitações delas, por exemplo não usar pontos negativos. Mas de qualquer forma temos podemos extrair hipóteses testáveis daqui.  O gráfico dela fica assim, o inverso do anterior...  <a href="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg"><img src="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg" alt="03" width="480" height="480" class="aligncenter size-full wp-image-1697" srcset="http://recologia.com.br/wp-content/uploads/2013/10/03.jpg 480w, http://recologia.com.br/wp-content/uploads/2013/10/03-150x150.jpg 150w, http://recologia.com.br/wp-content/uploads/2013/10/03-300x300.jpg 300w" sizes="(max-width: 480px) 100vw, 480px" /></a>   Agora partimos pro caso legal.  [latex] (\alpha_{12} > \alpha_{22}) \land  (\alpha_{21} > \alpha_{11})

Porque legal? Porque aqui a gente cria um equilíbrio interno instável, a exclusão de alguma das espécies vai ocorrer, mas qualquer uma pode vencer. Isso normalmente é chamado de controle do fundador, porque quem chega primeiro ta na vantagem, a decisão depende em parte das abundância iniciais que observamos. Aqui o equilíbrio cria uma saddle ou sela no “state-space” ou espaço de estado, espaço de possibilidades. Podemos dizer que de algumas direções, a sela funciona como um atrator, como na coexistência, mas de outro lado, ela age como um repulsor, como existe nos casos de exclusão competitiva. Mas para ver isso a gente tem que dar um passo além e usar matrizes jacobianas.

Mas antes vamos olhar o gráfico como ficou.

04

Veja para onde os vetores estão apontado, em cada parte do gráfico existe uma direção diferente, isso que estamos falando quando falamos de cela, é saddle em inglês, espero não ter traduzido errado.

Então podemos olhar agora nossos casos principais, coexistência, exclusão competitiva e controle do fundador, onde não temos um vencedor definido, depende das condições que começamos a observar o sistema.

05

Agora vale lembrar que usamos valores para as figuras ficarem “bonitinhas”. Mas em qualquer valor veremos essas coisas acontecendo. Saia inventado valores para as interações das espécies e veja o que acontece.

06

Bem pode parecer um bocado de trabalho, para chegar nessas conclusões não? Mas pense que no final vale a pena, agora a gente finalmente entende essa figura do livro de ecologia.

begonecology

Agora temos como fazer hipóteses subsidiadas em alguma teoria, perfeita não é, mas os pressupostos são simples, a salvo por alguns desvios, essas funções ajustam bem, e funcionam no mundo real. Sabemos que não é perfeito, e que em alguns casos elas falham, mas ainda assim ela nos trazem informação, já que sabemos como as coisas deveriam ser e se existe um desvio do que esperamos, podemos procurar o que esta acontecendo com alguma direção, já que se sabemos como o que acontece no mundo real é diferente desses casos, podemos ver também o que devemos mudar nesses modelos para tentar entender melhor o mundo. Ou seja, apesar do sacrifício de ter que entender todas esses gráficos e contas, no final ganhamos um melhor entendimento de como diabos o mundo funciona. E para entender essas coisas não existe muita mágica, se ficar batendo cabeça, cedo ou tarde a gente começa a entender.

Mas é isso, o script está no repositório do github, além de aqui em baixo, e logo vamos começar a tentar entender o que diabos é a tal da matriz jacobiana.

Referencia:
M Henry H Stevens 2011 A primer of ecology with R. Springer

###################################################
### Casos de equilibrio para duas populações 
###################################################
 
### Coexistencia estavel
 
 
layout(matrix(c(1:4),ncol=2,byrow=T))
 
a <- matrix(c(0.01,0.005,0.005,0.01),ncol=2,byrow=TRUE)
 
plot( (a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]),
     (a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2, pch=1, col=1, axes=FALSE,
     ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
     ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s')
 
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2)
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
box()
arrows(c(25,25,25), c(25,25,25), c(25,50,50), c(50,25,50),
       col=1, length=.12, lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(150,150,150), c(150,150,150), c(150,120, 120), c(120,150,120),
      col=1, length=.12,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(10, 10, 10), c(125,125,125), c(10,30,30), c(105,125,105),
       length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
arrows(c(125,125,125), c(10, 10,10), c(125,105,105), c(30,10,30), 
       length=0.12, col=1,  lwd=c(1,1,2), lty=c(2,1,1))
 
### Espécie 2 pode invadir, a espécie 1 não
 
a <- matrix(c(0.01, 0.02,0.005, 0.01), ncol=2, byrow=TRUE)
 
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2,  axes=FALSE,
      ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
      ylab=expression("N"[2]), xlab=expression("N"[1]))
 
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]), lwd=1,c(0, expression(1/alpha[22]),expression(1/alpha[12])))
box()
xs <- c(1, 1, .4, .4) * a[1,1]^-1
ys <- c(.2, .4, .4, .7)*a[2,2]^-1
 
arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4],lty=2:1, lwd=c(1,1), length=0.12)
points(0,1/a[2,2], cex=2)
 
### Espécie 1 pode invadir, espécie 2 não
 
a <- matrix(c(0.01, 0.005,0.02,0.01), ncol=2, byrow=TRUE)
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lwd=1, lty=2,  axes=FALSE,
      ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
      ylab=expression("N"[2]), xlab=expression("N"[1]))
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
 
ys <- c(1, 1, .4, .4) * a[1,1]^-1
xs <- c(.2, .4, .4, .7)*a[2,2]^-1
 
arrows(xs[1:3], ys[1:3], xs[2:4], ys[2:4],lty=1:2, lwd=c(1,1), length=0.12)
points(1/a[1,1],0, cex=2)
 
### Nenhuma especie pode invadir, com equilíbrio instável
 
a <- matrix(c(  0.01, 0.02,0.02, 0.01), ncol=2, byrow=TRUE)
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], lty=2,  axes=FALSE,
      ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
      ylab=expression("N"[2]), xlab=expression("N"[1]))
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
 
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
 
pxys <- matrix(c(.13,.1,.25,.1,
                 .1,.13,.1,.25,
                 .115,.115,.25,.25,
                 .7,.8,.4,.8,
                 .8,.7,.8,.4,
                 .75,.75,.4,.4,
                 .25,.4,.1, .6,
                 .4,.25,.6, .1),
              nc=4, byrow=T)
xys <- pxys
xys[,1:2] <-  a[1,1]^-1 * pxys[,1:2]
xys[,3:4] <-  pxys[,3:4] * a[2,2]^-1
arrows(xys[,1], xys[,2], xys[,3], xys[,4],lty=c(1,2,3,1,2,3,2,1), lwd=c(1), length=0.12)
n2star <- (a[1,1] - a[2,1])/(a[1,1]*a[2,2] - a[1,2]*a[2,1])
n1star <- (a[2,2] - a[1,2])/(a[1,1]*a[2,2] - a[1,2]*a[2,1])
points(n1star,n2star, cex=2)
 
#Interações ao acaso
par(mfrow=c(3,3))
 
for(i in c(1:9)) {
a <- matrix(runif(4,0,2),ncol=2,byrow=TRUE)
 
plot( (a[2,2]-a[1,2])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]),
     (a[1,1]-a[2,1])/(a[2,2]*a[1,1] - a[1,2]*a[2,1]), cex=2, pch=1, col=1, axes=FALSE,
     ylim=c(0, max(1/a[2,2], 1/a[1,2])), xlim=c(0, max(1/a[1,1], 1/a[2,1])),
     ylab=expression("N"[2]), xlab=expression("N"[1]), pty='s')
 
curve(1/a[2,2] - (a[2,1]/a[2,2])*x, 0, 1/a[2,1], add=TRUE, lty=2)
curve(1/a[1,2] - (a[1,1]/a[1,2])*x, 0, 1/a[1,1], add=T)
axis(1, c(0, 1/a[1,1], 1/a[2,1]),c(0, expression(1/alpha[11]),expression(1/alpha[21])))
axis(2, c(0, 1/a[2,2], 1/a[1,2]),c(0, expression(1/alpha[22]),expression(1/alpha[12])))
legend("topright",legend=paste(c("1,1=","1,2=","2,1=","2,2="),round(a,digits=2)),bty="n")
}