<jupyter_start><jupyter_text>Modelos MetapopulacionaisO conceito de metapopulações pode ser visto como uma maneira relativamente simples de tratar a heterogeneidade em populações sem abandonar o paradigma de populações homogêneas modeladas com Equações diferenciais ordinárias.![metapops](https://upload.wikimedia.org/wikipedia/commons/thumb/9/94/Metapopulation_%281%29.svg/2560px-Metapopulation_%281%29.svg.png) Exercício 1:Desenvolva um modelo metapopulacional com dois patches (duas sub-populações) A e B nos quais se propaga uma doença de transmissão direta modelada como SIR. Considere que existe um fluxo migratório entre as cidades.1. Escreva as equações das duas populações representando as migrações $A \rightarrow B$ como  $k_1$, e $B \rightarrow A$ como $k_2$.1. Utilizando uma parametrização que não permita equilíbrio endêmico, encontre os equilíbrios $S_A^* + S_B^*$ e $I_A^* + I_B^*$ .1. Calcule a Matriz Jacobiana do sistema e determine  a estabilidade dos equilíbrios.1. Crie uma função interativa que nos permita alterar os valores de $k_1$ e $k_2$. Realize simulações com $k_1 = k_2$, $k_1  k_2$, com $\beta_A \neq \beta_B$<jupyter_code>%display typeset<jupyter_output><empty_output><jupyter_text>Resposta ao item 1:$$ \dot{S}_A = -\beta_A S_A I_A + k_2 S_B - k_1 S_A$$$$\dot{I}_A = \beta_A S_A I_A + k_2 I_B - k_1 I_A -\gamma I_A$$$$\dot{S}_B = -\beta_B S_B I_B + k_1 S_A - k_2 S_B $$$$\dot{I}_B = \beta_B S_B I_B + k_1 I_A - k_2 I_B -\gamma I_B$$ 2) Vamos agora encontrar os equilíbrios do sistema usando o solver do sage:<jupyter_code>var('S_A I_A S_B I_B beta_A beta_B R_A R_B gamma k_1 k_2')
#beta_A = 1.2
#beta_B = 1.2
#gamma = 0.5
#k_1 = 0.5
#k_2 = 0.5
solve([-beta_A*S_A*I_A + k_2*S_B - k_1*S_A, 
       beta_A*S_A*I_A+ k_2*I_B - k_1*I_A - gamma*I_A,
       -beta_B*S_B*I_B - k_2*S_B + k_1*S_A,
       beta_B*S_B*I_B- k_2*I_B + k_1*I_A - gamma*I_B,
       gamma*I_A,
       gamma*I_B],[S_A, I_A, S_B, I_B, R_A, R_B])<jupyter_output><empty_output><jupyter_text>3) Calculamos a Matriz jacobiana para podermos realizar a análise de estabilidade do equilíbrio encontrado acima<jupyter_code>J = jacobian([-beta_A*S_A*I_A + k_2*S_B - k_1*S_A, 
       beta_A*S_A*I_A+ k_2*I_B - k_1*I_A - gamma*I_A,
       -beta_B*S_B*I_B - k_2*S_B + k_1*S_A,
       beta_B*S_B*I_B- k_2*I_B + k_1*I_A - gamma*I_B],[S_A, I_A, S_B, I_B])
J<jupyter_output><empty_output><jupyter_text>Agora Substituímos o valor do equilíbrio na matriz jacobiana:<jupyter_code>J_eq = J(S_A=1,I_A=0,S_B=k_1/k_2,I_B=0)
J_eq<jupyter_output><empty_output><jupyter_text>e também os valores dos parâmetros<jupyter_code>J_eq(gamma=1.0,k_1=.5,k_2=.5,beta_A=1.6,beta_B=1.6).eigenvalues()<jupyter_output><empty_output><jupyter_text>Como o maior autovalor é real e positivo o Equilíbrio é instável 4) Agora vamos criar uma simulação interativa para estudar o efeito da existência de patches na dinâmica do sistema.<jupyter_code>def model_2_pop(t, y, params):
    beta_A, beta_B, gamma, k_1, k_2, mu = params
    S_A, I_A, S_B, I_B, R_A, R_B = y
    return [-beta_A*S_A*I_A + k_2*S_B - k_1*S_A, 
       beta_A*S_A*I_A+ k_2*I_B - k_1*I_A - gamma*I_A,
       -beta_B*S_B*I_B - k_2*S_B + k_1*S_A,
       beta_B*S_B*I_B - k_2*I_B + k_1*I_A - gamma*I_B,
       gamma*I_A,
       gamma*I_B]<jupyter_output><empty_output><jupyter_text>Agora vamos criar uma versão do modelo anterior, com nascimentos e mortes, ou seja demografia:$$ \dot{S}_A = -\beta_A S_A I_A + k_2 S_B - k_1 S_A -\mu S_A +\mu(S_A+I_A+R_A)$$$$\dot{I}_A = \beta_A S_A I_A + k_2 I_B - k_1 I_A -\gamma I_A -\mu I_A$$$$\dot{R}_A = \gamma I_A-\mu R_A$$$$\dot{S}_B = -\beta_B S_B I_B + k_1 S_A - k_2 S_B -\mu S_B +\mu(S_B+I_B+R_B)$$$$\dot{I}_B = \beta_B S_B I_B + k_1 I_A - k_2 I_B -\gamma I_B -\mu I_B$$$$\dot{R}_B= \gamma I_B-\mu R_B$$<jupyter_code>def model_2_pop_demog(t, y, params):
    beta_A, beta_B, gamma, k_1, k_2, mu = params
    S_A, I_A, S_B, I_B, R_A, R_B = y
    return [-beta_A*S_A*I_A + k_2*S_B - k_1*S_A -mu*S_A +mu*(S_A+I_A+R_A), 
       beta_A*S_A*I_A+ k_2*I_B - k_1*I_A - gamma*I_A - mu*I_A,
       -beta_B*S_B*I_B - k_2*S_B + k_1*S_A - mu*S_B+mu*(S_B+I_B+R_B),
       beta_B*S_B*I_B - k_2*I_B + k_1*I_A - gamma*I_B - mu*I_B,
       gamma*I_A-mu*R_A,
       gamma*I_B-mu*R_B]<jupyter_output><empty_output><jupyter_text>Agora passamos à instanciação do nosso solver e à implementação de uma função de plotagem.<jupyter_code>T=ode_solver()
T.function = model_2_pop
T.algorithm = "rk8pd"#"rk4imp"
from itertools import cycle
labels = [r"$S_A$", "$I_A$", "$S_B$", "$I_B$","$R_A$", "$R_B$"]

def plot_sol(sol):
    plots=None
    c = cycle(['red','blue','green', 'black', 'yellow', 'orange', 'magenta', 'gray', 'pink', 'brown'])
    for i,co in zip(range(len(sol[0][1])),c):
#         co = c.next()
        if plots is None:
            plots = list_plot([(j[0],j[1][i]) for j in sol],color=co, plotjoined=True,      legend_label='%s'%labels[i], alpha=.8, gridlines=true)
        else:
            plots += list_plot([(j[0],j[1][i]) for j in sol],color=co, plotjoined=True,           legend_label='%s'%labels[i], alpha=.8, gridlines=true)
    show(plots)
@interact()
def simula(beta_A=(4.6,(.2,5,.1)), beta_B=(3.6,(.2,5,.1)), gamma=(1.2, 2,.1), k_1=(0.5,(0,1,.1)), k_2=(0.5,(0,1,.1)), mu=(.04,(0, 1,.1)), model=selector(['no demog', 'demog'], buttons=True)):
    if model == "demog":
        T.function = model_2_pop_demog
    else:
        T.function = model_2_pop
    inits = (.9,0.1, 1, 0, 0, 0)
    t_span = [0,30]
    T.ode_solve(t_span,inits,num_points=300, params=(beta_A, beta_B, gamma, k_1, k_2, mu))
    show(inits[0]*k_1/k_2)
    plot_sol(T.solution)<jupyter_output><empty_output><jupyter_text>Exercício 2:Modelos evolutivos em 2 patches.  Adapte o modelo de seleção natural sem mutação, para duas subpopulações, mantendo a generalidade para multiplos tipos.1. Escreva as equações em forma matricial2. Encontre os equilíbrios do sistema3. Escreva uma função interativa para estudar os parâmetros<jupyter_code>pretty_print(html('a) Vamos escrever o modelo evolutivo, começando com 4 tipos e 5 localidades, denotadas pelos vetores $X$ e $L$, respectivamente:'))
var('x_1 x_2 x_3 x_4 A B C D E')
X = Matrix(SR,[x_1, x_2, x_3, x_4])
L = Matrix(SR,[A,B,C,D,E])
show("X=", X,", L=",L)

D = Matrix(SR, X.ncols(),L.ncols(), lambda i,j: var('{}{}'.format(X[0,i],L[0,j])))
pretty_print(html("Seja $D$ a matriz com as densidades em cada uma das localidades:"))
show("D=",D)
pretty_print(html("Seja $M$ a matriz com as taxas de migração $k_{ij}$ entre cada uma das localidades:"))
M = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: var('k_{}{}'.format(L[0,i],L[0,j])))
MD = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: 0+ int(i==j)*sum(M[i]))
MLT = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: (1-int(j>=i))*var('k_{}{}'.format(L[0,i],L[0,j])))
MUT = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: (1-int(j<i))*var('k_{}{}'.format(L[0,i],L[0,j])))
M = M-(identity_matrix(5).elementwise_product(M))+MD  #MLT+MD+MUT
show("M=",M)
pretty_print(html(r"Logo, $F=D \times M$ é uma matriz em que cada elemento $f_{ij}$ representa a imigração do tipo $i$ em cada patch $j$ vindo de todos os outros patches:"))
F = D*M
show("F=", F)
pretty_print(html("Seja $f_i$ o vetor de fitness de cada tipo $i$:"))
f_i = Matrix(SR, [1/5, 1/3, 2/5, 1/2])
show('f=', f_i)
pretty_print(html(r"Nosso modelo fica então $\dot{x_{ij}}=f_i \times x_{ij} + F_{ij} - x_{ij}\sum_j M_{ij} -\phi x_i$:"))
def sel_model(t, D, params):
    D = D.reshape((4,5))
    f, M= params
    F=D*M
    E = M.sum(axis=1) # emigration
    X = D.sum(axis=1)
    phi = sum(np.dot(f,X))
    res = np.dot(f,X) + F- np.dot(X,E) - phi*D
    return res
Matrix(SR, L.ncols(),L.ncols(), lambda i,j: var('k_{}{}'.format(L[0,i],L[0,j])))*identity_matrix(5)
sum(M[i] for i in range(M.nrows()))
T2=ode_solver()
T2.function = sel_model
T2.algorithm = "rk8pd"
import numpy as np
inits = np.zeros(20)+0.25
t_span = [0,70]
T.ode_solve(t_span,list(inits),num_points=300, params=(f_i,D))
sel_model(0,inits,[f_i,np.array(D)])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17173858817660265, 'fraction_numerical': 0.029558493389872785, 'mean_word_length': 4.511340206185567, 'pattern_counts': {'":': 1, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '483420', 'license': 'BSD-Source-Code', 'n_tokens_mistral': 3881, 'n_tokens_neox': 3777, 'n_words': 810, 'path': 'Planilhas Sage/Aula 15 - Metapopulações.ipynb', 'repo_name': 'fccoelho/Modelagem-Matematica-IV', 'size': '8018'}