Stochastic SIR using SDEint python package

Question

I want to use the SDEint package to give a numerical solution (plot) of the following stochastic SIR model. Namely, a system of SDEs.

$$\begin{cases} dS = -\beta SIdt - \sigma SIdW \\ dI = (\beta SI - \gamma I)dt + \sigma SIdW \\ dR = \gamma Idt \end{cases}$$

Here is the code I have been trying so far. Below I will explain my issues.

import matplotlib.pyplot as plt
import numpy as np
import sdeint

tspan = np.linspace(0, 5, 1001)

y0 = np.array([900/1000, 100/1000, 0.,])
S = y0[0]
I = y0[1]
R = y0[2]

def f(y, t):

    f0 = -.80*S*I
    f1 = .80*S*I - 1.1*I
    f2 = 1.1*I
    return np.array([f0, f1, f2,])


def G(y,t):
    return np.array([[-0.9*S*I, 0, 0],
              [0.9*S*I, 0, 0],
              [0,0,0]])

result = sdeint.itoint(f, G, y0, tspan )

plt.plot(result)
plt.show()

$dS$ and $dI$ are effected by the same noise $dW$, previous I had used np.diag but that is for independent Wiener processes but they are not independent here, they should be the same. That is if I look at just the diffusion portion of the system I can write that as

$$ \left[ \begin{matrix} -\sigma SI \\ \sigma SI \\ 0 \end{matrix} \right]dW$$

Unsure of how to input this into I use the following matrix/array in my definition of G in my code

$$ \left[ \begin{matrix} -\sigma SI & 0 & 0 \\ \sigma SI & 0 & 0 \\ 0&0&0 \end{matrix} \right] $$

which when multiplied by a column vector of say three Wiener processes will yield what I want above.

The key issue I am noticing when I run my script, I get the plotted paths, and specifically for $I$, it dipping below $0$ which should not be the case, these numerical solutions should be remaining non-negative.

My second issue is the plotted path for $S$ and $I$ are coming out as reflections of each other, which should not be exactly the case as they are not symmetric in their respective definitions.

Any input, or recommendations for another package or method is appreciated.

Answer 1

This model is implemented using Julia's DifferentialEquations.jl in this tutorial. Here's a version of that code:

using DifferentialEquations
using StochasticDiffEq
using DiffEqCallbacks
using Random
using SparseArrays
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

A = zeros(3,2)

# Make `g` write the sparse matrix values
function sir_noise!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    ifrac = β*c*I/N*S
    #rfrac = γ*I
    du[1,1] = -sqrt(ifrac)
    du[2,1] = sqrt(ifrac)
    #du[2,2] = -sqrt(rfrac)
    #du[3,2] = sqrt(rfrac)
end;

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[2]
end;
function affect!(integrator)
  integrator.u[2] = 0.0
end;
cb = ContinuousCallback(condition,affect!);

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0]; # S,I,R
p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob_sde = SDEProblem(sir_ode!,sir_noise!,u0,tspan,p,noise_rate_prototype=A);
sol_sde = solve(prob_sde,LambaEM(),callback=cb);

df_sde = DataFrame(sol_sde(t)')
df_sde[!,:t] = t;

@df df_sde plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

That entire repository describes how to implement different varients of epidemic models, from delayed jump equations to SDEs so it's a nice resource.

Notice that there are two things here:

1. In the Chemical Langevin Equation (CLE), the noise term is not the reaction but the square root of the reaction. Your equation was missing the square root which is one thing that can more easily pull it negative. In fact, your equation is not guarenteed to have a positive solution for a similar reason to the fact that u' = -sqrt(u) does not have a positive solution.

2. This uses the ContinuousCallback mechanism for more readily impose positivity.

For an extended example which showcases how to do things like (GPU) parallelism, see the extended discussion in the QuantEcon stochastic differential equations notes.