Arduino

DSC04416

Desde que eu me interessei por evolução, eu comecei a procurar sobre como são os métodos utilizados para conseguir dados para trabalhar com genética e evolução, mais especificamente, como se faz PCR e depois sequência o resultado da PCR, porque esses são os dados que a gente vê nos problemas de bioinformática projeto Rosalind, aquele monte de letrinha que são sequências de DNA, RNA, etc.

Esse tipo de produto, máquinas que fazem PCR, custam uma fortuna, pelo menos eu não tenho condições de comprar, definitivamente um produto que não chega na mão de qualquer curioso. Inclusive é ruim pois, por exemplo no segundo grau ou mesmo na faculdade, quando os Azão azinho da vida poderiam fazer sentido e estudar ser mais legal, tudo que se vê é teoria e um papo sobre ervilha comumente, nada mais.

Fazer PCR depende de um termociclador, um aparelho caro que só grandes empresas produzem e que talvez o preço tenha a ver com a baixa procura, sei-la. Mas procurando sobre termocicladores, eu conheci um projeto chamado Open PCR.

ca5259b3fb380c31847aaa3cd487e6f0-original

Basicamente dois caras ficaram cansados do alto custo para se adquirir um termociclador, e queriam ter um. Viram como faz um termociclador e fizeram o deles em casa. Eles inclusive comentam que uma das motivações, era que eles queriam fazer experimentos com genética, com alunos, tipo ver alguns genes conhecidos relacionados a cor de cabelo, polegar opositor, etc.

A gente raspa a bochecha de todo mundo da sala, coleta material, pega uns primers prontos e outras coizinhas e tenta amplificar esses genes na galera e corre no gel, e isso torna a aula de genética na escola muito mais interessante. Mas tudo empacava nos custos, então eles reduziram os custos construindo o próprio termociclacdor, e agora eles tem uma pequena empresa que vende esses termocicladores por 600 dólares, com o intuito de abrir essa possibilidade a todos.

Mas não para por ai, eles deixam disponível todas as especificações do produto deles, logo, se você produzir as peças, que eles ensinam como, todas as especificações estão abertas, você pode fazer o seu próprio termociclador em casa, ou até melhorar o projeto, um hardware open-source.

Bem isso é algo que chama a atenção, pessoas comuns fazendo algo complexo em casa, como isso é possível?

Foi ai que eu vi que esses caras são na verdade parte de um movimento cultural chamado de Makers.

Esses caras, usaram uma coisa chamada Arduino, que é uma plaquinha com um microcontrolador que da para ligar no computador e programar ele, e ele tanto recebe informação de coisas como manda informação, pro computador ou coisas.

Como assim?

O Termociclador é uma máquina, que tem um sensor de temperatura, uma coisa que esquenta a amostra e outra que resfria. Os caras bolaram um lugar para colocar amostras, um pratinho que recebe amostras, e fizeram um programa que le a temperatura e faz os ciclos da PCR, primeiro esquenta a amostras, ai resfria, pelo tanto de vezes necessários de acordo com o programa, como com o Arduino o computador sabe a temperatura, o Arduino lê um sensor e manda a informação pro computador, o computador pode mandar informação para aquecer ou resfriar com um cooler de computador, e eu não sei o nome da coisinha que esquenta, tornando tudo isso possível.

Pode parecer complexo, mas se você ouvir a historia do porque os caras inventaram a Arduino, você pode mudar de idéia, aqui tem um documentário rapidinho sobre os caras que inventaram o Arduino.

Vale a pena ver, mas denovo a boa historia que nos anima, alguns caras queriam ensinar sobre eletrônica, e como funcionam as coisas. Mas você tem que ver tanta teoria antes de fazer um sistema que pisca uma luz, que desanima todo mundo, um artista por exemplo não quer saber muito de eletrônica para fazer luzes piscarem de forma bonita. Mas eles falaram não, e se a gente conseguir simplificar tudo isso, fazer um sisteminha que você pluga no computador, diz para o microcontrolador o que ele tem que fazer com um programinha e pimba, ta funcionando, assim as pessoas conseguem fazer coisas legais, e isso estimula todo mundo a querer aprender mais para fazer coisas mais legais.

Mas acredito que não imaginaram nas proporções que isso iria tomar, isso virou a diversão de quem queria aprender eletrônica. Mas além disso pessoas começar a refazer produtos que ja existem nas suas casas.

So que no meio científico, alguns caras, como os que fizeram o OpenPCR, tornaram mais barato e talvez mais democrático a ciência.

Mas não para aqui.

Um site chamado Instructables, onde as pessoas podem depositar seus projetos, um cara fez um termociclador que segundo ele da para montar a um custo de 85 dolares. Esse projeto aqui.

FMXHYHAH4AGJRQD.LARGE

Se até agora eu não fui claro, aqui o que aparece é uma oportunidade bem legal, que é uma forma muito legal de coletar dados.

No r-bloggers, eu ja tinha visto alguns caras, como aqui e aqui, que estão conectando o computador com o mundo através do Arduino para coletar dados, e os dados ja caem direto no R, muito legal.

Basicamente você compra um sensor como esse por exemplo.

sku_138531_1_small

Ai você coloca ele no arduino e faz o arduino mandar os dados para o seu computador, sem muito esforço é até possível você colocar um acessório no arduino para ele mandar os dados wireless para o seu computador, so deixar um pilha nele.

So que os tipos de sensores, de coisas que da para medir são imensas, temperatura, umidade, luz, movimento, vibração, coordenadas geográficas, da para medir de tudo com os muitos tipos de sensores ou combinações deles.

sku_142834_1

Agora se você pensa num experimento, como essas coisas são baratinhas e um sensor desse na maioria dos casos não passa de 10 dolares, deve ser muito útil, além de aumentar a quantidade de dados que você pode coletar e a qualidade das medias, ja que sempre tudo é medido por um único coletor, um computador, que você sempre pode verificar e testar a qualidade dos dados que ele entrega. Ele pode ficar coletando dados o tempo que você especificar, até 24 horas por dia.

Eu imagino que um aquário com câmeras registrando o movimento de peixes, ou dados de micro-clima em qualquer lugar, desde a mata até arvores, ou até características limnológicas da água, qual o micro-clima dentro de uma caverna e como isso influência no comportamento de morcegos, você deixa uma caixinha la coletando dados, somado a mais observações pessoais e conhecimento sobre quais são morceguinhos que estão la, deve dar para fazer muita coisa legal.

E grande parte dessas coisas já estão construídas, o trabalho é juntar vários projetos e fazer eles funcionarem juntos a nosso favor, coletando os dados que precisamos para o nosso experimento, e dar um jeito de gravar a informação, seja em memoria, mandar wireless para o conforto da onde você esta, além de manter o conforto do seu objeto de estudo, ja que você não vai estar la atrapalhando os bixinhos.

Além disso, assim como o R, existe uma comunidade gigante trabalhando nisso, então não é algo que você nunca encontrara ninguém para tirar duvidas ou ajudar a construir essas coisas, você não vai embarcar sozinho nisso.

Mas se ainda sim dúvidas forem colocadas sobre a viabilidade disso, até o LHC, o experimento mais caro que eu conheço, que deve juntar os físicos, engenheiros mais fodas do mundo, tem uns Arduinos pelos cantos coletando dados. Olhem o simbolo nessa caixinha, o mesmo da minha plaquinha la do começo do post.

lhcarduino

Então como não pareceu tão difícil assim, até eu fui olhar como é uma plaquinha dessa, comprei uma da china por uns trocados e fui ver como é a idéia disso.

Basicamente você faz um programinha na linguagem C, que não é o C propriamente dito, mas tudo funciona como no C. Fiz um programinha para piscar luzes em sequência. Meu programinha é assim.

/* Piscar luzes do pino 13 a 3 */
 
//inicializa variaveis que vou usar
int i;
int led[11];
 
// A rotina de setup, inicia assim que se aperta o reset da placa
void setup() {               
 
  //Atribui valor aos pinos
  for(i=0;i<=10;i++) {
    led[i] = i+3;
  }
 
  // Inicializa os pinos como saida
 
  for(i=0;i<=10;i++) {
    pinMode(led[i], OUTPUT);
  }
 
}
 
// A rotina de loop que vai rodar para sempre
void loop() {
  digitalWrite(led[0], LOW);
  for(i=1;i<=10;i++) {            // Loop com a luzes indo
    digitalWrite(led[i], HIGH);   // Acende a luz
    delay(500);               
    digitalWrite(led[i], LOW);    // Apaga a luz
  }   
  for(i=9;i>1;i--) {              // Loop com a luzes voltando
    digitalWrite(led[i], HIGH);   // Acende a luz
    delay(500);               
    digitalWrite(led[i], LOW);    // Apaga a luz
  }  
  digitalWrite(led[0], HIGH);     // Antes de recomeçar eu deixo a ultima ligada
}

E o resultado é:
DSC04415

Luzinhas piscando numa sequência que eu determinei, bem bobinho e inútil, mas a possibilidade de coletar dados e guardar pode ser bem inspiradora e animar o aprendizado. As pessoas vão tão longe com essas coisas, que até laser para matar mosquitos, e apenas os que transmitem malaria é possível fazer.
Mas realmente essa parece uma possibilidade bem legal e que pode render alguns frutos legais num futuro, mesmo não sendo algo muito sofisticado ou elegante, quem sabe.
Mas enquanto nada acontece o jeito é continuar estudando. Até a próxima.

Calculando derivadas usando o R

A palavra cálculo sempre anda bem longe da biologia. Aliás, muitas coisas são mal vistas pelas pessoas de modo geral.

Um exemplo comum é a fórmula da baskara, todo mundo ja deve ter visto algo assim pela internet, no facebook, twitter ou sei la.

fórmuladebáskara

Sem querer discutir como funciona o sistema educacional, mas eu mesmo me pergunto as vezes, para que diabos serve isso ou aquilo que a gente le nos livros. Mas ao mesmo tempo eu também sempre tenho dúvidas sobre como diabos é possível chegar a algumas conclusões, como por exemplo, como as pessoas medem a distância de planetas, se a gente só ve a luz dos planetas chegando na terra.

Se isso ainda é muito distante da vida cotidiana, pense se você realmente compreende como funciona a lampada fluorescente, a televisão de plasma ou como a geladeira gela as coisas.

O que fica faltando muitos vezes é a ligação de uma coisa com a outra, a gente fica com um monte de informação perdida sem conseguir entender a ligação entre elas, ou a utilidade de muitas coisas, como a fórmula de baskara, sendo que elas são a solução das outras dúvidas que temos e não temos resposta. E isso vale para qualquer coisa que a gente aprende, não só ciências ou matemática ou química ou sei la o que mais.

Mas o lance desse post é calcular derivadas com o R então chega desse mimimi senão vamos encalhar aqui.

A gente ja viu várias coisas de crescimento populacional que utilizavam a palavra derivada, derivada parcial, etc.

Segundo o wikipedia, derivada é a taxa de variação instantânea de uma função.

Mas o que isso quer dizer sem ter que fazer muita conta?

Bem primeiro vamos definir uma função simples no R.

f(x)=x^2

Uma função bem bobinha.
Dai no r a gente pode definir ela assim:

função <- function(x) { return(x^2) }

Ou seja a gente cria uma função, que para um determinado valor de x, ela retorna esse valor elevado ao quadrado.

Dai é fácil, é só sair usando ela para vários valores.

> função(1) [1] 1 > função(2) [1] 4 > função(3) [1] 9 > função(-2) [1] 4 > função(1.5) [1] 2.25

Um elevado ao quadrado é um mesmo, dois elevado ao quadrado é quatro, três elevado ao quadrado é nove, menos dois elevado ao quadrado é quatro também, porque um número negativo vezes outro numero negativo da um número positivo. Certo até aqui ta fácil.

Mas como a gente gosta de gráfico, e ta fácil, vamos fazer um gráfico, onde a gente define um valor para x no eixo x e no eixo y colocamos a resposta da função:

01

Agora pensa assim, quando o x era 1, o y é 1, quando o x era 2 o y era 4, então o aumento de 1 unidade de x aumentou 3 unidades em y, mas de 2 para 3 no eixo x, também aumentamos também 1 unidade em x, mas o y subiu de 4 para 9, um aumento de 5 unidades, muito mais.

Logo a taxa de aumento de 1 pra 2 é menor que de 2 para 3. O quão menor ou maior é mais ou menos o que a derivada mede , ou deriva, sei la se é essa a forma correta de falar. Eu nunca sou bom com definições. Mas isso é uma das cosias que interessa na derivada para o calculo.

Mas o legal é que tem como obter uma fórmula que diz como muda essa variação. Existe regras, teoremas, muito coisa para calcular derivadas de vários tipos de funções, como a gente so quer ter uma idéia que como é, a gente usa o r para calcular pra gente.

E como faz isso? É so usar a função deriv(), ou D(), da na mesma, de um ?deriv para ver mais sobre a função.

Mas a gente usa

dfunção <- deriv(~ x^2, "x")

O primeiro argumento que demos a função deriv é a formula de f(x), depois do “~”, e no segundo argumento quem é para ser derivado nessa formula, o x entre aspas. E essa função gera uma expressão com a derivada.

dfunção expression({ .value <- x^2 .grad <- array(0, c(length(.value), 1L), list(NULL, c("x"))) .grad[, "x"] <- 2 * x attr(.value, "gradient") <- .grad .value })

Note que um dos atributos é o resultado da derivada para um determinado número, esse aqui “.grad[, “x”] <- 2 * x". Essa é a derivada, se você como eu não sabe sair do [latex]f(x)=x^2[/latex] e chegar nesse [latex]\dot{f}(x)=2 \cdot x[/latex] (tem um monte de forma de escrever derivada, como eu to perdido eu vou usa esse pontinho ai, não sei usar o f linha no latex ainda), a gente pode consultar algum lugar, como o wolfram alpha para ver se ta certo. Olha aqui.
Ta certo, agora a gente pode atribuir valores a x e conseguir resultados

> x <- 1 > attr(eval(dfunção),"gradient") x [1,] 2 > x <- -2:2 > attr(eval(dfunção),"gradient") x [1,] -4 [2,] -2 [3,] 0 [4,] 2 [5,] 4

Agora o que isso quer dizer? O que esses números querem dizer? Essa é “inclinação da tangente à curva como a derivada de f(x)”, entre aspas porque eu to copiando e colando do wikipedia. Mas a gente pode copiar a formulinha para construir essa reta, essa formula aqui:

g(x)=\dot{f}(x) \cdot (x-a) + f(a)

Note logo de cara que isso é algo parecido com a+bx, uma coisa que multiplica um número e depois adiciona outra coisa, por isso que sai uma reta. E o “a” dessa formula ai é o ponto que estamos derivando, notem que cada ponto pode ter uma reta diferente.

Essa eu aprendi isso no primeiro livro (ou não) de cálculo que eu entendi, o Manga guide to calculus, esse foi legal, historinhas ensinando cálculo, mas é legal porque cálculo esta dentro do contexto da historinha, uma jovem repórter entendendo o mundo das notícias. Mas então dentro de um balãozinho ali tem um carinha falando isso ai.

Mas usando ela eu faço a reta para derivada de x=2

02

Certo a gente ta vendo uma reta, que encosta um pouquinho no gráfico da função.
Mas como para entender sempre é bom sair fazendo gráfico para um monte de valor, vamos fazer a derivada de vários valores, plotar a reta e ver o que acontece.

03

Olha ai, -1 a reta ta descendo, 0 ela ta retinha, 1 ela ta subindo, 2 ela ta subindo mais inclinada que o 1 e assim vai. Ou seja a derivada nos diz como vai mudando o valor de f(x) em relação ao x.

Mas f(x) é muito whatever, para que serve a função de f(x)=x^2? Ora, para aprender, porque até eu consigo calcular a derivada dessa dai.

Agora essa função aqui nos interessa.

N_{t}={{N_{0} \cdot e^{r \cdot t} }\over{1 + \alpha \cdot N_{0} \cdot (e^{r \cdot t} -1)}}

A gente ja viu nesse post aqui como se muda de crescimento discreto para crescimento contínuo. A gente só aplicou a mesma coisa para crescimento contínuo.
Notem que a parte de cima é a mesma coisa, só mudou que tem um termo denominador agora, que é relativo a capacidade suporte, que estamos tratando como {k}={1\over{\alpha}}. Essa mesma formula esta no wikipedia aqui, so que usando k ao invés de alpha, na parte de ecologia.

A gente pode então fixar os valores de r=0.5 , No=2 e alpha = 0.001 (k=1000) e escrever essa função no R assim:

flogistico <- function(t) { return((2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1))) }

Notem que eu ja substitui os valores de populacão inicial com 2, r com 0.5 e alpha com 0.001, assim só sobrou o t de tempo. E agora ela é um f(t), a outra é um f(x), mas pra gente t é t de tempo, pode ser dia, mês ou ano, algo que a gente mede. Mas é igual, so um x com algum sentido para a gente agora. O f(x) era uma resposta qualquer, mas aqui o f(t) é o tamanho da população naquele tempo.

A gente pode usar essa função para testar

> flogistico(0) [1] 2 > flogistico(1) [1] 3.29317 > flogistico(2) [1] 5.417945 > flogistico(3) [1] 8.901394

Legal que no tempo 0, a população é a população inicial, mas podemos calcular agora para qualquer tempo e podemos colocar num gráfico.

04

Esse gráfico a gente ve direto correto?

Como a gente colocou no gráfico agora, é uma função e da para derivar igualzinho a gente fez com o f(x)=x^2, o que valia para f(x) vale para a nossa função do tamanho populacional ao longo do tempo, as contas são iguais, podíamos não ter perguntas da relação de x com f(x), só que temos perguntas em relação a como o tempo afeta o tamanho populacional, se a população ta crescendo, diminuindo, vai extinguir, quanto podemos extrair dela sem extinguir ela, quanto tempo demora para ela voltar a crescer bastante e assim vai. Então vamos derivar ela

> dlogistico <- deriv(~ (2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1)), "t")

O comando é o mesmo, só derivamos usando t ao invés de x.
Podemos agora calcular o gradiente de variação como na outra função

> t=0 > attr(eval(dlogistico),"gradient") t [1,] 0.998 > t=5 > attr(eval(dlogistico),"gradient") t [1,] 11.63201 > t=30 > attr(eval(dlogistico),"gradient") t [1,] 0.07629933

Mas como esse monte de número tem pouco significado, vamos fazer um gráfico também com a inclinação da derivada para ver como fica nesse caso.

05

Olha ai, a inclinação é como está o crescimento da população, vamos fazer para vários momentos para ver se bate com a nossa intuição.

06

Opa, com a minha intuição bateu. Como a gente ja tinha comentado aqui, no começo o crescimento é alto, apesar da população pequena, e mais para frente o crescimento começa a decair. Agora os últimos gráficos do post de pertubação, onde um é com retas e o outro é com curvas deve fazer mais sentido, um é da derivada, a outro mostrando a função em si.

Podemos ainda criar uma sequência de valores de tempo para avaliar, avaliar somente a inclinação, que é uma informação bem interessante, lembrando que estamos até avaliando valores quebrados aqui, continuamente.

07

E olha ai que legal, ja falamos isso antes também aqui. O crescimento começa baixinho, ai vai aumentando, até chegar num pico, ai começa a decrescer e para la, até… ter uma pertubação, como foi dito aqui.

E tharam, se derivadas de f(x)=x^2 são inúteis para a biologia, essa ai é bem útil, e temos falado dela recorrentemente. E apesar de ser difícil fazer aquelas contas com muitas e muitas letrinhas, podemos pelo menos pegar o “feeling” do que quer dizer uma derivada para captar varias idéia legais que utilizam elas depois. E usar o R para fazer as contas. Muitas pessoas desenvolveram algoritmos complexos que fazem o trabalho pesado, é so a gente usar.
Além do mais eu sempre achei legal esses gráficos com esse risquinho reto encostando na curvinha.

Bem o script esta no repositório recologia do github além de aqui em baixo.

função <- function(x) { return(x^2) }
 
função(1)
função(2)
função(3)
função(-2)
função(1.5)
 
#figura 1
curve(função(x),-5,5,ylab="função(x)",xlab="Valor de x",frame=F)
 
dfunção <- deriv(~ x^2, "x")
 
dfunção
 
x <- 1
attr(eval(dfunção),"gradient")
 
x <- -2:2
attr(eval(dfunção),"gradient")
 
#figura 2
x <- 2
curve(função(x),-5,5,ylab="função(x)",xlab="Valor de x",frame=F)
curve(attr(eval(dfunção),"gradient")*(dx-x)+função(x),-5,5,ylab="função(x)",xlab="Valor de x",add=T,
      xname = "dx",lty=2,lwd=2,col="red")
 
#figura 3
curve(função(x),-5,5,ylab="função(x)",xlab="Valor de x",frame=F)
x <- -1
curve(attr(eval(dfunção),"gradient")*(dx-x)+função(x),-5,5,ylab="função(x)",xlab="Valor de x",add=T,
      xname = "dx",lty=2,lwd=2,col=2)
x <- 0
curve(attr(eval(dfunção),"gradient")*(dx-x)+função(x),-5,5,ylab="função(x)",xlab="Valor de x",add=T,
      xname = "dx",lty=2,lwd=2,col=3)
x <- 1
curve(attr(eval(dfunção),"gradient")*(dx-x)+função(x),-5,5,ylab="função(x)",xlab="Valor de x",add=T,
      xname = "dx",lty=2,lwd=2,col=4)
x <- 2
curve(attr(eval(dfunção),"gradient")*(dx-x)+função(x),-5,5,ylab="função(x)",xlab="Valor de x",add=T,
      xname = "dx",lty=2,lwd=2,col=5)
legend("center",lty=2,lwd=2,col=2:5,legend=paste("Derivada de", -1:2),bty="n")
 
#Nt=(No*exp(r*t))/(1+alpha*No(exp(r*t)-1))
#Vamos fixar o r=0.5 ,  No=2 e alpha = 0.01 (k=100)
 
flogistico <- function(t) { return((2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1))) }
 
flogistico(0)
flogistico(1)
flogistico(2)
flogistico(3)
 
#figura 4
curve((2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1)),0,30,ylab="Tamanho da população",xlab="Tempo",xname="t",frame=F)
 
dlogistico <- deriv(~ (2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1)), "t")
 
t=0
attr(eval(dlogistico),"gradient")
 
t=5
attr(eval(dlogistico),"gradient")
 
t=30
attr(eval(dlogistico),"gradient")
 
#figura 5
curve((2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1)),0,30,ylab="Tamanho da população",xlab="Tempo",xname="t",frame=F)
t <- 19
curve(attr(eval(dlogistico),"gradient")*(dl-t)+flogistico(t),0,30,add=T,xname = "dl",lty=2,lwd=2,col="red")
 
#figura 6
curve((2*exp(0.5*t))/(1+0.001*2*(exp(0.5*t)-1)),0,30,ylab="Tamanho da população",xlab="Tempo",xname="t",frame=F)
t <- 19
curve(attr(eval(dlogistico),"gradient")*(dl-t)+flogistico(t),0,30,add=T,xname = "dl",lty=2,lwd=2,col="red")
t <- 8.5
curve(attr(eval(dlogistico),"gradient")*(dl-t)+flogistico(t),0,30,add=T,xname = "dl",lty=2,lwd=2,col="blue")
legend("bottomright",col=c("red","blue"),lwd=2,lty=2,legend=c("Derivada de 19","Derivada de 8.5"),bty="n")
 
#figura 7
t<-seq(0,30,by=0.01)
plot(t,attr(eval(dlogistico),"gradient"),type="l",frame=F,ylab="O quanto a população esta mudando (inclinação)",
     xlab="Tempo")

Existe mais de uma maneira de se fazer a mesma coisa, um exemplo com o modelo de Cinética enzimática de Michaelis-Menten

Bem, tem um monte de gente que diz frases legais sobre o R. Essa é uma:

This is R. There is no if. Only how.

O autor é o “Simon ‘Yoda’ Blomberg, R-help (April 2005)”

Na minha tradução usando meu inglês nórdico fica…

No R não existe se é possível, apenas como é possível.

Certo, mas vamos olhar uma coisa eu acho bem legal. O modelo do Michaelis-Menten, eu particularmente nunca fui bom de química, além de ser ruim em muitas outras disciplinas, mas eu acho esse modelo legal porque foi um dos primeiros modelos não lineares que eu vi, que fez sentido, que eu entendi de verdade, eu acho.
Basicamente ele relaciona como a velocidade que uma reação acontece com a concentração de uma enzima. Se você tiver um monte de reagente numa solução, ela ocorre devagarinho, mas a gente pode adicionar uma enzina para aumentar a velocidade, lembrando que uma enzina não é consumida na reação. Mas se a gente adicionar so um pouquinho de enzina, a gente aumenta um pouquinho, se a gente adiciona bastante, a gente aumenta mais a velocidade da reação, mas até um limite. So que no começo adicionar enzimas, aumentar a concentração desta na solução, aumenta bem a velocidade da reação, mas conforme a gente coloca mais e mais, a solução satura, o efeito de aumentar a velocidade da reação satura, ou seja tem um limite, um teto de aumento de velocidade, acima disso não adianta mais adicionar enzina que não da em nada.
Essa idéia pode ser descrita na forma de uma equação, o modelo de Michaelis-Menten, não tenho certeza se meu exemplo é muito bom, mas a idéia é por ai.

A modelo é o seguinte:

 {Vel} = {{Vmax\cdot Con}\over{Km+Con}}

Certo, então a gente tem duas constantes, Vmax e Km, Vmax é a maior velocidade de reação possível, enquanto Km é a constante que diz o quanto o aumento da concentração de algo, uma enzima por exemplo, aumenta a velocidade da reação, até o limite de Vmax.

Vamos aquele velho esquema de de gerar alguns dados, onde Vmax=3.4 e Km=0.4, os mesmo valores do artigo do wikipedia sobre o modelo.
Com os dados em mão, vamos fazer um gráfico para dar uma olhada.

01

Certo, primeiro vamos olhar aquelas características que falamos, Vmax é o teto do modelo, a tendencia máxima da reação. O segundo parâmetro, Km, é o momento também onde a velocidade da reação é metade do Vmax, algumas características legais.
Outra coisa é que esse modelo não precisa só ser pensado em nível de reações químicas com enzimas, por exemplo, uma vez eu tentei usar o modelo de Michaelis-Menten para relacionar o número de interações de parasitas com o número de estudos deles, ou seja quanto mais estudos, em mais hospedeiros esses os parasitas eram descobertos, só que isso tinha um limite, chega uma hora que não interessa quantos estudos mais sejam feitos, a gente ja sabe quem uma determinada espécie de parasita ataca, logo há uma saturação, e mais estudos dificilmente vão descobrir hospedeiros novos, mas se ha poucos estudos para uma especie, é bem provável que novas espécies de hospedeiros sejam descobertas. Pensando desse jeito, o Vmax é o total de espécies de hospedeiros que uma espécie de parasita pode parasitar e Km seria o quanto um estudo contribui para a descoberta do total de hospedeiros possíveis de um parasita, sendo que essa contribuição satura, diminui conforme a gente acumula mais estudos.

Certo mas chega de justificar utilidades, temos dados, temos um modelo e queremos ajustar valores para os parâmetros, para as constantes dele.
De cara a gente vê que o modelo não é linear, não é uma linha, logo a gente precisa ajustar um modelo não linear. Vamos fazer isso usando um modelo não linear, no R usando o comando nls, que serve para ajustar modelos não lineares.
Bem ele tem várias coisas a se estudar, mas a principio vamos simplesmente ajustar o modelo, ele não é muito diferente de usar o comando lm. A gente diz a formula, fornece os dados e a função faz o resto. Mas como a função otimiza os valores das constantes, a gente tem que partir de algum ponto, ou seja tem que fornecer um ponto de inicio para começar a otimizar. Vamos fornecer valores quaisquer sem pensar muito, mas existem formas elegantes de conseguir valores iniciais.

nls.menten<-nls(Vel~(Vmax * Con)/(Km + Con),
                data=dados,start=list(Vmax=1,Km=1),algorithm="default")
summary(nls.menten)

Notem que o modelo que a gente usou na formula:

Vel~(Vmax * Con)/(Km + Con)

É exatamente o modelo la de cima, so que escrito do jeito do R, e usamos Con como nome da concentração, além disso fornecemos os dados, o start é uma lista com valores iniciais para todas as constantes, temos duas aqui Vmax e Km, e o algorithm é a estrategia de otimizar os valores, a gente ta usando o default, mas existe varias maneiras de otimizar valores, no nls por padrão vem com três, mas existe outras ainda, mas isso é outra historia.

Mas continuando, temos como resposta:

Formula: Vel ~ (Vmax * Con)/(Km + Con) Parameters: Estimate Std. Error t value Pr(>|t|) Vmax 3.48185 0.03615 96.32 <2e-16 *** Km 0.44281 0.02307 19.20 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.07808 on 28 degrees of freedom Number of iterations to convergence: 7 Achieved convergence tolerance: 3.189e-07

Olha ae, as constantes estão até com o nome original da formula, bem parecido com o que usamos para gerar os dados, que foi Vmax=3.4 e Km=0.04. O parâmetros convergiram legal para esses valores, e não deveria ser diferente, ja que geramos esses dados a partir exatamente desse modelo e esses valores de constantes. Se considerar os desvios, os valores originais estão dentro do intervalo de confiança de 95%.
Podemos colocar esse modelo no gráfico, junto aos dados e ver como ele fica.

02

Ta mas agora vamos olhar uma coisa interessante. Podemos transformar esse modelo, os dados para um modelo linear, uma reta.
Difícil pensar né, como uma curva dessas pode virar uma reta. Mas tudo o que a gente tem que fazer é inverter o valor de concentração e de Velocidade da reação.

Primeiro vamos fazer um gráfico com 1/Vel e 1/Con e ver como fica, para avaliar isso visualmente.

03

Olha ae que louco, agora é so fazer uma regressão ai para estimar a equação da reta.
Mas vamos fazer o seguinte, vamos usar um glm, que é modelo linear genearilizado em português, eu acho, ou generalized linear model em inglês.

Bem, ajustar isso no r também é bem simples, basta usar a função glm(), da seguinte forma.

glm.menten<-glm(Vel ~c(Con^-1), data = dados,family = gaussian(link = "inverse"))
summary(glm.menten)

Agora você deveria perguntar, ue mas vc não vai inverter as coisas? Eu estou invertendo a concentração aqui “c(Con^-1)”.
Primeira coisa, eles esta dentro do concatenar “c()” para nao enviar ele (^1) para a formula, o valor para ser invertido, e a formula continuar igual, a+bx.
Outra coisa, olhem que eu estou elevando ela a -1, se você não lembra de exponenciação, de uma revisada rápida no wikipedia, mas essa regra é a seguinte…

{a^{-1}} = {1\over{a}}

A ultima questão é, porque isso só foi feito para a concentração (con) e não para a resposta Vel?

A resposta esta aqui “family = gaussian(link = “inverse”)”, todo glm usa uma distribuição e uma função de ligação, uma das muitas possibilidades da função de ligação é a inversa, ou seja 1/valor, ou seja aqui eu defini que usaremos a distribuição Gaussiana, que é a mesma coisa que distribuição normal que até hoje eu acho meio confuso os muitos nomes dela, mas então, aqui falamos que queremos que a resposta seja invertida, a função de ligação é a inversa, ou seja a função vai inverter os dados de Vel, o glm ta fazendo isso para a gente.

Mas sabendo que ta tudo invertido, o resultado fica.

Call: glm(formula = Vel ~ c(Con^-1), family = gaussian(link = "inverse"), data = dados) Deviance Residuals: Min 1Q Median 3Q Max -0.14465 -0.05057 -0.02119 0.05772 0.13804 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.287204 0.002982 96.32 <2e-16 *** c(Con^-1) 0.127178 0.005528 23.01 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 0.006095741) Null deviance: 15.60591 on 29 degrees of freedom Residual deviance: 0.17068 on 28 degrees of freedom AIC: -63.938 Number of Fisher Scoring iterations: 4

Agora temos uma equação de reta, a+bx, intercepto e inclinação, ou como queira chamar. Mas com uma rápida batida de olho, a gente vê que esses valores não se parecem nem um pouco com Vmax ou Km que usamos.
Mas se a gente olhar o ajuste deles aos dados…

04

Fica muito bom. Mas o que diabos significa esses valores de intercepto e inclinação aqui?

Aqui a jogada é a seguinte, os valores foram invertidos certo? Então vamos relembrar o modelo original.

 {Vel} = {{Vmax \cdot Con} \over {Km + Con}}

So que Vel ta escrito na forma de  Vel^{-1} que é o mesmo que  1\over{Vel} .

Como estamos olhando para uma igualdade, temos que fazer a mesma coisa dos dois lados certo? Logo…

 {Vel^{-1}} = ({{Vmax \cdot Con} \over {Km + Con}})^{-1}

Que é o mesmo que:

 {1\over{Vel}} =  {{Km + Con} \over {Vmax \cdot Con}}

Certo, nenhuma algebra muito intensa aqui, a gente so inverteu os valores dos dois lados, so que do outro lado como ja era um fração, ela inverteu.

Mas agora é uma soma em cima, ou seja:

 {1\over{Vel}} =  {{Km} \over {Vmax \cdot Con}} + {{Con} \over {Vmax \cdot Con}}

Para facilitar vamos somente mudar as coisas de lugar:

 {1\over{Vel}} =  {{Con} \over {Vmax \cdot Con}} + {{Km} \over {Vmax \cdot Con}}

E vamos cortar aquele Con dividido por Con

 {1\over{Vel}} =  {1 \over {Vmax}} + {{Km} \over {Vmax \cdot Con}}

E do outro lado o que a gente pode ver? Bem é a razão de Km por Vmax que multiplica o inverso da Concentração.

 {1\over{Vel}} =  {1 \over {Vmax}} + {{Km} \over {Vmax}}\cdot{1\over{Con}}

Agora com o que isso te parace, imagine que 1/con é um número x qualquer, num é uma coisa + outra coisa que multiplica x?
Isso não é a+bx, ou intercepto + inclinação?
Vamos substituir ali o intercepto e inclinação por alpha e beta para a gente ver melhor.

 {1\over{Vel}} =  \alpha + \beta\cdot{1\over{Con}}

Fechou, é uma equação da reta, simples assim.
Agora a gente sabe o que fazer com os coeficientes la em cima.
So pensar pelo que substituimos aqui o alpha(\alpha) e o beta(\beta).

O alpha é:
\alpha={1 \over {Vmax}}

E o beta é:
\beta= {{Km} \over {Vmax}}

Então ali no modelo glm a gente tem que fazer o que para ver o Km e o Vmax?

Vmax vai ser:

Vmax={1 \over {\alpha}}

E Km vai ser:

Km = \beta \cdot {Vmax}

Que ja temos, mas mas podemos pensar também como

Km = \beta \cdot {\alpha^{-1}}

Agora no R é so a gente pegar os valores de coeficientes do modelo glm, com a função coef e fazer essas continhas, que para Vmax e Km fica:

coef(glm.menten)[1]^-1
coef(glm.menten)[2]*(coef(glm.menten)[1]^-1)

Isso da:

> coef(glm.menten)[1]^-1 (Intercept) 3.481847 > coef(glm.menten)[2]*(coef(glm.menten)[1]^-1) c(Con^-1) 0.4428143

Voala, quase como mágica, os mesmos valores estimados com o modelo não linear, usando nls.
a gente pode ainda agora ver os valores, usando o inverso, que entendemos o porque disso agora, e ver como fica os dois modelos ajustados com nls e glm aos dados.

05

E são a mesma coisa.
Muito legal né. Mas o que eu queria mostrar de verdade é como as vezes vemos uma pessoa trabalhando de um jeito, outra trabalhando de outro, mas as vezes as pessoas estão fazendo a mesma coisa. Nesse caso a gente ajustou o mesmo modelo, com uma função não linear, e depois com uma função linear que a gente conhece bem, o a+bx (conhecida tambem como \alpha +\beta \cdot x), só que é preciso pensar um pouco para visualizar onde foram parar os parâmetros que temos interesse no caso da função linear. Mas de qualquer forma ajustamos exatamente a mesma coisa. Fizemos a mesma coisa de duas formas completamente diferentes, para chegar ao mesmo resultado..

Eu acho um pouco legal pensar também, como algumas vezes, a gente pode estar advogando um modelo contra o outro, e eles não podem ser algo assim, simplesmente a mesma coisa escrita de outra forma.

Foi mais uma curiosidade, mas serviu para ver também como uma função de ligação, como simplesmente inverter uma função faz o milagre do enrretamento das coisas em alguns casos.

Eu vi essa curiosidade numa apostilinha de glm, já que nunca teria essa sacada na vida hehehe, no fim do texto eu deixo o link e a referência, mas acho que é apenas um curso de R sobre GLM, nada que de para citar, mas é legal a apostila.

Bem o script vai estar ainda no final aqui, mas eu vou deixar ele la no repositório Recologia do Github, junto com os outros, fica um lugar fácil para achar somente o script, como as vezes eu mesmo procuro coisas minhas para copiar ou lembrar comandos, fica mais fácil de achar tudo.

Bem é isso, até o próximo post.

Referência:

Ulrich Halekoh, Søren Højsgaard 2008 Generalized Linear Models (GLM) 32 pp.

#Gerando dados
set.seed(21)
m_menten <- function(C,Vmax,Km) { return( (Vmax * C)/( Km + C) ) }
con<-runif(30,0,4)
v<-m_menten(con,Vmax=3.4,Km=0.4)
v<-v+rnorm(30,0,0.08)
dados<-data.frame(Con=con,Vel=v)
 
#Figura 1, Caracteristicas do modelo
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,4.5),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
abline(h=3.4,lty=2,col="red")
text(4,3.5,"Vmax")
abline(h=0.5*(3.4),lty=2,col="red")
text(4,(0.5*(3.4))-0.1,expression("1/2 Vmax"))
lines(x=c(0.4,0.4),y=c(0,0.5*(3.4)),lty=2,col="red")
text(0.55,1,"Km")
 
#Ajustando um modelo não linear para os dados
nls.menten<-nls(Vel~(Vmax * Con)/(Km + Con),data=dados,start=list(Vmax=1,Km=1),algorithm="default")
summary(nls.menten)
 
#Figura 2 - ajuste modelo não linear
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,6),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
points(seq(0,6,0.01),predict(nls.menten,newdata=data.frame(Con=seq(0,6,0.01))),type="l",lty=2,col="blue")
legend("topleft",lty=2,col="blue",legend="Ajuste com nls()",bty="n")
 
#Figura 3.
plot(c(1/Vel)~c(1/Con),frame=F,ylim=c(0,1),xlim=c(0,6),ylab="1/Velocidade",xlab="1/Concentração",pch=19,data=dados)
 
#Ajustando um modelo Geral Linear
glm.menten<-glm(Vel ~c(Con^-1), data = dados,family = gaussian(link = "inverse"))
summary(glm.menten)
 
#Retornando os valores de Vmax e Km
coef(glm.menten)[1]^-1
coef(glm.menten)[2]*(coef(glm.menten)[1]^-1)
 
#Figura 4 - Ajuste com modelo linear
plot(c(1/Vel)~c(1/Con),frame=F,ylim=c(0,1),xlim=c(0,6),ylab="1/Velocidade",xlab="1/Concentração",pch=19,data=dados)
points(c(seq(0.1,6,0.01)),predict(glm.menten,newdata=data.frame(Con=c(1/seq(0.1,6,0.01)))),type="l",lty=2,col="green")
legend("topleft",lty=2,col="green",legend="Ajuste com glm()",bty="n")
 
 
#Figura 5 - Comparando os dois ajustes
plot(Vel~Con,frame=F,ylim=c(0,4),xlim=c(0,6),ylab="Velocidade",xlab="Concentração",pch=19,data=dados)
points(seq(0,6,0.01),predict(nls.menten,newdata=data.frame(Con=seq(0,6,0.01))),type="l",lty=2,col="blue")
points(c(seq(0.1,6,0.01))^-1,predict(glm.menten,newdata=data.frame(Con=c(1/seq(0.1,6,0.01))))^-1,type="l",lty=2,col="green")
legend("topleft",lty=2,col=c("blue","green"),legend=c("Ajuste com nls","Ajuste com glm"),bty="n")

Estabilidade e tempo de retorno a uma pertubação para populações.

Vamos ver mais uma coisa ainda sobre o crescimento logístico. Ja vimos varias características legais aqui e aqui.

Certo mas vamos olhar denovo aquele gráfico da dinâmica populacional ao longo do tempo segundo a equação de crescimento logístico, agora olhando por um pouco mais de tempo.

01

Certo então a tendencia é ela permanecer estável. Como se ouve-se um ponto um tamanho populacional que sempre que a população caminha no tempo, ela caminha para esse ponto. Nesse caso o valor de estabilidade é equivalente a capacidade suporte, ou k ou 1 \over {\alpha}, representações para todos os gostos, usa a letra que quiser.

Mas na verdade existe 2 pontos de atração aqui, um estável e um instável. O que é um ponto de atração instável, é por exemplo quando a população chega a zero, se ela chegar a zero, a população não vai mais crescer para a capacidade suporte, porque não ha mais nenhum indivíduo para ter filhotinhos e fazer a população crescer ir para a capacidade suporte. Mas nesse caso, qualquer pertubação, por exemplo colocar um indivíduo na população, uma migração, que coloniza essa população novamente, vai fazer ela começar a crescer denovo, um reintrodução de indivíduos num local onde a espécie sumiu, horas, aquele lugar ainda pode sustentar a espécie, então ela começa a subir denovo, então quando a população chegou a zero, ela ficara parada no zero a não ser que haja uma pertubação do tipo que reintroduz indivíduos ali.

O que nos poe devolta a essa figura:

02

O que estou dizendo é que, a taxa de crescimento, que no gráfico diz para que lado o tamanho populacional esta indo, tem dois pontos onde é exatamente zero, a e d, a é quando a população esta em zero, e fica ali estacionado até haver uma pertubação, por exemplo a reintrodução de indivíduos da espécie, e o ponto d, que é quando a população chegou a capacidade suporte, ai ela não cresce nem diminui. Se a população esta acima da capacidade suporte, ela diminui, olhe o ponto e, a taxa de crescimento fica negativa. Se ela esta menor, a população cresce, com velocidade que depende do N atual.
Lembrando que aqui estamos olhando uma população que tende a uma estabilidade, com taxa de crescimento de 1, ja vimos que isso é mais confuso quando a taxa de crescimento é muito grande aqui.

Mas agora a gente pode fazer a seguinte pergunta, quanto tempo demora, ou qual a taxa de retorno de uma população ao ser pertubada, pertubada uma certa quantidade, retirar ou reintroduzir indivíduos.

Podemos sem fazer contas, bolar um exemplo simples para pegar um “feeling” da coisa.
A pergunta é a seguinte, suponha que temos uma população como a la de cima, do primeiro gráfico, com o r de 1. Se a gente for la e tirar 20 indivíduos, quando ela ta parado na capacidade suporte de 100, quanto tempo demora para ela voltar para 100 denovo? Ou de forma mais simples, voltar para a capacidade suporte?
Talvez podemos tentar comparar mentalmente com uma população com crescimento mais lento, r=0.5, bem como a taxa de crescimento é mais devagar, com r=0.5, ela vai demorar mais para voltar a 100 correto? Logo o tempo de retorno depende da taxa de crescimento.
Graficamente é isso que estou falando.

03

Vejam como a população com r=0.5, verdinha, ela cresce mais devagar, porque a taxa de crescimento dela é menor, mas ela também, após ser pertubada, volta mais devagar.
Agora é importante notar que essa idéia funciona da mesma forma para o caso da população exceder a capacidade suporte do ambiente.
So que ao invés de estar subindo, ela esta descendo, assim:

04

Agora vem o problema, na verdade o tempo de retorno pode ser estabelecido usando calculo para calcular a inclinação da curva, ja que a inclinação da curva é uma derivada. Nesse caso, nos precisamos derivar a taxa de crescimento (dN/dt) com respeito a variável no eixo X, o tamanho da população N. Essa derivada é chamada de parcial, “partial derivative“, já que só queremos saber a inclinação em relação a taxa de crescimento e não ao tamanho populacional N.

Para achar a derivada parcial, primeiro nos olhamos a taxa de crescimento como:

  \dot N  = r \cdot N \cdot ( 1 - \alpha  \cdot N )

Esse ponto em cima do N é comum para derivadas em relação ao tempo, como taxa de crescimento, eu num faço idéia por que, por enquanto eu so acredito porque funciona.

Sua derivada parcial em relação ao tamanho populacional é:

{ \partial  \dot N \over \partial  N } = r - 2 \cdot r \cdot \alpha \cdot N

Como a gente ta interessado no equilíbrio e o equilíbrio é o momento em que  N  = {1 \over \alpha}, quando atingiu a capacidade suporte, a gente pode substituir isso ali em cima e ficamos com:

{ \partial  \dot N \over \partial  N } = r - 2 \cdot r \cdot \alpha \cdot {1 \over \alpha}

Ai como tem \alpha \cdot {1 \over \alpha}, que da 1.

{ \partial  \dot N \over \partial  N } = r - 2 \cdot r \cdot 1

r menos dois r da menos um r, qualquer coisa que multiplica 1 é ela mesma,

{ \partial  \dot N \over \partial  N } = - r

Taram, o valor de retorno é inversamente proporcional a taxa de crescimento r. Quanto maior a taxa, mais rápido o retorno, menor é esse tempo, o que já tinha pensado la em cima.

Isso quer dizer que perto da capacidade suporte, a taxa de crescimento é -rx, onde x é o tamanho da pertubação, o quão distante a população foi para em relação ao seu ponto de equilíbrio.

{ dx \over dt } = - r \cdot x

E podemos ver isso num gráfico, mudança da taxa de crescimento de acordo com a pertubação,

05

Mas o retorno em si sera exponencial, ainda uma curva, tanto numa pertubação quando adicionamos indivíduos como quando tiramos.

06

Bem isso pode parecer bem complicado, para mim parece, principalmente quando a gente inventa de enfiar calculo em biologia, mas responde a uma pergunta bem prática. O que acontece após uma pertubação. Isso gera várias outras respostas também muito úteis, como qual o impacto de uma pertubação x, sendo que podemos quantificar o tamanho dessa pertubação e já prever o que vai acontecer com a população, ou pelo menos esperar o que vai acontecer, ter um panorama na melhor das hipóteses ou pior das hipóteses.
Claro que isso, de uma forma bem simplista, sem levar outras coisas em consideração, como interação entre espécies.
Mas como visto nas teorias de rendimento máximo sustentável de uma população, valores de retorno, mudança na taxa de crescimento, essas coisas são muito importantes se queremos explorar o ambiente sem destruir ele e condenar a nos mesmos a ausência certos recursos.

O script esta no repositório do Recologia no Github além do fim do post aqui.
E sempre fica aquele alerta, se eu escrevi alguma coisa errado ali em cima, confira tudo antes de acreditar em algo, mas se achar o erro me avise num comentário ou telegrafo, sei la 🙂
Até a próxima.

###################################################
### Valores de Retorno
###################################################
#install.packages("primer")
rm(list=ls())
library(primer)
 
#Crescimento Logístico
dlogistic <- function(alpha=0.01, rd=1, N0=2, t=15) {
  N <- c(N0, numeric(t))
  for(i in 1:t) {
    N[i+1] <- N[i] + rd * N[i] * (1 - alpha * N[i])
  }
  return(N)
}
 
#Usando a função com parametros default que definimos
Nts <- dlogistic(,t=40)
 
#Observando num grafico
t <- 40
a <- 0.01
plot(0:t, Nts,frame=F,xlab="Tempo",ylab="Tamanho Populacional",type="b",
     xlim=c(0,45),ylim=c(0,120),pch=19,cex=0.5)
abline(h=1/a,lty=3,col="red")
 
 
###################################################
# Estabilidade
###################################################
pop.growth.rate <- expression( r * N * (1 - alpha * N) )
r <- 1
alpha <- 0.01
N <- 0:120
plot(N, eval(pop.growth.rate), type="l",frame=F,
     ylab="Taxa de crescimento populacional (dN/dt)", xlab="N")
abline(h=0)
legend('topright', "r=1", lty=1,bty="n")
N <- c(0, 10, 50, 100, 115)
points(N, eval(pop.growth.rate), cex=1.5)
text(N, eval(pop.growth.rate), letters[1:5],adj=c(.5,2))
arrows(20,2, 80,2, length=.1, lwd=3)
arrows(122,-2,109,-2, length=.1, lwd=3)
 
###################################################
# Exemplo de pertubação
###################################################
pop1 <- dlogistic(alpha=0.01,N0=2,rd=0.5,t=20)
pop1 <- c(pop1,dlogistic(alpha=0.01,N0=(pop1[20]-20),rd=0.5,t=20))
 
pop2 <- dlogistic(alpha=0.01,N0=2,rd=1,t=20)
pop2 <- c(pop2,dlogistic(alpha=0.01,N0=(pop2[20]-20),rd=1,t=20))
 
 
#Observando num grafico
plot(1:length(pop1), pop1,frame=F,xlab="Tempo",ylab="Tamanho Populacional",type="b",
     xlim=c(0,45),ylim=c(0,120),pch=19,cex=0.5,col="green")
points(1:length(pop2), pop2,type="b",pch=19,cex=0.5,col="blue")
abline(h=1/0.01,lty=3,col="red")
arrows(22,50, 22,75, length=.1, lwd=3)
text(22,45,"Pertubação")
legend('topright', c("r=0.5", "r=1"),col=c("green","blue"),pch=19,lty=1,bty="n")
 
##########
 
pop1 <- dlogistic(alpha=0.01,N0=2,rd=0.5,t=20)
pop1 <- c(pop1,dlogistic(alpha=0.01,N0=(pop1[20]+20),rd=0.5,t=20))
 
pop2 <- dlogistic(alpha=0.01,N0=2,rd=1,t=20)
pop2 <- c(pop2,dlogistic(alpha=0.01,N0=(pop2[20]+20),rd=1,t=20))
 
 
#Observando num grafico
plot(1:length(pop1), pop1,frame=F,xlab="Tempo",ylab="Tamanho Populacional",type="b",
     xlim=c(0,45),ylim=c(0,120),pch=19,cex=0.5,col="green")
points(1:length(pop2), pop2,type="b",pch=19,cex=0.5,col="blue")
abline(h=1/0.01,lty=3,col="red")
arrows(22,70, 22,95, length=.1, lwd=3)
text(22,65,"Pertubação")
legend('topright', c("r=0.5", "r=1"),col=c("green","blue"),pch=19,lty=1,bty="n")
 
###################################################
# Retorno
###################################################
mi <- 99.4
ma <- 100.6
N <- seq(mi,ma, length=30)
plot(N, eval(pop.growth.rate), type="l",frame=F,
     ylab="Taxa de crescimento populacional (dN/dt)", xlab="N")
abline(h=0)
r <- .5
lines(N, eval(pop.growth.rate), lty=2); r <- 1
segments(c(100-.3, 100+.3), c(-.05,-.05), c(100-.3,100+.3), c(.05,.05))
 
text(99.7,-.15, quote("N*"-"x"))
text(100.3,.15, quote("N*"+"x"))
points(100, 0)
text(100, -.1, "N*")
 
legend('topright', c("r=1", "r=0.5"), lty=1:2,bty="n")
###################################################
# Pertubação
###################################################
par( mgp = c( 2, .75, 0 ) )
x <- 1
plot(c(-1,8), 1.2*c(-x,x), type='n',
     ylab="N", xlab='', axes=FALSE)
axis(2, at=c(-x,0,x), labels=c(expression("N*"-"x"),"N*", expression("N*"+"x")) )
abline(h=0)
t <- 1
points(c(0, t), c(1, x*exp(-1*t) ), pch=c(19,1))
curve(1*exp(-1*x), 0, .98, add=T, lty=3)
arrows(0,-x-.2, 7,-x-.2)
mtext("Tempo", side=1)
text(0, 0, "Pertubação", srt=90, adj=c(-.1,-.1) )
arrows(0,0,0,.8, length=.1)
text(3, .7, expression("Recuperação na taxa de"=="e"^"-r"))
 
points(c(0+2, t+2), c(-1, -x*exp(-1*t) ), pch=c(19,1))
x2 <- seq(0,.98, length=30)
y <- 0 - 1 * exp(-1*x2)
lines(x2+2, y, lty=3)
text(2, 0, "Pertubação", srt=90, adj=c(1.02,-.1) )
arrows(2, 0, 2,-.8, length=.1)
text(5, -.7, expression("Recuperação na taxa de"=="e"^"-r"))

Problema – Rosalind – Inferring mRNA from Protein – MRNA

Ja houve um problema, esse aqui, onde tinhamos o proposito de, a partir de uma sequencia de RNA, prever qual proteína seria codificada.

Mas e o inverso disso? Isso é o que esse problema diz respeito, a partir de uma proteína, mais precisamente a sequência de aminoácidos, prever qual era o RNA original.

Aqui entra uma observação interessante também, já que uma sequência de RNA codifica uma proteína, mas uma proteína pode ser sintetizada por diferentes sequencias de RNA.

Essa é a pergunta especificamente.

Então uma sequência de aminoácidos, como no exemplo:

MA

Se olharmos esse problema anterior, vemos que quase sempre existem mais de uma possibilidade para cada aminoácido.

Para a sequencia acima, “M” é a abreviação de Metionina, em inglês Methionine, abreviado como “Met” ou a letra “M” sozinha. Ele é codificado pela sequencia “AUG” de RNA.

Até aqui é fácil, agora o “A” vem de Alanina, em inglês Alanine, abreviado como “Ala” ou “A” nesse caso. Ele pode ser codificado por GCU, GCC, GCA ou GCG. Isso nos da 4 possibilidades.

Logo

AUG GCU 1 AUG GCC 2 AUG GCA 3 AUG GCG 4

Todos podem ter gerado a sequencia MA de aminoácidos.
Mas pera ai, além dos aminoácidos, existe o codon que diz para o processo parar, que é UAA, UGA, UAG.

Ou seja, aumentamos em as possibilidades para:

AUG GCU UAA 1 AUG GCC UAA 2 AUG GCA UAA 3 AUG GCG UAA 4 AUG GCU UGA 5 AUG GCC UGA 6 AUG GCA UGA 7 AUG GCG UGA 8 AUG GCU UAG 9 AUG GCC UAG 10 AUG GCA UAG 11 AUG GCG UAG 12

Doze possibilidades, que é 1*4*3, o número de possibilidades do M vezes o número de possibilidades de A vezes o número de possibilidades de fim, que da 12.

Assim o problema fica fácil agora.
So temos que reorganizar as possibilidades por aminoácido, quantas possibilidades de codon cada um tem, por exemplo usando dicionários temos:

{'A': 4, 'M': 1, 'C': 2, 'E': 2, 'D': 2, 'G': 4, 'F': 2, 'I': 3, 'H': 2, 'K': 2, 'Stop': 3, 'L': 6, 'N': 2, 'Q': 2, 'P': 4, 'S': 6, 'R': 6, 'T': 4, 'W': 1, 'V': 4, 'Y': 2}

Agora é só entrar com um string, ai para cada letra, a gente pega a o número respectivo e multiplica, não esquecendo de multiplicar por 3, que é o número de possibilidades de fim de sequência de aminoácidos.

A tabela inicial a gente usa a mesma do ultimo exercício, apresentado la em cima.

Um último passo é que o resultado final, deve ser apresentado como o resto da divisão por um milhão, acho que para não ficar um número muito grande.
Não entendi exatamente se existe algo especial quanto ao modulo de 1 milhão.

Mas é isso.
Esse foi o post estreando o meu repositório do rosalind 🙂

RNA_CODON_TABLE = {
    'UUU': 'F',     'CUU': 'L',     'AUU': 'I',     'GUU': 'V',
    'UUC': 'F',     'CUC': 'L',     'AUC': 'I',     'GUC': 'V',
    'UUA': 'L',     'CUA': 'L',     'AUA': 'I',     'GUA': 'V',
    'UUG': 'L',     'CUG': 'L',     'AUG': 'M',     'GUG': 'V',
    'UCU': 'S',     'CCU': 'P',     'ACU': 'T',     'GCU': 'A',
    'UCC': 'S',     'CCC': 'P',     'ACC': 'T',     'GCC': 'A',
    'UCA': 'S',     'CCA': 'P',     'ACA': 'T',     'GCA': 'A',
    'UCG': 'S',     'CCG': 'P',     'ACG': 'T',     'GCG': 'A',
    'UAU': 'Y',     'CAU': 'H',     'AAU': 'N',     'GAU': 'D',
    'UAC': 'Y',     'CAC': 'H',     'AAC': 'N',     'GAC': 'D',
    'UAA': 'Stop',  'CAA': 'Q',     'AAA': 'K',     'GAA': 'E',
    'UAG': 'Stop',  'CAG': 'Q',     'AAG': 'K',     'GAG': 'E',
    'UGU': 'C',     'CGU': 'R',     'AGU': 'S',     'GGU': 'G',
    'UGC': 'C',     'CGC': 'R',     'AGC': 'S',     'GGC': 'G',
    'UGA': 'Stop',  'CGA': 'R',     'AGA': 'R',     'GGA': 'G',
    'UGG': 'W',     'CGG': 'R',     'AGG': 'R',     'GGG': 'G'
}
 
def freq_codon():
    frequencias = {}
    for k, v in RNA_CODON_TABLE.items():
        if not frequencias.has_key(v):
            frequencias[v] = 0
        frequencias[v] += 1
    return frequencias
 
 
def possibilidades_rna(seq_aminoacidos):
    f = freq_codon()
    n = f['Stop']
 
    for c in seq_aminoacidos:
        n *= f[c]
 
    return n
 
 
small_dataset = "MA"
#large_dataset = open('/home/augusto/Downloads/rosalind_mrna.txt').read().strip()
 
print possibilidades_rna(small_dataset) % 1000000