Le modèle SEIRD de l'Italie

La modélisation du covid-19 pour l'Italie

In [208]:
from IPython.display import HTML
HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
Le code python est caché par défault.
Pour l'activer ou le désactiver cliquez <a href="javascript:code_toggle()">ici</a>.''')
Out[208]:
Le code python est caché par défault. Pour l'activer ou le désactiver cliquez ici.
Le début du confinement en Italie était 9 mars 2020 et le début du déconfinement le 4 mai 2020. Il y a donc plusieurs $\beta$ en fonction de la propagation du virus, en autre la qualité du confinement, la prise de conscience de l'épidémie.
A partir du modèle SEIRD, nous avons cherché des paramètres qui permettent de se rapprocher au mieux du graphe ci-dessus. Les paramètres, pour une population totale de 60 431 283, nous semblant adaptés sont :
  • le nombre initial de personnes infectés : 14
  • la durée d'incubation : 3 jours
  • la durée de rétablissement : 17 jours
  • le taux de mortalité : 3%
Le taux de mortalité élevé en Italie peut s'expliquer par le fait que le popultion italienne soit plus vieille que la moyenne mondiale. En effet, 21,69% de la population a plus de 65 ans. L'Italie est le second pays le plus vieux au monde derrière le Japon.
Les graphes commencent le 31 janvier 2020, la date correspondant aux premiers cas resensés.
In [309]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import least_squares
import time
import datetime as dt

                    # Fonctions (Ne pas modifier)

def DerivSEIRD(y0, t, N, beta, sigma, gamma, mu):
    S, E, I, R, D = y0
    dSdt = -(beta/N)*I*S
    dEdt = (beta/N)*I*S - sigma*E
    dIdt = sigma*E - gamma*I - mu*I
    dRdt = gamma*I
    dDdt = mu*I
    return dSdt, dEdt, dIdt, dRdt, dDdt
    
    
def SEIRD(y0, n, N, beta, sigma, gamma, mu): 
    t = np.linspace(1, n, n)
    sol = odeint(DerivSEIRD, y0, t, args=(N, beta, sigma, gamma, mu))
    return sol.T


def SEIRD_phases(y0, Tn, N, Tbeta, sigma, gamma, mu):
    n = 0
    y = y0
    
    n0, beta0 = Tn[0], Tbeta[0]
    S, E, I, R, D = SEIRD(y, n0, N, beta0, sigma, gamma, mu)
    n = n+n0
    y = S[n-1], E[n-1], I[n-1], R[n-1], D[n-1]
    t = np.linspace(1, n, n)

    for i in range (1, len(Tn)):
        ni, betai = Tn[i], Tbeta[i]
        Si, Ei, Ii, Ri, Di = SEIRD(y, ni+1, N, betai, sigma, gamma, mu)[:,1:]
        S, E, I, R, D = np.concatenate((S,Si)), np.concatenate((E,Ei)), np.concatenate((I,Ii)), np.concatenate((R,Ri)), np.concatenate((D,Di))
        n = n+ni
        y = S[n-1], E[n-1], I[n-1], R[n-1], D[n-1]

    return S, E, I, R, D, t



def MoyMob(T, k):                            #Entier k : Nombre de valeurs prises
    l = len(T)                               #de chaque côté des éléments de T         
    L = np.zeros(l)
    
    if ((2*k+1)>l):
        return L
    
    else :
        for i in range(l):
            if ((i-k)<0):                           #Soit k < i
                m = i
                L[i] = np.mean( T[(i-m):(i+m+1)] )  #i-m = 0
            
            elif ((i+k)>(l-1)):                     #Soit k > l-1-i
                m = l-1-i
                L[i] = np.mean( T[(i-m):(i+m+1)] )  #i+m+1 = l
            
            else :
                L[i] = np.mean( T[(i-k):(i+k+1)] )
    
        return L



def Periode(n, p):
    q, r = n//p, n%p
    Tnp = [p]*q
    if (r!=0):
        Tnp.append(r)
    return Tnp




                    # Données : Epidémie



# dincub = Durée d'incubation : Temps avant d'être contaminant une fois contaminé (jours)
dincub = 3

# dret = Durée de rétablissement : Temps que l'on reste contaminant (jours)
dret = 17

# m = Taux de mortalité parmis les infectés en %
m = 3


# Pas besoin de modifier celles-ci car elles dépendent des précédentes

# sigma = Inverse de la durée d'incubation (jours^(-1)) : Taux d'invidus qui passe de E à I par jour
sigma = 1/dincub

# gamma = Inverse de la durée de rétablissement (jour^(-1)) : Taux d'individus qui passe de I à R par jour
gamma = 1/dret

# mu = Taux de mortalité journalier parmis les infectés (jours^(-1)) : Taux d'individus qui passe de I à D par jour
mu = (m/100)*gamma





                    # Données : Reelles (Modifier selon besoin)



# CodePays = Identifiant du pays que l'on veut étudier, sous la forme "CODE"
CodePays = "ITA"


# On définit un jour du début des observations.

# Pour commencer à une date précise, on la définit. Sinon on laisse 0, 0, 0.
# Exemple : pour le 13 Mars 2020, on remplace 0, 0, 0 par 13, 3, 2020
DateDeb = 0, 0, 0

# Pour commencer à un certain nombre de morts / infectés, on les définit.
# Il faudra par ailleurs remettre DateDeb à 0, 0, 0 pour que cela marche.
NbIDeb = 1
NbMDeb = 1000
    

# On définit la plage k utilisée pour lisser les données réelles
kI = 2
kM = 2


#On définit la période p : durée de chaque phase (tq beta fixe) pour la simulation (jours)
p = 6



# Ne pas modifier les commandes ci-dessous : dépendent de celles au dessus

data = pd.read_excel("https://www.ecdc.europa.eu/sites/default/files/documents/COVID-19-geographic-disbtribution-worldwide.xlsx")
data = data[data["countryterritoryCode"] == CodePays]
data = (np.array(data))[::-1]

i = 0
if (DateDeb == (0, 0, 0)):
    while (data[i][4]<NbIDeb):
        i = i+1
else :
    while ((data[i][1], data[i][2], data[i][3]) != (DateDeb)):
        i = i+1
    
data = data[i:,:]

IpJreel = data[:,4]
IpJliss = MoyMob(IpJreel, kI)

MpJreel = data[:,5]
MpJliss = MoyMob(MpJreel, kM)

LJours = data[:, 0]
n = len(LJours)

                    # Données : Population

# N = Population totale étudiée
N = data[0,9]

# D0 = Nombre initial de décédés
D0 = data[0,5]

# R0 = Nombre initial de rétablis
R0 = 0

# I0 = Nombre initial d'infectés
I0 = data[0,4]

# E0 = Nombre initial d'exposés = infectés non infectieux / contaminés non contaminants (en général 0)
E0 = 0

# S0 = Nombre initial de sains / susceptibles d'êtres infectés (le reste de la population)
S0 = N - E0 - I0 - R0 - D0

# y0 = Vecteur comportant ces données initiales
y0 = S0, E0, I0, R0, D0


            # Simulation et Graphes

def Reste(Tbeta):
    S, E, I, R, D, t = SEIRD_phases(y0, Tn, N, Tbeta, sigma, gamma, mu)
    MpJtheo = I*mu
    Res = MpJtheo - MpJliss
    return np.linalg.norm(Res, 2)


Tn = Periode(n, p)
Tbetai = [1]*len(Tn)

Tbeta = least_squares(Reste, Tbetai, bounds=(0, 1)).x

S, E, I, R, D, t = SEIRD_phases(y0, Tn, N, Tbeta, sigma, gamma, mu)

MpJtheo = I*mu
Ce premier graphe est l'évolution du $\beta$ qui permet de diviser l'épidémie en différentes phases. Nous pouvons remarquer son évolution particulièrement au cours du confinement. Le dernier $\beta$ calculé est :
In [311]:
print(Tbeta[-1])

0.2292269297388847
In [310]:
fig = plt.figure(facecolor='w', figsize=(8,6), dpi=80)
plt.plot(Tbeta, color="red", label = 'Beta mis à jour tous les 6 jours')
axes = plt.gca()
axes.xaxis.set_ticks(range(n))
axes.xaxis.set_ticklabels(['2020-01-31', '2020-02-06', '2020-02-12', '2020-02-18', '2020-02-24', '2020-03-01', '2020-03-07', '2020-03-13', '2020-03-19', '2020-03-25','2020-03-31', '2020-04-06', '2020-04-12', '2020-04-18', '2020-04-24','2020-04-30','2020-05-06', '2020-05-12', '2020-05-18', '2020-05-24','2020-05-30','2020-06-05','2020-06-11','2020-06-17','2020-06-23','2020-06-29','2020-07-05','2020-07-11','2020-07-17','2020-07-23','2020-07-29'], rotation = 90, color = 'black', fontsize = 8)
u = [6.33,6.33]
v = [0,1]
plt.plot(u, v, '--', color = "purple", lw = 1, label='Confinement')
a = [15.5,15.5]
b = [0,1]
plt.plot(a, b, '--', color = "purple", lw = 1)
plt.xlabel('Temps (/6jours)')
plt.ylabel('Beta')
legend = plt.legend()
plt.title("Evolution du Beta en fonction du temps",fontsize=18, color = "red")
plt.show()
Ensuite, les données les plus fiables, pour évaluer les paramètres ci-dessus, sont le nombre de morts du covid-19. Nous avons donc chercher à moyenner sur 2 jours les données réelles de morts. A partir de la courbe lisée, il est possible de trouver les paramètres correspondants pour approcher au mieux le nombre de morts théorique. En superposant les trois courbes, nous obtenons des courbes peu lisible, une preuve que le travail est cohérant.
Il est intéressant de remarquer que la courbe des personnes décédées par la Covid-19 est décalée d'une quinzaine de jours par rapport au moment où les personnes ont été infectées, d'où le pic pendant le confinement.
In [312]:
fig = plt.figure(facecolor='w', figsize=(10,8), dpi=80)
plt.suptitle('Nombre de morts du Covid-19 en Italie',fontsize=18, color="red")
plt.subplot(221)
plt.tick_params(axis = 'both', labelsize = 8)
plt.xticks(rotation = 'vertical')
plt.plot(LJours, MpJreel, 'blue', alpha=0.5, lw=2, label='Morts réels' )
plt.plot(LJours, MpJliss, 'black', alpha=0.5, lw=2, label='Morts lissés' )
plt.axvline(x="2020-03-09",color='purple',linestyle='--', label = 'Confinement')
plt.axvline(x="2020-05-04",color='purple',linestyle='--')
plt.xlabel('Temps (jours)')
plt.ylabel('Nombre (individus)')
plt.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = plt.legend()
legend.get_frame().set_alpha(0.5)

plt.subplot(222)
plt.tick_params(axis = 'both', labelsize = 8)
plt.xticks(rotation = 'vertical')
plt.plot(LJours, MpJliss, 'black', alpha=0.5, lw=2, label='Morts lissés' )
plt.plot(LJours, MpJtheo, 'green', alpha=0.5, lw=2, label='Morts théoriques')
plt.axvline(x="2020-03-09",color='purple',linestyle='--', label = 'Confinement')
plt.axvline(x="2020-05-04",color='purple',linestyle='--')
plt.xlabel('Temps (jours)')
plt.ylabel('Nombre (individus)')
plt.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = plt.legend()
legend.get_frame().set_alpha(0.5)
plt.show()

fig = plt.figure(facecolor='w', figsize=(22,20), dpi=80)
plt.subplot(223)
plt.tick_params(axis = 'both', labelsize = 8)
plt.xticks(rotation = 50)
plt.plot(LJours, MpJreel, 'blue', alpha=0.5, lw=2, label='Morts réels' )
plt.plot(LJours, MpJliss, 'black', alpha=0.5, lw=2, label='Morts lissés' )
plt.plot(LJours, MpJtheo, 'green', alpha=0.5, lw=2, label='Morts théoriques')
plt.axvline(x="2020-03-09",color='purple',linestyle='--', label = 'Confinement')
plt.axvline(x="2020-05-04",color='purple', linestyle='--')
plt.xlabel('Temps (jours)')
plt.ylabel('Nombre (individus)')
plt.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = plt.legend()
legend.get_frame().set_alpha(0.5)
plt.show()
Nous avons maintenant tous les paramètres nécessaires pour tracer le modèle SEIRD. Le choix de ne pas tracer la courbe des suceptibles s'explique par l'illisibilité du graphe, en effet nous pouvons l'observer :
In [313]:
t = LJours

fig = plt.figure(facecolor='w', figsize=(10,14), dpi=80)
plt.suptitle("Représentation du modèle SEIRD de l'Italie",fontsize=18, color="red")
plt.subplot(211)
plt.xticks(rotation = 50)
plt.tick_params(axis = 'both', labelsize = 8)
plt.plot(t, S, 'b', alpha=0.5, lw=2, label='Susceptibles')
plt.plot(t, E, 'g', alpha=0.5, lw=2, label='Exposés')
plt.plot(t, I, 'r', alpha=0.5, lw=2, label='Infectés')
plt.plot(t, R, 'y', alpha=0.5, lw=2, label='Rétablis avec immunité')
plt.plot(t, D, 'black', alpha=0.5, lw=2, label='Morts')
plt.axvline(x="2020-03-09",color='purple',linestyle='--', label = 'Confinement')
plt.axvline(x="2020-05-04",color='purple', linestyle='--')
plt.xlabel('Temps /Jours')
plt.ylabel('Population')
plt.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = plt.legend()
legend.get_frame().set_alpha(0.5)

plt.subplot(212)
plt.xticks(rotation = 50)
plt.tick_params(axis = 'both', labelsize = 8)
plt.plot(t, E, 'g', alpha=0.5, lw=2, label='Exposés')
plt.plot(t, I, 'r', alpha=0.5, lw=2, label='Infectés')
plt.plot(t, R, 'y', alpha=0.5, lw=2, label='Rétablis avec immunité')
plt.plot(t, D, 'black', alpha=0.5, lw=2, label='Morts')
plt.axvline(x="2020-03-09",color='purple',linestyle='--', label = 'Confinement')
plt.axvline(x="2020-05-04",color='purple', linestyle='--')
plt.xlabel('Temps /Jours')
plt.ylabel('Population')
plt.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = plt.legend()
legend.get_frame().set_alpha(0.5)

plt.show()
Par Léa Comin
In [ ]: