Aula 15 - Metapopulações.ipynb

Kernel: SageMath 9.7

Modelos Metapopulacionais

O 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.

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$.

2. 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^*$ .

3. Calcule a Matriz Jacobiana do sistema e determine a estabilidade dos equilíbrios.

4. 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$ e $k_1 > k_2$, com $\beta_A \neq \beta_B$

%display typeset

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:

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])

$\displaystyle \left[\left[S_{A} = r_{1}, I_{A} = 0, S_{B} = \frac{k_{1} r_{1}}{k_{2}}, I_{B} = 0, R_{A} = r_{2}, R_{B} = r_{3}\right]\right]$

3. Calculamos a Matriz jacobiana para podermos realizar a análise de estabilidade do equilíbrio encontrado acima

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

$\displaystyle \left(\begin{array}{rrrr} -I_{A} \beta_{A} - k_{1} & -S_{A} \beta_{A} & k_{2} & 0 \\ I_{A} \beta_{A} & S_{A} \beta_{A} - \gamma - k_{1} & 0 & k_{2} \\ k_{1} & 0 & -I_{B} \beta_{B} - k_{2} & -S_{B} \beta_{B} \\ 0 & k_{1} & I_{B} \beta_{B} & S_{B} \beta_{B} - \gamma - k_{2} \end{array}\right)$

Agora Substituímos o valor do equilíbrio na matriz jacobiana:

J_eq = J(S_A=1,I_A=0,S_B=k_1/k_2,I_B=0)
J_eq

$\displaystyle \left(\begin{array}{rrrr} -k_{1} & -\beta_{A} & k_{2} & 0 \\ 0 & \beta_{A} - \gamma - k_{1} & 0 & k_{2} \\ k_{1} & 0 & -k_{2} & -\frac{\beta_{B} k_{1}}{k_{2}} \\ 0 & k_{1} & 0 & -\gamma + \frac{\beta_{B} k_{1}}{k_{2}} - k_{2} \end{array}\right)$

e também os valores dos parâmetros

J_eq(gamma=1.0,k_1=.5,k_2=.5,beta_A=1.6,beta_B=1.6).eigenvalues()

$\displaystyle \left[\frac{3}{5}, -\frac{2}{5}, -1, 0\right]$

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.

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]

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$$

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]

Agora passamos à instanciação do nosso solver e à implementação de uma função de plotagem.

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)

Exercício 2:

Modelos evolutivos em 5 patches. Adapte o modelo de seleção natural sem mutação, para 4 subpopulações, mantendo a generalidade para multiplos tipos.

1. Escreva as equações em forma matricial

2. Encontre os equilíbrios do sistema

3. Escreva uma função interativa para estudar os parâmetros

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. como não faz sentido fala de migração de uma localidade "))
pretty_print(html("para si mesma, sua diagonal é zero. "))
M = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: var(f'k_{L[0,i]}{L[0,j]}'))
MD = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: 0 + int(i==j)*sum(M[i]))  # Diagonal de M
MLT = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: (1-int(j>=i))*var(f'k_{L[0,i]}{L[0,j]}')) # Triangulo inferior
MUT = Matrix(SR, L.ncols(),L.ncols(), lambda i,j: (1-int(j<=i))*var(f'k_{L[0,i]}{L[0,j]}')) # Triangulo superior
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:"))
Fi = D*M
show("Fi=", Fi)

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}- E_{ij} -\phi x_i$:"))

a) Vamos escrever o modelo evolutivo, começando com 4 tipos e 5 localidades, denotadas pelos vetores $X$ e $L$, respectivamente:

$\displaystyle \verb|X=| \left(\begin{array}{rrrr} x_{1} & x_{2} & x_{3} & x_{4} \end{array}\right) \verb|,|\verb| |\verb|L=| \left(\begin{array}{rrrrr} A & B & C & D & E \end{array}\right)$

Seja $D$ a matriz com as densidades em cada uma das localidades:

$\displaystyle \verb|D=| \left(\begin{array}{rrrrr} x_{\mathit{1A}} & x_{\mathit{1B}} & x_{\mathit{1C}} & x_{\mathit{1D}} & x_{\mathit{1E}} \\ x_{\mathit{2A}} & x_{\mathit{2B}} & x_{\mathit{2C}} & x_{\mathit{2D}} & x_{\mathit{2E}} \\ x_{\mathit{3A}} & x_{\mathit{3B}} & x_{\mathit{3C}} & x_{\mathit{3D}} & x_{\mathit{3E}} \\ x_{\mathit{4A}} & x_{\mathit{4B}} & x_{\mathit{4C}} & x_{\mathit{4D}} & x_{\mathit{4E}} \end{array}\right)$

Seja $M$ a matriz com as taxas de migração $k_{ij}$ entre cada uma das localidades. como não faz sentido fala de migração de uma localidade

para si mesma, sua diagonal é zero.

$\displaystyle \verb|M=| \left(\begin{array}{rrrrr} 0 & k_{\mathit{AB}} & k_{\mathit{AC}} & k_{\mathit{AD}} & k_{\mathit{AE}} \\ k_{\mathit{BA}} & 0 & k_{\mathit{BC}} & k_{\mathit{BD}} & k_{\mathit{BE}} \\ k_{\mathit{CA}} & k_{\mathit{CB}} & 0 & k_{\mathit{CD}} & k_{\mathit{CE}} \\ k_{\mathit{DA}} & k_{\mathit{DB}} & k_{\mathit{DC}} & 0 & k_{\mathit{DE}} \\ k_{\mathit{EA}} & k_{\mathit{EB}} & k_{\mathit{EC}} & k_{\mathit{ED}} & 0 \end{array}\right)$

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:

$\displaystyle \verb|Fi=| \left(\begin{array}{rrrrr} k_{\mathit{BA}} x_{\mathit{1B}} + k_{\mathit{CA}} x_{\mathit{1C}} + k_{\mathit{DA}} x_{\mathit{1D}} + k_{\mathit{EA}} x_{\mathit{1E}} & k_{\mathit{AB}} x_{\mathit{1A}} + k_{\mathit{CB}} x_{\mathit{1C}} + k_{\mathit{DB}} x_{\mathit{1D}} + k_{\mathit{EB}} x_{\mathit{1E}} & k_{\mathit{AC}} x_{\mathit{1A}} + k_{\mathit{BC}} x_{\mathit{1B}} + k_{\mathit{DC}} x_{\mathit{1D}} + k_{\mathit{EC}} x_{\mathit{1E}} & k_{\mathit{AD}} x_{\mathit{1A}} + k_{\mathit{BD}} x_{\mathit{1B}} + k_{\mathit{CD}} x_{\mathit{1C}} + k_{\mathit{ED}} x_{\mathit{1E}} & k_{\mathit{AE}} x_{\mathit{1A}} + k_{\mathit{BE}} x_{\mathit{1B}} + k_{\mathit{CE}} x_{\mathit{1C}} + k_{\mathit{DE}} x_{\mathit{1D}} \\ k_{\mathit{BA}} x_{\mathit{2B}} + k_{\mathit{CA}} x_{\mathit{2C}} + k_{\mathit{DA}} x_{\mathit{2D}} + k_{\mathit{EA}} x_{\mathit{2E}} & k_{\mathit{AB}} x_{\mathit{2A}} + k_{\mathit{CB}} x_{\mathit{2C}} + k_{\mathit{DB}} x_{\mathit{2D}} + k_{\mathit{EB}} x_{\mathit{2E}} & k_{\mathit{AC}} x_{\mathit{2A}} + k_{\mathit{BC}} x_{\mathit{2B}} + k_{\mathit{DC}} x_{\mathit{2D}} + k_{\mathit{EC}} x_{\mathit{2E}} & k_{\mathit{AD}} x_{\mathit{2A}} + k_{\mathit{BD}} x_{\mathit{2B}} + k_{\mathit{CD}} x_{\mathit{2C}} + k_{\mathit{ED}} x_{\mathit{2E}} & k_{\mathit{AE}} x_{\mathit{2A}} + k_{\mathit{BE}} x_{\mathit{2B}} + k_{\mathit{CE}} x_{\mathit{2C}} + k_{\mathit{DE}} x_{\mathit{2D}} \\ k_{\mathit{BA}} x_{\mathit{3B}} + k_{\mathit{CA}} x_{\mathit{3C}} + k_{\mathit{DA}} x_{\mathit{3D}} + k_{\mathit{EA}} x_{\mathit{3E}} & k_{\mathit{AB}} x_{\mathit{3A}} + k_{\mathit{CB}} x_{\mathit{3C}} + k_{\mathit{DB}} x_{\mathit{3D}} + k_{\mathit{EB}} x_{\mathit{3E}} & k_{\mathit{AC}} x_{\mathit{3A}} + k_{\mathit{BC}} x_{\mathit{3B}} + k_{\mathit{DC}} x_{\mathit{3D}} + k_{\mathit{EC}} x_{\mathit{3E}} & k_{\mathit{AD}} x_{\mathit{3A}} + k_{\mathit{BD}} x_{\mathit{3B}} + k_{\mathit{CD}} x_{\mathit{3C}} + k_{\mathit{ED}} x_{\mathit{3E}} & k_{\mathit{AE}} x_{\mathit{3A}} + k_{\mathit{BE}} x_{\mathit{3B}} + k_{\mathit{CE}} x_{\mathit{3C}} + k_{\mathit{DE}} x_{\mathit{3D}} \\ k_{\mathit{BA}} x_{\mathit{4B}} + k_{\mathit{CA}} x_{\mathit{4C}} + k_{\mathit{DA}} x_{\mathit{4D}} + k_{\mathit{EA}} x_{\mathit{4E}} & k_{\mathit{AB}} x_{\mathit{4A}} + k_{\mathit{CB}} x_{\mathit{4C}} + k_{\mathit{DB}} x_{\mathit{4D}} + k_{\mathit{EB}} x_{\mathit{4E}} & k_{\mathit{AC}} x_{\mathit{4A}} + k_{\mathit{BC}} x_{\mathit{4B}} + k_{\mathit{DC}} x_{\mathit{4D}} + k_{\mathit{EC}} x_{\mathit{4E}} & k_{\mathit{AD}} x_{\mathit{4A}} + k_{\mathit{BD}} x_{\mathit{4B}} + k_{\mathit{CD}} x_{\mathit{4C}} + k_{\mathit{ED}} x_{\mathit{4E}} & k_{\mathit{AE}} x_{\mathit{4A}} + k_{\mathit{BE}} x_{\mathit{4B}} + k_{\mathit{CE}} x_{\mathit{4C}} + k_{\mathit{DE}} x_{\mathit{4D}} \end{array}\right)$

Seja $f_i$ o vetor de fitness de cada tipo $i$:

$\displaystyle \verb|f=| \left(\begin{array}{rrrr} \frac{1}{5} & \frac{1}{3} & \frac{2}{5} & \frac{1}{2} \end{array}\right)$

Nosso modelo fica então $\dot{x_{ij}}=f_i \times x_{ij} + F_{ij}- E_{ij} -\phi x_i$:

Onde a emigração $E_{ij}$de cada tipo $i$ a partir de cada local $j$ é dada pelo efluxo total de cada localidade:

# Vamos somar as linhas de M para gerar os efluxo de cada localidade
of=matrix(sum(M.columns()))
of

$\displaystyle \left(\begin{array}{rrrrr} k_{\mathit{AB}} + k_{\mathit{AC}} + k_{\mathit{AD}} + k_{\mathit{AE}} & k_{\mathit{BA}} + k_{\mathit{BC}} + k_{\mathit{BD}} + k_{\mathit{BE}} & k_{\mathit{CA}} + k_{\mathit{CB}} + k_{\mathit{CD}} + k_{\mathit{CE}} & k_{\mathit{DA}} + k_{\mathit{DB}} + k_{\mathit{DC}} + k_{\mathit{DE}} & k_{\mathit{EA}} + k_{\mathit{EB}} + k_{\mathit{EC}} + k_{\mathit{ED}} \end{array}\right)$

Multiplicado pelas densidades relativas de cada tipo em cada localidade

for i in range(2):
    of=of.stack(of)
E = of.elementwise_product(D)
show("Ei=",E)

$\displaystyle \verb|Ei=| \left(\begin{array}{rrrrr} {\left(k_{\mathit{AB}} + k_{\mathit{AC}} + k_{\mathit{AD}} + k_{\mathit{AE}}\right)} x_{\mathit{1A}} & {\left(k_{\mathit{BA}} + k_{\mathit{BC}} + k_{\mathit{BD}} + k_{\mathit{BE}}\right)} x_{\mathit{1B}} & {\left(k_{\mathit{CA}} + k_{\mathit{CB}} + k_{\mathit{CD}} + k_{\mathit{CE}}\right)} x_{\mathit{1C}} & {\left(k_{\mathit{DA}} + k_{\mathit{DB}} + k_{\mathit{DC}} + k_{\mathit{DE}}\right)} x_{\mathit{1D}} & {\left(k_{\mathit{EA}} + k_{\mathit{EB}} + k_{\mathit{EC}} + k_{\mathit{ED}}\right)} x_{\mathit{1E}} \\ {\left(k_{\mathit{AB}} + k_{\mathit{AC}} + k_{\mathit{AD}} + k_{\mathit{AE}}\right)} x_{\mathit{2A}} & {\left(k_{\mathit{BA}} + k_{\mathit{BC}} + k_{\mathit{BD}} + k_{\mathit{BE}}\right)} x_{\mathit{2B}} & {\left(k_{\mathit{CA}} + k_{\mathit{CB}} + k_{\mathit{CD}} + k_{\mathit{CE}}\right)} x_{\mathit{2C}} & {\left(k_{\mathit{DA}} + k_{\mathit{DB}} + k_{\mathit{DC}} + k_{\mathit{DE}}\right)} x_{\mathit{2D}} & {\left(k_{\mathit{EA}} + k_{\mathit{EB}} + k_{\mathit{EC}} + k_{\mathit{ED}}\right)} x_{\mathit{2E}} \\ {\left(k_{\mathit{AB}} + k_{\mathit{AC}} + k_{\mathit{AD}} + k_{\mathit{AE}}\right)} x_{\mathit{3A}} & {\left(k_{\mathit{BA}} + k_{\mathit{BC}} + k_{\mathit{BD}} + k_{\mathit{BE}}\right)} x_{\mathit{3B}} & {\left(k_{\mathit{CA}} + k_{\mathit{CB}} + k_{\mathit{CD}} + k_{\mathit{CE}}\right)} x_{\mathit{3C}} & {\left(k_{\mathit{DA}} + k_{\mathit{DB}} + k_{\mathit{DC}} + k_{\mathit{DE}}\right)} x_{\mathit{3D}} & {\left(k_{\mathit{EA}} + k_{\mathit{EB}} + k_{\mathit{EC}} + k_{\mathit{ED}}\right)} x_{\mathit{3E}} \\ {\left(k_{\mathit{AB}} + k_{\mathit{AC}} + k_{\mathit{AD}} + k_{\mathit{AE}}\right)} x_{\mathit{4A}} & {\left(k_{\mathit{BA}} + k_{\mathit{BC}} + k_{\mathit{BD}} + k_{\mathit{BE}}\right)} x_{\mathit{4B}} & {\left(k_{\mathit{CA}} + k_{\mathit{CB}} + k_{\mathit{CD}} + k_{\mathit{CE}}\right)} x_{\mathit{4C}} & {\left(k_{\mathit{DA}} + k_{\mathit{DB}} + k_{\mathit{DC}} + k_{\mathit{DE}}\right)} x_{\mathit{4D}} & {\left(k_{\mathit{EA}} + k_{\mathit{EB}} + k_{\mathit{EC}} + k_{\mathit{ED}}\right)} x_{\mathit{4E}} \end{array}\right)$

def sel_model(t, D, params):
    D = np.array(D).reshape((4,5))
#     print(D)
    f, M= np.array(params[0]).reshape(4,1),np.array(params[1])
    F=D@M
    of = (M.sum(axis=1)).reshape(1,5)
    E = np.vstack([of,of,of,of])*D

    phi = sum(f*D, axis=0)
    res = f*D + F - E - phi*D 
    return list(res.ravel())

T2=ode_solver()
T2.function = sel_model
T2.algorithm = "rk8pd"

import numpy as np
import matplotlib.pyplot as plt
inits = np.random.random((4,5))
# Normalizing population on every locality
s = inits.sum(axis=0)
inits /= s.reshape(1,5)
print(inits.sum(axis=0))
inits = inits.ravel()
t_span = [0,70]
m = np.random.random((5,5)) # random migration rates
np.fill_diagonal(m,0)
m = np.around(m,2)
f = np.array(f_i)
T2.ode_solve(t_span,list(inits),num_points=300, params=(f,m))

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

import networkx as nx
def plot_network(mr):
    G = nx.from_numpy_matrix(mr, parallel_edges=True)
    G = nx.relabel_nodes(G, {0:'A',1:'B',2:'C',3:'D',4:'E'})

    elarge = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] > 0.5]
    esmall = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] <= 0.5]

    pos = nx.spring_layout(G)
    nx.draw_networkx_nodes(G, pos, node_size=700)
    nx.draw_networkx_edges(G, pos, edgelist=elarge, width=6)
    nx.draw_networkx_edges(
            G, pos, edgelist=esmall, width=6, alpha=0.5, edge_color="b", style="dashed"
        )

    # node labels
    nx.draw_networkx_labels(G, pos, font_size=20, font_family="sans-serif")
    # edge weight labels
    edge_labels = nx.get_edge_attributes(G, "weight")
    nx.draw_networkx_edge_labels(G, pos, edge_labels)

    ax = plt.gca()
    ax.margins(0.08)
    plt.axis("off")
    plt.tight_layout()
plot_network(m)

def plot_matrix_pop(sol):
    series = []
    ts = []
    for t, s in sol:
        ts.append(t)
        series.append(s)
    plt.plot(ts, np.array(series))
    plt.grid()
    plt.xlabel('time')
    plt.ylabel('Abundance')
    plt.legend(['$x_{1A}$','$x_{1B}$','$x_{1C}$','$x_{1D}$','$x_{1E}$',
                '$x_{2A}$','$x_{2B}$','$x_{2C}$','$x_{2D}$','$x_{2E}$',
                '$x_{3A}$','$x_{3B}$','$x_{3C}$','$x_{3D}$','$x_{3E}$',
                '$x_{4A}$','$x_{4B}$','$x_{4C}$','$x_{4D}$','$x_{4E}$'
               ], fancybox=True, loc='lower center', bbox_to_anchor=(0.5, 1.05),
          ncol=4,)

plot_matrix_pop(T2.solution)

inits = sel_model(0,inits,[f_i,m])
# inits

Epigrass: Modelos metapopulacionais geo-referenciados

Para criar e simular modelos metapopulacionais mais complexos a biblioteca Epigrass pode ajudar.