Aula 11 - Dinâmica Evolutiva 2.ipynb

Kernel: SageMath 9.7

Espaços de Sequências

Genomas são sequèncias de 4 nucleotídeos, denotados pelas letras A (adenina), C (citosina), T (timina) e G (guanina). Toda a informação codificada nesta sequência (genótipo) é transmitida de geração para geração, quase que perfeitamente. Durante a reprodução, dos organismos, processos geradores de variabilidade tais como mutação e recombinação podem introduzir "erros" no processo de cópia do genoma. Genomas de Comprimento L permitem $4^L$ sequencias diferentes. Modelos computacionais de evolução de sequências frequentemente simplificam esta representação para sequências binárias cuja manipulação é mais natural em computadores.

A figura abaixo mostra um exemplo de um espaço de sequências binárias, particularmente para uma sequencia de comprimento 3, $L=3$. Para sequências deste tamanho, o espaço consiste em reticulado dem 3 dimensões no qual cada vértice representa uma sequência possível. Cada dimensão no espaço de sequências assim definido, contém tantos pontos quantos valores possíveis para cada locus da sequencia. Logo, para sequências de DNA, são 4 pontos e para sequências binárias, 2. Para representar o genoma humano precisaríamos de cerca de 3 bilhões de dimensões.

A distância entre duas sequências no espaço de sequências como o ilustrado acima é o caminho mais curto, passando pelas arestas do reticulado (e não pelas diagonais). Assim, temos que a distância entre $000$ e $110$ é $2$ e não $\sqrt{2}$ (pela distância Euclidiana) . Esta distância é chamada de distância de Hamming, ou distância de Manhattan.

A Equação da Quasi-Espécie

Com o conceito de sequência como uma série de elementos que podem sofrer mutação ou outros tipos de error de replicação durante a reprodução, Podemos introduzir o conceito de quasi-espécie que é o conjunto de variantes de sequências em uma população que podem pretencer a uma mesma espécie biológica. Se atribuirmos um fitness $f_i$ ao genoma i de uma população infinita de de indivíduos de genoma com comprimento L, podemos escrever a equação da quasi-espécie:

$$\dot{x}_i=\sum_{j=0}^n x_j f_j q_{ji}−\phi x_i$$

Onde $Q$ é a matrix de mutação, como vimos anteriormente. Caso não haja erros na reprodução, a matriz de mutação é a matriz identidade e a equação da Quasi-espécie se reduz à equação da seleção, vista anteriormente. Caso existam erros, e algum $f_i$ seja maior do que os demais, o equilíbrio não maximiza o fitness médio da população.

Na equação acima, é comum combinar o vetor $\overrightarrow{f}$ com a matrix $Q$, para obter a matrix de mutação-seleção, $W$:

$$W=[w_{ij}]=\left[f_j q_{ji}\right].$$

Assim em notação vetorial,

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

cujo equilíbrio é dado por:

$$\overrightarrow{x} W=\phi \overrightarrow{x}$$

Exercício 1:

Mostre por meio de simulação que o fitness médio da população é menor que o fitness da sequência com maior $f$, quando a matrix $Q$ é diferente da identidade.

%display typeset

# Criando uma matriz de mutação
import numpy as np
from itertools import cycle
from scipy import stats as st
n = 6
Q = np.eye(n) + st.uniform(0,1).rvs((n,n))
s = Q.sum(axis=1)
s.shape = n,1
Q /= np.ones((n,n)) * s
show(matrix(Q))
pretty_print(html("Como $Q$ é uma matriz estocástica, todas as linhas devem somar 1. Então se somarmos as colunas de Q obtemos um vetor de uns:"))
Q.sum(axis=1)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 0.34439944715935816 & 0.20942515335657008 & 0.07291730805590858 & 0.1964095742930142 & 0.017276038753472552 & 0.1595724783816765 \\ 0.1407283370408197 & 0.37340428909713097 & 0.022240364812021057 & 0.07048105490363249 & 0.19293026292396273 & 0.20021569122243316 \\ 0.06705345268863451 & 0.1432300367550686 & 0.28781588896814325 & 0.08983729874236686 & 0.21453492240056 & 0.1975284004452268 \\ 0.25075160820209336 & 0.04213624017259127 & 0.19506922694823764 & 0.37827341769137557 & 0.045366648239725464 & 0.08840285874597655 \\ 0.15026592836664554 & 0.19330645863575716 & 0.03678244550457313 & 0.010865434117902788 & 0.4904660643772306 & 0.11831366899789092 \\ 0.06526320113240203 & 0.21076650018949555 & 0.1884910234519407 & 0.10516628465575893 & 0.06488921102911993 & 0.3654237795412828 \end{array}\right)$$

Como é uma matriz estocástica, todas as linhas devem somar 1. Então se somarmos as colunas de Q obtemos um vetor de uns:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|[1.|\phantom{\verb!x!}\verb|1.|\phantom{\verb!x!}\verb|1.|\phantom{\verb!x!}\verb|1.|\phantom{\verb!x!}\verb|1.|\phantom{\verb!x!}\verb|1.]|$$

a = np.array([[.1,.3,.5]])
a.shape = 3,1
a*np.ones((3,3))

f = st.uniform(0,1).rvs((n,1))
show("f=", matrix(f))
show("Q=",matrix(Q))
W = f * Q 
show("=")
show("W=f*Q=",matrix(W))
pretty_print(html("Conferindo, a soma das colunas deve nos retornar o vetor $f$:"))
W.sum(axis=1)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|f=| \left(\begin{array}{r} 0.3046252457817722 \\ 0.0311532001864665 \\ 0.4471922503706508 \\ 0.6687989051053951 \\ 0.030430053756782205 \\ 0.015151634667324831 \end{array}\right)$$

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|Q=| \left(\begin{array}{rrrrrr} 0.4194730635143879 & 0.006333386415575917 & 0.1624474612693166 & 0.1229559005560127 & 0.15853004880042393 & 0.1302601394442829 \\ 0.020515525149425746 & 0.44808079565502495 & 0.24860132268935106 & 0.08252068154937059 & 0.16919375344999887 & 0.031087921506828757 \\ 0.14051986330189542 & 0.19888725201498988 & 0.3779941778106052 & 0.10764613717452759 & 0.04228283953648077 & 0.13266973016150102 \\ 0.20786135450074783 & 0.012939002503978237 & 0.142594713908697 & 0.3378683665455292 & 0.17196634313786613 & 0.12677021940318176 \\ 0.04373568028380345 & 0.13368967114773106 & 0.10440126647994409 & 0.14329610587865407 & 0.37422610144213186 & 0.2006511747677355 \\ 0.20969100347783437 & 0.05301439205303421 & 0.04229763726082219 & 0.224762759545287 & 0.0416396597028503 & 0.42859454796017193 \end{array}\right)$$

=

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|W=f*Q=| \left(\begin{array}{rrrrrr} 0.12778208507190336 & 0.001929309393475751 & 0.04948559781579049 & 0.037455471427194506 & 0.04829225507962548 & 0.0396805269937826 \\ 0.0006391242619105483 & 0.013959150726752181 & 0.007744726772361711 & 0.0025707833118311946 & 0.005270926871527472 & 0.0009684882420833936 \\ 0.06283939389175085 & 0.08894083779861808 & 0.16903606700212845 & 0.048138518326784764 & 0.018908558164379963 & 0.05932887518698865 \\ 0.13901744630382454 & 0.00865359070781661 & 0.09536718853595359 & 0.22596599361539824 & 0.11501090200558353 & 0.08478378393681867 \\ 0.0013308791021255787 & 0.004068183879751991 & 0.003176936151260843 & 0.004360508205024998 & 0.011387720384075104 & 0.0061058260345436925 \\ 0.003177161477720887 & 0.0008032547004979032 & 0.000640878347067004 & 0.0034055232194499653 & 0.0006309089114893155 & 0.006493908011099756 \end{array}\right)$$

Conferindo, a soma das colunas deve nos retornar o vetor :

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|[0.30462525|\phantom{\verb!x!}\verb|0.0311532|\phantom{\verb!xx!}\verb|0.44719225|\phantom{\verb!x!}\verb|0.66879891|\phantom{\verb!x!}\verb|0.03043005|\phantom{\verb!x!}\verb|0.01515163]|$$

def quasispecies(t, y):
    """quasispecies model"""
    y = np.array(y)
    y.shape = 1,n
    d = np.dot(y,W) - sum(f*y)*y
    return d[0]

c = cycle(['red','blue','green', 'black', 'yellow', 'orange', 'magenta', 'gray', 'pink', 'brown'])
def plot_sol(sol):
    #fitness medio
    plots = list_plot([(j[0],sum(f*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)

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

plot_sol(T.solution)

Calculando a matriz de mutação para mutações pontuais

Considere o conjunto de todas as sequências de tamanho L. A distância de Hamming $h_{ij}$ conta o número de posições diferentes entre duas sequências $i$ e $j$. Por exemplo, a distância de Hamming entre a sequência 1010 e a sequência 1100 é 2. seja $u$ a probabilidade de que uma mutação ocorra em uma posição da sequência. Logo, $(1-u)$ é a probabilidade da sequência ser copiada corretamente. Nós podemos escrever a probabilidade da replicação da sequência $i$ gerar a sequência $j$ como

$$q_{ij}=u^{h_{ij}} (1-u)^{L-h_{ij}}$$

ou seja, uma mutação precisa ocorrer em tantas posições quantos sejam as posições diferentes entre $i$ e $j$, que é precisamente o que nos dá a distância de Hamming. E nenhuma mutação deve ocorrer nas posições remanescentes. Este cálculo assume independencia entre as mutações nas várias posições e não considera outros tipos de mutações como inserções e deleções.

Teoria de Jogos Evolutivos

A teoria dos jogos, inventada por John von Neumann, versa estratégias racionais para maximizar os retornos em um jogo cujas regras são conhecidas. Jonh Maynard Smith e George Price em um artigo publicado na revista Nature em 1973, combinaram a teoria dos jogos com conceitos de evolução gerando o conceito de jogos evolutivos.

Na teoria dos jogos evolutivos uma população infinitamente grande de jogadores, com estratégias fixas competem em rodadas do jogo. Neste contexto, o retorno esperado do cada encontro conforme explicitado na matrix de retornos (Payoff Matrix) é interpretada como o fitness da estratégia. Estratégia com maior retorno vêem sua frequência crescer, enquanto que as demais podem ser excluídas competitivamente.

A Equação do Replicador

A equação do Replicador é uma extensão da equação da Quasi-espécie, pois ela permite que o fitness de cada tipo seja definida em função da estrutura da população ao invés de ser constante para cada tipo. Na sua forma mais simples, sem mutação, pode escrita da seguinte forma:

$$\dot{x_i} = x_i [ f_i(\overrightarrow{x}) - \phi(x)], \quad \phi(\overrightarrow{x}) = \sum_{j=1}^{n}{x_j f_j(\overrightarrow{x})}$$

A função de fitness $f_i(x)$ é definida, matricialmente, como $(Ax)_i$ que é o payoff esperado de $x_i$, se assumirmos que o fitness depende linearmente da distribuição de frequências dos tipos. $A$ é a chamada matriz de payoff e contém toda a informação do fitness dos tipos. O Fitness médio da população pode ser escrito como $x^T Ax$.

Duas estratégias $A$ e $B$

seja $x_A$ a frequência de $A$ e $x_B$ a frequência de $B$. O vetor $\overrightarrow{x} = (x_A, x_B)$ define a composição da população. Seja $f_A(\overrightarrow{x})$ o fitness de $A$ e $f_B(\overrightarrow{x})$ o fitness de $B$. Então a dinâmica de seleção pode ser escrita como

$$\dot{x}_A = x_A [f_A(\overrightarrow{x}) - \phi]$$$$\dot{x}_B = x_B [f_B(\overrightarrow{x}) - \phi]$$

O fitness médio fica sendo $\phi=x_A f_A(\overrightarrow{x}) + x_B f_B(\overrightarrow{x})$.

como $x_A +x_B = 1$ sempre, podemos introduzir a variável $x$, tal que $x_A=x$ e $x_B=1-x$. As funções de fitness passam a ser apenas $f_A(x)$ e $f_B(x)$. e o fitness médio fica sendo

$$\phi=x f_A(x) (1-x) f_B(x)$$$$=x f_A(x) + f_B - x f_B(x)$$

então podemos escrever o modelo de seleção como

$$\dot{x}= x [f_A(x) - x f_A(x)+f_B(x) -x f_B(x)]$$$$\dot{x}=x[(1-x)f_A(x) + (1-x)f_B(x)]$$$$\dot{x}=x (1-x)[f_A(x)+f_B(x)]$$

Payoff matrix:

-	A	B
A	$a$	$b$
B	$c$	$d$

Agora, como o fitness de cada tipo depende da estrutura da população, definimos o fitness de cada tipo em função das interações possíveis. Logo a matriz de payoff acima pode ser interpretada como:

se $A$ encontra com $A$, ambos tem payoff $a$, se $B$ se encontra com $B$ ambos tem payoff $d$, se $A$ se encontra com $B$, $A$ tem payoff $b$ e $B$ tem payoff $c$. A Matriz de Payoff tem sua origem na teoria de jogos, na qual seus valores correspondem aos ganhos de dois jogadores. Na teoria dos jogos evolutivos, os ganhos são interpretados como fitness. Neste caso

$$f_A = a x_A +b x_B$$$$f_B = c x_A + d x_B$$

se nós plugarmos estas funções lineares de fitness no modelo de seleção, temos

$$\dot{x} =x(1-x)[(a-b-c+d)x +b -d]$$

var('x a b c d')
pretty_print(html("Estrutura da população:"))
X = matrix([[x],[(1-x)]])
show(X)
pretty_print(html("Matriz de payoff:"))
A = matrix([[a,b],[c,d]])
show(A)

Estrutura da população:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{r} x \\ -x + 1 \end{array}\right)$$

Matriz de payoff:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} a & b \\ c & d \end{array}\right)$$

Calculando o fitness médio de ambas as estratégias:

(A*X)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{r} -b {\left(x - 1\right)} + a x \\ -d {\left(x - 1\right)} + c x \end{array}\right)$$

var('a b c d')
A = matrix([[a,b],[c,d]])
f = X[0]*((A*X)[0] - (X.transpose()*(A*X))[0])
pretty_print(html("A e B são bi-estáveis se $a>c$ e $b<d$"))
show(f.expand())
pretty_print(html("Fazendo $a=3,b=1,c=2,d=2$"))
show(f(a=3,b=1,c=2,d=2).expand())
pf = plot(f(a=3,b=1,c=2,d=2),(x,0,1), legend_label=r"$a>c,\, b<d$: Bi-estável")

A e B são bi-estáveis se e

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-a x^{3} + b x^{3} + c x^{3} - d x^{3} + a x^{2} - 2 \, b x^{2} - c x^{2} + 2 \, d x^{2} + b x - d x$$

Fazendo

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-2 \, x^{3} + 3 \, x^{2} - x$$

var('a b c d')
A = matrix([[a,b],[c,d]])
g = X[0]*((A*X)[0] - (X.transpose()*(A*X))[0])
pretty_print(html("A e B coexistem se $a<c$ e $b>d$, conforme curva vermelha abaixo."))
show(g.expand())
pretty_print(html("Fazendo $a=1,b=2,c=2,d=1$"))
show(g(a=1,b=2,c=2,d=1).expand())
pg = plot(g(a=1,b=2,c=2,d=1),(x,0,1),color="red", legend_label=r"$a<c,\, b>d$: coexistência")
show(pf+pg)

A e B coexistem se e , conforme curva vermelha abaixo.

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-a x^{3} + b x^{3} + c x^{3} - d x^{3} + a x^{2} - 2 \, b x^{2} - c x^{2} + 2 \, d x^{2} + b x - d x$$

Fazendo

$$\renewcommand{\Bold}[1]{\mathbf{#1}}2 \, x^{3} - 3 \, x^{2} + x$$

h = x*(1-x)*((a-b-c+d)*x +b-d)
show(h(a=3,b=3,c=2,d=2).expand())
plot(h(a=3,b=3,c=2,d=2),(x,0,1),color="green",legend_label=r"$a>c,\, b>d$: $A$ domina")

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-x^{2} + x$$

Há cinco possibilidades para a dinâmica da seleção entre dois tipos: A domina B, B domina A, A e B são Bi-estáveis, A e B coexistem em um equilíbrio estável, e A e B são variantes Neutras um do outro. Esta dinâmica é denominada seleção dependente da frequência.

Exercício 2:

Construa o Gráfico das outras 3 situações.

O Equilíbrio de Nash

Nó já vimos que os tipos dos modelos de seleção podem ser interpretados como as estratégias de um jogo. O equilíbrio de Nash, é um equilíbrio no qual nenhum dos jogadores pode desviar de sua estratégia e obter um melhor payoff.

no caso da matriz de payoff para duas estratégias, A e B, as seguintes condições levam a um equilíbrio de Nash:

1. A é um equilíbrio e Nash estrito se $a>c$

2. A é um equilíbrio de Nash se $a\geq c$

3. B é um equilíbrio de Nash estrito, se $d>b$

4. B é um equilíbrio de Nash se $d\geq b$ Consideremos o seguinte exemplo de jogo:

matrix([[3,0],[5,1]])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 3 & 0 \\ 5 & 1 \end{array}\right)$$

Se ambos os jogadores escolherem a estratégia A: a1,1, então um dos jogadores pode melhorar o seu retorno trocando para a estratégia B. Se ambos jogarem B nenhum jogador consegue melhores condições mudando de estratégia. Logo a estratégia B é um equilíbrio de Nash. Note que apesar de B dominar A, o Payoff de jogar no equilíbrio de Nash é menor do que jogar a estratégia dominada A. Este é um exemplo do famoso dilema do prisioneiro.

Considere este outro jogo:

matrix([[5,0],[3,1]])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 5 & 0 \\ 3 & 1 \end{array}\right)$$

Neste caso ambas estratégias são esquilíbrios de Nash.

Estratégia Evolutivamente Estável (ESS)

Este conceito foi introduzido pelo biólogo John Maynard Smith sem que ele tivesse conhecimento dos conceito de equilíbrio de Nash na teoria de jogos.

Imagine uma grande população de jogadores do tipo A. Uma quantidade infinitesimal $\epsilon$ de um mutante tipo B, é introduzida na população. Se temos a matriz de payoff abaixo, e as frequências de B e A são, respectivamente, $\epsilon$ e $1-\epsilon$, que condições seletivas impedem a invasão de B em A?

var('a b c d')
matrix([[a,b],[c,d]])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} a & b \\ c & d \end{array}\right)$$

Para que B não invada, o fitness de A tem que ser maior que o fitness de B, logo:

$$a(1-\epsilon) + b\epsilon > c(1-\epsilon) + d \epsilon$$

Cancelando os termos com $\epsilon$, chegamos à seguinte conclusão:

$$a>c$$

Contudo, se $a=c$, então temos que

$$b>d$$

Em suma, A é uma ESS se:

1. $a>c$

2. $a=c$ e $b>d$

Esta conclusão no entanto só vale para populações infinitamente grandes e quantidade infinitesimalmente pequenas da estratégia invasora.

Mais de duas estratégias

Seja o payoff da estratégia S_i versus a estratégia S_j denotado por $E(S_i,S_j)$

1. A estratégia $S_k$ é um equilíbrio de Nash estrito se: $$E(S_k,S_k)>E(S_i,S_k)\,\,\,\,\,\ \forall i$$

2. A estratégia $S_k$ é um equilíbrio de Nash se: $$E(S_k,S_k) \geq E(S_i,S_k)\,\,\,\,\,\ \forall i$$

3. A estratégia $S_k$ é ESS se $\forall i \neq k$ temos uma das duas situações abaixo: $$E(S_k,S_k)>E(S_i,S_k)$$ ou $$E(S_k,S_k)=E(S_i,S_k) \,\,\, e \,\,\, E(S_k,S_i)>E(S_i,S_i)$$ Note que a ESS garante que a seleção se oporá a qualquer invasor potencial. O mesmo é verdadeiro para um equilíbrio de Nash estrito, mas não para um esquilíbrio de Nash. Se $E(S_k,S_k)=E(S_j,S_k)$ e E(S_k,S_j)<E(S_j,S_j) então $S_k$ é um equilíbrio de Nash mas a seleção favorecerá a invasão de $S_k$ por $S_j$. Portanto, faz sentido adicionar uma quarta definição:

4. A estratégia $S_k$ é ESS contra invasões por meio de seleção (ESS fraco) se $\forall i \neq k$ temos uma das duas situações $$E(S_k,S_k)>E(S_i,S_k)$$ ou $$E(S_k,S_k)=E(S_i,S_k)\,\,\, e \,\,\, E(S_k,S_i)\geq E(S_i,S_i)$$ Logo,

Nash estrito $\Rightarrow$ ESS $\Rightarrow$ ESS fraco $\Rightarrow$ Nash.

Três estratégias: Pedra Papel e Tesoura

Este é um jogo interessante no simplex $S_3$ no qual temos três estratégias $Pe$ (pedra), $Pa$ (papel) e $T$ (tesoura). A pedra domina a tesoura que domina o papel, que por sua vez domina a pedra. Este jogo apresenta uma estratégia cíclica de dominação, que pode ser descrita por uma matriz $3x3$ como a apresentada abaixo (ordem das linhas: $Pe$, $T$ e $Pa$).

matrix([[4,2,1],[3,1,3],[5,0,2]])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 4 & 2 & 1 \\ 3 & 1 & 3 \\ 5 & 0 & 2 \end{array}\right)$$

Na matriz acima a ordem das linhas e colunas é a da dominância: $Pe \rightarrow T \rightarrow Pa$. Antes de analisarmos este modelo convém observar que a dinâmica da equação replicante não se altera se adicionarmos uma constante a cada coluna da matriz. Desta forma, podemos simplificar a nossa matriz de payoff, subtraindo o elemento da diagonal de cada coluna para obter: uma matriz de payoff com a diagonal 0.

matrix([[0,1,-1],[-1,0,1],[1,-1,0]])

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 0 & 1 & -1 \\ -1 & 0 & 1 \\ 1 & -1 & 0 \end{array}\right)$$

Vamos simular a dinâmica deste jogo:

html("Seja $X$ o nosso vetor de frequências (pedra, tesoura e papel):")
X = matrix([[0.2],[0.5],[0.3]])
show(X)
pretty_print("Seja a matriz de Payoff")
A = matrix([[0,1,-1],[-1,0,1],[1,-1,0]])
show(A)
pretty_print(html("Podemos então definir o nosso sistema de equações diferenciais"))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{r} 0.200000000000000 \\ 0.500000000000000 \\ 0.300000000000000 \end{array}\right)$$

Seja a matriz de Payoff

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 0 & 1 & -1 \\ -1 & 0 & 1 \\ 1 & -1 & 0 \end{array}\right)$$

Podemos então definir o nosso sistema de equações diferenciais

import numpy as np
A = np.matrix([[0,1,-1],[-1,0,1],[1,-1,0]])
def fun(t,X):
    X = np.matrix(X).T
    return [float(X[0]*((A*X)[0] - (X.T*(A*X))[0])),
            float(X[1]*((A*X)[1] - (X.T*(A*X))[0])),
            float(X[2]*((A*X)[2] - (X.T*(A*X))[0])),
    ]
r = fun(0,[.56,.22,.22])
show(r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left[7.212008767965017 \times 10^{-19}, -0.07480000000000002, 0.07480000000000002\right]$$

T = ode_solver()
T.function = fun
inits = [1/3,1/3,1/3]
inits= [.56,.22,.22]
tspan = [0,100]
T.ode_solve(tspan, inits, num_points=500)

def plot_sol(sol, labels=['Pedra','Tesoura','Papel']):
    a=list_plot([(i[0],i[1][0]) for i in sol],color='red', plotjoined=True, legend_label=labels[0], alpha=.8)
    c=list_plot([(i[0],i[1][1]) for i in sol],color='blue', plotjoined=True, legend_label=labels[1], alpha=.8, axes_labels=["t","freq"], gridlines=True)
    r = list_plot([(i[0],i[1][2]) for i in sol],color='green', plotjoined=True, legend_label=labels[2], alpha=.8, axes_labels=["t","freq"], gridlines=True)
    a.legend()
    c.legend()
    r.legend()
    show(c+a+r)
plot_sol(T.solution)

np.linalg.det(A)

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

point3d([(i[1][0],i[1][1],i[1][2]) for i in T.solution], size=5)

inits = [.56,.22,.22]
A = np.matrix([[0,1,-1.1],[-1.1,0,1],[1,-1.1,0]])
show(matrix(A))
pretty_print(html("Podemos calcular o determinante da matriz de Payoffs:"))
show(np.linalg.det(A))
pretty_print(html("Como o determinante é positivo, temos um único equilíbrio interno, globalmente estável:"))
pa = A[2,0]*A[0,1]*A[1,2]
pb = A[1,0]*A[2,1]*A[0,2]
#print pa, pb
show("a1a2a3>b1b2b3 ",pa > pb)
T.ode_solve([0,350], inits, num_points=500)
plot_sol(T.solution)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 0.0 & 1.0 & -1.1 \\ -1.1 & 0.0 & 1.0 \\ 1.0 & -1.1 & 0.0 \end{array}\right)$$

Podemos calcular o determinante da matriz de Payoffs:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-0.33100000000000024$$

Como o determinante é positivo, temos um único equilíbrio interno, globalmente estável:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|a1a2a3>b1b2b3| \verb|True|$$

p3 = point3d(inits, size=5, color="green")
l3 = line3d([(i[1][0],i[1][1],i[1][2]) for i in T.solution], size=10)
show(p3+l3)

inits = [.56,.22,.22]
del(A)
A = np.matrix([[0,1.1,-1],[-1,0,1.1],[1.1,-1,0]])
show(matrix(A))
pretty_print(html("Podemos calcular o determinante da matriz de Payoffs:"))
pa = A[2,0]*A[0,1]*A[1,2]
pb = A[1,0]*A[2,1]*A[0,2]
#print pa, pb
show("pa>pb: ",pa > pb)
tspan=[0,350]
show(np.linalg.det(A))
pretty_print(html("Como o determinante é positivo, temos um único equilíbrio interno, globalmente estável:"))
#inits = [1/3,1/3,1/3]
T.ode_solve(tspan, inits, num_points=500)
plot_sol(T.solution)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 0.0 & 1.1 & -1.0 \\ -1.0 & 0.0 & 1.1 \\ 1.1 & -1.0 & 0.0 \end{array}\right)$$

Podemos calcular o determinante da matriz de Payoffs:

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\verb|pa>pb:| \verb|True|$$

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

Como o determinante é positivo, temos um único equilíbrio interno, globalmente estável:

l=line3d([(i[1][0],i[1][1],i[1][2]) for i in T.solution], size=2)
p=point3d([(i[1][0],i[1][1],i[1][2]) for i in T.solution], size=2,color='red')
show(l+p)

Replicante Mutante

Se adicionarmos mutação à equação do replicante acima temos a equação do replicante mutante:

$$\dot{x_i} = \sum_{j=1}^{n}{x_j f_j(x) Q_{ji}} - \phi(x)x_i,$$

Para uma discussão mais aprofundada sobre este modelo, veja este artigo.