So I would like to explicitly evaluate the R_0 of an epidemiological model written in DifferentialEquations.jl using ModelingToolkit.jl.
Here is a simple SIRD model:
using DifferentialEquations, ModelingToolkit
# Model parameters
β = 0.05 # infection rate
λ_R = 0.55 # inverse of transition time from infected to recovered
λ_D = 0.83 # inverse of transition time from infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, λ_R, λ_D]...)
# regional contact matrix and regional population
## regional contact matrix
regional_all_contact_matrix = [3.45536 0.485314 0.506389 0.123002 ; 0.597721 2.11738 0.911374 0.323385 ; 0.906231 1.35041 1.60756 0.67411 ; 0.237902 0.432631 0.726488 0.979258] # 4x4 contact matrix
## regional population stratified by age
N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.
# Initial conditions
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...)
# Time
𝒯 = (1.0,20);
function SIRD!(du,u,p,t)
# Parameters to be calibrated
β, λ_R, λ_D = p
# initializa this parameter (death probability stratified by age, taken from literature)
δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
δ = vcat(repeat([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))
C = regional_all_contact_matrix
# State variables
S = @view u[4*0+1:4*1]
I = @view u[4*1+1:4*2]
R = @view u[4*2+1:4*3]
D = @view u[4*3+1:4*4]
D_tot = @view u[4*4+1:4*5]
# Differentials
dS = @view du[4*0+1:4*1]
dI = @view du[4*1+1:4*2]
dR = @view du[4*2+1:4*3]
dD = @view du[4*3+1:4*4]
dD_tot = @view du[4*4+1:4*5]
# Force of infection
Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]]
# System of equations
@. dS = -Λ*S
@. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
@. dR = λ_R*(1-δ)*I
@. dD = λ_D*δ*I
@. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]
end;
# create problem and modeltoolkitize it
problem = ODEProblem(SIRD!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem);
# The theory states that R₀ = eigmax(J_f*J_V⁻¹), in the early stage limit.
# If i want to calculate J_V, i just have to take all but the first 4 columns and 4 rows of the jacobian of sys and evaluate it at S = [1,1,1,1] and at the above parameter values.
#evaluate jacobian
J = ModelingToolkit.calculate_jacobian(sys);
#strip the first 4 columns and rows
J_V = J[5:end,5:end];
But then if I try to evaluate it at specific parameter values, I don’t get any substitution (the \alpha_i are still there):
subs_params = [param =>p for (param,p) in zip(sys.ps, 𝒫)]
J_V = substitute.(J_V, (subs_params,))
16×16 Array{Expression,2}:
-1 * (0.99997α₂ + 3.0e-5α₃) + 4.77782e-6 * x₁(t) * α₁ … Constant(0) Constant(0)
8.26486e-7 * x₂(t) * α₁ Constant(0) Constant(0)
1.25307e-6 * x₃(t) * α₁ Constant(0) Constant(0)
3.28954e-7 * x₄(t) * α₁ Constant(0) Constant(0)
0.99997α₂ Constant(0) Constant(0)
Constant(0) … Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
Constant(0) … Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ … Constant(0) Constant(0)
I also can’t substitute variable values (here I try to set S = [1,1,1,1]):
subs_var = [var => 0 for var in sys.states[1:4] ]
J_V = substitute.(J_V, (subs_var,))
16×16 Array{Expression,2}:
-1 * (0.99997α₂ + 3.0e-5α₃) + 4.77782e-6 * x₁(t) * α₁ … Constant(0) Constant(0)
8.26486e-7 * x₂(t) * α₁ Constant(0) Constant(0)
1.25307e-6 * x₃(t) * α₁ Constant(0) Constant(0)
3.28954e-7 * x₄(t) * α₁ Constant(0) Constant(0)
0.99997α₂ Constant(0) Constant(0)
Constant(0) … Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
Constant(0) … Constant(0) Constant(0)
Constant(0) Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ Constant(0) Constant(0)
3.0e-5α₃ … Constant(0) Constant(0)
So how may I substitute variable and parameter values in the jacobian?
Is there an alternative way?
Thank you very much