# Substituting state and parameter values in jacobian of DifferentialEquations.jl system using ModelingToolkit

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)]) for i in 1:size(C)]

# System of equations
@. dS = -Λ*S
@. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
@. dR = λ_R*(1-δ)*I
@. dD = λ_D*δ*I
@. dD_tot = dD+dD+dD+dD

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

2 Likes

Huh, your code works fine when I run it.

julia> J_V
(-1 * ((0.99997 * α₂) + (3.0e-5 * α₃))) + (4.777823254167543e-6 * x₁(t) * α₁)

julia> subs_params = [param =>p for (param,p) in zip(sys.ps, 𝒫)];

julia> sJ_V = substitute.(J_V, (subs_params,));

julia> sJ_V
-0.5500084000000001 + (4.777823254167543e-6 * x₁(t) * 0.05)


Is everything updated?

julia> versioninfo()
Julia Version 1.5.3
Commit 788b2c77c1 (2020-11-09 13:37 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
JULIA_REVISE_INCLUDE = 1

(jl_nkocE5) pkg> st
Status /tmp/jl_nkocE5/Project.toml
[0c46a032] DifferentialEquations v6.15.0
[961ee093] ModelingToolkit v4.0.5
[91a5bcdd] Plots v1.9.1

(jl_nkocE5) pkg>

4 Likes

Thank you very much, that was indeed the case.I write here my experience for others who may face the same problem.
I saw that I had ModelingToolkit 3.x.y installed, so I tried to do

pkg> update ModelingToolkit


and it just went up to 3.20.z or something. So I created a new julia environment ( went to the project folder and did pkg> activate .) and then installed there everything I need. I noticed that it installs ModelingToolkit v4.0.5 , but when it then installs JuliaVariables.jl (required by other packages) it downgrades ModelingToolkit to 3.x.y . So for me the solution was to create an ad-hoc environment with the minimum packages required to do what I need.

Anyway, now I’m facing another related issue.

The model is the same as before:

using DifferentialEquations, ModelingToolkit, LinearAlgebra

# 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)]) for i in 1:size(C)]

# System of equations
@. dS = -Λ*S
@. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
@. dR = λ_R*(1-δ)*I
@. dD = λ_D*δ*I
@. dD_tot = dD+dD+dD+dD

end;

# create problem and modeltoolkitize it
problem = ODEProblem(SIRD!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem);



The theory states that:

R₀ = eigmax(J_F * J_V^{-1})

in the early stage limit.
I would like to implement the next-generation matrix method to be able to extend this to more complex models .

So let’s evaluate J_F:

nic = 1 # number of infective compartments
nac = 4 # number of age classes

# define F-system
F = [eq.rhs for eq in sys.eqs[nac+1:nac+nic*nac]];
# substitute paremeter values
sF = substitute.(F, ([param => p for (param,p) in zip(sys.ps,𝒫)],));

# substitute S = [1,1,1,1] (early stage limit)
sFs = substitute.(sF, ([var => 1.0 for var in sys.states[1:nac]],));

# evaluate jacobian
J_F = ModelingToolkit.jacobian(sFs,sys.states[nac+1:nac+nic*nac])

4×4 Array{Num,2}:
-0.550008     2.77592e-8   1.9023e-8    4.35582e-9
4.13243e-8  -0.550011     3.42366e-8   1.14519e-8
6.26536e-8   7.72413e-8  -0.550582     2.3872e-8
1.64477e-8   2.47458e-8   2.72912e-8  -0.567125


Evaluate J_V:

# define V-system
V = [eq.rhs for eq in sys.eqs[nac+1:nac+nic*nac]];
# substitute paremeter values
sV = substitute.(V, ([param => p for (param,p) in zip(sys.ps,𝒫)],));

# substitute S = [0,0,0,0] (by definition of the V-system)
sVs = substitute.(sV, ([var => 0.0 for var in sys.states[1:nac]],));

# evaluate jacobian
J_V = ModelingToolkit.jacobian(sVs,sys.states[nac+1:nac+nic*nac])

4×4 Array{Num,2}:
-0.550008   0          0          0
0         -0.550011   0          0
0          0         -0.550582   0
0          0          0         -0.567125


Compute R_0:

eigmax(J_F*inv(J_V))

MethodError: no method matching eigvals!(::Array{Num,2}; permute=true, scale=true)
Closest candidates are:
eigvals!(!Matched::SymTridiagonal{var"#s828",V} where V<:AbstractArray{var"#s828",1} where var"#s828"<:Union{Float32, Float64}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\tridiag.jl:292 got unsupported keyword arguments "permute", "scale"
eigvals!(!Matched::SymTridiagonal{var"#s828",V} where V<:AbstractArray{var"#s828",1} where var"#s828"<:Union{Float32, Float64}, !Matched::UnitRange) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\tridiag.jl:295 got unsupported keyword arguments "permute", "scale"
eigvals!(!Matched::SymTridiagonal{var"#s828",V} where V<:AbstractArray{var"#s828",1} where var"#s828"<:Union{Float32, Float64}, !Matched::Real, !Matched::Real) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\tridiag.jl:300 got unsupported keyword arguments "permute", "scale"
...

Stacktrace:
 eigvals(::Array{Num,2}; kws::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol,Symbol},NamedTuple{(:permute, :scale),Tuple{Bool,Bool}}}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:326
 eigmax(::Array{Num,2}; permute::Bool, scale::Bool) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 eigmax(::Array{Num,2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 top-level scope at In:1


Converting to Array{Float64,2} does not work:

eigmax(Array{Float64,2}J_F*inv(J_V))

MethodError: no method matching *(::Type{Array{Float64,2}}, ::Array{Num,2})
Closest candidates are:
*(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:538
*(!Matched::ChainRulesCore.One, ::Any) at C:\Users\claud\.julia\packages\ChainRulesCore\LsNyE\src\differential_arithmetic.jl:98
*(!Matched::LightGraphs.LinAlg.Noop, ::Any) at C:\Users\claud\.julia\packages\LightGraphs\HsNig\src\linalg\graphmatrices.jl:226
...

Stacktrace:
 top-level scope at In:1


So how may I get eigmax to work?

Curiously, if I do the same with a more complex SEIHRD model, I get a different error:

eigmax(J_F*inv(J_V))

MethodError: no method matching floatmin(::Type{Num})
Closest candidates are:
floatmin() at float.jl:791
floatmin(!Matched::Type{Float16}) at float.jl:738
floatmin(!Matched::Type{Float32}) at float.jl:739
...

Stacktrace:
 safemin(::Type{T} where T) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\util.jl:6
 _scale!(::Array{Num,2}) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\util.jl:15
 gschur!(::Array{Num,2}; wantZ::Bool, scale::Bool, kwargs::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:permute,),Tuple{Bool}}}) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\GenericSchur.jl:619
 eigvals!(::Array{Num,2}; kwargs::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol,Symbol},NamedTuple{(:permute, :scale),Tuple{Bool,Bool}}}) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\GenericSchur.jl:15
 eigvals(::Array{Num,2}; kws::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol,Symbol},NamedTuple{(:permute, :scale),Tuple{Bool,Bool}}}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:326
 eigmax(::Array{Num,2}; permute::Bool, scale::Bool) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 eigmax(::Array{Num,2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 top-level scope at In:1