I am creating the simplest and lowest level SIR simulation with LightGraphs and ModelingToolkit (I am a beginner). For now, I have a simple network of 10 nodes and 7 edges. And it is a directed graph.

Here is the top part of my code up until the line that crashes:

```
using ModelingToolkit
using DiffEqBase
using DiffEqJump
using LightGraphs
nedges = 7
nverts = 10
dg = DiGraph(nverts, nedges)
# Allow for different beta on each edge
@parameters t β[1:nedges] γ[1:nedges] ;
@variables S[1:nverts](t);
@variables I[1:nverts](t);
@variables R[1:nverts](t);
rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
for (i,e) ∈ enumerate(edges(dg))]
rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]]) # note: src to src, yet there is no edge
for (i,e) ∈ enumerate(edges(dg))]
rxs = vcat(rxsS, rxsI)
rs = ReactionSystem(rxs, t, u, k)
js = convert(JumpSystem, rs)
S0 = Array{Float64,1}(undef, nverts)
I0 = Array{Float64,1}(undef, nverts)
R0 = Array{Float64,1}(undef, nverts)
fill!(S0, 1.);
fill!(I0, 0.);
fill!(R0, 0.);
S0[1] = 0.; # One person is infected
I0[1] = 1.;
R0[1] = 0.;
# Two column vectors
γ = Array{Float64,1}(undef, nedges)'
β = Array{Float64,1}(undef, nedges)'
fill!(γ, 0.25)
fill!(β, 0.50)
# Row vector
pmap = hcat(β,γ)
u0 = vcat(S0,I0,R0)
tspan = (0.0,10.0)
dprob = DiscreteProblem(js, u0, tspan, pmap)
```

`DiscreteProblem`

throws an error:

```
ERROR: MethodError: Cannot `convert` an object of type
Float64 to an object of type
Variable
Closest candidates are:
convert(::Type{Variable}, ::Operation) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/ModelingToolkit.jl:72
convert(::Type{T}, ::T) where T at essentials.jl:171
Variable(::Any) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:50
...
Stacktrace:
[1] varmap_to_vars(::Array{Float64,1}, ::Array{Variable,1}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:251
[2] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Array{Float64,1}, ::Tuple{Float64,Float64}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
[3] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Array{Float64,1}, ::Tuple{Float64,Float64}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
[4] top-level scope at none:0
```

which I assume is due to the way I combine my reactions:

```
rxs = vcat(rxsS, rxsI)
```

Alternatively, it may have to do with my unconventional reactions equations, which I am not sure are allowed:

```
rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
for (i,e) ∈ enumerate(edges(dg))]
rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]]) # note: src to src, yet there is no edge
for (i,e) ∈ enumerate(edges(dg))]
```

Notice that rxsS and rxsI have reactions defined as

```
I(src) + S(dst) --> I(dst)
I(src) --> I(src)
```

which is not the usual way of writing reaction equations on a network, but the ones above reflect what is really happening. The infected person at a node eventually recovers at that same node. There is interaction with another node.

Any help would be greatly appreciated.