Criar modelo para simular isolamento parcial da população (ex idosa) - "quarentena vertical"

[vitoventurieri]

Olá pensei em criar um modelo para simular o distanciamento social dos suscetíveis, mas apenas de parte da população, como por exemplo, para simular o distanciamento social das pessoas com mais de 55 anos

Pensando nas ODEs seria algo assim acho (me corrijam se estiver errado):

#V = suscetíveis idosos; J= suscetíveis jovens; ômega= (fator atenuador de contatos para o grupo idoso, que será menor ou igual a 1)
        dVdt = -beta * ômega * V * I / N
        dJdt = -beta * J * I / N
        dEdt = - dJdt -dVdt  -  alpha*E
        dIdt = alpha*E - gamma*I
        dRdt = gamma * I

Info para as condiçoes iniciais para o Brasil, segundo a piramide etária do IBGE, temos:

Pessoas com 55 anos ou mais	28.866.818
Pessoas com até 55 anos	155.109.809

Ou seja, teríamos as seguintes frações de distribuição da populacao no inicio em cada compartimento:

Pessoas com 55 anos ou mais	0,1569048
Pessoas com até 55 anos	0,8430952

[giulianonetto]

Galera, comecei um draft dessa implementação:

Parâmetros:

Sg0 = 155_109_809
Sr0 = 28_866_818
E0 = 10
I0 = 10
R0 = 0
N = Sg0 + Sr0 + E0 + I0 + R0

_params = {
    't_max': 240, 
    'total_population': N,
    'r0_dist': make_lognormal_from_interval(3.5, 4, 0.95),
    'gamma_inv_dist': make_lognormal_from_interval(10, 12, 0.95),
    'alpha_inv_dist': make_lognormal_from_interval(5, 7, 0.95),
    'omega': 0.1,
    'init_conditions': (Sg0, Sr0, E0, I0, R0)
}

Setup inicial e sampling

size = 500
t_space = np.arange(0, _params['t_max'])
N = _params['total_population']
Sg, Sr, E, I, R = [np.zeros((_params['t_max'], size)) for _ in range(5)]
Sg[0, ], Sr[0, ], E[0, ], I[0, ], R[0, ] = _params['init_conditions']

r0 = _params['r0_dist'].rvs(size)
gamma = 1/_params['gamma_inv_dist'].rvs(size)
alpha = 1/_params['alpha_inv_dist'].rvs(size)
beta_general = r0*gamma
beta_risk = beta_general * _params['omega']

SEIR

def transition_probability(_lambda):
    return expon(scale=_lambda).cdf(1)

for t in t_space[1:]:
    SgE = npr.binomial(Sg[t-1, ].astype(int),
                      transition_probability(_lambda=1/(beta_general*I[t-1, ]/N)))    
    SrE = npr.binomial(Sr[t-1, ].astype(int),
                      transition_probability(_lambda=1/(beta_risk*I[t-1, ]/N)))
    EI = npr.binomial(E[t-1, ].astype(int),
                      transition_probability(_lambda=1/alpha))
    IR = npr.binomial(I[t-1, ].astype(int),
                      transition_probability(_lambda=1/gamma))

    dSg =  0 - SgE
    dSr = 0 - SrE
    dE = SgE + SrE - EI
    dI = EI - IR
    dR = IR - 0

    Sg[t, ] = Sg[t-1, ] + dSg
    Sr[t, ] = Sr[t-1, ] + dSr
    E[t, ] = E[t-1, ] + dE
    I[t, ] = I[t-1, ] + dI
    R[t, ] = R[t-1, ] + dR

Resultado parcial

azul = I amarelo = E vermelho = S geral verde = S grupo de risco

no caso setei uns parametros com pouca incerteza e com um cenário meio "bad", mas pra mim tá demorando muito o pico das curvas I e E.

Sugestões são bem-vindas. Depois podemos separar os recuperados em "curados" e mortes também.

[vitoventurieri]

Tenta colocar nas condições iniciais o numero de casos que temos + ou - agora, deve demorar bem menos o pico. Se não for isso, tenho que olhar com mais calma