Aula 10 - Dinâmica Evolutiva.ipynb

Kernel: SageMath 9.5

Introdução à Evolução

Seleção natural

Seleção natural é o fenômeno relacionado com a diferença na capacidade de diferentes indivíduos de formar descendentes, ou seja, crescer. Esta diferença decorre do efeito do ambiente sobre os indivíduos, cada qual com seu grau de adaptação. Esta capacidade de reprodução também representada pela taxa de reprodução, em modelos evolutivos, Corresponde ao grau de adaptação do indívíduo ao seu ambiente e é chamada de Fitness.

Para modelar os efeitos da seleção natural, precisamos de pelo menos dois "tipos" de indivíduos. Vamos chamá-los de A e B. Denotaremos o fitness de A por $a$ e o de B por $b$. Seja $x(t)$ o número de indivúduos do tipo A no tempo $t$ e $y(t)$ o número de indivíduos do tipo B no tempo $t$. Sejam ainda $x_0$ e $y_0$ os números de A e B no tempo $t=0$.

Cujas soluções são:

%display typeset

var("x x_0 y_0 a b y t")
x  = function('x')(t)
y  = function('y')(t)
xdot = diff(x,t)==a*x
ydot = diff(y,t)==b*y
show(xdot)
show(ydot)
solx = desolve(xdot,x, ics=[x_0], ivar=t)
soly = desolve(ydot,y, ics=[y_0], ivar=t)
show(solx.subs(_C=x_0))
show(soly.subs(_C=y_0))
xdot/ydot

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\partial}{\partial t}x\left(t\right) = a x\left(t\right)$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\partial}{\partial t}y\left(t\right) = b y\left(t\right)$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}x_{0} e^{\left(a t\right)}$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}y_{0} e^{\left(b t\right)}$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\frac{\partial}{\partial t}x\left(t\right)}{\frac{\partial}{\partial t}y\left(t\right)} = \frac{a x\left(t\right)}{b y\left(t\right)}$$

Seja $t=τ$ o tempo até que a população A dobre de tamanho, é fácil encontrar que $τ=log_2 /a$. De maneira similar, B leva $log_2 /b$ para dobrar de tamanho.

var('tau')
f = x_0*exp(a*tau)==2*x_0
solve(f,tau)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\tau = \frac{\log\left(2\right)}{a}\right]$$

Seja a razão entre as duas populações: $ρ(t)=x(t)/y(t)$. Logo temos a seguinte equação diferencial

$$\dot{ρ}=\frac{\dot{x}y−x\dot{y}}{y^2}=(a−b)ρ,$$

onde $\dot{\rho}=\frac{d\rho}{dt}$. Abaixo temos a derivação de $\dot{ρ}$ usando o sage:

var('rho x y t')
x = function('x')(t)
y = function('y')(t)
rho = function('rho')(t) == x/y
drho = diff(rho,t)
drho.simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\partial}{\partial t}\rho\left(t\right) = \frac{y\left(t\right) \frac{\partial}{\partial t}x\left(t\right) - x\left(t\right) \frac{\partial}{\partial t}y\left(t\right)}{y\left(t\right)^{2}}$$

cuja solução é $ρ(t)=ρ_0e^{(a−b)t}$. Se $a>b,$ $ρ$ tende a infinito, em cujo caso A vencerá a competição com B, e se $a<b$, $ρ$ tenderá a zero, dando a vantagem a B.

Agora imaginemos o caso em a população total é constante, por exemplo em decorrência da existência de uma capacidade de suporte ambiental. Neste caso, para simplificar vamos considera que tanto $x(t)$ quanto $y(t)$ representam a abundância relativa de seus tipos respectivos, ou como também poderíamos chamar, a sua frequência, logo $x+y=1$. Temos então o seguinte sistema de equações: $$\begin{align} \dot{x}&=(a−ϕ)x\\ \dot{y}&=(b−ϕ)y \end{align}$$ Sendo $ϕ$ o fitness médio da população, $ϕ=ax+by$, o crescimento/decrescimento de cada população, passa a ser relativo ao fitness médio, conforme as equações acima.

Podemos simplicar o sistema acima, substituindo $y$ por $1−x$:

$$\dot{x}=x(1−x)(a−b)$$

A equação acima apresenta dois equilíbrios: $x=0$ e $x=1$. No primeiro, todos os indivíduos são do tipo B e no segundo todos os indivíduos são do tipo A. Após o sistema atingir qualquer dos equilíbrios, nada mais pode ocorrer pois já não haverão indivíduos do outro tipo. Qual dos equilíbrios irá absorver a dinâmica do sistema depende dos valores $a$ e $b$. Se $a<b$ o sistema convergirá para o equilíbrio $x=0$, e no caso contrário, para $x=1$. É a "sobrevivência do mais apto".

var("a b x y")
c=(a-(a*x+b*y))*x
s=x*(1-x)*(a-b)
show(c.subs(y=1-x))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{\left(b {\left(x - 1\right)} - a x + a\right)} x$$

Sobrevivência do mais Apto

Este modelo pode ser extendido para descrever a seleção entre $n$ tipos diferentes: Seja $x_i(t)$ a frequência do tipo $i$, $(i=1…n)$. Seja $f_i$ o fitness do tipo $i$. Assim o fitness médio da população é dado por

$$ϕ=\sum_{i=1}^n x_if_i$$

A dinâmica da seleção continua como antes:

$$\dot{x}_i=x_i(f_i−ϕ)$$

Da mesma forma a frequencia de um dado tipo aumenta apenas se a sua fitness (aptidão) for maior que a da média da população, caso contrário diminuirá. A população total continua constante: $\sum_{i=1}^n x_i=1$ e $\sum^n_{i=1} \dot{x}_i=0$.

Exercício:

Implementar uma simulação de competição com 10 tipos para visualizar o efeito das extinções sobre o fitness médio da população e a dinâmica dos outros tipos

import numpy as np
n = 10
fis = np.array([random() for i in range(n)])
def evolucao(t,y):
    phi = sum(np.array(y)*fis)
    d = y*(fis-phi)
    return d

T = ode_solver()
T.algorithm = "rk8pd"
T.function = evolucao
T.ode_solve(y_0=[1/n]*n,t_span=[0,50], num_points=500)

from itertools import cycle

def plot_sol(sol):
    #fitness medio
    c = cycle(['red','blue','green', 'black', 'yellow', 'orange', 'magenta', 'gray', 'pink', 'brown'])
    plots = list_plot([(j[0],sum(fis*np.array(j[1]))) for j in sol], plotjoined=True, linestyle='-.', legend_label=r"$\phi$")
    for i in range(len(sol[0][1])):
        co = c.__next__()
        plots += list_plot([(j[0],j[1][i]) for j in sol],color=co, plotjoined=True, legend_label='$x_{%s}$'%i, alpha=.8, gridlines=true)
    show(plots)

show("Initial conditions: ", np.array(T.solution[0][1]))
show("fi = ",fis)
show("Fitness medio: ",sum(fis*np.array(T.solution[0][1])))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Initial|\verb| |\verb|conditions:| \verb|[0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1|\verb| |\verb|0.1]|$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|fi|\verb| |\verb|=| \begin{array}{l} \verb|[0.89603651|\verb| |\verb|0.15176991|\verb| |\verb|0.10013761|\verb| |\verb|0.83357806|\verb| |\verb|0.81809195|\verb| |\verb|0.63081744|\\ \verb| |\verb|0.3875324|\verb| |\verb|0.96596628|\verb| |\verb|0.75777258|\verb| |\verb|0.01459507]| \end{array}$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Fitness|\verb| |\verb|medio:| 0.5556297821898809$$

plot_sol(T.solution)

O Simplex

o conjunto de pontos com a propriedade $\sum_{i=1}^n x_i=1$ é chamado simplex $S_n$. Em modelos Evolutivos como os que estamos estudando, as população evoluem no simplex. O Simplex $S_n$ é uma estrutura $(n−1)$-dimensional embutida em um espaço Euclidiano $n$-dimensional. Então para dois tipos, temos um segmento de reta, $S_2$, para três , um triângulo, $S_3$ para 4 um tetraedro, $S_4$ e assim por diante. o simplex $S_n$ tem $n$ faces e cada uma delas é um simplex $S_{n−1}$. Os vértices de um simplex no contexto de nosso modelo evolutivo, são os pontos onde a penas um tipo está presente, e os demais foram extintos. Os pontos do simplex, que não pertencem a nenhuma das faces ou são vértices, são pontos internos, que representam estados onde há co-existência entre todos os $n$ tipos.

Relaxando a premissa de taxas de crescimento lineares

No modelo de dois tipos apresentado acima, assumimos que a taxas de crescimento eram funções lineares da frequência dos tipos. Podemos relaxar esta premissa re-escrevendo as equações da seguinte forma: $$\begin{align} \dot{x}&=ax^c−ϕx\\ \dot{y}&=by^c−ϕy \end{align}$$

Neste novo modelo se $c=1$ temos nosso modelo antigo, mas se $c<1$ temos um crescimento sub-exponencial, ou seja na ausência da limitação por densidade, $ϕ$, a curva de crescimento seria mais lenta do que uma exponencial. Quando $c>1$, o crescimento torna-se super-exponencial e a curva de crescimento hiperbólica.

var('x y a  b c phi t')
x = function('x')(t)
f = diff(x,t) == a*x^c -phi*x
pretty_print(html('$c>1$'))
show(desolve(f(c=2),x, ics=[0], ivar=t, contrib_ode=True))
pretty_print(html('$c<1$'))
show(desolve(f(c=0.9),x, ics=[0], ivar=t, contrib_ode=True))
pretty_print(html('$c=1$'))
show(desolve(f(c=1),x, ics=[0], ivar=t, contrib_ode=True))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{\log\left(a x\left(t\right) - \phi\right) - \log\left(x\left(t\right)\right)}{\phi} = C + t$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{10 \, \log\left(\phi x\left(t\right)^{\frac{1}{10}} - a\right)}{\phi} = C + t$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}C e^{\left({\left(a - \phi\right)} t\right)}$$

Para manter a população constante agora, $x+y=1$, fazemos $ϕ=ax^c+by^c$. Então, de forma similar a como fizemos anteriormente, podemos simplificar o sistema acima a

$$\dot{x}=x(1−x)f(x)$$

onde

$$f(x)=ax^{c−1}−b(1−x)^{c−1}$$

a=.5
b=.6

@interact
def variando_c(c=(0.2,1.8,.2)):
    phi(x) = a*x^(c-1)-b*(1-x)^(c-1)
    P = plot((x*(1-x)*phi).subs(a=.5, b=.6, c=c),(x,0,1))
    show(P)

Exercício 1:

Encontrar o equilíbrio interno ao simplex acima, $x^∗$, como uma expresão dos parâmetros do modelo, decrevendo a sua estabilidade, quando $c>1$ e quando $c<1$.

s = solve(x*(1-x)*phi,x)
s

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = 0, x = 1\right]$$

Então para $c≠1$, existe apenas um outro ponto fixo entre 0 e 1:

$$x^∗=\frac{1}{1+\sqrt[c-1]{a/b}}$$

Invasibilidade

Invasão significa que uma quantidade infinitesimal de um dos tipos consegue crescer em uma população quase que completamente constituída pelo outro tipo.

Exercício 2:

Simule a dinâmica de um modelo com dois tipos e mostre o efeito do valor de c sobre a invasibilidade de um dos tipos.

Exercício 3:

Construa o diagrama de bifurcação do sistema $\dot{x}=x(1−x)f(x)$, onde $f(x)=ax^{c−1}−b(1−x)^{c−1}$. Encontre o ponto de bifurcação e identifique o tipo de bifurcação.

Mutação

O processo natural de mutação, é a principal fonte de geração de variabilidade, que é a matéria prima sobre a qual a seleção natural atua. Mutações são erros que ocorrem durante a replicação do material genético de um indivíduo durante a reprodução.

Vamos explorar a forma mais simples de inserir o conceito de mutação no modelo evolutivo que estamos desenvolvendo. Chamemos de $u_1$ a taxa de mutação de A para B, ou seja, $u_1$ é a probabilidade de que a reprodução de um indivíduo do tipo A gere um indivíduo do tipo B. De maneira similar, vamos chamar de $u_2$ a taxa de mutação de B para A. Isto nos leva às seguintes equações: $$\begin{align} \dot{x}&=x(1−u_1)+yu_2−ϕx\\ \dot{y}&=xu_1+y(1−u_2)−ϕy \end{align}$$ Se assumirmos que A e B têm o mesmo fitness, ou seja, ($a=b=1$), o fitness médio da população, é constante e dado por $ϕ=ax+by=1$, dado que $x+y=1$. Com isso, o sistema acima se reduz à seguinte equação diferencial:

$$\dot{x}=u_2−x(u_1+u_2).$$

Que apresenta o seguinte equilíbrio para a frequência do tipo A:

var('x u_1 u_2')
solve(u_2-x*(u_1+u_2),x)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = \frac{u_{2}}{u_{1} + u_{2}}\right]$$

Concluímos que neste caso de igual fitness a mutação leva à co-existência dos tipos, e as suas frequências no equilíbrio dependerá das taxas de mutação.

Exercício 4:

Encontre a razão entre as frequências de A e B no equilíbrio.

Matriz de Mutação

Podemos extender a dinâmica de mutação para $n$ tipos. Para isso precisamos introduzir a matriz de mutação, $Q=[q_{ij}]$. A probabilidade e do tipo $i$ mutar para o tipo $j$ é dada por $q_{ij}$. Como cada tipo gera, na reprodução, um outro indivíduo do seu tipo ou de outro tipo, temos que $\sum^n_{j=1}q_{ij}=1$. Logo $Q$ é uma matriz estocástica $n\times n$.

Matrizes estocásticas têm as seguintes propriedades:

1. Todos os elementos são números no intervalo [0,1] (probabilidades)

2. São quadradas

3. A soma de cada linha é 1

Matrizes estocásticas têm sempre um autovalor igual a 1 , e nenhum autovalor pode ter um valor absoluto maior do que 1.

O dinâmica de mutação fica então:

$$\dot{x}_i=\sum_{j=1}^n x_jq_{ji}−ϕx_i$$

para $i=1,\ldots,n.$

Ou, em notação vetorial:

$$\dot{\overrightarrow{x}}=\overrightarrow{x} Q−\phi\overrightarrow{x}$$

Novamente o fitness médio é $ϕ=1$. O equilíbrio é dado pelo autovetor à esquerda, associado com o autovalor 1:

$$\overrightarrow{x^∗}Q=\overrightarrow{x^∗}$$

Exercício 5:

Simule um modelo de seleção com mutação, com matrizes de mutação simétricas e assimétricas.

import numpy as np
n = 10
# fis = np.ones(n)/n#
fis = np.array([random() for i in range(n)])
Q =np.random.random((n,n))
qs = Q.sum(axis=1)
for i,r in enumerate(Q):
    Q[i,:] = r/qs[i]
print(Q.sum(axis=1))
def evolucaoQ(t,y):
    y = np.array(y).reshape((1,n))
    phi = sum(y*fis)
    d = y@Q - phi*y
    return d[0]

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

T = ode_solver()
T.algorithm = "rk8pd"
T.function = evolucaoQ
T.ode_solve(y_0=[1/n]*n,t_span=[0,20], num_points=500)

plot_sol(T.solution)

Q.sum(axis=0)/Q.sum(axis=1)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{l} \verb|[1.10908426|\verb| |\verb|0.90846055|\verb| |\verb|0.96265006|\verb| |\verb|1.10370426|\verb| |\verb|1.14345312|\verb| |\verb|1.04049078|\\ \verb| |\verb|1.10628897|\verb| |\verb|0.70135343|\verb| |\verb|0.86568679|\verb| |\verb|1.05882779]| \end{array}$$

fis

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{l} \verb|[0.35442058|\verb| |\verb|0.56695803|\verb| |\verb|0.09375394|\verb| |\verb|0.20363271|\verb| |\verb|0.45522536|\verb| |\verb|0.31394509|\\ \verb| |\verb|0.59047749|\verb| |\verb|0.57565622|\verb| |\verb|0.10331312|\verb| |\verb|0.77935505]| \end{array}$$

sum(fis)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}4.036737577547483$$