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)[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鈥檛 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鈥檛 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]
(-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[1]
-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鈥檓 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)[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^{-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鈥檚 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:
 [1] 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
 [2] 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
 [3] eigmax(::Array{Num,2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 [4] top-level scope at In[6]:1
 [5] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

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:
 [1] top-level scope at In[7]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

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:
 [1] safemin(::Type{T} where T) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\util.jl:6
 [2] _scale!(::Array{Num,2}) at C:\Users\claud\.julia\packages\GenericSchur\YzsRq\src\util.jl:15
 [3] 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
 [4] 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
 [5] 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
 [6] 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
 [7] eigmax(::Array{Num,2}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:374
 [8] top-level scope at In[85]:1
 [9] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

Any help is greatly apprecciated.

1 Like