Calcoli scientifici con Julia/Modello SIR

Il modello SIR è un modello compartimentale semplice per la dinamica di una malattia infettiva. La popolazione è divisa in tre classi:

• $S(t)$: suscettibili
• $I(t)$: infetti
• $R(t)$: rimossi (guariti o immuni)

Le equazioni sono:

• $\frac{dS}{dt}=-\beta SI$
• $\frac{dI}{dt}=\beta SI-\gamma I$
• $\frac{dR}{dt}=\gamma I$

Dove $\beta$ è il tasso di trasmissione e $\gamma$ il tasso di guarigione.

La curva dei suscettibili S(t) tende a a decrescere, mentre la curva degli infetti I(t) cresce in modo accelerato (fase esponenziale) raggiungendo un picco ben definito per poi decrescere perché il pool di suscettibili si è ridotto sufficientemente tendendo a 0, la curva dei rimossi (guariti o immuni) R(t) cresce raggiungendo il massimo.

Il tasso di trasmissione misura la velocità con cui un suscettibile entra in contatto con un infetto e contrae la malattia. Più grande è $\beta$, più rapidamente si diffonde l’infezione.

Il tasso di guarigione misura la velocità con cui gli infetti guariscono (o vengono rimossi dall’epidemia perché immuni o deceduti).Più grande è $\gamma$, più rapidamente si estingue l’infezione.

Attraverso il parametro $R_0=\frac{\beta }{\gamma }$ si stabilisce se l'epidemia tende a crescere o a estinguersi:

• Se $R_0>1$ l’epidemia tende a crescere.
• Se $R_0<1$ l’epidemia tende a spegnersi.

Esempio realistico:

• Influenza stagionale: $\beta =0,4$, $\gamma =0,2$ $R_0=2$
• COVID-19 (ceppo iniziale): $\beta =0,5-0,6$ , $\gamma =0,1$ $R_0=5-6$

Esempio in Julia

Se l'epidemia dura 160 giorni, le condizioni iniziali del sistema di equazioni differenziali u0 danno il 99% di suscettibili, l'1% di infetti e 0% di guariti o immuni con tasso di trasmissione 0,5 e tasso di guarigione 0,1 allora R0=5 quindi l'epidemia cresce raggiungendo un picco come si vede dal seguente grafico tracciato da Julia, ottenuto risolvendo il sistema di equazioni differenziali:

using DifferentialEquations
using Plots

function SIR!(du, u, p, t)
S, I, R = u
β, γ = p
du\[1] = -β \* S \* I
du\[2] = β \* S \* I - γ \* I
du\[3] = γ \* I
end

β = 0.5
γ = 0.1
p = (β, γ)

u0 = \[0.99, 0.01, 0.0]
tspan = (0.0, 160.0)

prob = ODEProblem(SIR!, u0, tspan, p)
sol = solve(prob, Tsit5())

plot(sol, xlabel="Giorni", ylabel="Frazione della popolazione",
label=\["Sani" "Infetti" "Guariti"], lw=2)