Aula 7: Equação Diferencial Ordinária (ODE)

MSP4092: Biomatemática: aspectos quantitativos da vida

suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(deSolve, warn.conflicts=FALSE))

Adicionar

• 5.2. THE SIR MODEL OF INFECTIOUS DISEASE, p. 360:
  • Explorations of mathematical models in biology with MATLAB - Shahin - 2014.pdf

Material

• HTML de Rmarkdown em RPubs
• Prof. José Siqueira: ResearchGate

Pensamento

“Significação das teorias físicas: Existem muitas pessoas que se espantam ao ver quão efêmeras são as teorias científicas. Essas pessoas veem teorias serem sucessivamente abandonadas depois de alguns anos de prosperidade; veem ruínas se acumulando sobre ruínas; preveem que teorias hoje na moda sucumbirão, por sua vez, dentro em pouco, e concluem que elas são inteiramente vãs. É o que chamam a falência da ciência. […] O objetivo das teorias matemáticas não é nos revelar a verdadeira natureza das coisas – o que não seria uma pretensão razoável – e sim coordenar as leis físicas que a experiência nos revela mas que, sem o auxílio da matemática, não poderíamos sequer enunciar. Pouco nos importa que o éter exista realmente: é um problema para os metafísicos (filósofos). O importante para nós é que tudo se passa como se ele existisse, e essa é uma hipótese cômoda para a explicação de fenômenos. Afinal, temos outras razões para crer na existência de objetos materiais? Essa também é uma hipótese cômoda e que nunca deixará de o ser, ao passo que um dia, certamente, o éter será rejeitado, por inútil.”

POINCARÉ, H.-J. A ciência e a hipótese. Brasília: Universidade de Brasília, 1984, p. 127 e 157.

Conteúdo

1. Número real (Capítulo 1)
2. Função e relação (Capítulo 3)
3. Funções potência, periódica, exponencial e logarítmica (Capítulos 4, 5 e 6)
4. Método gráfico (Capítulo 7)
5. Série e limite (Capítulo 8)
6. Derivada e integral (Capítulos 9 e 10)

7. Equação diferencial ordinária (ODE) (Capítulo 11)

8. Função de duas ou mais variáveis independentes (Capítulo 12)
9. Probabilidade (Capítulo 13)
10. Matriz, vetor e número complexo (Capítulos 14 e 15)

Introdução

• Differential Equations: Paul’s Online Notes

O primeiro texto científico que resolveu uma equação diferencial foi escrito por Isaac Newton, em sua obra intitulada “Mathematical Principles of Natural Philosophy” (ou “Princípios Matemáticos da Filosofia Natural”, em tradução livre). O livro foi publicado em 1687 e é comumente referido como “Principia”.

No “Principia”, Newton formulou as leis do movimento e a lei da gravitação universal, além de desenvolver os fundamentos do cálculo diferencial e integral. Na segunda parte do livro, Newton introduziu uma série de métodos matemáticos para resolver problemas físicos, incluindo a resolução de equações diferenciais.

O trabalho de Newton foi um marco importante no desenvolvimento do cálculo e da física, estabelecendo as bases para a resolução de equações diferenciais e a compreensão dos fenômenos naturais. Desde então, equações diferenciais têm sido uma ferramenta essencial em muitas áreas da ciência e da engenharia.

Uma equação diferencial é uma equação que envolve uma função desconhecida e uma ou mais de suas derivadas. Ela descreve a relação entre a função desconhecida e suas derivadas em termos das variáveis independentes.

Existem diferentes tipos de equações diferenciais, como equações diferenciais ordinárias (EDO, ODE) e equações diferenciais parciais (EDP, PDE). As EDO envolvem apenas derivadas em relação a uma única variável independente, enquanto as EDP envolvem derivadas parciais em relação a duas ou mais variáveis independentes.

A solução de uma equação diferencial depende de sua forma específica. Existem métodos analíticos e numéricos para resolver equações diferenciais. Os métodos analíticos geralmente envolvem técnicas matemáticas, como integração, fator integrante, transformadas de Laplace e Fourier, entre outros.

A solução particular de uma equação diferencial é a solução específica que satisfaz as condições iniciais ou de contorno especificadas no problema. As condições iniciais são valores dados da função desconhecida e suas derivadas em um ponto inicial, enquanto as condições de contorno são valores dados em pontos específicos do domínio.

Cada equação diferencial terá sua própria forma e método de solução específico.

Movimento sob força constante pela segunda lei de Newton

• Newton’s laws of motion: Wikipedia

A segunda lei de Newton afirma que a força resultante que atua sobre um corpo é igual ao produto de sua massa pela aceleração. Se \(x(t)\) representa a posição do corpo no instante \(t\), então sua aceleração instantânea é dada por \(x^{\prime\prime}(t)\)

A segunda lei de Newton pode ser escrita em forma vetorial como

\[ \vec{F}=m\vec{a} \]

No movimento unidimensional ao longo do eixo \(x\), a aceleração é a segunda derivada da posição em relação ao tempo:

\[ \vec{a}(t)=x^{\prime\prime}(t)\,\vec{i} \] Sendo que o vetor \(\vec{i}\) é o vetor unitário associado ao eixo horizontal \(x\). Ele possui módulo igual a \(1\) e aponta no sentido positivo do eixo \(x\).

Assim,

\[ \vec{F}(t)=m x^{\prime\prime}(t)\,\vec{i} \]

Tomando o módulo dos dois lados,

\[ |\vec{F}(t)|= \left|m x^{\prime\prime}(t)\,\vec{i}\right| \]

Como \(m>0\) e \(|\vec{i}|=1\), resulta

\[ |\vec{F}(t)|=m|x^{\prime\prime}(t)| \]

Portanto, no caso unidimensional, quando se trabalha apenas com o módulo da força e da aceleração, a segunda lei pode ser escrita como

\[ F=m|x^{\prime\prime}(t)| \]

Se o sinal do movimento for preservado, usa-se a forma escalar com sinal:

\[ F=m x^{\prime\prime}(t) \]

Logo, a derivada entra na aceleração porque a aceleração é a taxa de variação da velocidade, e a velocidade é a taxa de variação da posição:

\[ v(t)=x^{\prime}(t), \qquad a(t)=v^{\prime}(t)=x^{\prime\prime}(t) \] A unidade de medida da força é Newton que corresponde a kg m s^-2.

Diagrama de força sobre bloco

No caso particular em que a força é constante, isto é, \(F=k>0\), com \(k\) constante, obtemos a equação diferencial ordinária linear

\[ m x^{\prime\prime}(t)=k \]

Supondo as condições de contorno

\[ x(0)=0 \qquad \text{e} \qquad x(\tau)=d>0 \]

temos o sistema de equações

\[ \begin{cases} m x^{\prime\prime}(t)=k\\ x(0)=0\\ x(\tau)=d \end{cases} \]

Reescrevendo:

\[ x^{\prime\prime}(t)=\dfrac{k}{m} \]

Problema: Qual é a função \(x(t)\) que derivada duas vezes resulta um valor constante?

Resposta: \(x(t)=c_1 t^2 + c_2 t + c_3\): \(x^{\prime\prime}(t)=2c_1\)

D[c_1 t^2 + c_2 t + c_3, t, t]

A solução desse sistema é

\[ x(t)=\dfrac{k}{2m}t(t-\tau)+\dfrac{d}{\tau}t \] ou

\[ x(t)=\frac{k}{2m}t^2+\left(\dfrac{d}{\tau}-\frac{\tau k}{2m}\right)t \] A solução da EDO tem a seguinte estrutura:

\[x(t)=\dfrac{d}{2} t^2+ v_0 t+x_0\] Sendo que \(a\) é a aceleração constante, \(\nu\) é uma velocidade inicial e \(x_0=x(t=0)=0\) é a posição inicial.

A trajetória do corpo é quadrática no tempo.

A trajetória é aproximadamente linear se \(k\) ;e pequeno; ou \(m\) ;e grande; ou o intervalo de observação \(t\) é pequeno; ou a velocidade inicial \(v_0\) domina o efeito acelerativo.

Se \(k=0\), a solução é exatamente linear:

\[ x(t)=\dfrac{d}{\tau}t+x_0 \]

A velocidade é aproximadamente constante quando o termo dependente de \(t\) é pequeno em relação à velocidade inicial \(\frac{d}{\tau} - \frac{\tau k}{2m}\).

A primeira derivada representa a velocidade instantânea:

\[ x'(t) = \frac{k}{m}t + \left( \frac{d}{\tau} - \frac{\tau k}{2m} \right) = v(t) \]

A velocidade varia linearmente com o tempo devido ao termo \(\frac{k}{m}t\), cuja inclinação é igual à aceleração constante \(\frac{k}{m}\). A velocidade será aproximadamente constante quando a contribuição temporal for pequena em relação à velocidade inicial \(\frac{d}{\tau} - \frac{\tau k}{2m}\).

Fisicamente, isso ocorre quando a força constante \(k\) é pequena, quando a massa \(m\) é grande ou quando o intervalo de tempo observado é suficientemente curto. Nessas condições, a velocidade varia muito pouco ao longo do tempo e o movimento pode ser aproximado por um movimento uniforme.

A segunda derivada de \(x(t)\) representa a aceleração instantânea:

\[ x^{\prime\prime}(t)=\dfrac{k}{m}=a \] O movimento torna-se aproximadamente uniforme, com velocidade praticamente constante e posição aproximadamente linear no tempo se a velocidade variar muito pouco ao longo do tempo. Isso acontece quando a força constante \(k\) aplicada ao corpo é muito pequena ou quando a massa do corpo \(m\) é muito grande em comparação com essa força, i.e, \(k\ll m\). Nessas condições, o corpo praticamente mantém sua velocidade constante, caracterizando um movimento aproximadamente uniforme. Em termos físicos, isso significa que o efeito da força sobre o movimento é desprezível durante o intervalo de tempo observado.

Verificação da solução na equação da força:

\[ m x^{\prime\prime}(t)=m \dfrac{k}{m}=k \]

Além disso,

\[ x(0)=0 \]

e

\[ x(\tau)=d \]

Portanto, a função

\[ x(t)=\frac{k}{2m}t^2+\left(\dfrac{d}{\tau}-\frac{\tau k}{2m}\right)t \]

satisfaz tanto a equação diferencial como as condições de contorno.

m x''(t)= k, x(0)=0, x(tau)=d

m <- 1
k <- 2
a <- 3

x <- function(t, m, k, a) {
  (k/(2*m))*t*(t - 1) + a*t
}

curve(x(x, m = m, k = k, a = a),
      from = 0, to = 1,
      lwd = 2,
      xlab = "t",
      ylab = "x(t)",
      main = "Solução de mx''(t)=k com x(0)=0 e x(1)=3")

abline(h = 0, lty = 2)
abline(v = c(0, 1), lty = 3)
points(c(0, 1), c(0, a), pch = 19)
text(0, 0, labels = "(0,0)", pos = 4)
text(1, a, labels = paste0("(1,", a, ")"), pos = 2)

Nessa equação, \(m\) representa a massa do objeto, \(x\) representa a posição do objeto em relação ao tempo \(t\), e \(F\) representa a força aplicada sobre o objeto. A derivada segunda de \(x\) em relação a \(t\), denotada como \(x^{\prime\prime}(t)\), representa a aceleração do objeto.

A solução particular dessa equação diferencial depende da forma específica da força \(F\) e das condições iniciais ou de contorno do problema. Por exemplo, se a força \(F\) for uma força constante, como a gravidade, a solução particular pode ser obtida integrando a equação diferencial.

No entanto, é importante destacar que a equação do movimento de Newton é um exemplo específico de uma equação diferencial para o movimento de um objeto sob uma força constante. Existem várias outras equações diferenciais que descrevem o movimento de objetos em diferentes situações, dependendo das forças envolvidas e das condições específicas do problema.

Segunda Lei de Newton e Conservação da Quantidade de Movimento

Em formulações modernas, a segunda lei de Newton pode ser vista como consequência da dinâmica da quantidade de movimento.

Historicamente, Isaac Newton formulou a segunda lei originalmente em termos de “mudança da quantidade de movimento”, e não diretamente como

\[ \vec{F}=m\vec{a} \]

A forma mais geral da segunda lei é

\[ \vec{F}=\frac{d\vec{Q}}{dt} \]

em que \(\vec{Q}=m\vec{v}\) é a quantidade de movimento linear.

Se a massa é constante, então

\[ \vec{F} = \frac{d}{dt}(m\vec{v}) = m\frac{d\vec{v}}{dt} = m\vec{a} \]

Assim, \(\vec{F}=m\vec{a}\) é um caso particular da equação geral da quantidade de movimento.

A conservação da quantidade de movimento surge quando a força externa resultante sobre o sistema é nula:

\[ \vec{F}_{\mathrm{ext}}=0 \]

Nesse caso,

\[ \frac{d\vec{Q}_{\mathrm{total}}}{dt}=0 \]

o que implica

\[ \vec{Q}_{\mathrm{total}}=\text{constante} \]

Portanto, existe uma hierarquia conceitual:

\[ \vec{F}=\frac{d\vec{Q}}{dt} \quad\Longrightarrow\quad \vec{F}=m\vec{a} \quad\Longrightarrow\quad \vec{Q}_{\mathrm{total}}=\text{constante} \ \text{quando}\ \vec{F}_{\mathrm{ext}}=0 \]

Em mecânica moderna, a conservação da quantidade de movimento é frequentemente considerada mais fundamental do que a forma escolar

\[ \vec{F}=m\vec{a} \]

Na mecânica lagrangiana e hamiltoniana, por exemplo, a conservação do momento linear decorre da invariância espacial do sistema, resultado formalizado pelo teorema de Noether.

Segunda lei de Newton e movimento harmônico simples

• Harmonic oscillator: Wikipedia

O movimento harmônico simples pode ser obtido a partir da combinação da segunda lei de Newton com a lei de Hooke.

A segunda lei de Newton afirma que a o vetor da força resultante aplicada a um corpo é igual ao produto de sua massa pelo vetor da aceleração. Se \(x(t)\) representa o deslocamento da massa no instante \(t\).

A equação segunda lei de Newton pode ser escrita como

\[ F=m x^{\prime\prime}(t) \]

No caso de uma massa presa a uma mola ideal, a força elástica obedece à lei de Hooke:

\[ F = -k x(t) \]

em que \(k>0\) é a constante elástica da mola (grau de rigidez da mola, e seu inverso é o grau de complacência). O sinal negativo indica que a força é restauradora, isto é, atua sempre no sentido oposto ao deslocamento, puxando a massa de volta para a posição de equilíbrio.

Sistema bloco-mola em equilíbrio

Igualando as duas expressões da força, obtemos

\[ m x^{\prime\prime}(t) = -k x(t) \]

Reescrevendo:

\[ x^{\prime\prime}(t) = -\dfrac{k}{m} x(t) \]

Problema: Qual é a função \(x(t)\) que derivada duas vezes é igual ao valor negativo da função?

Resposta: \(x(t)=\sin(t)\): \(x^{\prime\prime}(t)=-\sin(t)\)

D[sin(t), t, t]

Reescrevendo novamente,

\[ x^{\prime\prime}(t) + \frac{k}{m}x(t)=0 \]

Essa é a equação diferencial ordinária linear homogênea do movimento harmônico simples.

As condições de contorno são:

\[ x(0)=0 \qquad \text{e} \qquad x(\tau)=d \]

A solução da EDO é:

\[ \begin{align} x(t)&=d\csc\left(\tau\sqrt{\frac{k}{m}}\right)\sin\left(t\sqrt{\frac{k}{m}}\right)\\ x(t)&= d\, \frac{ \sin\left(t\sqrt{\dfrac{k}{m}}\right) }{ \sin\left(\tau\sqrt{\dfrac{k}{m}}\right) } \end{align} \]

A primeira derivada em relação ao tempo representa a velocidade instantânea:

\[ v(t)=x'(t) \]

Derivando,

\[ v(t)= d\,\sqrt{\frac{k}{m}} \frac{ \cos\left(\sqrt{\dfrac{k}{m}}\,t\right) }{ \sin\left(\sqrt{\dfrac{k}{m}}\,\tau\right) } \]

Portanto,

\[ v(0)= \frac{d\, \sqrt{\dfrac{k}{m}} }{ \sin\left(\tau\sqrt{\dfrac{k}{m}}\right) } \] A segunda derivada em relação ao tempo representa a aceleração instantânea:

\[ a(t)=x^{\prime\prime}(t)=v'(t) \]

Derivando,

\[ a(t)= - d\frac{\,k}{m} \frac{ \sin\left(t\sqrt{\dfrac{k}{m}}\right) }{ \sin\left(\tau\sqrt{\dfrac{k}{m}}\right) } \]

Portanto,

\[ a(0)=0 \]

m x''(t)= -k x(t), x(0)=0, x(tau)=d

m <- 1
a <- 2
k_vals <- c(100, 200)

csc <- function(z) {
  1/sin(z)
}

x_fun <- function(t, k) {
  omega <- sqrt(k/m)
  a * csc(omega) * sin(omega * t)
}

titulo <- bquote(
  atop(
    m*x^{"''"} == -k*x ~~ "," ~~ x(0) == 0 ~~ "," ~~ x(1) == a,
    paste("m = ", .(m), ",  k = ", .(k_vals[1]), " e ", .(k_vals[2]), ",  a = ", .(a))
  )
)

tt <- seq(0, 1, length.out = 500)
yy1 <- x_fun(tt, k_vals[1])
yy2 <- x_fun(tt, k_vals[2])

curve(x_fun(x, k_vals[1]),
      from = 0, to = 1,
      lwd = 2,
      xlab = "t",
      ylab = "x(t)",
      main = titulo,
      ylim = range(c(yy1, yy2)))

curve(x_fun(x, k_vals[2]),
      from = 0, to = 1,
      add = TRUE,
      lwd = 2,
      lty = 2)

abline(h = 0, lty = 2)
abline(v = c(0, 1), lty = 3)
points(c(0, 1), c(0, a), pch = 19)
text(0, 0, labels = "(0,0)", pos = 4)
text(1, a, labels = paste0("(1,", a, ")"), pos = 2)

legend("topright",
       legend = paste("k =", k_vals),
       lwd = 2,
       lty = 1:2,
       bty = "n")

Introdução à equação diferencial ordinária

• Ordinary differential equation: Wikipedia

Quando a derivada \(y^{\prime} = f^{\prime}(t)\) de uma função desconhecida \(y = f(t)\) é dada, geralmente temos que encontrar a antiderivada. Algumas vezes a derivada \(y^{\prime}\) não é dada como uma função de \(t\), mas está envolvida em uma equação que contém também a função desconhecida \(y = f(t)\). Como exemplo, considere a equação

\[ y^{\prime} + ay = bt+ c \]

com coeficientes reais conhecidos \(a\), \(b\), \(c\).

Essa equação é chamada de equação diferencial, pois contém não apenas a função desconhecida, mas também sua derivada. O problema consiste em encontrar uma função adequada que satisfaça a equação diferencial.

As equações diferenciais ocorrem frequentemente na análise de sistemas fisiológicos e de sistemas ecológicos. Podemos falar brevemente de análise de sistemas. Quando uma quantidade varia em uma parte de um sistema, sua taxa de variação geralmente depende das quantidades em outras partes.

Além disso, qualquer mudança de uma quantidade pode influenciar indiretamente a própria quantidade, um fenômeno que é chamado de feedback.

O estudo dos sistemas de realimentação teve origem na engenharia, mas sua aplicação às ciências da vida acaba sendo mais frutífera.

A variável independente geralmente é o tempo. Portanto, denotamos por \(t\) na maior parte deste capítulo. Há poucas exceções.

Exemplos onde a variável independente não é o tempo são dados no Exemplo 11.3.5 e na Seção 11.6 sobre alometria.

Variáveis dependentes são indicadas por \(x = x(t)\), \(y = y(t)\), \(m = m(t)\), \(N = N(t)\), \(Q = Q(t)\) etc.

Os capítulos 9 e 10 são pré-requisitos para a compreensão das equações diferenciais.

Como o escopo deste livro é limitado, não nos aprofundaremos na teoria. Um estudo abrangente de equações diferenciais seria uma tarefa que duraria vários anos. Felizmente, existem listas de equações diferenciais e suas soluções disponíveis.

Existem equações diferenciais cujas soluções não podem ser escritas de forma gerenciável. Eles são resolvidos por computadores, aplicando métodos de análise numérica ou por simulação de computador.

Interpretação geométrica

y' = y/2

slope field y' = y/2

vector field y' = y/2

y'(t) = y(t)/2, y(0)=1

y'(t) = y(t)/2, y(0)=k

Figura 11.2a. O campo de inclinação (slope field) de y’=ay para a=1/2. Também são mostradas duas curvas integrais.

y' = y - t^2

slope field y' = y - t^2

vector field y' = y - t^2

y'(t) = y(t) - t^2, y(0)=1

y'(t) = y(t) - t^2, y(0)=a

Fig. 11.1. Campo de inclinação dado pela equação diferencial y’ = y - t². O campo de inclinação se assemelha à imagem de um fluido ou gás em movimento. As linhas de corrente são as curvas integrais.

Classificação de equação diferencial: DSolve do WolframAlpha

https://reference.wolfram.com/language/tutorial/DSolveOrdinaryDifferentialEquations.html

Equação diferencial ordinária linear geral

Uma equação diferencial linear é uma equação diferencial na qual a função desconhecida e suas derivadas aparecem de forma linear. Isso significa que a função e suas derivadas são multiplicadas por coeficientes (que podem ser funções de variáveis independentes), mas nunca são elevadas a uma potência diferente de 1 nem aparecem multiplicadas entre si.

Uma equação diferencial linear de ordem \(n\) pode ser escrita na forma geral:

\[ a_n(x) \dfrac{d^n y}{dx^n} + a_{n-1}(x) \dfrac{d^{n-1} y}{dx^{n-1}} + \dots + a_1(x) \dfrac{dy}{dx} + a_0(x) y = g(x) \]

Sendo que:

• \(y = y(x)\) é a função desconhecida.
• \(a_i(x)\) são funções conhecidas das variáveis independentes (geralmente \(x\)).
• \(g(x)\) é o termo não homogêneo ou fonte.
• A equação é chamada homogênea se \(g(x) = 0\) e não homogênea se \(g(x) \neq 0\).

Exemplo de equação diferencial linear de segunda ordem:

\[ (1-x)y^{\prime\prime} + x^2 y^{\prime} + \ln(x) y = \cos(x) \]

A característica principal da linearidade é que nenhuma potência da função \(y\) ou de suas derivadas aparece, nem produtos entre elas.

Equação diferencial ordinária linear mais simples

\[ \begin{cases} y^{\prime}(t) = a y(t)\\ y(0)=y_0 \end{cases} \] Solução: \(y(t) =y_0e^{at}\)

A equação diferencial

\[ y^{\prime}(t)=a y(t) \]

é a equação diferencial ordinária linear de primeira ordem mais simples e um dos modelos fundamentais do cálculo diferencial. Ela expressa a ideia de que a taxa de variação de uma quantidade é proporcional ao seu próprio valor no instante considerado.

Fig. 11.2b. Para y positivo, a solução de y’ = ay é uma função crescente se a > 0 e uma função decrescente se a < 0.

Esse tipo de relação aparece naturalmente em diversos contextos, como crescimento populacional, juros compostos contínuos, decaimento radioativo, eliminação de fármacos e crescimento de microrganismos. Em todos esses casos, a hipótese central é a mesma: quanto maior o valor atual de \(y(t)\), maior, em módulo, tende a ser sua variação instantânea.

Se \(a>0\), o modelo descreve crescimento exponencial; se \(a<0\), descreve decaimento exponencial; e, se \(a=0\), a solução é constante. Assim, o parâmetro \(a\) controla o comportamento dinâmico da solução.

O problema de valor inicial

\[ \begin{align} y^{\prime}(t)& = a y(t), \quad y(0)=y_0 \end{align} \]

tem solução

\[ y(t)=y_0 e^{at} \]

Essa função satisfaz a equação diferencial porque sua derivada é

\[ y^{\prime}(t)=a y_0 e^{at}=a y(t) \]

e satisfaz também a condição inicial, pois

\[ y(0)=y_0 e^0=y_0 \]

Portanto, essa equação fornece o modelo matemático mais elementar para processos cuja velocidade de crescimento ou de decaimento é proporcional à quantidade existente.

a <- -0.5
y0 <- 1

y <- function(t) {
  y0 * exp(a*t)
}

curve(y,
      from = 0, to = 10,
      lwd = 2,
      xlab = "t",
      ylab = "y(t)",
      main = bquote(atop(y^{"'"} == .(a)*y ~~ "," ~~ y(0) == .(y0),
                         y(t) == .(y0)*e^{.(a)*t})))

y'(t) = a y(t), y(0) = b

DSolve[y'[t] == a y[t], y, t]

Exemplo 11.3.1: Crescimento da célula

Suponha que uma célula tenha massa \(m_0\).

Em um ambiente ideal, a célula cresce. Assim, sua massa é uma função do tempo e podemos escrever \(m = m(t)\) com \(m = m_0\) em \(t = 0\).

Suponha que as substâncias químicas passem rapidamente pela parede celular e que o crescimento seja determinado apenas pela velocidade do metabolismo interno a célula. Como o rendimento do metabolismo depende da massa das moléculas participantes, é razoável esperar que a taxa de crescimento seja proporcional à massa em cada instante de tempo, ou seja,

\[ \dfrac{dm}{dt} \propto m \]

ou

\[ \dfrac{dm}{dt} =a m \]

com uma certa constante positiva \(a\).

Claro, há uma limitação: se a massa \(m\) da célula atingir um certo tamanho, a célula se dividirá em vez de continuar a crescer. Assim, adicionamos uma restrição, \(m < m_1\).

Portanto, a solução geral segue:

\[ m(t)=ce^{at} \]

Pela nossa suposição de que \(m = m_0\) no instante de tempo \(t = 0\), podemos determinar a constante \(c\). Obtemos \(c=m_0\).

Portanto, a solução particular:

\[ m(t)=m_0e^{at} \]

Com nossas suposições, fomos um pouco além da experiência. Apresentamos alguns argumentos teóricos. Costuma-se dizer que estamos construindo um modelo. Se nosso modelo é ou não biologicamente significativo, só pode ser testado por experimentos. Aqui e em modelos subsequentes, compartilhamos a visão de G. F. Gause (Gause, 1934, p. 10):

Gause, G. F. (1934). The Struggle for Existence. USA: Williams and Wilkins Company, p. 10

“Não há dúvida de que [crescimento] é um problema biológico, e que deve ser resolvido por experimentação e não na mesa de um matemático. Mas, para penetrar mais profundamente na natureza desses fenômenos, devemos combinar o método experimental com a teoria matemática, uma possibilidade que foi criada por pesquisadores [brilhantes]. A combinação do método experimental com a teoria quantitativa é em geral uma das ferramentas mais poderosas nas mãos da ciência contemporânea.”

Vale a pena discutir o modelo de crescimento acima sob diferentes aspectos. Como \(dm/dt\) foi considerado proporcional a \(m\), podemos introduzir a taxa de crescimento específico ou relativo definida por

\[ \dfrac{dm/dt}{m}=a \]

Este é o quociente entre a taxa de crescimento absoluto \(dm/dt\) e a massa \(m\).

Nossa equação diferencial afirma: Em cada instante de tempo, a taxa de crescimento relativa permanece constante.

A noção de uma taxa de crescimento relativo precisa de alguma ilustração.

Suponha que uma planta que atingiu a massa \(m = 300\) g cresça 12 g nas próximas 24 horas. Então a taxa média de crescimento é de 12 g/24 horas = 0.5 g/h. Assumindo que a taxa de crescimento não flutua, podemos considerar 0.5 g/h como uma boa aproximação da taxa de crescimento instantânea \(dm/dt\).

Podemos perguntar: esta taxa de crescimento é grande ou pequena?

A resposta depende muito da massa atual da planta. Para uma planta de massa \(m = 10\) g apenas, nossa taxa de crescimento seria tremenda, enquanto para uma grande árvore de massa (viva) \(m = 1000\) kg, o mesmo crescimento deve ser chamado de minúsculo. Portanto, temos que relacionar 0.5 g/h com a massa atual, no nosso caso com 300 g. O quociente é

\[ \dfrac{0.5\;g/h}{300\; g}=0.0017 \; h^{-1} \]

Este quociente é chamado de taxa de crescimento específico. Com a mesma taxa específica de crescimento, a árvore de massa (viva) de 1000 kg ganharia 1.7 kg/hora.

A taxa de crescimento específico é um conceito importante. Há duas etapas envolvidas. Primeiro, ao formar \(dm/dt\), relacionamos o aumento da massa com o tempo, o que nos dá alguma medida da velocidade de crescimento.

Em segundo lugar, relacionamos a velocidade de crescimento com a massa presente.

Vamos finalmente considerar outro aspecto da equação diferencial.

Com o aumento de \(m\), a taxa de crescimento \(dm/dt\) também aumenta. Essa taxa de crescimento, por sua vez, determina os valores futuros de \(m\). Assim, temos um exemplo simples de um mecanismo de feedback com um único loop:

Exemplo 11.3.2: Processo de nascimento

Seja \(N\) o número de indivíduos em uma população animal ou vegetal. Este número é dependente do tempo, de modo que podemos escrever \(N = N(t)\).

A rigor, \(N(t)\) assume apenas valores inteiros e é uma função descontínua de \(t\). Entretanto, \(N(t)\) pode ser aproximado por uma função contínua e diferenciável assim que o número de indivíduos for grande o suficiente.

Em alguns microrganismos, a reprodução ocorre por simples divisão celular. A reprodução por simples divisão celular, também conhecida como fissão binária, ocorre principalmente em organismos procariontes e alguns organismos unicelulares. Os principais grupos de microorganismos que utilizam esse método de reprodução incluem:

1. Bactérias: As bactérias se reproduzem quase exclusivamente por fissão binária, onde a célula se divide em duas células-filhas geneticamente idênticas.
2. Arqueas: Semelhante às bactérias, as arqueas também se reproduzem por fissão binária.
3. Alguns Protistas Unicelulares: Alguns protistas, como as amebas e alguns tipos de algas unicelulares, podem se reproduzir por divisão celular.

Esse processo de divisão envolve a duplicação do material genético seguido pela divisão do citoplasma, resultando em duas células-clone independentes.

Em indivíduos multicelulares, distinguimos entre reprodução vegetativa e sexual. Incluiremos todas essas possibilidades em nosso estudo.

Assumimos que a proporção de indivíduos reprodutivos (viabilidade) permanece constante na população crescente. Além disso, assumimos fertilidade constante. Então a taxa de natalidade é proporcional ao número \(N(t)\) de indivíduos. Se por fim excluirmos a morte, a emigração e a imigração, a taxa de crescimento coincide com a taxa de natalidade. Por isso:

\[ \dfrac{dN}{dt}=\lambda N \]

sendo que \(\lambda>0\) é uma constante chamada de taxa de natalidade específica.

A solução é:

\[ N(t)=N_0e^{\lambda t} \]

sendo que \(N_0\) denota o tamanho da população em \(t = 0\).

N0 <- 100
lambda <- 0.2

t <- seq(0, 20, by = 0.1)
N <- N0 * exp(lambda * t)

plot(t, N,
     type = "l",
     lwd = 2,
     xlab = "t",
     ylab = "N(t)",
     main = expression(N(t) == N[0] * e^{lambda * t}))

eta <- lambda * t

plot(t, eta,
     type = "l",
     lwd = 2,
     xlab = "t",
     ylab = expression(eta),
     main = expression(eta == lambda * t))

abline(h = 0, lty = 2)
abline(v = 0, lty = 2)

Esse processo de nascimento acaba sendo bastante realista em uma grande população que cresce em condições ideais, ou seja, quando todos os fatores inibidores do crescimento estão ausentes.

Em uma população pequena, não podemos esperar que a ocorrência de nascimentos seja distribuída uniformemente ao longo do tempo. Em vez disso, enfrentamos flutuações aleatórias. Então o processo tem que ser modificado à luz da teoria da probabilidade.

n'(t) = lambda n(t), n(0)=a

plot 10 e^(t), t from 0 to 5 axes label "t" "y"

Exemplo 11.3.3: Processo de nascimento e morte

Consideremos uma população animal ou vegetal nas condições descritas no exemplo anterior. Agora vamos estender o modelo permitindo a morte. A mudança líquida no tamanho da população pode ser positiva ou negativa.

Tratamos \(N = N(t)\) como uma função contínua e diferenciável do tempo, embora isso signifique apenas uma aproximação da realidade.

Da mesma forma, assumimos um grande número de nascimentos e mortes, de modo que o número de nascimentos \(B = B(t)\) e de mortes \(D = D(t)\) também podem ser considerados como funções diferenciáveis.

A taxa de variação líquida é igual à taxa de nascimento menos a taxa de mortalidade.

A taxa \(dN/dt\) pode ser positiva ou negativa conforme predominem as ocorrências de nascimento ou de óbito.

Dentro de um intervalo de tempo de comprimento \(dt\), obtemos

\[ \text{variação líquida} = \text{número de nascimentos - número de mortes} \]

ou, em uma notação conveniente, após a divisão por \(dt\), obtemos a taxa de variação:

\[ \dfrac{dN}{dt}=\dfrac{dB}{dt}-\dfrac{dD}{dt} \]

No Exemplo 11.3.2, estabelecemos hipóteses tais que a taxa de natalidade torna-se proporcional ao número de indivíduos \(N(t)\).

Sob suposições correspondentes sobre a morte, a taxa de mortalidade também se torna proporcional a \(N(t)\).

Portanto,

\[ \dfrac{dB}{dt}=\lambda N \\ \dfrac{dD}{dt}=\mu N \]

sendo que \(\lambda\) denota a taxa de natalidade específica e \(\mu\), a taxa de mortalidade específica.

Portanto,

\[ \dfrac{dN}{dt}=\lambda N - \mu N\\ \dfrac{dN}{dt}=(\lambda - \mu) N \]

A solução desta equação diferencial é:

\[ N(t)=N_0 \;e^{(\lambda - \mu)t} \]

sendo que \(N_0\) representa o tamanho da população no tempo \(t = 0\).

Quando prevalece a taxa de natalidade, ou seja, quando \(\lambda > \mu\), o tamanho da população aumenta exponencialmente. Temos uma erupção.

Quando, em vez disso, \(\lambda< \mu\), o tamanho da população diminui e a população é extinta.

Somente para \(\lambda =\mu\) a população permanecerá estável.

Esse modelo de processo de nascimento e morte não leva em conta flutuações aleatórias. É, portanto, chamado de não estocástico ou determinístico.

Exemplo 11.3.5: Cinética química

O pentóxido de nitrogênio gasoso se decompõe conforme indicado pela equação

\[ 2N_2O_5 \to 4NO_2+O_2 \]

Estamos interessados na velocidade dessa reação quando a temperatura é mantida constante.

Seja \(C=[N_2O_5]\) a concentração de pentóxido de nitrogênio gasoso medida em mol por litro.

A concentração \(C = C(t)\) é uma função decrescente do tempo de forma que a derivada \(dC/dt\) é negativa.

Essa derivada é chamada de taxa de reação.

A velocidade da reação depende da concentração \(C = [N_2O_5]\).

Intuitivamente, esperamos que quanto maior a concentração, mais freqüentemente ocorrerão colisões de duas moléculas de \(N_2O_5\) com o possível surgimento das novas ligações \(NO_2\) e \(O_2\).

Pode-se teorizar que sob temperatura constante a taxa de reação é proporcional a \(C\), que é,

\[ \dfrac{dC}{dt}=-kC \]

sendo que \(k\) denota uma constante positiva.

A solução deste equação diferencial é

\[ C(t)=C_0\; e^{-kt} \]

sendo \(C_0\) a concentração de \(N_2O_5\) no tempo \(t=0\).

Os fatos experimentais estão de acordo com este modelo. A concentração \(C\) tenderá assintoticamente para zero. Nunca chegará a zero exatamente.

c'(t) = -k c(t), c(0)=a

plot 10 e^t, t from 0 to 5 axes label "t" "C"

Equação diferencial ordinária linear

\[ \begin{cases} y^{\prime}(t)= ay(t) + b\\ y(0)=y_0 \end{cases} \] Solução: \(y(t)=\left(y_0+\dfrac{b}{a}\right)e^{at}-\dfrac{b}{a}\)

A equação diferencial

\[ y^{\prime}(t)=a y(t)+b \]

é uma equação diferencial ordinária linear de primeira ordem com coeficientes constantes. Ela generaliza o modelo mais simples

\[ y^{\prime}(t)=a y(t) \]

ao incluir o termo constante \(b\), que representa uma contribuição externa fixa à taxa de variação de \(y(t)\).

Esse modelo aparece em várias situações aplicadas. Em termos conceituais, a variação instantânea de \(y(t)\) passa a depender de duas componentes: uma parte proporcional ao próprio estado atual, dada por \(a y(t)\), e uma parte independente de \(y(t)\), dada por \(b\). Assim, o modelo pode descrever processos de crescimento ou decaimento com entrada, produção, reposição ou remoção constante.

Quando \(a\neq 0\), o problema de valor inicial

\[ \begin{align} y^{\prime}(t)&= a y(t) + b, \quad y(0)=y_0 \end{align} \]

tem solução

\[ y(t)=\left(y_0+\dfrac{b}{a}\right)e^{at}-\dfrac{b}{a} \]

Essa solução mostra que o comportamento de \(y(t)\) resulta da combinação de uma componente exponencial com um nível de equilíbrio constante. De fato, se existir um valor constante \(y^{\ast}\) tal que \(y^{\prime}(t)=0\), então

\[ a y^{\ast}+b=0 \]

isto é,

\[ y^{\ast}=-\dfrac{b}{a} \]

Portanto, a solução pode ser interpretada como uma trajetória que evolui em torno desse valor de equilíbrio.

Além disso, essa função satisfaz a equação diferencial, pois sua derivada é

\[ y^{\prime}(t)=a\left(y_0+\dfrac{b}{a}\right)e^{at} \]

e, substituindo \(y(t)\) no lado direito,

\[ a y(t)+b = a\left[\left(y_0+\dfrac{b}{a}\right)e^{at}-\dfrac{b}{a}\right]+b = a\left(y_0+\dfrac{b}{a}\right)e^{at} \]

Logo,

\[ y^{\prime}(t)=a y(t)+b \]

Além disso,

\[ y(0)=\left(y_0+\dfrac{b}{a}\right)e^0-\dfrac{b}{a}=y_0 \]

Portanto, a expressão apresentada é de fato a solução do problema de valor inicial.

a <- -0.5
b <- 2
y0 <- 1

y <- function(t) {
  (y0 + b/a) * exp(a*t) - b/a
}

curve(y,
      from = 0, to = 10,
      lwd = 2,
      xlab = "t",
      ylab = "y(t)",
      main = bquote(atop(y^{"'"} == .(a)*y + .(b) ~~ "," ~~ y(0) == .(y0),
                         y(t) == (.(y0) + .(b)/.(a))*e^{.(a)*t} - .(b)/.(a))))

abline(h = -b/a, lty = 2)

y'(t)= ay(t) + b, y(0)=c

DSolve[y'[t] == a y[t] + b, y[t], t]

c_1 = y(0) + b/a

y(t) = (y(0) + b/a) e^(a t) - b/a

• Mathematica

DSolve[{y'[t] == a y[t] + b, y[0] == c}, y[t], t]

\[ y(t) = \dfrac{b \left(e^{a t}-1\right)}{a}+c e^{a t} \]

Exemplo 11.4.1: Crescimento limitado

Nenhum organismo e nenhuma população cresce indefinidamente.

Existem limitações impostas pela escassez de alimentos ou abrigo, pela falta de espaço, por condições físicas intoleráveis ou por algum mecanismo de controle.

Assuma que existe um limite superior fixo para o tamanho \(y\) de um indivíduo, um tecido, uma população ou uma cultura. O tamanho pode ser um volume, uma massa, um diâmetro, uma quantidade etc.

Denotamos o limite superior por \(B\).

Então \(y = y(t)\) pode se aproximar de \(B\) assintoticamente.

Isso implica que a taxa de crescimento \(dy/dt\) tende a 0 à medida que \(B-y\) se torna cada vez menor.

Uma formulação matemática plausível de tal modelo é dada pela equação diferencial:

\[ \dfrac{dy}{dt} = k(B - y) \]

sendo que \(k\) é uma constante positiva que determina a rapidez com que \(dy/dt\) tende a 0.

Se \(y\) é pequeno em relação a \(B\), então temos aproximadamente \(y' \approx kB = \text{constante}\), ou seja, o tamanho de \(y\) aumenta aproximadamente em função linear do tempo.

No entanto, se \(y\) estiver próximo de \(B\), então \(B - y\) é uma pequena quantidade positiva, assim como a taxa de crescimento \(dy/dt\).

Esta equação diferencial pode ser reescrita na forma:

\[ y^{\prime}=kB-ky=ay+b \]

com \(a = - k\) e \(b = kB\).

Assim, obtemos a solução particular, se \(y(0)=0\):

\[ y(t)=B\left(1-e^{-kt}\right) \]

Este modelo foi proposto por E. A. Mitscherlich em 1939. Ele se ajusta muito bem a alguns dados experimentais na agricultura.

y'(t)=k (B-y(t)), y(0)=0

plot 2 (1 - e^(-t)), t from 0 to 5 axes label "t" "y"

limit B (1 - e^(-k t)) as t -> infinity

Exemplo 11.4.2: Processo de nascimento e imigração

No Exemplo 11.3.2, estudamos um processo de nascimento. Várias suposições foram feitas sobre a população. Mantemos essas suposições com a exceção de que agora permitimos a imigração de indivíduos a uma taxa constante. A taxa é medida em número de indivíduos por unidade de tempo e denotada por \(\nu>0\). Então:

\[ \dfrac{dN}{dt}=\lambda N+\nu \]

Portanto, a solução particular é:

\[ N(t)=\left(N_0+\dfrac{\nu}{\lambda}\right)e^{\lambda t}-\dfrac{\nu}{\lambda} \]

Exemplo 11.4.3: Lei de Newton do resfriamento

Seja \(T = T(t)\) a temperatura do corpo no instante de tempo \(t\), \(T_0\) sua temperatura no instante \(t = 0\), e \(T_s\) a temperatura constante do meio ambiente.

Considere um corpo sem aquecimento interno cuja temperatura do corpo é superior à do meio ambiente. O corpo irá então resfriar-se. Desejamos saber como a temperatura do corpo decresce em função do tempo.

A derivada \(\dfrac{dT}{dt}\) é chamada de taxa de resfriamento. Como \(T\) decresce, essa taxa é negativa. A derivada depende da diferença \(T - T_s\). Sob condições favoráveis, a taxa de resfriamento é proporcional a \(T - T_s\), ou seja,

\[ \dfrac{dT}{dt} = -k(T - T_s) \tag{11.4.8} \]

sendo que \(k\) é uma constante positiva determinada pelas condições físicas da troca de calor. Como o lado direito pode ser escrito como \((-k)T + kT_s\), a equação diferencial está na forma da equação (11.4.1). A partir de (11.4.2), obtemos a solução geral:

\[ T(t) = ce^{-kt} + T_s \tag{11.4.9} \]

Por fim, aplicamos a condição inicial \(T = T_0\) no instante \(t = 0\). Isso leva a \(c = T_0 - T_s\) e, portanto:

\[ T(t) = T_s + (T_0 - T_s)e^{-kt} \tag{11.4.10} \]

À medida que \(t\) tende ao infinito, o segundo termo tende a zero, e \(T\) aproxima-se de \(T_s\) assintoticamente. A equação (11.4.8) é conhecida como lei do resfriamento de Newton.

Esta e a aplicação a seguir podem ser consideradas ilustrações da análise compartimental (cf. Exemplo 11.3.6).

y'(t) = - k (y(t)-c), y(0)=a

Exemplo 11.4.4: Difusão na célula

Suponhamos que uma célula de volume constante esteja suspensa em um líquido homogêneo que contém um soluto de concentração \(c_0\), constante tanto no espaço quanto no tempo. Seja \(c = c(t)\) a concentração do soluto dentro da célula no instante de tempo \(t\), e suponha-se que o soluto esteja quase uniformemente distribuído ao longo da célula em todos os momentos, de modo que \(c = c(t)\) dependa apenas do tempo.

Fig. 11.4. Difusão de moléculas através da parede celular. Na figura, assume-se que c₀ > c(t) e, portanto, que mais moléculas entram na célula do que saem dela.

Pela difusão, moléculas do soluto entrarão na célula a partir do líquido circundante, mas também haverá moléculas do soluto que sairão da célula. Assim, há um fluxo de moléculas através da membrana celular em ambas as direções (ver Fig. 11.4). O fluxo líquido é do líquido para o interior da célula se \(c_0\) for maior que \(c(t)\), e vice-versa. Estamos interessados em determinar a função \(c(t)\).

Seja \(m = m(t)\) a massa de soluto na célula, \(A\) a área da membrana celular, e \(V\) o volume da célula. Pela definição de concentração, temos

\[ m(t) = V \, c(t). \tag{11.4.11} \]

A derivada \(\dfrac{dm}{dt}\) representa a taxa de variação de \(m\), e pode ser chamada de taxa de fluxo líquido em nosso problema. A primeira lei de Fick afirma que \(\dfrac{dm}{dt}\) é proporcional à área da membrana e à diferença de concentração nos dois lados da membrana. Assim, temos

\[ \dfrac{dm}{dt} = kA(c_0 - c) \tag{11.4.12} \]

Se \(c < c_0\), isto é, se o soluto tem menor concentração dentro da célula do que fora, \(m\) aumentará. Portanto, \(k\) é uma constante positiva. Essa constante é determinada pela estrutura e espessura da membrana, e é chamada de permeabilidade da membrana para o soluto em questão.

Com base na equação (11.4.11), podemos substituir \(\dfrac{dm}{dt}\) por \(V \dfrac{dc}{dt}\) em nossa equação diferencial. Obtemos, então:

\[ \dfrac{dc}{dt} = \dfrac{kA}{V} (c_0 - c) \tag{11.4.13} \]

Integramos essa equação utilizando a solução explícita da equação (11.4.2), obtendo:

\[ c(t) = K \exp\left( -\dfrac{kA}{V}t \right) + c_0 \tag{11.4.14} \]

sendo que \(K\) é a constante de integração. Quando \(t\) tende ao infinito, \(c(t)\) aproxima-se assintoticamente de \(c_0\). A constante \(K\) pode ser determinada a partir de uma condição inicial, por exemplo, \(c = c^{\ast}\) no instante \(t = 0\). Deixamos ao leitor a tarefa de discutir os dois casos: \(c^{\ast} > c_0\) e \(c^{\ast} < c_0\).

y'(t) = k (A/V) (a - y(t))

Esta aplicação foi adaptada de Thrall et al. (1967, CA 10).

Ressalta-se que este modelo é uma aproximação grosseira da realidade. A difusão através de membranas celulares é um processo complexo que não pode ser tratado adequadamente neste contexto.

Equação diferencial ordinária não-linear

\[ \begin{cases} y^{\prime}(t) = a y^2(t) + b y(t) + c\\ y(0) = y_0 \end{cases} \] Solução: \(y(t) = A + \dfrac{B-A}{1+\dfrac{B-y_0}{y_0-A}e^{a(B-A)t}}\)

Sendo que \(A\) e \(B\) são raízes reais de \(a y^2 + b y + c=0\).

A equação diferencial

\[ y^{\prime}(t)=a(y(t))^2+b y(t)+c \]

é uma equação diferencial ordinária não-linear de primeira ordem do tipo de Riccati com coeficientes constantes. Diferentemente dos modelos lineares, a taxa de variação de \(y(t)\) passa a depender quadraticamente da própria variável, o que pode produzir comportamentos mais ricos, como crescimento acelerado, saturação, instabilidade ou aproximação a equilíbrios.

Uma forma útil de analisar essa equação é fatorar o polinômio quadrático

\[ a y^2+b y+c=a(y-A)(y-B) \]

em que \(A\) e \(B\) são as raízes da equação algébrica

\[ a y^2+b y+c=0 \]

Nessa forma, a equação diferencial pode ser escrita como

\[ y^{\prime}(t)=a(y(t)-A)(y(t)-B) \]

Assim, os valores \(A\) e \(B\) são soluções de equilíbrio, pois, se \(y(t)=A\) ou \(y(t)=B\), então \(y^{\prime}(t)=0\).

Quando \(A\neq B\) e \(y_0\neq A\), a solução do problema de valor inicial

\[ \begin{align} y^{\prime}(t) & = a y^2(t) + b y(t) + c, \quad y(0) = y_0 \end{align} \]

pode ser escrita como

\[ y(t)=A+\dfrac{B-A}{1+\dfrac{B-y_0}{y_0-A}e^{a(B-A)t}} \]

Essa expressão mostra como a solução evolui em função dos equilíbrios \(A\) e \(B\), da condição inicial \(y_0\) e dos parâmetros do modelo. Dependendo dos sinais de \(a\) e de \(B-A\), a trajetória pode convergir para um equilíbrio, afastar-se dele ou até explodir em tempo finito.

A solução tem forma logística generalizada. Ela descreve uma transição entre dois níveis de equilíbrio, \(A\) e \(B\).

A função logística clássica pode ser escrita como

\[ L(t)=\frac{K}{1+Ce^{-rt}} \]

ou, de forma mais geral,

\[ L(t)=A+\frac{B-A}{1+Ce^{-rt}} \]

Portanto, a solução apresentada tem estrutura logística, com

\[ C=\frac{B-y_0}{y_0-A} \qquad \text{e} \qquad r=-a(B-A) \]

Mais precisamente, trata-se de uma função logística deslocada e escalonada, pois a trajetória não varia necessariamente entre \(0\) e \(K\), mas entre os equilíbrios \(A\) e \(B\).

Em resumo,

\[ y^{\prime}=ay^2+by+c \]

não é, em geral, chamada equação logística. No entanto, quando o polinômio quadrático pode ser fatorado como

\[ y^{\prime}=a(y-A)(y-B) \]

sua solução pode ser escrita em forma logística.

Equivalência entre as formas da equação logística generalizada:

Considere a equação diferencial

\[ y^{\prime} = a(y-A)(y-B) \]

onde \(A\) e \(B\) são pontos de equilíbrio da dinâmica.

Expandindo o produto:

\[ (y-A)(y-B) = y^2-(A+B)y+AB \]

Logo,

\[ y^{\prime} = ay^2-a(A+B)y+aAB \]

Portanto, a equação também pode ser escrita na forma quadrática geral

\[ y^{\prime} = ay^2+by+c \]

com

\[ b=-a(A+B) \]

e

\[ c=aAB \]

Assim, as duas EDOs são equivalentes.

Solução usando os pontos de equilíbrio:

A solução da EDO

\[ y^{\prime} = a(y-A)(y-B) \]

com condição inicial

\[ y(0)=d \]

é

\[ y(t) = A+ \frac{B-A}{ 1+ \dfrac{B-d}{d-A}e^{a(B-A)t} } \]

Esta forma evidencia diretamente:

• os equilíbrios \(A\) e \(B\)
• o valor inicial \(d\)
• a assíntota da solução
• o comportamento logístico generalizado

Relação com a solução exponencial simbólica:

Sistemas algébricos computacionais frequentemente devolvem soluções na forma

\[ y(t)= \frac{ Be^{aAt+Ac_1} - Ae^{aBt+Bc_1} }{ e^{aAt+Ac_1} - e^{aBt+Bc_1} } \]

onde \(c_1\) é uma constante arbitrária.

Esta expressão é matematicamente equivalente à forma anterior.

Impondo

\[ y(0)=d \]

obtém-se uma expressão para \(c_1\) em função de \(d\).

Após simplificação algébrica, a solução reduz-se à forma racional

\[ y(t) = A+ \frac{B-A}{ 1+ \dfrac{B-d}{d-A}e^{a(B-A)t} } \]

que é muito mais adequada para interpretação matemática e gráfica.

Interpretação dinâmica:

Os valores

\[ y=A \qquad \text{e} \qquad y=B \]

são soluções estacionárias da EDO, pois

\[ y^{\prime}=0 \]

quando

\[ y=A \qquad \text{ou} \qquad y=B \]

Dependendo do sinal de \(a\):

• uma das soluções estacionárias será estável
• a outra será instável

A solução geral aproxima-se assintoticamente de um desses equilíbrios.

A <- 10  
B <- 100  
y0 <- 20  
a <- -0.01  

y <- function(t) {
  
  A + (B - A) /
    (1 + ((B - y0) / (y0 - A)) * exp(a * (B - A) * t))
  
}

curve(
  y,
  from = 0, to = 11,
  lwd = 2,
  xlab = "t",
  ylab = "y(t)",
  ylim = c(A - 5, B + 5),
  main = bquote(
    y(t) == .(A) + frac(
      .(B) - .(A),
      1 + frac(.(B) - .(y0), .(y0) - .(A)) *
        e^{.(a) * (.(B) - .(A)) * t}
    )
  )
)

abline(h = A, lty = 2)
abline(h = B, lty = 2)

text(
  10,
  A,
  labels = expression(A == 10),
  pos = 3
)

text(
  10,
  B,
  labels = expression(B == 100),
  pos = 1
)

points(0, y0, pch = 19)

text(
  0,
  y0,
  labels = paste0("(0,", y0, ")"),
  pos = 4
)

y'(t)= a (y(t))^2 + b y(t) + c

y'(t)= a (y(t) -A) (y(t) - B)

solve for k: y_0=A+(B-A)/(1+k)

A + (B-A)/(1+ke^(a(B-A)t)) /. k = (y_0 - B)/(A - y_0)

DSolve[y'[t] == a (y[t])^2 + b y[t] + c, y[t], t]

a y^2 + b y + c=0

Exemplo 11.5.1. Crescimento limitado II

Voltamos ao estudo do crescimento das populações.

Seja \(y = y(t)\) o número de indivíduos em uma população no instante de tempo \(t\).

A equação diferencial \(y^{\prime} = ay\) com \(a> 0\) fornece um crescimento exponencial irrestrito, enquanto \(y^{\prime} = a(B - y)\) com \(a> 0\), \(B > 1\) resulta em crescimento que é quase linear no início e se nivela mais tarde.

Para obter um modelo de crescimento que seja biologicamente mais significativo, podemos combinar as duas abordagens, ou seja, assumir que \(y^{\prime}\) é proporcional a \(y\) assim como a \(B - y\). Essa ideia leva à equação diferencial

\[ \dfrac{dy}{dt}=\lambda y(B-y) \]

sendo que \(\lambda>0\).

Portanto, a solução particular é, sendo \(y(0)=y_0\), \(A=0\) e \(a=-\lambda\):

\[ \begin{align} y(t)&= A + \dfrac{B-A}{1+\dfrac{B-y_0}{y_0-A}e^{a(B-A)t}} \\\\ y(t)&= \dfrac{B}{1+\dfrac{B-y_0}{y_0}e^{-\lambda B t}} \end{align} \]

Sendo que \(A\) e \(B\) são raízes reais de \(\lambda y(B-y)=-\lambda y^2 + \lambda B y=0\).

y'(t) = lambda y(t) (B - y(t))

Em nosso modelo, \(y\) nunca pode exceder \(B\).

Portanto, o denominador da solução deve ser maior que \(1\) e \((B-y_0)/y_0\) deve ser restrito a valores positivos.

A quantidade \(y\) aumenta monotonicamente desde que a equação diferencial implica que \(dy/dt> 0\).

Para \(t\to-\infty\), \(y\) tende a 0, e para \(t\to+\infty\), \(y\) tende a \(B\).

O crescimento começa lentamente, então se torna mais rápido e, finalmente, desacelera. O crescimento é mais rápido na vizinhança do ponto de inflexão. Para obter a sua localização temos que igualar a segunda derivada de \(y = y(t)\) a 0.

\[ \dfrac{d^2y}{dt^2}=\lambda (B-2y)\dfrac{dy}{dt} \]

Esta expressão só pode ser nula se \(B - 2y = 0\) ou \(y=B/2\), ou seja, o ponto de inflexão está a meio caminho entre as linhas \(y = 0\) e \(y = B\).

Para obter a abcissa, fazemos \(y = B/2\) na equação da solução particular e resolvemos em relação a \(t\):

\[ t=\dfrac{\ln\left(\dfrac{B-y_0}{y_0}\right)}{B\lambda} \]

Esta abcissa particular é positiva ou negativa dependendo se \((B-y_0)/y_0>1\) ou \((B-y_0)/y_0< 1\).

y’(t)= lambda y(t) (B - y(t)), y(0)=1

solve for t: D[(B e^((B t)/3))/(-1 + B + e^((B t)/3)), {t,2}]=0

A Fig. 11.5 mostra um gráfico da função da solução particular em uma aplicação particular. A curva é em forma de S ou sigmoide.

A fórmula da solução particular é geralmente conhecida como a função logística.

Ela foi introduzido na dinâmica populacional por Verhulst em 1838.

O gráfico da função logística é chamado de curva logística.

Existem inúmeros dados experimentais de crescimento, especialmente para protozoários e bactérias, para os quais o ajuste de uma curva logística foi bastante bem-sucedido. Mas em algumas outras populações o ajuste foi ruim e a previsão foi enganosa.

Fig. 11.5. Crescimento de uma população de Drosophila sob condições experimentais controladas.

Exemplo 11.5.2: Propagação de infecção

Como uma doença infecciosa se espalha em uma comunidade de indivíduos suscetíveis?

Claro, isso depende de muitas circunstâncias. Para simplificar, faremos algumas suposições. Em uma população de indivíduos igualmente suscetíveis, introduzimos um único infeccioso. Pelo contato entre os indivíduos, a doença se espalhará, i.e., o número de infectantes aumentará. No início, o número de infecciosos aumentará lentamente, depois o processo será acelerado e, finalmente, estabilizado quando a maioria dos indivíduos se tornarem infectados. Além disso, assumimos que um indivíduo que já foi infectado permanecerá infectado durante o processo e que ninguém será removido.

Seja \(x = x(t)\) o número de susceptíveis, \(y = y(t)\) o número de infectados no instante \(t\) e \(n\) o tamanho total da população na qual um infectado foi introduzido. Assumimos agora que o processo começou no tempo \(t=0\) com o único infectado. Assim, a qualquer momento,

\[ x+y=n+1 \]

Tratamos \(x\) e \(y\) como variáveis reais. A taxa na qual o número de infectados aumenta é então \(dy/dt\). Quanto mais infectados e susceptíveis estiverem presentes, mais frequentemente ocorrerão contatos que levarão à infecção. Portanto, é plausível assumir que \(dy/dt\) é proporcional a \(y\) e também a \(x\). Assim, obtemos a equação diferencial:

\[ \begin{align} \dfrac{dy}{dt}&=\beta yx\\ \dfrac{dy}{dt}&=\beta y(n+1-y) \end{align} \]

sendo que \(\beta>0\) é chamada taxa de infecção específica.

A solução geral é:

\[ y(t)=\dfrac{n+1}{1+ke^{-\beta(n+1)t}} \]

Se \(y(0) = 1 = (n + 1)/(1 + k)\), então \(k = n\).

Portanto, a solução particular é:

\[ y(t)=\dfrac{n+1}{1+ne^{-\beta(n+1)t}} \]

Sob nossas suposições simplificadoras, a propagação da doença segue uma lei logística.

O ponto de inflexão é:

\[ \begin{align} t_i&=\dfrac{1}{\beta}\dfrac{\ln(n)}{n+1}\\ y(t_i)&=\dfrac{n+1}{2} \end{align} \]

y'(t)=beta y(t) (n +1 - y(t)), y(0) = 1

plot (100+1)/(1+100 e^(-0.005 (100+1) t)), t from 0 to 20 axes label "t" "y"

Solve[D[(n+1)/(1+ne^(- beta (n+1) t)),{t,2}]==0,t]

(n+1)/(1+ne^(-beta (n+1) t)) /. t = ln(n)/(beta (n+1))

Exemplo 11.5.3: Cinética química

Consideremos a reação entre fluoreto de n-amila e etóxido de sódio:

\[ n\text{--} \mathrm{C}_5\mathrm{H}_{11}\mathrm{F} + \mathrm{NaOC}_2\mathrm{H}_5 \rightarrow \mathrm{NaF} + n\text{--} \mathrm{C}_5\mathrm{H}_{11}\mathrm{OC}_2\mathrm{H}_5 \]

Seja \(A\) a concentração inicial de fluoreto de n-amila e \(B\) a concentração inicial de etóxido de sódio, e suponha-se que \(A \neq B\). Como durante a reação uma molécula do primeiro composto reage com exatamente uma molécula do segundo composto, \(A\) e \(B\) diminuem na mesma quantidade \(x = x(t)\) a qualquer instante \(t\). Assim, as concentrações restantes dos componentes são \(A - x\) e \(B - x\).

Chamamos \(\dfrac{dx}{dt}\) de taxa de reação. Como a reação exige colisões entre moléculas de fluoreto de n-amila com moléculas de etóxido de sódio, é plausível supor que \(\dfrac{dx}{dt}\) seja proporcional ao número de moléculas presentes de ambos os componentes. Isso equivale a dizer que \(\dfrac{dx}{dt}\) é proporcional a \(A - x\) e a \(B - x\). Assim, obtemos a equação diferencial:

\[ \dfrac{dx}{dt} = r(A - x)(B - x) \tag{11.5.19} \]

com uma taxa de reação específica \(r > 0\). Para que \(r\) permaneça constante, é necessário assumir que a temperatura é mantida constante durante a reação. A equação (11.5.19) está na forma da equação (11.5.4) e a solução geral está apresentada em (11.5.9). Para satisfazer a condição inicial \(x(0)=0\), temos \(k=-B/A\). Após rearranjar os termos, a solução particular pode ser escrita da forma:

\[ x(t) = A \left( 1 - \dfrac{A - B}{A - B e^{-(A - B)rt}} \right) \tag{11.5.20} \]

Quando \(t \rightarrow \infty\), \(x\) tende a \(A\), se \(A < B\). Deve-se observar que essa função não é uma função logística.

plot 1 + 1/(1 - 2 e^x), x:0..4, axes labels "t" "x"

Há boa concordância entre o modelo e o experimento. Para uma discussão completa deste exemplo, ver Latham (1964, p. 103).

Exemplo 11.5.4. Autocatálise

Outro tipo de reação química leva a uma equação diferencial relacionada. Considere o processo pelo qual o tripsinogênio é convertido em tripsina (uma enzima). A reação inicia-se apenas na presença de alguma quantidade de tripsina, ou seja, o produto da reação atua como catalisador.

Seja \(y_0\) a concentração inicial de tripsina no instante \(t = 0\), e seja \(y = y(t)\) a concentração adicional adquirida pela reação no instante \(t\), de modo que a concentração total seja \(y_0 + y(t)\). Seja \(B\) a concentração inicial de tripsinogênio. Como cada molécula de tripsinogênio gera uma molécula de tripsina, \(B\) diminui na mesma proporção em que a concentração de tripsina aumenta, isto é, em \(y\). Portanto, no instante \(t\), as concentrações de tripsinogênio e de tripsina são, respectivamente, \(B - y\) e \(y_0 + y\).

É razoável assumir que a taxa de reação \(\dfrac{dy}{dt}\) seja proporcional a \(y_0 + y\) e a \(B - y\). Assim, obtemos a equação diferencial

\[ \dfrac{dy}{dt} = r(y_0 + y)(B - y), \tag{11.5.21} \]

em que \(r\) é uma constante. Para comparar essa equação com a equação (11.5.4), tomamos \(r = -a\) e \(y_0 = -A\), de forma que se obtenha uma concordância formal. A solução geral, conforme indicado em (11.5.9), é:

\[ y(t) = \dfrac{B + y_0}{1 + ke^{-r(B + y_0)t}}-y_0 \tag{11.5.22} \]

Com a condição inicial \(t = 0\), \(y = 0\), conseguimos determinar \(k\). Assim, obtemos uma solução particular com \(k = \dfrac{B}{y_0}\).

A solução particular é:

\[ \begin{align} y(t) &= \dfrac{B + y_0}{1 + \dfrac{B}{y_0}e^{-r(B + y_0)t}}-y_0 \\ y(t) &= y_0\left(\dfrac{y_0+B}{y_0 + Be^{-r(y_0+B)t}}-1\right) \end{align} \]

Pela comparação com (11.5.11), é fácil verificar que o gráfico da equação (11.5.22) é uma curva logística. Conforme mostrado na Fig. 11.6, a concordância com os dados experimentais é satisfatória.

Fig. 11.6. Ativação autocatalítica do tripsinogênio cristalino. A curva ajustada aos pontos observados é uma curva logística. A figura foi reproduzida de Northrop, Kunitz e Herriot (1948, p. 126). Ver também B. Stevens (1965, p. 83).

Exemplo: Modelo geral para o crescimento ontogênico

Crescimento celular:

O crescimento de uma célula depende do fluxo de nutrientes através da sua superfície.

Seja \(W = W(t)\) a massa da célula. Suponha que, por um tempo limitado, a taxa de crescimento \(\dfrac{dW}{dt}\) seja proporcional à área da superfície, i.e., \(\dfrac{dW}{dt} \propto S\).

Se a forma da célula em crescimento não se altera, então a área da superfície, \(S\), é proporcional ao quadrado de uma dimensão linear (como o diâmetro, conforme definido na Seção 4.2), e, portanto, \(S\propto W^{2/3}\). Assim, temos a seguinte equação diferencial:

\[ \dfrac{dW}{dt} = kW^{2/3} \]

sendo que \(k > 0\) é uma constante.

A solução desta equação diferencial é:

\[ W(t) = \dfrac{1}{27}\left(3\sqrt[3]{W_0}+kt\right)^3 \]

x'(t) = kx(t)^(2/3), x(0)=a

A solução indica que a massa da célula aumenta com o cubo de uma função linear do tempo. Isso reflete o fato de que, se o crescimento da célula é impulsionado pela absorção de nutrientes pela superfície, e essa área aumenta com a dimensão linear ao quadrado (e, portanto, com \(W^{2/3}\)), o crescimento resultante do volume (ou massa) segue uma lei cúbica. Este tipo de modelo foi discutido por von Bertalanffy (1951), que propôs equações de crescimento base.

Crescimento animal:

Conforme Bertalanffy (1957), o crescimento animal pode ser considerado como resultado de uma contraposição entre síntese e destruição, do anabolismo e catabolismo dos materiais constituintes do corpo. Haverá crescimento enquanto a construção prevalecer sobre a degradação; o organismo atinge um estado estacionário se e quando ambos os processos forem iguais. Podemos expressar isso em uma fórmula geral:

\[ \dfrac{dW}{dt} = \eta W^m - \kappa W^n \]

Em palavras: A mudança da massa corporal \(W\) é dada pela diferença entre os processos de construção e degradação; \(\eta\) e \(\kappa\) são constantes do anabolismo e catabolismo, respectivamente, e os expoentes \(m\) e \(n\) indicam que estes são proporcionais a alguma potência da massa corporal \(W\).

A solução para \(n = 1\) é (Bertalanffy, 1941b):

\[ W(t) = \left( \dfrac{\eta}{\kappa} - \left( \dfrac{\eta}{\kappa} - W_0^{1-m} \right) e^{-(1-m)\kappa t} \right)^{\frac{1}{1-m}} \]

com \(W_0\) sendo a massa corporal no tempo \(t = 0\).

O expoente \(n\) em nossa equação básica denota a dependência do anabolismo na massa corporal. Se inserirmos para esse valor o que é experimentalmente encontrado para a dependência do metabolismo basal no tamanho, as leis de crescimento para o organismo em questão seguem automaticamente. Assim, podemos prever o tipo de crescimento de um animal a partir de seu tipo metabólico, e essa previsão tem se mostrado correta em um grande número de casos, muitas vezes de maneira bastante inesperada.

Em um primeiro tipo, a taxa respiratória é proporcional à potência \(\dfrac{2}{3}\) da massa corporal, de acordo com a regra da superfície. Consequentemente, a lei do crescimento assume a seguinte forma:

\[ \dfrac{dW}{dt} = \eta W^{2/3} - \kappa W \]

Para o caso específico onde \(m = \dfrac{2}{3}\) e \(n = 1\), a solução da equação diferencial é:

\[ W(t) = \left( \dfrac{\eta}{\kappa} - \left( \dfrac{\eta}{\kappa} - \sqrt[3]{W_0} \right) e^{-\frac{1}{3}\kappa t} \right)^{3} \]

Sendo que:

• \(W_0\) é massa corporal inicial no tempo \(t = 0\).
• \(\eta\) é a constante do anabolismo.
• \(\kappa\) é a constante do catabolismo.

O limite de \(W(t)\) quando \(t\) tende para infinito é:

\[ \lim_{t \to \infty} W(t) = \left( \dfrac{\eta}{\kappa} \right)^{3} \]

Este é o valor assintótico superior para a massa corporal \(W(t)\) quando \(t\) tende para infinito.

Conforme West et al. (2001) e Savage et al. (2007), o expoente \(\dfrac{3}{4}\) em vez de \(\dfrac{2}{3}\) é bem apoiado por dados sobre mamíferos, aves, peixes, moluscos e plantas.

\[ \dfrac{dm}{dt}=a\; m^{3/4}\left(1-\left(\dfrac{m}{M}\right)^{1/4}\right) \]

Sendo que \(m=m(t)\) é massa corportal no instante \(t\), \(m_0\) é a massa corporal ao nascer em \(t=0\) e \(M\) é a massa corporal máxima.

Solução:

\[ \begin{align} \sqrt[4]{\dfrac{m(t)}{M}} &= 1-\left(1-\sqrt[4]{\dfrac{m_0}{M}}\right)\exp\left(-\dfrac{a}{4\sqrt[4]{M}}t\right) \\ 1-\sqrt[4]{\dfrac{m(t)}{M}} &= \left(1-\sqrt[4]{\dfrac{m_0}{M}}\right)\exp\left(-\dfrac{a}{4\sqrt[4]{M}}t\right) \\ \\ R(t) &= 1-\sqrt[4]{\dfrac{m(t)}{M}} \\ R(0) &= 1-\sqrt[4]{\dfrac{m_0}{M}} \\ \\ R(t) &= R(0)\exp\left(-\dfrac{a}{4\sqrt[4]{M}}t\right) \\ \dfrac{R(t)}{R(0)} &= \exp\left(-\dfrac{a}{4\sqrt[4]{M}}t\right) \\ \ln\left(\dfrac{R(t)}{R(0)}\right) &= \ln\left(\exp\left(-\dfrac{a}{4\sqrt[4]{M}}t\right)\right) \\ \ln\left(\dfrac{R(t)}{R(0)}\right) &= -\dfrac{a}{4\sqrt[4]{M}}t \\ \end{align} \]

\(R(t)\) é a proporção do poder metabólico devotado ao crescimento do animal.

Se \(y=\ln\left(\dfrac{R(t)}{R(0)}\right)\) e \(x=-\dfrac{a}{4\sqrt[4]{M}}t\), a inclinação é unitário em valor absoluto da regressão linear \(y=\alpha+\beta x+\varepsilon\).

\[ \begin{align} r(t)&=1-R(t) \quad \text{: massa normalizada adimensional}\\ \tau(t)&=\dfrac{a}{4\sqrt[4]{M}}t-R(0) \quad \text{: tempo adimensional}\\\\ r(t)&=1-e^{-\tau(t)}\quad \text{: curva de crescimento de qualquer animal} \end{align} \]

Lei alométrica

Allometry: Wikipedia

É fácil observar que, na criança, a cabeça cresce mais lentamente que o corpo, e os olhos crescem ainda mais lentamente que a cabeça. Esse padrão diferencial de crescimento é estudado pela alometria.

Seja \(x = x(t)\) o tamanho (comprimento, volume ou massa) de um órgão e \(y= y(t)\) o tamanho de outro órgão ou parte no mesmo indivíduo no instante de tempo \(t\).

As taxas de crescimento dessas duas partes são \(dx/dt\) e \(dy/dt\).

Quando afirmamos que uma parte cresce mais lentamente que outra, não estamos comparando apenas os incrementos absolutos por unidade de tempo, isto é, não basta observar se \(dx/dt \ne dy/dt\). Isso porque uma estrutura menor tende, naturalmente, a apresentar menor ganho absoluto do que uma estrutura maior. A comparação relevante é entre as taxas relativas de crescimento, isto é,

\[ \frac{1}{x}\frac{dx}{dt} \quad \text{e} \quad \frac{1}{y}\frac{dy}{dt} \]

Assim, dizemos que \(x\) cresce mais lentamente que \(y\) quando sua taxa relativa de crescimento é menor:

\[ \frac{1}{x}\frac{dx}{dt} < \frac{1}{y}\frac{dy}{dt} \]

No contexto da alometria, portanto, o interesse não está apenas em quanto cada parte aumenta, mas em como o crescimento de uma parte se relaciona proporcionalmente ao crescimento de outra.

Numerosos dados empíricos apoiam a afirmação de que as taxas específicas de crescimento (elasticidades instantâneas) são aproximadamente proporcionais, i.e., que a equação funciona com precisão satisfatória:

\[ \begin{align} \eta(y(t))&\propto \eta(x(t)) \\ \eta(y(t))&=k\; \eta(x(t)) \\ t\dfrac{dy/dt}{y}&=kt\dfrac{dx/dt}{x} \\ \dfrac{dy/y}{dt/t}&=k\dfrac{dx/x}{dt/t} \\ \dfrac{dy/y}{dt}&=k\dfrac{dx/x}{dt} \\ \dfrac{1}{y}\dfrac{dy}{dt}&=k\dfrac{1}{x}\dfrac{dx}{dt} \end{align} \]

A constante \(k>0\) depende apenas da natureza dos dois órgãos ou partes em consideração.

A equação anterior é chamada de lei alométrica.

A equação da lei alométrica pode ser simplificada eliminando o tempo:

\[ \begin{align} \dfrac{dy}{dx}&= k \;\dfrac{y}{x}\\ x&, \; y, \; k>0 \\ \\ x\dfrac{dy/dx}{y}&= k \\ \\ \eta(y(x))&= k \\\\ \ln(y)&=k\ln(x)+C\\ y&= c x^k \end{align} \]

A lei alométrica pode ser expressa pela equação diferencial, pela função de potência ou pela equação linear com logaritmo. As três expressões são equivalentes.

Observe que a lei alométrica não se refere primariamente a velocidade de crescimento, pois o tempo é eliminado. Um indivíduo pode crescer em função do tempo seguindo uma lei exponencial, logística ou qualquer outra. Isso não afeta a relação alométrica.

A lei alométrica foi aplicada com sucesso não apenas ao crescimento relativo de partes de um corpo, mas também ao metabolismo, a problemas de dose-resposta, a diferenças raciais e à história evolutiva.

Existem também tentativas de justificar teoricamente a lei alométrica.

Forame magno: Em anatomia, foramen magnum ou forame magno é a grande abertura através do osso occipital localizada no centro da fossa posterior do neurocrânio. É o maior dos forâmens do crânio e comunica a cavidade craniana com o canal vertebral.

Fig. 11.7. Relação entre volume endocraniano e área do forame magno em cinco ordens de mamíferos. As letras referem-se às retas ajustadas aos dados de cada grupo. I insetívoros, R roedores, P primatas prossímios, M macacos do novo e velho mundo combinados, C carnívoros fissípedes (lobo, leão, urso-pardo, guaxinim etc.), A artiodáctilos (porco, javali, camelo, lhama, ruminantes).

dy/dx = k y/x

Sistema de equações diferenciais ordinárias lineares

Em sistemas fisiológicos, assim como em populações, geralmente há duas ou até mais funções do tempo em consideração.

Sejam \(x = x(t)\) e \(y = y(t)\) duas dessas funções do tempo.

Assumimos algum tipo de interação entre as quantidades \(x\) e \(y\).

Assim, as taxas de variação \(dx/dt\) e \(dy/dt\) podem depender das quantidades \(x\) e \(y\).

Um caso relativamente simples consiste nas duas equações simultâneas:

\[ \begin{align} x^{\prime}&=ax+by \\ y^{\prime}&= cx+dy \end{align} \]

sendo que \(a\), \(b\), \(c\), \(d\) são constantes reais dadas.

As funções desconhecidas \(x = x(t)\) e \(y = y(t)\) e suas derivadas aparecem na primeira potência.

Portanto, o conjunto das duas equações diferenciais lineares é dito ser um sistema de equações diferenciais lineares.

Vale a pena estudar o conteúdo lógico de tal sistema.

As grandezas \(x\) e \(y\) interagem de tal forma que determinam totalmente as taxas \(dx/dt\) e \(dy/dt\).

Com o passar do tempo, \(x\) e \(y\) aumentarão ou diminuirão de acordo com as taxas \(dx/dt\) e \(dy/dt\).

Assim, as taxas determinarão os valores futuros de \(x\) e \(y\).

A informação originalmente fornecida para \(dx/dt\) e para \(dy/dt\) é então realimentada para \(x\) e \(y\).

Obtemos dois loops de feedback, conforme indicado na seguinte apresentação gráfica:

O sistema é resolvido pelo método das funções-teste.

Tentamos as seguintes funções especiais:

\[ \begin{align} x&=Ae^{\lambda t}\\ y&=Be^{\lambda t} \end{align} \]

sendo que \(A\), \(B\), \(\lambda\) são algumas constantes a serem determinadas posteriormente.

Para evitar soluções triviais, podemos assumir que \(A \ne 0\) e \(B \ne 0\).

A equação característica é:

\[ \lambda^2+\left(a+d\right)\lambda+(ad-bc)=0 \]

lambda^2-(a+d) lambda+(ad-bc)=0

determinant {{a-lambda, b}, {c, d-lambda}}=0

Portanto, se \(\lambda_1 \ne \lambda_2\):

\[ x(t)=A_1e^{\lambda_1 t}+A_2e^{\lambda_2 t} \\ y(t)=B_1e^{\lambda_1 t}+B_2e^{\lambda_2 t} \]

Exemplo: Teórico

\[ x^{\prime}=3x-2y \\ y^{\prime}= 2x-2y \]

\[ x(t)=ke^{- t}-2me^{2 t} \\ y(t)=2ke^{- t}-me^{2 t} \]

As soluções são explosivas (tendem para infinito).

eigenvalues ((3, -2), (2, -2))

{x'(t) = 3 x(t) - 2 y(t), y' = 2 x(t) - 2 y(t)}

eigenvalues ((a, b), (c, d))

{x'(t) = a x(t) + b y(t), y' = c x(t) + d y(t)}

eigenvalues ((1, 1), (1, 1))

{x'(t) = x(t) + y(t), y' = x(t) + y(t)}

eigenvalues ((1, 1), (1, -1))

{x'(t) = x(t) + y(t), y' = x(t) - y(t)}

eigenvalues ((1, -1), (1, 1))

{x'(t) = x(t) - y(t), y' = x(t) + y(t)}

eigenvalues ((1, -1), (1, -1))

{x'(t) = x(t) - y(t), y' = x(t) - y(t)}

Exemplo 11.7.2: Passagem de alimentos em ruminantes

Há animais ruminantes e semi-ruminantes. Os ruminantes, como veados, ovelhas, cabras e bois, possuem um estômago complexo, dividido em compartimentos especializados, que permite a fermentação do alimento e a digestão eficiente de material vegetal rico em celulose.

Os semi-ruminantes, por sua vez, também realizam fermentação digestiva, mas apresentam um sistema gástrico menos compartimentalizado que o dos ruminantes típicos. Exemplos incluem camelos, lhamas e alpacas.

Alimentos recém-comidos, mas não mastigados, são passados para um compartimento de armazenamento chamado rúmen.

Mais tarde, quando mastigado, o alimento passa pelo abomaso para o duodeno, onde é processado posteriormente. A partir daí, ele entra lentamente nos intestinos. Para obter uma descrição matemática da passagem do alimento pelo trato digestivo, Blaxter, Graham e Wainman (1956) propuseram o seguinte modelo:

Seja \(r = r(t)\) a quantidade de alimento no rúmen no instante \(t\). Em \(t = 0\) esta quantidade é uma quantidade conhecida, digamos \(r_0\).

Por \(u = u(t)\) denotamos a quantidade de alimento no abomaso no instante \(t\). Em \(t = 0\) temos \(u= 0\). A taxa de diminuição de \(r= r(t)\) é provavelmente proporcional a \(r\).

Também é razoável supor que \(du/dt\) consiste em dois termos: uma taxa de aumento igual à taxa de diminuição de \(r\) e uma taxa de diminuição proporcional a \(u\). Assim obtemos as equações diferenciais lineares:

\[ \begin{align} \dfrac{dr}{dt}&=-k_1r \\ \dfrac{du}{dt}&=k_1r-k_2u \end{align} \]

sendo que \(k_1\) e \(k_2\) são constantes positivas. Eles podem ser chamados de taxas específicas de digestão. Assumimos que \(k_1 \ne k_2\). A primeira equação contém um sinal de menos, pois \(r\) diminui.

Além disso, denotamos a quantidade total de alimento que entrou no duodeno no tempo \(t\) por \(v= v(t)\) e a quantidade de fezes por \(w= w(t)\).

Como o duodeno recebe exatamente a mesma quantidade que sai do abomaso, obtemos a equação diferencial

\[ \dfrac{dv}{dt}=k_2u \]

sendo que \(v= 0\) no tempo \(t = 0\).

Finalmente, assumimos que as fezes deixam o animal com um atraso de tempo constante que denotamos por \(\tau\).

Desprezando a perda de matéria que entra nos vasos sanguíneos, obtemos \(w(t)\) no tempo \(t - \tau\):

\[ w(t) = r_0 - r (t - \tau) - u (t - \tau), \; t> \tau \]

O processo geral pode ser resumido graficamente por:

A equação característica do sistema \(\{r^{\prime}, u^{\prime}\}\) é:

\[ \lambda^2+\left(k_1+k_2\right)\lambda+k_1k_2=0 \]

de modo que \(\lambda_1=-k_1\) e \(\lambda_2=-k_2\) onde por suposição \(\lambda_1 \ne \lambda_2\).

A solução geral pode ser escrita na forma:

\[ r(t)=A_1e^{-k_1t}+A_2e^{-k_2t} \\ u(t)=B_1e^{-k_1t}+B_2e^{-k_2t} \]

r'(t)= -k_1 r(t), u'(t) = k_1 r(t) - k_2 u(t)

A solução particular é, se \(r(0)=r_0, u(0)=0\):

\[ \begin{align} r(t)&=r_0\;e^{-k_1t} \\ u(t)&=r_0\dfrac{k_1}{k_2-k_1}\left(e^{-k_1t}-e^{-k_2t}\right)\\ w(t) &= r_0 - r (t - \tau) - u (t - \tau), \quad t>\tau \end{align} \]

Modelamos o trato digestivo como um sistema de três compartimentos que pode ser interpretado como um exemplo de análise de compartimento.

Fig. 11.8. A passagem do alimento pelo trato digestivo dos ruminantes. A quantidade de alimento no rúmen e no abomaso é denotada por r(t), u(t), respectivamente. w(t) é a quantidade de fezes.

Plot[{100 E^(-0.1 t), 100 (0.1 / (0.2-0.1)) (E^(-0.1 t) - E^(-0.2 t)), 100 - 100 E^(-0.1 (t-10)) - 100 (0.1 / (0.2-0.1)) (E^(-0.1 (t-10)) - E^(-0.2 (t-10)))}, {t, 0, 50}]

Modelo compartimental de doença infecciosa

Compartmental models (epidemiology): Wikipedia

Modelos compartimentais são modelos matemáticos usados para descrever a dinâmica de populações cujos indivíduos são agrupados em categorias, chamadas compartimentos, segundo alguma característica de interesse. Em epidemiologia, por exemplo, é comum classificar os indivíduos como suscetíveis, infectados e recuperados. Em ecologia, os compartimentos podem representar classes etárias ou estágios de desenvolvimento. Em farmacocinética, podem representar regiões do organismo entre as quais uma substância é distribuída.

A ideia central é simples: em vez de acompanhar cada indivíduo separadamente, acompanha-se o número, ou a proporção, de indivíduos em cada compartimento ao longo do tempo. O foco passa a ser o fluxo entre compartimentos.

Se \(X_1(t), X_2(t), \dots, X_k(t)\) denotam os compartimentos no instante \(t\), então um modelo compartimental descreve como essas quantidades variam no tempo por meio de equações diferenciais ordinárias, como

\[ \frac{dX_j}{dt} = \text{entradas no compartimento } j - \text{saídas do compartimento } j\\ \qquad j=1,\dots,k \]

Cada equação expressa um balanço: a taxa de variação de um compartimento é dada pela diferença entre o que entra e o que sai. Essa estrutura torna os modelos compartimentais intuitivos e úteis para representar processos de transmissão, crescimento, decaimento, migração e troca entre classes.

Um exemplo clássico é o modelo epidemiológico SI, com dois compartimentos: \[ S(t) = \text{suscetíveis}, \qquad I(t) = \text{infectados} \]

Parâmetros epidemiológicos do modelo SIS

Sem dinâmica vital, uma forma simples do modelo é \[ \begin{align} \frac{dS}{dt} &= -\beta S I\\ \frac{dI}{dt} &= \beta S I \end{align} \] sendo que \(\beta>0\) é o parâmetro de transmissão. Nesse caso, indivíduos saem do compartimento dos suscetíveis e entram no dos infectados.

Uma propriedade importante de muitos modelos compartimentais é a conservação da população total. Isso significa que os indivíduos apenas mudam de compartimento, sem entrada ou saída da população.

No modelo SI sem dinâmica vital, há dois compartimentos: suscetíveis \(S(t)\) e infectados \(I(t)\). A população total é

\[ N(t)=S(t)+I(t) \]

Um modelo SI simples é

\[ \frac{dS}{dt}=-\beta\frac{S I}{N} \]

\[ \frac{dI}{dt}=\beta\frac{S I}{N} \]

A primeira equação diz que os suscetíveis diminuem à taxa \(\beta SI/N\). A segunda equação diz que os infectados aumentam exatamente na mesma taxa. Portanto, o que sai de \(S\) entra em \(I\).

Para verificar isso matematicamente, derivamos \(N(t)\):

\[ \frac{dN}{dt} = \frac{d}{dt}\{S(t)+I(t)\} \]

Logo,

\[ \frac{dN}{dt} = \frac{dS}{dt} + \frac{dI}{dt} \]

Substituindo as duas equações do modelo, obtemos

\[ \frac{dN}{dt} = -\beta\frac{S I}{N} + \beta\frac{S I}{N} \]

Os dois termos são iguais em módulo e opostos em sinal. Portanto, eles se cancelam:

\[ -\beta\frac{S I}{N} + \beta\frac{S I}{N} = 0 \]

Assim,

\[ \frac{dN}{dt}=0 \]

O zero não aparece por hipótese adicional. Ele aparece porque a perda de suscetíveis é exatamente igual ao ganho de infectados.

Como a derivada de \(N(t)\) é igual a zero, \(N(t)\) não varia com o tempo. Portanto,

\[ N(t)=N(0) \]

Isso não significa que \(S(t)\) e \(I(t)\) sejam constantes. Ao contrário, \(S(t)\) diminui e \(I(t)\) aumenta. A população total permanece constante porque cada indivíduo que sai do compartimento \(S\) entra no compartimento \(I\). Assim, a perda em \(S\) é exatamente igual ao ganho em \(I\).

Assim, modelos compartimentais combinam simplicidade conceitual e utilidade prática. Eles permitem transformar um processo real em um sistema matemático capaz de descrever, explicar e prever a evolução temporal de populações ou substâncias distribuídas entre compartimentos.

Conforme Naziazeno et al. (2024),

“Neste trabalho, discute-se a importância da modelagem matemática para compreender e prevenir pandemias, com foco em um modelo compartimental do tipo Suscetível-Infectado-Suscetível-Infectado-Recuperado. Trata-se de uma extensão do modelo SIR clássico, que incorpora a possibilidade de reinfecção, estágios intermediários de suscetibilidade e a competição entre diferentes variantes virais. O estudo propõe e analisa esse modelo, investigando propriedades como a existência de equilíbrios em torno de pontos críticos e o cálculo do número básico de reprodução. Além disso, o modelo é aplicado à cidade de Manaus em um período crítico da pandemia de COVID-19, permitindo avaliar sua utilidade na compreensão da dinâmica da doença e na formulação de estratégias de controle. Os resultados indicam que o modelo amplia o entendimento da epidemiologia da COVID-19 ao fornecer subsídios para o planejamento de intervenções mais direcionadas e ajustadas às características das variantes circulantes. Também contribui, em termos didáticos, para a discussão de sistemas dinâmicos e para a compreensão do papel da modelagem matemática em saúde pública.”

Unidade e interpretação dos parâmetros \(\beta\), \(\gamma\) e \(\mu\)

A unidade de \(\beta\) depende de como o termo de transmissão é escrito. A regra prática é verificar a unidade do lado esquerdo da equação e exigir que o lado direito tenha a mesma unidade.

Se \(S\) e \(I\) são números de indivíduos e o termo de transmissão é

\[ \beta SI \]

então \(\beta SI\) deve ter unidade de indivíduos por unidade de tempo. Como \(SI\) tem unidade de indivíduos ao quadrado, segue que

\[ [\beta]=\frac{1}{\text{indivíduo}\cdot T} \]

Nesse caso, \(\beta\) é uma taxa de transmissão por indivíduo por unidade de tempo.

Se \(S\) e \(I\) são proporções, então são adimensionais. No modelo

\[ \frac{dI}{dt}=\beta SI \]

o produto \(SI\) é adimensional, e a derivada tem unidade de inverso do tempo. Portanto,

\[ [\beta]=T^{-1} \]

Outra forma comum é usar números absolutos, mas padronizar pela população total:

\[ \beta\frac{SI}{N} \]

Nesse caso, \(SI/N\) tem unidade de indivíduo. Logo, para que o termo represente indivíduos por unidade de tempo, temos novamente

\[ [\beta]=T^{-1} \]

Portanto, a unidade de \(\beta\) não é fixa. Ela depende da forma usada para escrever a incidência.

No caso normalizado, em que

\[ S(t)+I(t)=1 \]

as variáveis \(S(t)\) e \(I(t)\) representam proporções. Assim, \(\beta\), \(\gamma\) e \(\mu\) têm a mesma unidade matemática:

\[ [\beta]=[\gamma]=[\mu]=T^{-1} \]

Apesar disso, os significados epidemiológicos são diferentes. O parâmetro \(\beta\) mede a intensidade da transmissão, \(\gamma\) mede a taxa de recuperação e \(\mu\) mede a taxa de natalidade ou mortalidade natural.

Por exemplo, em um modelo SIS normalizado,

\[ \frac{dI}{dt}=\beta SI-\gamma I-\mu I \]

o termo \(\beta SI\) representa a entrada de novos infectados, o termo \(\gamma I\) representa a saída por recuperação e o termo \(\mu I\) representa a saída por mortalidade natural.

Suponha tempo medido em dias e

\[ \beta=0.8\;\text{dia}^{-1} \]

\[ \gamma=0.3\;\text{dia}^{-1} \]

\[ \mu=0.01\;\text{dia}^{-1} \]

Se, em certo instante,

\[ S=0.9 \]

\[ I=0.1 \]

então a transmissão instantânea é

\[ \beta SI=0.8\times 0.9\times 0.1=0.072\;\text{dia}^{-1} \]

A recuperação é

\[ \gamma I=0.3\times 0.1=0.03\;\text{dia}^{-1} \]

A mortalidade natural entre infectados é

\[ \mu I=0.01\times 0.1=0.001\;\text{dia}^{-1} \]

Assim,

\[ \frac{dI}{dt}=0.072-0.03-0.001=0.041\;\text{dia}^{-1} \]

Isso significa que, naquele instante, a proporção de infectados está aumentando a uma taxa de 0.041 por dia, isto é, 4.1 pontos percentuais por dia.

Os inversos de \(\gamma\) e \(\mu\) têm interpretação direta como escalas médias de tempo. Por exemplo,

\[ \frac{1}{\gamma}=\frac{1}{0.3}\approx 3.33\;\text{dias} \]

indica uma escala média de recuperação de aproximadamente 3.33 dias.

Também,

\[ \frac{1}{\mu}=\frac{1}{0.01}=100\;\text{dias} \]

indica que a dinâmica vital atua em escala mais lenta.

Para \(\beta\), a interpretação do inverso é menos direta. Como a transmissão efetiva depende de \(S\) e \(I\), a quantidade relevante é

\[ \beta SI \]

No exemplo,

\[ \frac{1}{\beta SI}=\frac{1}{0.072}\approx 13.89\;\text{dias} \]

Essa é a escala temporal associada à transmissão efetiva naquele estado da população. Portanto, \(\beta\) mede a intensidade potencial de transmissão, mas a velocidade real da infecção depende também dos valores atuais de \(S\) e \(I\).

Por que taxas iguais a \(\mu\) para nascimento e morte?

Usa-se o mesmo parâmetro \(\mu\) para nascimento e morte porque o modelo supõe população total constante.

Os nascimentos entram como

\[ \mu N \]

e as mortes naturais saem como

\[ -\mu S \qquad \text{e} \qquad -\mu I \]

Assim,

\[ \frac{dN}{dt}=\mu N-\mu(S+I) \]

Como

\[ N=S+I \]

segue que

\[ \frac{dN}{dt}=0 \]

Portanto, \(\mu\) representa nascimento e morte apenas por hipótese de simplificação, impondo que a taxa per capita de natalidade seja igual à taxa per capita de mortalidade.

Se essas taxas fossem diferentes, o correto seria usar parâmetros distintos.

Modelos SI e SIS com e sem dinâmica vital

Modelos compartimentais descrevem a dinâmica de infecções ao dividir a população em classes epidemiológicas. Nos modelos SI e SIS, a diferença central é simples: no modelo SI o indivíduo infectado não retorna ao estado suscetível; no modelo SIS, retorna. A presença ou ausência de dinâmica vital altera a reposição de suscetíveis e a possibilidade de persistência endêmica.

Esses modelos não descrevem o patógeno em si, mas a dinâmica populacional da infecção causada por ele. Em humanos, são úteis para doenças respiratórias, dermatológicas, sexualmente transmissíveis ou hospitalares, dependendo da presença ou ausência de imunidade e da duração da infecção. Em animais não humanos, são aplicáveis a infecções em rebanhos, animais de companhia, vetores e populações silvestres, sobretudo quando o interesse é estudar persistência, surtos, reinfecção e manutenção endêmica.

A escolha entre SI e SIS depende da biologia do hospedeiro e do patógeno. Se a infecção é praticamente permanente no horizonte estudado, SI pode ser suficiente. Se o hospedeiro se recupera e volta a ser suscetível, SIS é mais apropriado. A inclusão de dinâmica vital é necessária quando nascimento e morte têm impacto relevante na escala temporal da infecção.

Em termos epidemiológicos, o ponto central é este: SI descreve acúmulo irreversível de infectados; SIS descreve circulação contínua entre suscetibilidade e infecção. Sem dinâmica vital, a população é fechada; com dinâmica vital, a população se renova. Essa combinação determina se a doença desaparece, cresce ou se estabiliza em equilíbrio endêmico.

Embora os modelos SI e SIS sejam frequentemente resolvidos por meio de uma única equação diferencial, eles são, em sua formulação original, sistemas de duas equações diferenciais ordinárias não lineares.

Por exemplo, no modelo SI sem dinâmica vital,

\[ \begin{align} \frac{dS}{dt} &= -\beta SI \\ \frac{dI}{dt} &= \beta SI \end{align} \]

há duas variáveis de estado, \(S(t)\) e \(I(t)\), e duas equações diferenciais acopladas. Além disso, o termo \(SI\) torna o sistema não linear, pois envolve o produto entre as variáveis dependentes.

De modo análogo, no modelo SIS sem dinâmica vital,

\[ \begin{align} \frac{dS}{dt} &= -\beta SI + \gamma I \\ \frac{dI}{dt} &= \beta SI - \gamma I \end{align} \]

também temos um sistema de duas EDOs não lineares.

Portanto, do ponto de vista matemático, SI e SIS são sistemas dinâmicos bidimensionais não lineares.

No entanto, em muitos casos existe uma relação de conservação da população total. Se não há nascimentos, mortes nem migração, então

\[ S(t)+I(t)=N \]

ou, na forma normalizada,

\[ S(t)+I(t)=1 \]

Essa relação permite eliminar uma das variáveis. Por exemplo, se

\[ S(t)+I(t)=1 \]

então

\[ S(t)=1-I(t) \]

e o sistema pode ser reduzido a uma única EDO em \(I(t)\), ou equivalentemente em \(S(t)\).

No caso do modelo SI sem dinâmica vital, a substituição leva a

\[ \frac{dI}{dt}=\beta(1-I)I \]

que é uma equação diferencial não linear de primeira ordem.

Assim, há duas descrições matemáticas equivalentes:

1. o modelo original como sistema de duas EDOs não lineares;
2. o modelo reduzido como uma única EDO não linear obtida pela relação de conservação.

Isso não muda a natureza do problema. O modelo continua sendo originalmente um sistema não linear com duas variáveis de estado, mas sua estrutura permite uma redução de dimensão, o que facilita a resolução analítica e a interpretação do comportamento dinâmico.

SI sem dinâmica vital

No modelo compartimental SI sem dinâmica vital, a população é dividida em dois compartimentos: suscetíveis, \(S(t)\), e infectados, \(I(t)\). Como não há nascimentos, mortes nem migração, a população total permanece constante ao longo do tempo.

Assim, \[ N(t)=S(t)+I(t) \]

Além disso, pelas equações do modelo, \[ \begin{align} \frac{dS}{dt} &= -\beta S I\\ \frac{dI}{dt} &= \beta S I \end{align} \]

SI sem dinâmica vital

Somando as duas equações, obtemos \[ \frac{dN}{dt}=\frac{dS}{dt}+\frac{dI}{dt}=0 \]

Portanto, \[ N(t)=N(0) \] isto é, a população total é constante.

Muitos livros simplificam o modelo normalizando a população total inicial para 1. Nesse caso, \[ N(0)=1 \]

Como \(N(t)\) é constante, segue que \[ N(t)=1 \quad \text{para todo } t \]

Isso não significa que exista apenas um indivíduo na população. Significa apenas que \(S(t)\) e \(I(t)\) passam a representar proporções da população total, e não números absolutos. Assim,

\[ S(t)+I(t)=1 \]

Logo, \[ I(t)=1-S(t) \] e a equação para os suscetíveis pode ser escrita como \[ \frac{dS}{dt}=-\beta S(1-S) \]

A solução é uma função logística

\[ I(t)=\frac{I_0e^{\beta t}}{1+I_0\left(e^{\beta t}-1\right)} \] e

\[ S(t)=\frac{1-I_0}{1+I_0\left(e^{\beta t}-1\right)} \]

Portanto, a escrita \(N(t)=1\) faz sentido no modelo SI sem dinâmica vital, desde que se entenda que a população foi normalizada. Em rigor, primeiro se impõe \(N(0)=1\); depois, como \(N(t)\) é constante, conclui-se que \(N(t)=1\) para todo \(t\).

y'(t) = \beta (1 - y(t)) y(t)

\[ y(t) = \dfrac{e^{β t}}{c_1 + e^{β t}} \]

y0 = 1/(c + 1) for c

\[ c_1 = \dfrac{1}{y_0} - 1\\y_0 \ne 0 \] y(t) = e^(β t)/(1/y0 - 1 + e^(β t))

\[ \begin{align} y(t)&=\frac{y_0 e^{\beta t}}{1+y_0\left(e^{\beta t}-1\right)}\\ I(t)&=\frac{I_0e^{\beta t}}{1+I_0\left(e^{\beta t}-1\right)} \end{align} \]

I0 <- 0.01
t  <- seq(0, 10, by = 0.01)

beta1 <- 0.4
beta2 <- 0.8
beta3 <- 1.2

S1 <- (1 - I0) / (1 + I0 * (exp(beta1 * t) - 1))
S2 <- (1 - I0) / (1 + I0 * (exp(beta2 * t) - 1))
S3 <- (1 - I0) / (1 + I0 * (exp(beta3 * t) - 1))

titulo <- substitute(
  S(t) == frac(1 - I[0], 1 + I[0] * (e^(beta * t) - 1)) ~ "," ~ I[0] == i0,
  list(i0 = I0)
)

plot(t, S1, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo)

lines(t, S2, lwd = 2, lty = 2)
lines(t, S3, lwd = 2, lty = 3)

legend("right",
       legend = c(expression(beta == 0.4),
                  expression(beta == 0.8),
                  expression(beta == 1.2)),
       lwd = 2,
       lty = c(1, 2, 3),
       bty = "n")

SI com dinâmica vital

No modelo compartimental SI com dinâmica vital, a população é dividida em dois compartimentos: suscetíveis, \(S(t)\), e infectados, \(I(t)\). Diferentemente do modelo SI sem dinâmica vital, aqui são considerados nascimentos e mortes naturais.

Assumindo taxa per capita de natalidade e mortalidade iguais a \(\mu > 0\), e taxa de transmissão \(\beta > 0\), o modelo pode ser escrito como

\[ \begin{align} \frac{dS}{dt} &= \mu N - \beta S I - \mu S \\ \frac{dI}{dt} &= \beta S I - \mu I \end{align} \]

sendo que \[ N(t)=S(t)+I(t) \]

SI com dinâmica vital

Somando as duas equações, obtemos

\[ \begin{align} \frac{dN}{dt} &= \frac{dS}{dt}+\frac{dI}{dt} \\ &= \mu N - \beta S I - \mu S + \beta S I - \mu I \\ &= \mu N - \mu(S+I) \\ &= \mu N - \mu N \\ \frac{dN}{dt}&= 0 \end{align} \]

Portanto, a população total permanece constante. Se normalizarmos a população total para 1, então

\[ N(t)=S(t)+I(t)=1 \]

e, portanto,

\[ S(t)=1-I(t) \]

Substituindo em \(\dfrac{dI}{dt}\), segue que

\[ \begin{align} \frac{dI}{dt} &= \beta S I - \mu I \\ &= \beta (1-I)I - \mu I \\ &= (\beta-\mu)I - \beta I^2 \\ \frac{dI}{dt}&= I \left[(\beta-\mu)-\beta I\right] \end{align} \]

Se \(\beta-\mu \neq 0\), a solução com condição inicial \(I(0)=I_0\) é uma função logística

\[ I(t)= \frac{1-\frac{\mu}{\beta}}{1+\left(\dfrac{1-\frac{\mu}{\beta}}{I_0}-1\right)e^{-(\beta-\mu)t}} \]

ou, equivalentemente

\[ I(t)=\frac{(\beta-\mu)I_0e^{(\beta-\mu)t}}{\beta-\mu+\beta I_0\left(e^{(\beta-\mu)t}-1\right)} \]

e

\[ S(t)=\frac{\beta-\mu+I_0\left(\mu e^{(\beta-\mu)t}-\beta\right)}{\beta-\mu+\beta I_0\left(e^{(\beta-\mu)t}-1\right)} \]

Do ponto de vista qualitativo, o comportamento do sistema depende da comparação entre \(\beta\) e \(\mu\).

Se \[ \beta > \mu \] a infecção pode persistir na população.

Se \[ \beta \leq \mu \] a mortalidade natural domina a transmissão, e a fração de infectados tende a zero.

Os pontos de equilíbrio são obtidos impondo

\[ \frac{dS}{dt}=0 \qquad \text{e} \qquad \frac{dI}{dt}=0 \]

Para os infectados, os equilíbrios são

\[ I^*=0 \qquad \text{e} \qquad I^*=1-\frac{\mu}{\beta} \]

Para os suscetíveis, os equilíbrios são

\[ S^*=1 \qquad \text{e} \qquad S^*=\frac{\mu}{\beta} \]

Interpretação:

se \[ \beta \leq \mu \] o único equilíbrio biologicamente viável é \[ I^*=0 \qquad \text{e} \qquad S^*=1 \]

se \[ \beta > \mu \] existe equilíbrio endêmico, dado por

\[ S^*=\frac{\mu}{\beta} \qquad \text{e} \qquad I^*=1-\frac{\mu}{\beta} \]

Portanto, no modelo SI com dinâmica vital, a persistência da infecção depende de a taxa de transmissão superar a taxa de mortalidade natural.

Assim, o modelo SI com dinâmica vital preserva a estrutura compartimental do modelo SI simples, mas incorpora reposição de suscetíveis por nascimentos e remoção de indivíduos por mortalidade natural, o que altera de modo importante a dinâmica de longo prazo.

t  <- seq(0, 10, by = 0.01)
S0 <- 0.99
I0 <- 0.01

sol_SI <- function(beta, mu, t, I0) {
  r <- beta - mu
  K <- r / beta
  
  I <- K / (1 + ((K - I0) / I0) * exp(-r * t))
  S <- 1 - I
  
  list(S = S, I = I)
}

# Caso 1: beta <= mu
beta1 <- 0.4
mu1   <- 0.8

out1 <- sol_SI(beta1, mu1, t, I0)

S1e <- 1
I1e <- 0

# Caso 2: beta > mu
beta2 <- 1.2
mu2   <- 0.4

out2 <- sol_SI(beta2, mu2, t, I0)

S2e <- mu2 / beta2
I2e <- 1 - mu2 / beta2

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 4, 1))

titulo1 <- substitute(
  beta <= mu ~ "," ~ beta == b ~ "," ~ mu == m,
  list(b = beta1, m = mu1)
)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo1)

lines(t, out1$I, lwd = 2, lty = 2)
abline(h = S1e, lwd = 1.5, lty = 3)
abline(h = I1e, lwd = 1.5, lty = 4)

legend("right",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

titulo2 <- substitute(
  beta > mu ~ "," ~ beta == b ~ "," ~ mu == m,
  list(b = beta2, m = mu2)
)

plot(t, out2$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo2)

lines(t, out2$I, lwd = 2, lty = 2)
abline(h = S2e, lwd = 1.5, lty = 3)
abline(h = I2e, lwd = 1.5, lty = 4)

legend("right",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

par(oldpar)

t  <- seq(0, 10, by = 0.01)
S0 <- 0.99
I0 <- 0.01

sol_SI <- function(beta, mu, t, I0) {
  S <- ((beta - mu) + I0 * (mu * exp((beta - mu) * t) - beta)) /
    ((beta - mu) + beta * I0 * (exp((beta - mu) * t) - 1))
  
  I <- 1 - S
  
  list(S = S, I = I, Se = mu / beta, Ie = 1 - mu / beta)
}

beta <- c(1.2, 1.5, 2.0)
mu   <- c(0.4, 0.5, 0.8)

out1 <- sol_SI(beta[1], mu[1], t, I0)
out2 <- sol_SI(beta[2], mu[2], t, I0)
out3 <- sol_SI(beta[3], mu[3], t, I0)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção de suscetíveis",
     main = expression(
       S(t) == frac((beta - mu) + I[0] * (mu * e^((beta - mu) * t) - beta),
                    (beta - mu) + beta * I[0] * (e^((beta - mu) * t) - 1))
     ))

lines(t, out2$S, lwd = 2, lty = 2)
lines(t, out3$S, lwd = 2, lty = 3)

abline(h = out1$Se, lwd = 1.5, lty = 4)
abline(h = out2$Se, lwd = 1.5, lty = 5)
abline(h = out3$Se, lwd = 1.5, lty = 6)

legend("topright",
       legend = c(expression(beta == 1.2 ~ "," ~ mu == 0.4),
                  expression(beta == 1.5 ~ "," ~ mu == 0.5),
                  expression(beta == 2.0 ~ "," ~ mu == 0.8),
                  expression(S[1]^"*"),
                  expression(S[2]^"*"),
                  expression(S[3]^"*")),
       lwd = c(2, 2, 2, 1.5, 1.5, 1.5),
       lty = c(1, 2, 3, 4, 5, 6),
       bty = "n",
       cex = 0.9)

SIS sem dinâmica vital

No modelo compartimental SIS sem dinâmica vital, a população é dividida em dois compartimentos: suscetíveis, \(S(t)\), e infectados, \(I(t)\). Um indivíduo suscetível pode tornar-se infectado após contato com um infectado, e um indivíduo infectado pode recuperar-se e retornar ao compartimento dos suscetíveis. Não há nascimentos, mortes nem migração, de modo que a população total permanece constante.

O modelo é dado por

\[ \begin{align} \frac{dS}{dt} &= -\beta S I + \gamma I \\ \frac{dI}{dt} &= \beta S I - \gamma I \end{align} \]

sendo que \(\beta > 0\) é a taxa de transmissão e \(\gamma > 0\) é a taxa de recuperação.

SIS sem dinâmica vital

Somando as equações, obtemos

\[ \begin{align} \frac{d}{dt}(S+I) &= \frac{dS}{dt}+\frac{dI}{dt} \\ &= -\beta S I + \gamma I + \beta S I - \gamma I \\ \frac{d}{dt}(S+I)&= 0 \end{align} \]

Logo,

\[ S(t)+I(t)=N \]

Se a população for normalizada, então

\[ S(t)+I(t)=1 \]

e, portanto,

\[ S(t)=1-I(t) \]

Substituindo em \(\dfrac{dI}{dt}\), segue que

\[ \begin{align} \frac{dI}{dt} &= \beta S I - \gamma I \\ &= \beta (1-I)I - \gamma I \\ \frac{dI}{dt}&= (\beta-\gamma)I - \beta I^2 \\ \frac{dI}{dt}&= I(\beta-\gamma - \beta I) \end{align} \]

A solução analítica para os infectados é uma função logística

\[ I(t)=\frac{\beta-\gamma I_0 e^{(\beta-\gamma)t}}{\beta-\gamma+\beta I_0\left(e^{(\beta-\gamma)t}-1\right)} \]

e

\[ S(t)=\frac{\beta-\gamma+I_0\left(\gamma e^{(\beta-\gamma)t}-\beta\right)}{\beta-\gamma+\beta I_0\left(e^{(\beta-\gamma)t}-1\right)} \]

Os pontos de equilíbrio dos infectados são obtidos impondo

\[ \frac{dI}{dt}=0 \]

Logo,

\[ I^*=0 \qquad \text{ou} \qquad I^*=1-\frac{\gamma}{\beta} \]

Como

\[ S^*=1-I^* \]

os pontos de equilíbrio dos suscetíveis são

\[ S^*=1 \qquad \text{ou} \qquad S^*=\frac{\gamma}{\beta} \]

Interpretação:

se

\[ \beta \leq \gamma \]

o único equilíbrio biologicamente viável é

\[ I^*=0 \qquad \text{e} \qquad S^*=1 \]

se

\[ \beta > \gamma \]

existe equilíbrio endêmico, dado por

\[ I^*=1-\frac{\gamma}{\beta} \qquad \text{e} \qquad S^*=\frac{\gamma}{\beta} \]

Portanto, no modelo SIS sem dinâmica vital, a infecção pode persistir na população porque os indivíduos recuperados retornam continuamente ao compartimento dos suscetíveis.

t  <- seq(0, 10, by = 0.01)
S0 <- 0.99
I0 <- 0.01

sol_SIS <- function(beta, gamma, t, S0, I0) {
  I <- ((beta - gamma) * I0 * exp((beta - gamma) * t)) /
    (beta * S0 - gamma + beta * exp((beta - gamma) * t) * I0)
  
  S <- 1 - I
  
  list(S = S, I = I)
}

# Caso 1: beta <= gamma
beta1  <- 0.4
gamma1 <- 0.8

out1 <- sol_SIS(beta1, gamma1, t, S0, I0)

S1e <- 1
I1e <- 0

# Caso 2: beta > gamma
beta2  <- 1.2
gamma2 <- 0.4

out2 <- sol_SIS(beta2, gamma2, t, S0, I0)

S2e <- gamma2 / beta2
I2e <- 1 - gamma2 / beta2

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 4, 1))

titulo1 <- substitute(
  beta <= gamma ~ "," ~ beta == b ~ "," ~ gamma == g,
  list(b = beta1, g = gamma1)
)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo1)

lines(t, out1$I, lwd = 2, lty = 2)
abline(h = S1e, lwd = 1.5, lty = 3)
abline(h = I1e, lwd = 1.5, lty = 4)

legend("right",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

titulo2 <- substitute(
  beta > gamma ~ "," ~ beta == b ~ "," ~ gamma == g,
  list(b = beta2, g = gamma2)
)

plot(t, out2$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo2)

lines(t, out2$I, lwd = 2, lty = 2)
abline(h = S2e, lwd = 1.5, lty = 3)
abline(h = I2e, lwd = 1.5, lty = 4)

legend("right",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

par(oldpar)

t  <- seq(0, 10, by = 0.01)
I0 <- 0.01
S0 <- 1 - I0

sol_SIS <- function(beta, gamma, t, I0) {
  expo <- exp((beta - gamma) * t)
  
  S <- ((beta - gamma) + I0 * (gamma * expo - beta)) /
    ((beta - gamma) + beta * I0 * (expo - 1))
  
  list(
    S  = S,
    Se = gamma / beta
  )
}

beta  <- c(1.2, 1.5, 2.0)
gamma <- c(0.4, 0.5, 0.8)

out1 <- sol_SIS(beta[1], gamma[1], t, I0)
out2 <- sol_SIS(beta[2], gamma[2], t, I0)
out3 <- sol_SIS(beta[3], gamma[3], t, I0)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção de suscetíveis",
     main = expression(
       S(t) == frac((beta - gamma) + I[0] * (gamma * e^((beta - gamma) * t) - beta),
                    (beta - gamma) + beta * I[0] * (e^((beta - gamma) * t) - 1))
     ))

lines(t, out2$S, lwd = 2, lty = 2)
lines(t, out3$S, lwd = 2, lty = 3)

abline(h = out1$Se, lwd = 1.5, lty = 4)
abline(h = out2$Se, lwd = 1.5, lty = 5)
abline(h = out3$Se, lwd = 1.5, lty = 6)

legend("topright",
       legend = c(expression(beta == 1.2 ~ "," ~ gamma == 0.4),
                  expression(beta == 1.5 ~ "," ~ gamma == 0.5),
                  expression(beta == 2.0 ~ "," ~ gamma == 0.8),
                  expression(S[1]^"*"),
                  expression(S[2]^"*"),
                  expression(S[3]^"*")),
       lwd = c(2, 2, 2, 1.5, 1.5, 1.5),
       lty = c(1, 2, 3, 4, 5, 6),
       bty = "n",
       cex = 0.9)

SIS com dinâmica vital

No modelo compartimental SIS com dinâmica vital, a população é dividida em dois compartimentos: suscetíveis, \(S(t)\), e infectados, \(I(t)\). Um indivíduo suscetível pode tornar-se infectado após contato com um infectado, e um indivíduo infectado pode recuperar-se e retornar ao compartimento dos suscetíveis. Além disso, o modelo inclui nascimentos e mortes naturais.

Assumindo taxa per capita de natalidade e mortalidade iguais a \(\mu > 0\), taxa de transmissão \(\beta > 0\) e taxa de recuperação \(\gamma > 0\), o modelo é dado por

\[ \begin{align} \frac{dS}{dt} &= \mu N - \beta S I + \gamma I - \mu S \\ \frac{dI}{dt} &= \beta S I - \gamma I - \mu I \end{align} \]

sendo que

\[ N(t)=S(t)+I(t) \]

SIS com dinâmica vital

Somando as duas equações, obtemos

\[ \begin{align} \frac{dN}{dt} &= \frac{dS}{dt}+\frac{dI}{dt} \\ &= \mu N - \beta S I + \gamma I - \mu S + \beta S I - \gamma I - \mu I \\ &= \mu N - \mu(S+I) \\ \frac{dN}{dt}&= 0 \end{align} \]

Logo, a população total permanece constante. Se a população for normalizada, então

\[ S(t)+I(t)=1 \]

e, portanto,

\[ S(t)=1-I(t) \]

Substituindo em \(\dfrac{dI}{dt}\), segue que

\[ \begin{align} \frac{dI}{dt} &= \beta S I - \gamma I - \mu I \\ &= \beta (1-I)I - (\gamma+\mu)I \\ \frac{dI}{dt}&= \left[\beta-(\gamma+\mu)\right]I - \beta I^2 \end{align} \]

A solução analítica para os infectados é uma função logística

\[ I(t)=\frac{(\beta-\gamma-\mu)I_0e^{(\beta-\gamma-\mu)t}}{\beta-\gamma-\mu+\beta I_0\left(e^{(\beta-\gamma-\mu)t}-1\right)} \] e

\[ S(t)=\frac{\beta-\gamma-\mu+I_0\left((\gamma+\mu)e^{(\beta-\gamma-\mu)t}-\beta\right)}{\beta-\gamma-\mu+\beta I_0\left(e^{(\beta-\gamma-\mu)t}-1\right)} \]

Os pontos de equilíbrio dos infectados são obtidos impondo

\[ \frac{dI}{dt}=0 \]

Logo,

\[ I^*=0 \qquad \text{ou} \qquad I^*=1-\frac{\gamma+\mu}{\beta} \]

Como

\[ S^*=1-I^* \]

os pontos de equilíbrio dos suscetíveis são

\[ S^*=1 \qquad \text{ou} \qquad S^*=\frac{\gamma+\mu}{\beta} \]

Interpretação:

se

\[ \beta \leq \gamma+\mu \]

o único equilíbrio biologicamente viável é

\[ I^*=0 \qquad \text{e} \qquad S^*=1 \]

se

\[ \beta > \gamma+\mu \]

existe equilíbrio endêmico, dado por

\[ I^*=1-\frac{\gamma+\mu}{\beta} \qquad \text{e} \qquad S^*=\frac{\gamma+\mu}{\beta} \]

Portanto, no modelo SIS com dinâmica vital, a persistência da infecção depende de a taxa de transmissão superar a soma da taxa de recuperação com a taxa de mortalidade natural.

t  <- seq(0, 10, by = 0.01)
I0 <- 0.01

sol_SIS <- function(beta, gamma, mu, t, I0) {
  expo <- exp((beta - gamma - mu) * t)
  
  den <- beta - gamma - mu + beta * I0 * (expo - 1)
  
  I <- ((beta - gamma - mu) * I0 * expo) / den
  
  S <- (beta - gamma - mu + I0 * ((gamma + mu) * expo - beta)) / den
  
  list(S = S, I = I)
}

# Caso 1: beta <= gamma + mu
beta1  <- 0.6
gamma1 <- 0.4
mu1    <- 0.3

out1 <- sol_SIS(beta1, gamma1, mu1, t, I0)

S1e <- 1
I1e <- 0

# Caso 2: beta > gamma + mu
beta2  <- 1.4
gamma2 <- 0.4
mu2    <- 0.3

out2 <- sol_SIS(beta2, gamma2, mu2, t, I0)

S2e <- (gamma2 + mu2) / beta2
I2e <- 1 - (gamma2 + mu2) / beta2

oldpar <- par(mfrow = c(1, 2), mar = c(4, 4, 4, 1))

titulo1 <- substitute(
  beta <= gamma + mu ~ "," ~ beta == b ~ "," ~ gamma == g ~ "," ~ mu == m,
  list(b = beta1, g = gamma1, m = mu1)
)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo1)

lines(t, out1$I, lwd = 2, lty = 2)
abline(h = S1e, lwd = 1.5, lty = 3)
abline(h = I1e, lwd = 1.5, lty = 4)

legend("right",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

titulo2 <- substitute(
  beta > gamma + mu ~ "," ~ beta == b ~ "," ~ gamma == g ~ "," ~ mu == m,
  list(b = beta2, g = gamma2, m = mu2)
)

plot(t, out2$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção",
     main = titulo2)

lines(t, out2$I, lwd = 2, lty = 2)
abline(h = S2e, lwd = 1.5, lty = 3)
abline(h = I2e, lwd = 1.5, lty = 4)

legend("topright",
       legend = c("S(t)", "I(t)", expression(S^"*"), expression(I^"*")),
       lwd = c(2, 2, 1.5, 1.5),
       lty = c(1, 2, 3, 4),
       bty = "n",
       cex = 0.9)

par(oldpar)

t  <- seq(0, 10, by = 0.01)
I0 <- 0.01

sol_SIS <- function(beta, gamma, mu, t, I0) {
  expo <- exp((beta - gamma - mu) * t)
  
  S <- (beta - gamma - mu + I0 * ((gamma + mu) * expo - beta)) /
    (beta - gamma - mu + beta * I0 * (expo - 1))
  
  list(
    S  = S,
    Se = (gamma + mu) / beta
  )
}

beta  <- c(1.4, 1.8, 2.2)
gamma <- c(0.4, 0.5, 0.6)
mu    <- c(0.2, 0.4, 0.7)

out1 <- sol_SIS(beta[1], gamma[1], mu[1], t, I0)
out2 <- sol_SIS(beta[2], gamma[2], mu[2], t, I0)
out3 <- sol_SIS(beta[3], gamma[3], mu[3], t, I0)

plot(t, out1$S, type = "l", lwd = 2,
     ylim = c(0, 1),
     xlab = "Tempo",
     ylab = "Proporção de suscetíveis",
     main = expression(
       S(t) == frac(beta - gamma - mu + I[0] * ((gamma + mu) * e^((beta - gamma - mu) * t) - beta),
                    beta - gamma - mu + beta * I[0] * (e^((beta - gamma - mu) * t) - 1))
     ))

lines(t, out2$S, lwd = 2, lty = 2)
lines(t, out3$S, lwd = 2, lty = 3)

abline(h = out1$Se, lwd = 1.5, lty = 4)
abline(h = out2$Se, lwd = 1.5, lty = 5)
abline(h = out3$Se, lwd = 1.5, lty = 6)

legend("bottomright",
       legend = c(expression(beta == 1.4 ~ "," ~ gamma == 0.4 ~ "," ~ mu == 0.2),
                  expression(beta == 1.8 ~ "," ~ gamma == 0.5 ~ "," ~ mu == 0.4),
                  expression(beta == 2.2 ~ "," ~ gamma == 0.6 ~ "," ~ mu == 0.7),
                  expression(S[1]^"*"),
                  expression(S[2]^"*"),
                  expression(S[3]^"*")),
       lwd = c(2, 2, 2, 1.5, 1.5, 1.5),
       lty = c(1, 2, 3, 4, 5, 6),
       bty = "n",
       cex = 0.9)

Os equilíbrios são

\[ S_1^*=\frac{0.4+0.2}{1.4}\approx 0.429 \qquad S_2^*=\frac{0.5+0.4}{1.8}=0.5 \qquad S_3^*=\frac{0.6+0.7}{2.2}\approx 0.591 \]

Modelos compartimentais farmacocinéticos

A farmacocinética é a área das ciências farmacêuticas que se ocupa do estudo do que acontece ao medicamento dentro do organismo, procurando responder a questões do tipo:

• Como varia com o tempo a quantidade de fármaco no organismo?
• Como varia com o tempo a quantidade de fármaco no sangue?
• Como se distribui o fármaco pelo organismo?

Para encontrar respostas a este tipo de questões, aplica-se o princípio de conservação da massa aos fármacos no organismo. Para simplificar o problema, este é descrito através de modelos farmacocinéticos, nomeadamente modelos compartimentais e fisiológicos.

Além da farmacocinética, a análise baseada em modelos farmacocinéticos também é utilizada na toxicologia. A motivação é semelhante: determinar como varia a concentração do agente tóxico e quais órgãos são atingidos, permitindo compreender o mecanismo de toxicidade.

Modelos compartimentais

Num modelo compartimental, o organismo é representado de forma simplificada como uma combinação de compartimentos que comunicam entre si.

Um compartimento não representa necessariamente um órgão ou região anatômica, mas sim um grupo de tecidos com afinidade semelhante para o fármaco ou irrigação sanguínea semelhante.

Assume-se ainda que o fármaco se encontra uniformemente distribuído dentro de cada compartimento.

Para determinar a quantidade de medicamento em cada compartimento, aplica-se o princípio de conservação da massa.

Para o medicamento \(A\) e para o compartimento \(i\):

A velocidade de desaparecimento do fármaco do compartimento, nos modelos farmacocinéticos, engloba os termos de consumo e saída da equação de conservação da massa e é normalmente proporcional à quantidade de fármaco existente no compartimento.

O símbolo

\[ \dot m \]

representa a derivada temporal de \(m\), isto é,

\[ \dot m = \frac{dm}{dt} \]

Em farmacocinética, \(\dot m\) representa a taxa de variação da massa de fármaco por unidade de tempo.

Por outro lado, a forma de administração do medicamento condiciona o termo de entrada no compartimento sanguíneo:

• Administração intravenosa:

\[ \dot m_{A,\text{entrada sangue}}=0 \]

e

\[ m_{A,\text{sangue}}(0)=D \] Sendo que \(D\) representa a dose do fármaco administrada.

• Perfusão contínua:

\[ \dot m_{A,\text{entrada sangue}}=R \]

e

\[ m_{A,\text{sangue}}(0)=0 \]

• Administração extravascular:

Assume-se geralmente que a velocidade de absorção é proporcional à quantidade disponível no local de absorção:

\[ \dot m_{A,\text{entrada sangue}} = k_a m_{A,\text{LAb}} \]

e

\[ m_{A,\text{LAb}}(0)=F D \]

onde

\[ \begin{align} k_a & = \text{constante de absorção} \\ m_{A,\text{LAb}} & = \text{massa de fármaco no local de absorção} \\ F & = \text{fração da dose absorvida} \end{align} \]

Modelo monocompartimental

O modelo compartimental mais simples considera todo o organismo como um único compartimento, Figura 11.5

Neste modelo, assume-se que o fármaco se distribui uniforme e instantaneamente por todos os tecidos do organismo e é eliminado a uma velocidade proporcional à quantidade existente em cada instante

No caso de uma administração intravenosa em dose única,

\[ D(t>0)=0 \]

e

\[ m_{\text{organismo}}(0)=D \]

Aplicando o princípio de conservação da massa, obtém-se a equação diferencial

\[ \frac{dm}{dt} = -k_{el}m \]

Separando variáveis:

\[ \frac{dm}{m} = -k_{el}\,dt \]

Integrando:

\[ \int \frac{dm}{m} = \int -k_{el}\,dt \]

Logo,

\[ \ln(m) = -k_{el}t + C \]

Como

\[ m(0)=D \]

segue que

\[ C=\ln(D) \]

Portanto,

\[ \ln\left(\frac{m}{D}\right) = -k_{el}t \]

e assim

\[ m(t) = De^{-k_{el}t} \]

Se \(v\) representa o volume de distribuição, então a concentração plasmática é

\[ c(t) = \frac{m(t)}{v} = \frac{D}{v}e^{-k_{el}t} \]

Quando a distribuição do fármaco não é instantânea, é conveniente considerar modelos com dois ou mais compartimentos, levando a sistemas de equações diferenciais.

Exemplo: Modelo Compartimental da Lidocaína

• Barreira, 2014, p. 414-7, Exemplo 11.13

A lidocaína (2-(dietilamino)-N-(2,6-dimetilfenil)-acetamida) é um fármaco do grupo dos antiarrítmicos da classe I e dos anestésicos locais utilizado no tratamento da arritmia cardíaca e da dor localizada.

A sua farmacocinética pode ser descrita por um modelo compartimental com dois compartimentos:

• compartimento 1: sangue
• compartimento 2: coração

As constantes do modelo são

\[ k_{12}=0.0923\ \text{min}^{-1} \qquad k_{21}=0.0441\ \text{min}^{-1} \qquad k_{el}=0.0294\ \text{min}^{-1} \]

Considere a administração intravenosa de \(D=5\ \text{mg}\) na corrente sanguínea.

Determine a massa de fármaco no coração uma hora após a administração de 5 mg de lidocaína na corrente sanguínea.

Solução:

Aplicando o princípio de conservação da massa ao compartimento “coração”, obtém-se

\[ k_{12}m_1 = k_{21}m_2 + \frac{dm_2}{dt} \]

Esta equação não pode ser resolvida isoladamente porque nela aparece a massa do compartimento “sangue”. Assim, é necessário escrever também a equação de conservação para o compartimento “sangue”:

\[ k_{21}m_2 = k_{12}m_1 + k_{el}m_1 + \frac{dm_1}{dt} \]

Sabemos ainda que

\[ m_1(0)=5 \]

e

\[ m_2(0)=0 \]

Logo, o problema conduz ao seguinte sistema de equações diferenciais:

\[ \begin{cases} \dfrac{dm_2}{dt} = 0.0923m_1-0.0441m_2 \\ \dfrac{dm_1}{dt} = -0.1217m_1+0.0441m_2 \\ m_1(0)=5 \\ m_2(0)=0 \end{cases} \]

Em notação matricial:

\[ \begin{bmatrix} \dfrac{dm_1}{dt} \\ \dfrac{dm_2}{dt} \end{bmatrix} = \begin{bmatrix} -0.1217 & 0.0441 \\ 0.0923 & -0.0441 \end{bmatrix} \begin{bmatrix} m_1 \\ m_2 \end{bmatrix} \] A solução geral do sistema é

\[ \begin{bmatrix} m_1 \\ m_2 \end{bmatrix} = 3.093 \begin{bmatrix} 0.3885 \\ 1 \end{bmatrix} e^{-0.00823t} - 3.093 \begin{bmatrix} -1.228 \\ 1 \end{bmatrix} e^{-0.1576t} \]

As equações que permitem calcular a massa em cada compartimento em função do tempo são

\[ m_1(t) = 1.202e^{-0.00823t} + 3.798e^{-0.1576t} \]

e

\[ m_2(t) = 3.093 \left( e^{-0.00823t} - e^{-0.1576t} \right) \]

Assim, ao fim de uma hora,

\[ m_2(60) = 3.093 \left( e^{-0.00823\cdot 60} - e^{-0.1576\cdot 60} \right) \approx 1.887\text{ mg} \]

m1 <- function(t) {
  
  1.202 * exp(-0.00823 * t) +
    3.798 * exp(-0.1576 * t)
  
}

m2 <- function(t) {
  
  3.093 * (
    exp(-0.00823 * t) -
      exp(-0.1576 * t)
  )
  
}

t <- seq(0, 120, length.out = 2000)

plot(
  t,
  m1(t),
  type = "l",
  lwd  = 3,
  lty  = 1,
  xlab = "t (min)",
  ylab = "massa (mg)",
  main = "Modelo compartimental da lidocaína",
  ylim = c(
    0,
    max(
      m1(t),
      m2(t)
    )
  )
)

lines(
  t,
  m2(t),
  lwd = 3,
  lty = 2
)

legend(
  "topright",
  legend = c(
    expression(m[1](t) ~ "(sangue)"),
    expression(m[2](t) ~ "(coração)")
  ),
  lty = c(1, 2),
  lwd = 3,
  bty = "n"
)

abline(
  v = 60,
  lty = 3
)

points(
  60,
  m2(60),
  pch = 19
)

text(
  60,
  m2(60),
  labels = expression(m[2](60)),
  pos = 3
)

Modelos de base fisiológica

A principal desvantagem dos modelos compartimentais reside no facto de os compartimentos não corresponderem necessariamente a entidades anatômicas ou fisiológicas reais.

Além disso, estes modelos não permitem calcular diretamente a concentração do princípio ativo ou agente tóxico em todos os órgãos ou tecidos, pois para isso seria necessário conhecer características anatomofisiológicas específicas, tais como:

• dimensão do órgão
• débito sanguíneo
• permeabilidade das membranas
• afinidade do princípio ativo pelo tecido

Reunindo estas características fisiológicas e físico-químicas do princípio ativo, é possível construir um modelo farmacocinético fisiológico que inclua tantos compartimentos quanto os órgãos ou tecidos pelos quais o princípio ativo se distribui.

Exemplo: Modelo PBPK de Exposição por Inalação

• Barreira, 2014, p. 425-30, Exemplo 11.15

Determinar como variam as concentrações no sangue, fígado e resto do corpo quando o organismo passa a inalar ar com concentração de \(0.34 \ \text{mg/L}\) de uma dada substância.

Considere que o seguinte modelo \(PBPK\) descreve bem o comportamento desta substância.

Utilize os parâmetros fisiológicos da tabela abaixo.

Parâmetro fisiológico	Abreviatura	Valor
Débito cardíaco	\(Q_C\)	5.64 L/hr
Taxa de ventilação alveolar	\(Q_P\)	4.5 L/h
Fluxo sanguíneo hepático	\(Q_F\)	2.11 L/h
Fluxo sanguíneo para o resto do corpo	\(Q_{RC}\)	3.53 L/h
Volume do fígado	\(V_F\)	0.012 L
Volume resto do corpo	\(V_{RC}\)	0.261 L
Coeficiente de partição fígado-sangue	\(P_F\)	2.7
Coeficiente de partição resto do corpo-sangue	\(P_{RC}\)	6.34
Coeficiente de partição ar-sangue	\(P_S\)	40
Velocidade máxima do metabolismo	\(V_{\max}\)	3.6 mg/h
Constante de Michaelis-Menten	\(K_M\)	0.36 mg/L

Solução:

Equações do modelo PBPK:

Como a substância é inalada, as concentrações no sangue arterial \((C_A)\) e no sangue venoso \((C_V)\) são dadas por

\[ C_A = \frac{ Q_P C_{in}+Q_C C_V }{ Q_C+\dfrac{Q_P}{P_S} } \]

e

\[ C_V = \frac{ Q_F C_{VF}+Q_{RC} C_{VRC} }{ Q_C } \]

Para os restantes tecidos, a aplicação do princípio de conservação da massa conduz às seguintes equações

Fígado:

\[ \frac{dM_F}{dt} = Q_F C_A - Q_F C_{VF} - \frac{ V_{\max}C_{VF} }{ K_M+C_{VF} } \]

\[ C_F = \frac{M_F}{V_F} \]

\[ C_{VF} = \frac{C_F}{P_F} \]

Resto do corpo:

\[ \frac{dM_{RC}}{dt} = Q_{RC}C_A - Q_{RC}C_{VRC} \]

\[ C_{RC} = \frac{M_{RC}}{V_{RC}} \]

\[ C_{VRC} = \frac{C_{RC}}{P_{RC}} \]

Condições iniciais:

\[ C_V=C_A=C_F=C_{RC}=0 \]

library(deSolve)

pbpk <- function(t, y, p) {
  
  MF  <- y["MF"]
  MRC <- y["MRC"]
  
  Cin <- ifelse(t <= 1, p["Cin"], 0)
  
  CF  <- MF / p["VF"]
  CRC <- MRC / p["VRC"]
  
  CVF  <- CF / p["PF"]
  CVRC <- CRC / p["PRC"]
  
  CV <- (p["QF"] * CVF + p["QRC"] * CVRC) / p["QC"]
  
  CA <- (p["QP"] * Cin + p["QC"] * CV) /
    (p["QC"] + p["QP"] / p["PS"])
  
  dMF <- p["QF"] * CA -
    p["QF"] * CVF -
    (p["Vmax"] * CVF) / (p["KM"] + CVF)
  
  dMRC <- p["QRC"] * CA -
    p["QRC"] * CVRC
  
  list(c(dMF, dMRC))
}

param <- c(
  QC   = 5.64,
  QP   = 4.5,
  QF   = 2.11,
  QRC  = 3.53,
  VF   = 0.012,
  VRC  = 0.261,
  PF   = 2.7,
  PRC  = 6.34,
  PS   = 40,
  Vmax = 3.6,
  KM   = 0.36,
  Cin  = 0.34
)

y0 <- c(MF = 0, MRC = 0)

tempo <- seq(0, 10, by = 0.01)

sol <- ode(
  y = y0,
  times = tempo,
  func = pbpk,
  parms = param
)

sol <- as.data.frame(sol)
names(sol) <- c("time", "MF", "MRC")

sol$CF <- sol$MF / param["VF"]
sol$CRC <- sol$MRC / param["VRC"]

sol$CVF <- sol$CF / param["PF"]
sol$CVRC <- sol$CRC / param["PRC"]

sol$CV <- (param["QF"] * sol$CVF + param["QRC"] * sol$CVRC) /
  param["QC"]

sol$Cin <- ifelse(sol$time <= 1, param["Cin"], 0)

sol$CA <- (param["QP"] * sol$Cin + param["QC"] * sol$CV) /
  (param["QC"] + param["QP"] / param["PS"])

plot(
  sol$time,
  sol$CV,
  type = "l",
  lwd = 3,
  lty = 1,
  xlab = "t (h)",
  ylab = "C (mg/L)",
  main = "Modelo PBPK por inalação",
  ylim = c(0, 3)
)

lines(sol$time, sol$CA,  lwd = 3, lty = 2)
lines(sol$time, sol$CF,  lwd = 3, lty = 3)
lines(sol$time, sol$CRC, lwd = 3, lty = 4)

abline(v = 1, lty = 3)

legend(
  "topright",
  legend = c(
    expression(C[V] ~ "sangue venoso"),
    expression(C[A] ~ "sangue arterial"),
    expression(C[F] ~ "fígado"),
    expression(C[RC] ~ "resto do corpo")
  ),
  lty = 1:4,
  lwd = 3,
  bty = "n"
)

Sistema de equações diferenciais não-lineares: modelo presa-predador de Lotka-Volterra

Equação de Lotka-Volterra: Wikipedia

O sistema

\[ \begin{align} x^{\prime} &= ax-bxy\\ y^{\prime} &= cxy-dy \end{align} \]

com \(a>0\), \(b>0\), \(c>0\) e \(d>0\), é conhecido como modelo presa-predador de Lotka-Volterra.

Neste modelo, \(x(t)\) representa a população de presas no instante \(t\), enquanto \(y(t)\) representa a população de predadores. O termo \(xy\) representa a interação entre as duas populações. Essa interação é proporcional ao produto das populações porque, quanto maior o número de presas e de predadores, maior tende a ser o número de encontros entre eles.

A primeira equação pode ser escrita como

\[ x^{\prime}=x(a-by) \]

O termo \(ax\) representa o crescimento natural das presas na ausência de predadores. O termo \(-bxy\) representa a redução da população de presas pela ação dos predadores. Portanto, se houver poucos predadores, a população de presas tende a crescer. Se houver muitos predadores, a população de presas tende a diminuir.

A segunda equação pode ser escrita como

\[ y^{\prime}=y(cx-d) \]

O termo \(cxy\) representa o crescimento da população de predadores devido à disponibilidade de presas. O termo \(-dy\) representa a mortalidade natural dos predadores na ausência de alimento. Portanto, se houver poucas presas, a população de predadores tende a diminuir. Se houver muitas presas, a população de predadores tende a aumentar.

O equilíbrio biologicamente relevante ocorre quando as duas populações permanecem constantes. Isso exige

\[ x^{\prime}=0 \qquad \text{e} \qquad y^{\prime}=0 \]

O equilíbrio positivo é

\[ (x^*,y^*)= \left(\frac{d}{c},\frac{a}{b}\right) \]

No exemplo da figura,

\[ a=1,\qquad b=0.002,\qquad c=0.00001,\qquad d=0.08 \]

Assim,

\[ x^*=\frac{d}{c} = \frac{0.08}{0.00001} = 8000 \]

e

\[ y^*=\frac{a}{b} = \frac{1}{0.002} = 500 \]

Logo, o ponto de equilíbrio positivo é

\[ (x^*,y^*)=(8000,500) \]

A equação

\[ a\ln(y)-by=cx-d\ln(x)+k \]

descreve as órbitas do sistema no plano de fase \((x,y)\). Cada valor de \(k\) define uma curva diferente. A figura mostra uma dessas curvas fechadas.

Fig. 11.11. Gráfico logarítmico duplo de pontos (x, y) que satisfazem a equação da solução com as constantes particulares a = 1, b = 0.002, c = 0.00001, d = 0.08, k = 5.8. O conjunto solução é uma curva fechada. Se considerarmos x e y como funções do tempo, o ponto (x, y) move-se ao longo da curva fechada no sentido anti-horário, descrevendo assim ciclo por ciclo.

Plot[ln(y) - 0.002 y = 0.00001 x - 0.08 ln(x) + 5.8, {x,1000, 22000}, {y,300,700}]

Esse gráfico não representa \(x(t)\) ou \(y(t)\) separadamente em função do tempo. Ele representa a trajetória conjunta das duas populações. Cada ponto da curva corresponde a um par \((x,y)\), isto é, a uma quantidade de presas e a uma quantidade de predadores.

A interpretação biológica do ciclo é simples. Quando há poucos predadores, a população de presas cresce. Com mais presas disponíveis, a população de predadores passa a crescer. Com muitos predadores, a população de presas diminui. Com poucas presas, a população de predadores também diminui. Com menos predadores, as presas voltam a crescer, reiniciando o ciclo.

Portanto, o pico da população de presas ocorre antes do pico da população de predadores. Isso é biologicamente esperado, pois primeiro aumenta a disponibilidade de alimento e, somente depois, aumenta a população que depende desse alimento.

No exemplo, o equilíbrio ocorre em

\[ x=8000 \qquad \text{e} \qquad y=500. \]

Assim, quando \(y<500\), a população de presas aumenta. Quando \(y>500\), a população de presas diminui. Quando \(x<8000\), a população de predadores diminui. Quando \(x>8000\), a população de predadores aumenta.

O modelo é útil porque mostra que uma interação biológica simples pode gerar oscilações populacionais. A oscilação não precisa ser imposta externamente; ela surge da própria relação entre presa e predador.

Apesar disso, o modelo é idealizado. Ele não considera limitação ambiental das presas, competição entre indivíduos da mesma espécie, saturação alimentar dos predadores, migração, sazonalidade ou variações aleatórias do ambiente. Portanto, deve ser entendido como um modelo conceitual, não como uma descrição completa de um ecossistema real.

A mensagem central do modelo é que muitas presas favorecem o crescimento dos predadores; muitos predadores reduzem as presas; poucas presas reduzem os predadores; poucos predadores permitem novo crescimento das presas. Essa retroalimentação produz ciclos populacionais no plano de fase.

lotka_volterra <- function(t, estado, parametros) {
  x <- estado["x"]
  y <- estado["y"]
  
  alpha <- parametros["alpha"]
  beta  <- parametros["beta"]
  gamma <- parametros["gamma"]
  delta <- parametros["delta"]
  
  dx <- alpha * x - beta * x * y
  dy <- gamma * x * y - delta * y
  
  list(c(dx, dy))
}

parametros <- c(
  alpha = 1,
  beta  = 1,
  gamma = 1,
  delta = 1
)

estado_inicial <- c(
  x = 1.5,
  y = 1.0
)

tempo <- seq(0, 20, by = 0.01)

solucao <- deSolve::ode(
  y = estado_inicial,
  times = tempo,
  func = lotka_volterra,
  parms = parametros
)

solucao <- as.data.frame(solucao)

x_eq <- parametros["delta"] / parametros["gamma"]
y_eq <- parametros["alpha"] / parametros["beta"]

plot(
  solucao$x,
  solucao$y,
  type = "l",
  lwd = 2,
  xlab = "Presas",
  ylab = "Predadores",
  main = "Modelo presa-predador: gráfico de fase"
)

points(
  estado_inicial["x"],
  estado_inicial["y"],
  pch = 19
)

text(
  estado_inicial["x"],
  estado_inicial["y"],
  labels = "Início",
  pos = 4
)

points(
  x_eq,
  y_eq,
  pch = 19
)

text(
  x_eq,
  y_eq,
  labels = "Equilíbrio",
  pos = 4
)

indices <- seq(1, nrow(solucao) - 20, by = 150)

arrows(
  solucao$x[indices],
  solucao$y[indices],
  solucao$x[indices + 20],
  solucao$y[indices + 20],
  length = 0.08
)

mtext(
  "Parâmetros: alpha = 1, beta = 1, gamma = 1, delta = 1; \ncondição inicial: x(0) = 1.5, y(0) = 1.0",
  side = 3,
  line = 0.3,
  cex = 0.8
)

lotka_volterra <- function(t, estado, parametros) {
  x <- estado["x"]
  y <- estado["y"]
  
  a <- parametros["a"]
  b <- parametros["b"]
  c <- parametros["c"]
  d <- parametros["d"]
  
  dx <- a * x - b * x * y
  dy <- c * x * y - d * y
  
  list(c(dx, dy))
}

parametros <- c(
  a = 1,
  b = 0.002,
  c = 0.00001,
  d = 0.08
)

valores_k <- c(5.84, 5.8, 5.7)

x0 <- 5000

obter_y0 <- function(k, x0, parametros) {
  a <- parametros["a"]
  b <- parametros["b"]
  c <- parametros["c"]
  d <- parametros["d"]
  
  f_y <- function(y) {
    a * log(y) - b * y - c * x0 + d * log(x0) - k
  }
  
  y_critico <- a / b
  
  uniroot(
    f_y,
    interval = c(y_critico, 5000)
  )$root
}

tempo <- seq(0, 120, by = 0.01)

solucoes <- list()
y0_valores <- numeric(length(valores_k))

for (i in seq_along(valores_k)) {
  k <- valores_k[i]
  
  y0 <- obter_y0(
    k = k,
    x0 = x0,
    parametros = parametros
  )
  
  y0_valores[i] <- y0
  
  estado_inicial <- c(
    x = x0,
    y = y0
  )
  
  solucao <- deSolve::ode(
    y = estado_inicial,
    times = tempo,
    func = lotka_volterra,
    parms = parametros
  )
  
  solucoes[[i]] <- as.data.frame(solucao)
}

x_eq <- parametros["d"] / parametros["c"]
y_eq <- parametros["a"] / parametros["b"]

x_lim <- range(
  unlist(
    lapply(
      solucoes,
      function(z) z$x
    )
  )
)

y_lim <- range(
  unlist(
    lapply(
      solucoes,
      function(z) z$y
    )
  )
)

plot(
  NA,
  NA,
  xlim = x_lim,
  ylim = y_lim,
  xlab = "Presas",
  ylab = "Predadores",
  main = "Modelo presa-predador: gráfico de fase"
)

for (i in seq_along(solucoes)) {
  solucao <- solucoes[[i]]
  
  lines(
    solucao$x,
    solucao$y,
    lwd = 2,
    lty = i
  )
  
  indices <- seq(
    1,
    nrow(solucao) - 80,
    by = 2500
  )
  
  arrows(
    solucao$x[indices],
    solucao$y[indices],
    solucao$x[indices + 80],
    solucao$y[indices + 80],
    length = 0.08,
    lty = i
  )
  
  points(
    x0,
    y0_valores[i],
    pch = 19
  )
}

points(
  x_eq,
  y_eq,
  pch = 19
)

text(
  x_eq,
  y_eq,
  labels = "Equilíbrio",
  pos = 4
)

legend(
  "topright",
  legend = paste(
    "k =",
    valores_k,
    "| x(0) =",
    x0,
    "| y(0) =",
    round(y0_valores, 2)
  ),
  lwd = 2,
  lty = seq_along(valores_k),
  bty = "n",
  cex = 0.8
)

mtext(
  "Parâmetros: a = 1, b = 0.002, c = 0.00001, d = 0.08",
  side = 3,
  line = 0.3,
  cex = 0.8
)

Populações em função do tempo no modelo presa-predador

O sistema presa-predador de Lotka-Volterra é dado por

\[ \begin{align} x^{\prime} &= ax-bxy\\ y^{\prime} &= cxy-dy \end{align} \]

em que \(x(t)\) representa a população de presas no tempo \(t\) e \(y(t)\) representa a população de predadores no tempo \(t\).

O gráfico no plano de fase mostra a trajetória no plano \((x,y)\). Nesse caso, o eixo horizontal é a população de presas e o eixo vertical é a população de predadores. Cada ponto do gráfico representa simultaneamente uma quantidade de presas e uma quantidade de predadores.

Já o gráfico em função do tempo é diferente. Nele, o eixo horizontal representa o tempo \(t\), enquanto o eixo vertical representa as populações \(x(t)\) e \(y(t)\). Portanto, em vez de observar a trajetória no plano \((x,y)\), observamos como cada população varia ao longo do tempo.

A primeira equação pode ser escrita como

\[ x^{\prime}=x(a-by) \]

Assim, a população de presas aumenta quando

\[ a-by>0 \]

isto é,

\[ y<\frac{a}{b} \]

A população de presas diminui quando

\[ y>\frac{a}{b} \]

Portanto, o crescimento ou a redução das presas depende diretamente da quantidade de predadores.

A segunda equação pode ser escrita como

\[ y^{\prime}=y(cx-d) \]

Assim, a população de predadores aumenta quando

\[ cx-d>0 \]

isto é,

\[ x>\frac{d}{c} \]

A população de predadores diminui quando

\[ x<\frac{d}{c} \]

Portanto, o crescimento ou a redução dos predadores depende diretamente da quantidade de presas.

O ponto de equilíbrio positivo é

\[ (x^*,y^*)= \left( \frac{d}{c}, \frac{a}{b} \right) \]

Nesse ponto, temos

\[ x^{\prime}=0 \qquad \text{e} \qquad y^{\prime}=0 \]

Isso significa que, se o sistema estiver exatamente nesse ponto, as duas populações permanecem constantes.

Para obter as funções \(x(t)\) e \(y(t)\), não basta usar a curva fechada do plano de fase. A equação da órbita,

\[ a\ln(y)-by=cx-d\ln(x)+k \]

descreve apenas a relação entre \(x\) e \(y\) ao longo da trajetória. Ela informa por onde o sistema passa, mas não informa em que instante \(t\) o sistema passa por cada ponto.

Para obter as populações como funções do tempo, é necessário resolver o problema de valor inicial

\[ \begin{align} x^{\prime} &= ax-bxy\\ y^{\prime} &= cxy-dy\\ x(0) &= x_0\\ y(0) &= y_0 \end{align} \]

A condição inicial \((x_0,y_0)\) define o ponto de partida do sistema. A partir desse ponto, o sistema determina a evolução temporal das duas populações.

Em geral, as funções \(x(t)\) e \(y(t)\) não são obtidas por fórmulas elementares simples. Por isso, usa-se solução numérica. O método numérico calcula aproximações sucessivas

\[ x(t_0),y(t_0),\quad x(t_1),y(t_1),\quad x(t_2),y(t_2),\ldots \]

para uma sequência de tempos

\[ t_0<t_1<t_2<\cdots<t_n \]

Assim, em vez de uma expressão explícita para \(x(t)\) e \(y(t)\), obtemos uma tabela aproximada com os valores das populações ao longo do tempo.

No R, essa solução pode ser obtida pela função ode() do pacote deSolve. Essa função recebe o sistema de equações diferenciais, os parâmetros, a condição inicial e a sequência de tempos. O resultado é uma tabela com as colunas time, x e y.

No exemplo da figura, os parâmetros são

\[ a=1,\qquad b=0.002,\qquad c=0.00001,\qquad d=0.08 \]

O ponto de equilíbrio positivo é

\[ x^*=\frac{d}{c} = \frac{0.08}{0.00001} = 8000 \]

e

\[ y^*=\frac{a}{b} = \frac{1}{0.002} = 500 \]

Logo,

\[ (x^*,y^*)=(8000,500) \]

O valor de \(k\) define uma órbita fechada no plano de fase. Para obter o gráfico em função do tempo associado a essa órbita, escolhemos um valor inicial \(x_0\) e calculamos \(y_0\) de modo que o ponto inicial satisfaça

\[ a\ln(y_0)-by_0=cx_0-d\ln(x_0)+k \]

Por exemplo, usando

\[ k=5.8 \]

e escolhendo

\[ x_0=5000 \]

podemos calcular \(y_0\) numericamente resolvendo

\[ a\ln(y_0)-by_0-cx_0+d\ln(x_0)-k=0 \]

Esse cálculo fornece um ponto inicial pertencente à curva fechada desejada. Depois disso, o sistema diferencial é resolvido numericamente para produzir as sequências \(x(t)\) e \(y(t)\).

O gráfico de \(x(t)\) e \(y(t)\) contra o tempo mostra oscilações periódicas. Primeiro, a população de presas aumenta. Com mais presas disponíveis, há mais alimento para os predadores, e a população de predadores começa a aumentar. Depois, com muitos predadores, a população de presas diminui. Com poucas presas, falta alimento para os predadores, e a população de predadores também diminui. Com menos predadores, as presas voltam a crescer. O ciclo então se repete.

Essa dinâmica produz uma defasagem entre as duas populações. O pico das presas ocorre antes do pico dos predadores. Isso é biologicamente esperado, pois o aumento dos predadores depende do aumento prévio da disponibilidade de presas.

No gráfico em função do tempo, essa defasagem aparece como duas curvas oscilatórias deslocadas uma em relação à outra. A curva das presas sobe primeiro; a curva dos predadores sobe depois. Em seguida, a curva das presas cai primeiro; a curva dos predadores cai depois.

Portanto, há duas representações complementares. O gráfico de fase mostra a órbita fechada no plano \((x,y)\). O gráfico temporal mostra como \(x(t)\) e \(y(t)\) variam ao longo do tempo.

A interpretação biológica é simples: muitas presas favorecem o crescimento dos predadores; muitos predadores reduzem as presas; poucas presas reduzem os predadores; poucos predadores permitem novo crescimento das presas. Essa retroalimentação gera oscilações populacionais.

É importante notar que essas oscilações são consequência do próprio modelo. Elas não são impostas externamente. Surgem da interação não-linear entre as populações, representada pelo produto \(xy\).

Apesar disso, o modelo é idealizado. Ele não inclui capacidade de suporte ambiental, competição intraespecífica, saturação alimentar, migração, sazonalidade ou variações aleatórias. Por isso, o gráfico periódico deve ser interpretado como comportamento teórico de um modelo simples, não como descrição completa de um ecossistema real.

lotka_volterra <- function(t, estado, parametros) {
  x <- estado["x"]
  y <- estado["y"]
  
  a <- parametros["a"]
  b <- parametros["b"]
  c <- parametros["c"]
  d <- parametros["d"]
  
  dx <- a * x - b * x * y
  dy <- c * x * y - d * y
  
  list(c(dx, dy))
}

parametros <- c(
  a = 1,
  b = 0.002,
  c = 0.00001,
  d = 0.08
)

k <- 5.8

x0 <- 5000

f_y <- function(y) {
  parametros["a"] * log(y) -
    parametros["b"] * y -
    parametros["c"] * x0 +
    parametros["d"] * log(x0) -
    k
}

y0 <- uniroot(f_y, interval = c(500, 1000))$root

estado_inicial <- c(
  x = x0,
  y = y0
)

tempo <- seq(0, 100, by = 0.01)

solucao <- deSolve::ode(
  y = estado_inicial,
  times = tempo,
  func = lotka_volterra,
  parms = parametros
)

solucao <- as.data.frame(solucao)

plot(
  solucao$time,
  solucao$x,
  type = "l",
  lwd = 2,
  xlab = "Tempo",
  ylab = "População",
  ylim = range(solucao$x, solucao$y),
  main = "Modelo presa-predador: populações em função do tempo"
)

lines(
  solucao$time,
  solucao$y,
  lwd = 2,
  lty = 2
)

legend(
  "topright",
  legend = c("Presas: x(t)", "Predadores: y(t)"),
  lwd = 2,
  lty = c(1, 2),
  bty = "n"
)

mtext(
  paste(
    "Parâmetros: a = 1, b = 0.002, c = 0.00001, d = 0.08, k = 5.8;",
    "x(0) = 5000, y(0) =",
    round(y0, 4)
  ),
  side = 3,
  line = 0.3,
  cex = 0.75
)

Manual do WolframAlpha

• Ordinary Differential Equations (ODEs)

• Symbolic Differential Equation Solving

• An elementary introduction to the Wolfram language

Equação diferencial ordinária no YouTube

• Brilliant: This is why you’re learning differential equations
• Physics with Elliot: Physics Students Need to Know These 5 Methods for Differential Equations
• Dr. Trefor Bazett: The MATH of Pandemics | Intro to the SIR Model
• Brian Sullivan: SIR Model For Disease Spread- 1. Introduction
• Jeffrey Chasnov: Derivation of the basic reproductive ratio from the SIR model | Differential Equations for Engineers
• The Efficient Engineer: Entendendo a Equação de Bernoulli

Referências

• Allen, LJS (2007) An Introduction to Mathematical Biology. NJ: Pearson.
• Barreira, S (2014) Matemática aplicada às ciências farmacêuticas com Excel. Vol. 2. Lisboa: Escolar Editora.
• Bates, JH & Sobel, BE (2003) The conceptual basis of mathematics in cardiology: (II). Calculus and differential equations. Coronary artery disease 14(2): 135–148.
• Batschelet, E (1979) Introduction to mathematics for life scientists. 3^rd ed. NY: Springer.
• Batschelet, E (1978) Introdução à matemática para biocientistas. Tradução da 2ª ed. São Paulo: EDUSP e Rio de Janeiro: Interciência.
• Bertalanffy, L von (1957) Quantitative Laws in Metabolism and Growth. The Quarterly Review of Biology 32(3): 217-231.
• Edelstein-Keshet, L (2005) Mathematical Models in Biology. USA: Society for Industrial and Applied Mathematics (SIAM).
• Garfinkel, A et al. (2017) Modeling Life: The Mathematics of Biological Systems. USA: Springer.
• Giordano, FR et al. (2015) A first course in mathematical modeling. 5th ed. OH: Thomson.
• Jones, DS, Plank, MJ & Sleeman, BD (2009) Differential Equations and Mathematical Biology. 2nd ed. UK: CRC.
• Kreider, DL; Kuller, RG & Ostberg, DR (1972) Equações diferenciais. Tradução: Elza F. Gomide. SP: Blucher.
• Naziazeno, R et al. (2024) Explorando dinâmicas epidêmicas: estudo de um modelo compartimental e suas aplicações. Rev. Bras. Ensino Fís. 46: e20240128. https://doi.org/10.1590/1806-9126-RBEF-2024-0128
• Poularikas, AD (1999) The Handbook of Formulas and Tables for Signal Processing. Boca Raton: CRC/Springer/IEEE.
• Savage, VM et al. (2007) Scaling of number, size, and metabolic rate of cells with body size in mammals. Proceedings of the National Academy of Sciences of the United States of America 104(11): 4718–23.
• Siqueira, JO (2012) Fundamentos de métodos quantitativos: aplicados em Administração, Economia, Contabilidade e Atuária usando WolframAlpha e SCILAB. SP: Saraiva. Soluções dos exercícios em https://www.researchgate.net/publication/326533772_Fundamentos_de_metodos_quantitativos_-_Solucoes.
• West, GB, Brown, JH & Enquist, BJ (2001) A general model for ontogenetic growth. Nature 413(6856): 628–631. doi:10.1038/35098076.
• Zwillinger, D (1997) Handbook of differential equations. 3rd ed. USA: Academic.