Continuação do post anterior que comecei a escrever como um tutorial para os meus alunos que iniciam projetos de iniciação científica e ainda não estão familiarizados com esse processo de calibrar os parâmetros de modelos matemáticos. Caso não tenha lido volte uma casa (leia o post anterior aqui) e depois continue.
Para dar continuidade, vamos pegar um exemplo simples de problema considerando apenas a variação no tempo. Vamos supor que existe um determinado fenômeno que temos interesse em estudar. O primeiro passo é identificar se já existe um modelo matemático na literatura que represente esse fenômeno.
Um exemplo interessante e suficientemente simples é o modelo matemático utilizado para representar epidemias que foi proposto por Kermack e MacEndrick chamado SIR (talvez você tenha visto uma variação dele sendo usada para previsões durante a pandemia Covid-19). As letras S, I e R significam respectivamente:
Suscetíveis
Infectados
Recuperados
Esse modelo representa, portanto, três populações de indivíduos e suas dinâmicas variando com relação ao tempo. Para entender o modelo, vamos pensar que inicialmente, quando não existe ainda a epidemia e surge um novo agente, todas as pessoas podem ser suscetíveis a se contaminar. As pessoas se tornam infectadas ao entrarem em contato com esse determinado agente. Com o passar do tempo, vamos supor que os indivíduos infectados se recuperam. (Existem variações desse modelo para considerar outros tipos de epidemias em que a pessoa possa ter um tempo de incubação, pode não se recuperar etc. Mas vamos considerar apenas essa 3 populações para simplificar o exemplo).
Para chegar nas equações do modelo, vamos precisar de um termo que representa a taxa em que os indivíduos se tornam infectados: um parâmetro (beta) multiplicado a variável S que representa a velocidade em que um individuo deixa de ser suscetível e passa a ser infectado e essa taxa é proporcional ao valor da população S em determinado momento. E com isso a população S depende do tempo, o que representamos como S(t).
beta * S(t)
Como esta quantidade de indivíduos está “saindo” de uma população e “aparecendo” em outra, esse termo aparece com um sinal negativo (subtraindo) na equação dos Suscetíveis e o mesmo termo aparece com sinal positivo (somando) na equação dos Infectados. A notação dX/dt é adotada para representar equações diferenciais em que uma variável X(t) depende do tempo.
dS/dt = – beta * S(t) // equação dos Suscetíveis
dI/dt = + beta * S(t) // equação dos Infectados
Considerando apenas essas duas equações, se escolhemos valores para S(0) e I(0), que significa escolher as condições iniciais no tempo igual a zero, podemos resolver esse sistema. Simplificando bastante, da pra imaginar que a população de Suscetíveis vai sempre diminuir e a população de Infectados vai sempre aumentar, certo? Consegue imaginar isso? Olha o gráfico abaixo:
Bom, mas o indivíduo infectado pode se recuperar também né? Isso ocorre em uma determinada velocidade e é representado por outro termo em que um parâmetro diferente (gamma) é multiplicado a variável I (infectados) e aparece com sinal negativo na equação de Infectados e sinal positivo na equação R (recuperados), representando que a taxa de recuperação é proporcional, ou depende da quantidade de pessoas infectadas. Aí finalmente chegamos ao modelo SIR:
| dS/dt = – beta * S(t) | // equação dos Suscetíveis |
| dI/dt = + beta * S(t) – gamma * I(t) | // equação dos Infectados |
| dR/dt = + gamma * I(t) | // equação dos Recuperados |
Esse é o modelo epidemiológico mais simplificado e usaremos esse para fins didáticos. Consideraremos apenas as três variáveis (S, I e R) e apenas dois parâmetros (beta e gamma).
Para utilizar esse modelo para representar uma epidemia precisamos saber qual a população inicial de indivíduos suscetíveis, infectados e recuperados, o que chamamos de “condição inicial”. E precisamos saber também quais são as taxas de infecção e recuperação.
Vamos supor que temos as seguintes condições iniciais (valor das variáveis no início da simulação, no tempo igual a zero) para cada variável do modelo:
S(0) = 10000 // população de 10000 indivíduos que podem se contaminar
I(0) = 0 // ninguém foi infectado
R(0) = 0 // como não tem ninguém infectado, também não tem ninguém recuperado ainda
E vamos supor também os seguintes valores para os parâmetros:
beta = 0.1 (taxa de infecção)
gamma = 0.01 (taxa de recuperação)
Com isso, precisamos ainda definir qual o algoritmo numérico que será utilizado para resolver o sistema de 3 equações, qual o passo de tempo que utilizaremos e quanto tempo de simulação desejamos representar. Algumas horas? 1 semana? 1 mês? 6 meses? etc.
Aqui vou assumir que quem está lendo já entende alguma coisa sobre algoritmos numéricos e vamos utilizar a função de integração solve_ivp() da biblioteca scipy.integrate (lembrando que estamos usando a linguagem de programação Python).
Uma observação antes de continuar é que passos de tempo menores permitem uma melhor aproximação do resultado. Imagine que estamos avançando em direção a uma solução. Se dermos passos muito grandes e por algum motivo erramos um pouco a direção, fica difícil retomar o caminho sem dar um salto grande para corrigir. Mas se estamos caminhando dando passos bem pequenos, podemos fazer ajustes na trajetória sem precisar dar grandes saltos de forma que quem esta olhando de fora talvez nem perceba que chegamos a sair ligeiramente da rota correta. No caso da função solve_ivp() ela utiliza como padrão o método Runge-Kutta e o passo de tempo é determinado de acordo com o método.
Vamos supor que simularemos o que acontece no período de 100 dias, e que queremos dar passos de um dia. Nesse caso, calcularemos um valor para cada dia a partir da resolução do sistema considerando as condições iniciais acima (link para o código-fonte no final da postagem). O gráfico fica assim:
Talvez tenha visto gráficos parecidos com esse durante a pandemia… a linha laranja está representando os infectados e durante a pandemia modelos matemáticos semelhantes a esse foram bastante utilizados para prever os picos e propor soluções para diminuir ou “achatar” essa curva. Mais a frente vou deixar um link para o código no Google Colab para que possa criar uma copia e alterar os parâmetros e fazer novas simulações. Que mudanças no modelo fazem achatar a curva?
Nesse post vimos um exemplo de sistema de equações diferenciais e um exemplo de resolução desse sistema. Esse modelo representa as dinâmicas de como uma doença infecciosa se espalha. Mas como sabemos quais valores de parâmetros devemos utilizar no modelo? Esse é o assunto do próximo post.
Acesse aqui o código-fonte utilizado nesse exemplo.
Caso tenha dúvidas ou sugestões pode deixar um comentário!
Até o próximo post! 🙂
Barbaridades 