Aula 16 - Teoria Evolutiva de Grafos.ipynb

Kernel: SageMath 9.7

Teoria Evolutiva de Grafos

Neste notebook vamos explorar modelos evolutivos em que os invidíduos são representados pelos nós de um grafo, ao invés de pertencerem a uma população homogênea. A arestas do grafo representam as interações competitivas. Se existe uma aresta entre os nós $i$ e $j$, significa que a prole de $i$ pode vir a substituir $j$.

As principais questões que a teoria evolutiva de grafos visa responder são:

1. Será que a topologia dos grafos pode influenciar a taxa de evolução?

2. Podemos encontrar uma topologia que reduza a probabilidade de fixação de um mutante vantajoso?

3. Será que existem grafos capazes de anular completamente o efeito da seleção?

4. É possivel caracterizar os grafos que apresentam dinâmicas evolutivas similaers a populações não estruturadas, ou seja homogêneas?

5. etc.

Seja uma população com $N$ indivíduos $i$, $i=1,2,\ldots , N$. A probabilidade da prole de $i$ substituir $j$ é denotada por $w_{ij}$. Logo, o processo é determinado por uma matriz $W=[w_{ij}]$.

A matriz $W$ é estocástica pois é composta apenas por probabilidades e, como a prole de cada nó tem que necessáriamente ir para algum lugar, $\sum_{j=1}^N w_{ij} = 1$.

A matriz $W$ define um grafo direcionado ponderado.

Esta idéia deriva do Processo de Moran, um processo estocástico muito usado em biologia para simular evolução em populações finitas. No processo de Moran, a cada passo de tempo um indivíduo é selecionado aleatóriamente para se reproduzir e outro para morrer (de maneira a manter a população constante). A probabilidade de seleção de um tipo para reprodução, é proporcional ao seu fitness. em um grafo completo a probabilidade de fixação de um tipo é dada por:

$$\rho = \frac{1-1/r}{1-1/r^n}$$

Onde $r$ é o fitness relativo (razão entre o fitness do indivíduo e o fitness médio da população) do tipo em questão.

Uma população não estruturada, pode ser representada por um grafo completo. Seu processo evolutivo é equivalente ao do processo de Moran.

var('r n')
n=50
plot((1-r^-1)/(1-r^-n), (r,0,2))

Ciclo Direcionado

Qual a probabilidade de fixação de um mutante que surja em uma posição aleatória de um ciclo direcionado?

%display typeset

g = digraphs.Circuit(5)
g

g.adjacency_matrix()

Inicialmente, todos os indivíduos são do tipo $A$. Depois de algum tempo, um mutante $B$ surge com fitness relativo $r$. Este indivíduo dá origem a uma linhagem que tem dois destinos possíveis: a extinção ou a dominação total.

Seja $m$ o número de indivíduos $B$. Para reduzir m em 1, o indivíduo $A$ imediatamente antecedente ao cluster de $B$ deve ser selecionado para reprodução. Logo a probabilidade de passar de $m$ para $m-1$ é dada por

$$P_{m,m-1}=\frac{1}{N-m+rm}$$

Para aumentar $m$ em 1, o indivíduo B no final do cluster, tem que ser escolhido para reprodução. logo a probabilidade de passar de $m$ para $m+1$ é dada por

$$P_{m,m+1}=\frac{r}{N-m+rm}$$

A razão destas duas probabilidades é

$$\gamma_m = \frac{P_{m,m-1}}{P_{m,m+1}} = \frac{1}{r}$$

Pela teoria do processo de Moran, a probabilidade de fixação de um processo de nascimento e morte é dada por

$$\rho=\frac{1}{1+\sum_{k=1}^{N-1}\prod_{m=1}^k \gamma_m}$$

daí obtemos

$$\rho=\frac{1-1/r}{1-1/r^N}$$

A probabilidade de fixação do ciclo direcionado é a mesma do processo de Moran. Pode-se mostrar que o mesmo é verdadeiro para o ciclo não direcionado.

Grafos de Linha e Estrela

Nestes dois tipos de grafo , pode-se verificar graficamente que a probabilidade de fixação $\rho=1/N$, ou seja maior do que no processo de Moran, ou nos grafos apresentados acima. Estes grafos são supressores de seleção, pois a probabilidade de fixação não depende do fitness.

Exercício 1:

Calcule a distribuição dos tempos até a fixação de grafos em linha e em estrela. Considere uma taxa de nascimento constante. Compare as distribuições, Discuta.

import numpy as np
import pylab as P
from numpy.random import multinomial as mnrvs
r=0.9
g = np.array([1, 0, 0, 0, 0])

N = len(g)
passos = []
for i in range(1000):
    step = 1
    while sum(g)<5:
        fitness = np.array([r if i==1 else 1 for i in g])
        m = sum(g)
        n = len(g)-m
        repr = list(mnrvs(1, fitness/fitness.sum(), 1)[0]).index(1)
        #print("vertice reproduzindo no passo %s: %s"%(step,repr))
        if repr < 4:
            g[repr+1] = g[repr]
        #print("m=%s"%sum(g))
        step +=1
    passos.append(step)
    g = np.array([1, 0, 0, 0, 0])
P.clf()
P.hist(passos, bins=30)
P.title(u"Grafo linear - tempo até fixação");

from numpy.random import randint
r=1.2
g = np.array([1, 0, 0, 0, 0])
fitness = np.array([r,1,1,1,1])
N = len(g)
Spassos = []
for i in range(1000):
    step = 1
    while sum(g)<5:
        fitness = np.array([r if i==1 else 1 for i in g])
        m = sum(g)
        n = len(g)-m
        repr = list(mnrvs(1, fitness/fitness.sum(), 1)[0]).index(1)
        #print("vertice reproduzindo no passo %s: %s"%(step,repr))
        if repr == 0:
            destino  = randint(1,5,1)[0]
            g[destino] = g[0]
        #print(i,"m=%s"%sum(g),g)
        step +=1
    Spassos.append(step)
    #print(g,Spassos)
    g = np.array([1, 0, 0, 0, 0])
P.clf()
P.hist(passos,bins=30, color='b', label='Linear', alpha=0.5)
P.hist(Spassos,bins=30, color='r', label='Estrela', alpha=0.5)
P.legend();

O ciclo bidirecional

No caso do ciclo bidirecional é fácil ver que:

$P_{m,m-1}=\frac{1}{N-m+rm}$ e $P_{m,m+1}=\frac{r}{N-m+rm}$

Logo, apresenta a mesma prbabilidade de fixação do ciclo direcionado discutido acima.

Equilibrando Seleção e Deriva

Se um grafo $G$ tem a mesma probabilidade de fixação do que um processo de Moran, dizemos que ele é $\rho$-equivalente ao processo de Moran. Em tal grafo seleção e deriva estão equilibradas.

Mutante mais adaptado, $r>1$: Se $\rho_G>\rho_M$, então o grafo $G$ favorece a seleção em detrimento da deriva, ou seja, aumenta a probabilidade de fixação de um mutante mais adaptado. Logo, $G$ é um amplificador de seleção.

Se $\rho_G<\rho_M$, então o grafo $G$ favorece a deriva em detrimento da seleção, ou seja, reduz a probabilidade de fixação de um mutante mais adaptado. Logo, $G$ é um supressor de seleção.

Mutante menos adaptado, $r<1$ Se $\rho_G>\rho_M$, o grafo $G$ é supressor de seleção.

Se $\rho_G=1/N$ para qualquer $r$, então o grafo apresenta o máximo de supressão de seleção, eleiminando-a completamente.

Exercício 2:

Usando a biblioteca de grafos do Sage, gere grafos aleatórios usando os geradores descritos aqui. e implemente uma função para calcular o $\rho_G$.

grafo=digraphs.RandomDirectedGNP(8,.2)
grafo

Exercício 2:

1. Implemente uma função para simular a reprodução em grafos de acordo com as regras descritas no início desta planilha.

2. Utilizando os grafos do exercício anterior, construa um gráfico do tempo até a fixação em cada grafo em função do tamanho do grafo. Escolha com parcimônia o tamanho dos grafos. Como o processo de Moran é estocástico você terá que calcular multiplas vezes a simulação para cada par (grafo, tamanho).

np.random.randint(8)

import random
def SimulaMoran(grafo,r):
    n = grafo.order()
    a = np.random.randint(8)
    mutantes=[a]
    nummut=1
    yield mutantes,nummut
    nãomutantes=list(range(a))+list(range(a+1,n))
    while (nummut!=0)&(nummut!=n):
        x=random.uniform(0,n-nummut+nummut*r)
        if x<n-nummut:
            a=nãomutantes[floor(x)]
            vizinhos=grafo.neighbors_out(a)
            if vizinhos==[]:
                continue
            b=random.choice(grafo.neighbors_out(a))
            if b in mutantes:
                mutantes.remove(b)
                nãomutantes.append(b)
                nummut-=1
            yield mutantes,nummut
        else:
            a=mutantes[floor((x-n+nummut)/r)]
            vizinhos=grafo.neighbors_out(a)
            if vizinhos==[]:
                continue
            b=random.choice(grafo.neighbors_out(a))
            if b in nãomutantes:
                nãomutantes.remove(b)
                mutantes.append(b)
                nummut+=1
        yield mutantes,nummut

P=[]
i=0
for mutantes,nummut in SimulaMoran(grafo,1.5):
 P.append((i,nummut))
 i+=1
line(P)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [9], in <cell line: 3>()
      1 P=[]
      2 i=Integer(0)
----> 3 for mutantes,nummut in SimulaMoran(grafo,RealNumber('1.5')):
      4  P.append((i,nummut))
      5  i+=Integer(1)
Input In [7], in SimulaMoran(grafo, r)
     21     yield mutantes,nummut
     22 else:
---> 23     a=mutantes[floor((x-n+nummut)/r)]
     24     vizinhos=grafo.neighbors_out(a)
     25     if vizinhos==[]:
File /ext/sage/9.7/src/sage/functions/other.py:578, in Function_floor.__call__(self, x, **kwds)
    575     return r"\left \lfloor %s \right \rfloor"%latex(x)
    577 #FIXME: this should be moved to _eval_
--> 578 def __call__(self, x, **kwds):
    579     """
    580     Allows an object of this class to behave like a function. If
    581     ``floor`` is an instance of this class, we can do ``floor(n)`` to
   (...)
    593         -3
    594     """
    595     return _eval_floor_ceil(self, x, "floor", **kwds)
File src/cysignals/signals.pyx:310, in cysignals.signals.python_check_interrupt()
KeyboardInterrupt: 

O teorema Isotérmico

A temperatura de um vértice é a soma de todos os pesos chegando a ele. O teorema isotérmico diz que:

Se todos os vértices têm a mesma temperatura,então a probabilidade de fixação é equivalente à do processo de Moran.

Como nestes grafos evolutivos, as arestas representam probabilidade de que a prole de $i$ substitua a de $j$, para o grafo ser isotérmico, a matriz deve ser duplamente estocástica.

Exercício 4:

Prove que para um grafo populacional ser $\rho$-equivalente ao processo de Moran, ele precisa ser isotérmico, ou seja, sua matriz $W=[w_{ij}]$ é duplamente estocástica.

Amplificando e suprimindo a Seleção

Um vértice raiz é um vértice sem nenhuma aresta levando a ele. Uma raiz tem temperatura zero. Se um grafo tem apenas uma raiz, sua probabilidade de fixação é de $1/N$. A única forma de um mutante dominar toda a população é se ele for introduzido no vértice raiz.

Se um grafo possui multiplas raízes, então nenhum mutante será capaz de se fixar. mutantes que surjam em uma das raízes darão origem a linhagens que jamais serão extintas. Logo grafos com multiplas raízes permitem a coexistência de múltiplas linhagens.

Construindo um supressor de seleção

Uma maneira simples de construir um grasfo supressor de seleção é dividir a populações em duas subpopulações de tamanho $n_1$ e $n_2$ tal que $n_1+n_2=N$. A primeira subpopulação forma um grafo completo e a segunda formará um grafo com a única restrição de que todos os nós desta segunda subpopulação possam ser alcançados a partir da primeira. Assim a probabilidade de fixação fica sendo:

$$\rho_G=\frac{1-1/r}{1-1/r^{n_1}}$$

Logo, neste caso $1/N < \rho_G < \rho_M(N)$

n1=5
n2 =6
g = digraphs.Complete(n1)
h = digraphs.Path(n2)
G = g+h
G.add_edge([n1+1,n1+3])
G.add_edge([1,n1])
G

Estimando o $\rho_G$ de um grafo qualquer:

Aqui vamos usar o método de monte-carlo para estimar $\rho_G$, ou seja, vamos realizar a dinâmica evolutiva um grande número de vezes em um grafo e obter assim a probabilidade de fixação.

def passo_evolutivo(Gr, As, r=1):
    """
    Passo evolutivo.
    Gr: grafo
    As: Vertices do tipo A
    r: fitness relativo rb/ra
    """
    a = len(As)
    Bs = [i for i in Gr.vertices() if i not in As]
    b = len(Bs)
    # Seleciona quem vai se reproduzir e descobre a vizinhanca
    if np.random.random() <= 1/(a+r*b):
        nbs = Gr.neighbors_out(As[randrange(a)])
        reprodutor = 'A'
    else:
        nbs = Gr.neighbors_out(Bs[randrange(b)])
        reprodutor = 'B'
    n = len(nbs)
    if n != 0: # Caso existam vizinhos
        m = nbs[randrange(n)] # m recebera a prole
        if reprodutor == 'A':
            if m in As:
                pass
            else:
                As.append(m)
        elif reprodutor == 'B':
            if m in As:
                As.remove(m)
            else:
                pass
    return As
    
def estima_rho(Gr, r=1, n_iter_1=1000, n_iter_2=1000):
    l = len(Gr.vertices())
    s = len(Gr.sources())
    if s == 1: # se o Grafo possui uma unica raiz
        rho = 1./l
    elif s > 1: #Grafo possui mais de uma raiz
        rho = 0
    else: # grafo não tem raízes
        rhos = 0
        for i in range(n_iter_1):
            As = list(range(l)) # Inicialmente todos são do tipo A
            As.pop(randrange(l)) # um mutante entra em uma posição aleatoria
            j = 0
            while len(As)>0 and len(As)<l and j<n_iter_2:
                As = passo_evolutivo(Gr,As,r)
                j += 1
            if len(As) == 0: # Ocorreu a fixação em até n_iter_2 passos
                rhos += 1
        rho = float(rhos/n_iter_1)
    return rho

O $\rho_G$ esperado é, neste caso:

$\frac{1-1/r}{1-1/r^{n1}}$

var('rho_M')
r=10.7
rho_G = (1-1/r)/(1-1/r^n1)
show(rho_M,": ", (1-1/r)/(1-1/r^(n1+n2)))
show("Esperado(rho_G): ", rho_G)
print("Estimado: ")
estima_rho(G,r,100,700000)

Estimado: 

rho_G*(n1/(n1+n2))