1D PDE: SIS diffusion model - Julia PDE and UDE

Source

\begin{align} \frac{\partial S}{\partial t} &= d_S S_{xx} - \beta(x)\frac{SI}{S+I} + \gamma(x)I \\ \frac{\partial I}{\partial t} &= d_I I_{xx} + \beta(x)\frac{SI}{S+I} - \gamma(x)I \end{align}

(1)

where $x \in (0, 1)$

Solve the steady-state problem $\frac{\partial S}{\partial t} = \frac{\partial I}{\partial t} = 0$

The boundary condition: $\frac{\partial S}{\partial x} = \frac{\partial I}{\partial x} = 0$ for x = 0, 1

The conservative relationship: $\int^{1}_{0} (S(x) + I(x) ) dx = 1$

Notations:

$x$ : location
$t$ : time
$S(x, t)$ : the density of susceptible populations
$I(x, t)$ : the density of infected populations
$d_S$ / $d_I$ : the diffusion coefficients for susceptible and infected individuals
$\beta(x)$ : transmission rates
$\gamma(x)$ : recovery rates

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots

Setup parameters, variables, and differential operators

function sis_model(;dx = 0.01, order = 2)
    @independent_variables t x
    @parameters dS=0.5 dI=0.1 brn=3 ϵ=0.1
    @variables S(..) I(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dxx = Differential(x)^2
    # Helper functions
    γ(x) = x + 1
    ratio(x, brn, ϵ) = brn + ϵ * sinpi(2x)

    # 1D PDE for disease spreading
    eqs = [
        Dt(S(t, x)) ~ dS * Dxx(S(t, x)) - ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) + γ(x) * I(t, x),
        Dt(I(t, x)) ~ dI * Dxx(I(t, x)) + ratio(x, brn, ϵ) * γ(x) * S(t, x) * I(t, x) / (S(t, x) + I(t, x)) - γ(x) * I(t, x)
    ]

    # Boundary conditions (including initial conditions)
    bcs = [
        S(0, x) ~ 0.9 + 0.1 * sinpi(2x),
        I(0, x) ~ 0.1 + 0.1 * cospi(2x),
        Dx(S(t, 0)) ~ 0.0,
        Dx(S(t, 1)) ~ 0.0,
        Dx(I(t, 0)) ~ 0.0,
        Dx(I(t, 1)) ~ 0.0
    ]

    # Space and time domains
    domains = [
        t ∈ Interval(0.0, 10.0),
        x ∈ Interval(0.0, 1.0)
    ]

    # Build the PDE system
    @named pdesys = PDESystem(eqs, bcs, domains,
        [t, x], ## Independent variables
        [S(t, x), I(t, x)],  ## Dependent variables
        [dS, dI, brn, ϵ],    ## parameters
    )
    # Finite difference method (FDM) converts the PDE system into an ODE problem
    discretization = MOLFiniteDifference([x => dx], t, approx_order=order)
    prob = discretize(pdesys, discretization)
    return (; prob, t, x)
end

sis_model (generic function with 1 method)

There are 202 ODEs in the system.

@time "Build problem" prob, t, x = sis_model()

Build problem: 113.901841 seconds (257.08 M allocations: 12.024 GiB, 2.14% gc time, 96.10% compilation time: 3% of which was recompilation)

(prob = SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ModelingToolkitBase.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x13acd881, 0x1aa1c4a0, 0x02b29da7, 0x359015f8, 0x7050ea78), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1401282876548370056, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xc6e98927, 0xf57adf3c, 0xd1d4ae53, 0x1a8c5e6e, 0x9b5ef3db), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Union{Nothing, SciMLBase.OverrideInitData}, Union{Nothing, SciMLBase.ODENLStepData}}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}(SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x13acd881, 0x1aa1c4a0, 0x02b29da7, 0x359015f8, 0x7050ea78), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1401282876548370056, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xc6e98927, 0xf57adf3c, 0xd1d4ae53, 0x1a8c5e6e, 0x9b5ef3db), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}, Nothing, ModelingToolkitBase.System, Union{Nothing, SciMLBase.OverrideInitData}, Union{Nothing, SciMLBase.ODENLStepData}}(ModelingToolkitBase.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x13acd881, 0x1aa1c4a0, 0x02b29da7, 0x359015f8, 0x7050ea78), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1401282876548370056, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xc6e98927, 0xf57adf3c, 0xd1d4ae53, 0x1a8c5e6e, 0x9b5ef3db), Nothing}}(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0x13acd881, 0x1aa1c4a0, 0x02b29da7, 0x359015f8, 0x7050ea78), Nothing}(nothing), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__argₛᵧₘ1401282876548370056, :___mtkunknowns___, :___mtkparameters___, :__argₛᵧₘ1249670312406203976), ModelingToolkitBase.var"#_RGF_ModTag", ModelingToolkitBase.var"#_RGF_ModTag", (0xc6e98927, 0xf57adf3c, 0xd1d4ae53, 0x1a8c5e6e, 0x9b5ef3db), Nothing}(nothing)), LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, ModelingToolkitBase.ObservedFunctionCache{ModelingToolkitBase.System, Nothing}(Model pdesys:
Equations (198):
  198 standard: see equations(pdesys)
Unknowns (198): see unknowns(pdesys)
  (I(t))[100]
  (I(t))[99]
  (I(t))[98]
  (I(t))[97]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ
Observed (6): see observed(pdesys), Dict{Any, Any}(), false, false, ModelingToolkitBase, false, true, nothing), nothing, Model pdesys:
Equations (198):
  198 standard: see equations(pdesys)
Unknowns (198): see unknowns(pdesys)
  (I(t))[100]
  (I(t))[99]
  (I(t))[98]
  (I(t))[97]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ
Observed (6): see observed(pdesys), nothing, nothing), [0.19980267284282716, 0.1992114701314478, 0.19822872507286887, 0.1968583161128631, 0.19510565162951538, 0.19297764858882513, 0.190482705246602, 0.18763066800438638, 0.18443279255020154, 0.18090169943749476  …  0.9587785252292473, 0.9535826794978997, 0.9481753674101716, 0.9425779291565073, 0.9368124552684678, 0.9309016994374948, 0.9248689887164855, 0.9187381314585725, 0.9125333233564304, 0.9062790519529313], (0.0, 10.0), ModelingToolkitBase.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}([0.5, 0.1, 3.0, 0.1], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], (), (), (), ()), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}(), MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, ModelingToolkitBase.PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}(MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}(PDEBase.VariableMap(Any[I(t, x), S(t, x)], SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}[x], Symbolics.Num[dS, dI, brn, ϵ], t, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Tuple{Float64, Float64}}(t => (0.0, 10.0), x => (0.0, 1.0)), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, ReadOnlyArrays.ReadOnlyVector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.SmallVec{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}}}}(I => [t, x], S => [t, x]), SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}[S, I], Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Int64}(x => 1), Dict{Int64, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}(1 => x), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Symbolics.Num}()), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}}(I(t, x) => [(I(t))[1], (I(t))[2], (I(t))[3], (I(t))[4], (I(t))[5], (I(t))[6], (I(t))[7], (I(t))[8], (I(t))[9], (I(t))[10]  …  (I(t))[92], (I(t))[93], (I(t))[94], (I(t))[95], (I(t))[96], (I(t))[97], (I(t))[98], (I(t))[99], (I(t))[100], (I(t))[101]], S(t, x) => [(S(t))[1], (S(t))[2], (S(t))[3], (S(t))[4], (S(t))[5], (S(t))[6], (S(t))[7], (S(t))[8], (S(t))[9], (S(t))[10]  …  (S(t))[92], (S(t))[93], (S(t))[94], (S(t))[95], (S(t))[96], (S(t))[97], (S(t))[98], (S(t))[99], (S(t))[100], (S(t))[101]]), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(x => 0.0:0.01:1.0), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}(x => 0.0:0.01:1.0), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Float64}(x => 0.01), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}(I(t, x) => CartesianIndices((101,)), S(t, x) => CartesianIndices((101,))), Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}(I(t, x) => CartesianIndices((101,)), S(t, x) => CartesianIndices((101,))), nothing), MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Symbolics.Num, Float64}(x => 0.01), t, 2, MethodOfLines.UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}()), PDESystem
Equations: Symbolics.Equation[Differential(t, 1)(S(t, x)) - ((-(1 + x)*(brn + sinpi(2x)*ϵ)*I(t, x)*S(t, x)) / (I(t, x) + S(t, x))) - dS*Differential(x, 2)(S(t, x)) - (1 + x)*I(t, x) ~ 0, -(((1 + x)*(brn + sinpi(2x)*ϵ)*I(t, x)*S(t, x)) / (I(t, x) + S(t, x))) + Differential(t, 1)(I(t, x)) - dI*Differential(x, 2)(I(t, x)) + (1 + x)*I(t, x) ~ 0]
Boundary Conditions: Symbolics.Equation[S(0, x) ~ 0.9 + 0.1sinpi(2x), I(0, x) ~ 0.1 + 0.1cospi(2x), Differential(x, 1)(S(t, 0)) ~ 0.0, Differential(x, 1)(S(t, 1)) ~ 0.0, Differential(x, 1)(I(t, 0)) ~ 0.0, Differential(x, 1)(I(t, 1)) ~ 0.0]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0 .. 10.0), Symbolics.VarDomainPairing(x, 0.0 .. 1.0)]
Dependent Variables: Symbolics.Num[I(t, x), S(t, x)]
Independent Variables: Symbolics.Num[t, x]
Parameters: Symbolics.Num[dS, dI, brn, ϵ]
Default Parameter ValuesModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}}(), false, Base.RefValue{Any}(Model pdesys:
Equations (202):
  202 standard: see equations(pdesys)
Unknowns (202): see unknowns(pdesys)
  (I(t))[1]
  (I(t))[2]
  (I(t))[3]
  (I(t))[4]
  ⋮
Parameters (4): see parameters(pdesys)
  dS
  dI
  brn
  ϵ), nothing, Pair{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}[(I(t))[1] => 0.1 + 0.1cospi(0.0), (I(t))[2] => 0.1 + 0.1cospi(0.02), (I(t))[3] => 0.1 + 0.1cospi(0.04), (I(t))[4] => 0.1 + 0.1cospi(0.06), (I(t))[5] => 0.1 + 0.1cospi(0.08), (I(t))[6] => 0.1 + 0.1cospi(0.1), (I(t))[7] => 0.1 + 0.1cospi(0.12), (I(t))[8] => 0.1 + 0.1cospi(0.14), (I(t))[9] => 0.1 + 0.1cospi(0.16), (I(t))[10] => 0.1 + 0.1cospi(0.18)  …  (S(t))[92] => 0.9 + 0.1sinpi(1.82), (S(t))[93] => 0.9 + 0.1sinpi(1.84), (S(t))[94] => 0.9 + 0.1sinpi(1.86), (S(t))[95] => 0.9 + 0.1sinpi(1.88), (S(t))[96] => 0.9 + 0.1sinpi(1.9), (S(t))[97] => 0.9 + 0.1sinpi(1.92), (S(t))[98] => 0.9 + 0.1sinpi(1.94), (S(t))[99] => 0.9 + 0.1sinpi(1.96), (S(t))[100] => 0.9 + 0.1sinpi(1.98), (S(t))[101] => 0.9 + 0.1sinpi(2.0)])), t = t, x = x)

Solving time-dependent SIS epidemic model¶

KenCarp47 and FBDF are good at solving reaction-diffusion problems.

@time sol = solve(prob, FBDF(), saveat=0.2)

  7.078958 seconds (18.38 M allocations: 838.811 MiB, 1.51% gc time, 98.34% compilation time: <1% of which was recompilation)

retcode: Success
Interpolation: Dict{Symbolics.Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}
t: 51-element Vector{Float64}:
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2
  1.4
  1.6
  1.8
  ⋮
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0ivs: 2-element Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymbolicUtils.SymReal}}:
 t
 xdomain:([0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8  …  8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0], 0.0:0.01:1.0)
u: Dict{Symbolics.Num, Matrix{Float64}} with 2 entries:
  I(t, x) => [0.2 0.199803 … 0.199803 0.2; 0.204072 0.203977 … 0.245383 0.24554…
  S(t, x) => [0.904194 0.906279 … 0.893721 0.895806; 0.879621 0.879606 … 0.7900…

Grid points

discrete_x = sol[x]
discrete_t = sol[t]

51-element Vector{Float64}:
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2
  1.4
  1.6
  1.8
  ⋮
  8.4
  8.6
  8.8
  9.0
  9.2
  9.4
  9.6
  9.8
 10.0

Results (Matrices)

@variables S(..) I(..)
sol
S_solution = sol[S(t, x)]
I_solution = sol[I(t, x)]

51×101 Matrix{Float64}:
 0.2       0.199803  0.199211  0.198229  …  0.199211  0.199803  0.2
 0.204072  0.203977  0.203691  0.203211     0.2449    0.245383  0.245544
 0.239567  0.239532  0.239427  0.239248     0.320235  0.320608  0.320732
 0.300266  0.300276  0.300306  0.300352     0.413704  0.413992  0.414088
 0.376036  0.376073  0.376184  0.376363     0.503742  0.503967  0.504042
 0.452112  0.452158  0.452297  0.452525  …  0.573952  0.574128  0.574187
 0.516857  0.516901  0.517035  0.517254     0.620401  0.620539  0.620584
 0.566217  0.566254  0.566367  0.566552     0.647111  0.647219  0.647254
 0.601153  0.601183  0.601272  0.601419     0.660266  0.660352  0.66038
 0.624956  0.624979  0.625049  0.625163     0.665592  0.665662  0.665685
 ⋮                                       ⋱                      ⋮
 0.67636   0.676364  0.676374  0.676391     0.656473  0.656498  0.656507
 0.676357  0.676361  0.676372  0.676388     0.656475  0.6565    0.656509
 0.676354  0.676358  0.676369  0.676385     0.656478  0.656503  0.656511
 0.676351  0.676354  0.676365  0.676382  …  0.65648   0.656506  0.656514
 0.676347  0.67635   0.676361  0.676378     0.656483  0.656509  0.656517
 0.676343  0.676346  0.676357  0.676374     0.656487  0.656512  0.656521
 0.676338  0.676342  0.676352  0.676369     0.656491  0.656516  0.656524
 0.676332  0.676336  0.676347  0.676363     0.656495  0.656521  0.656529
 0.676325  0.676329  0.67634   0.676357  …  0.656501  0.656526  0.656534

Visualize the solution

surface(discrete_x, discrete_t, S_solution, xlabel="Location", ylabel="Time", title="Susceptible")

surface(discrete_x, discrete_t, I_solution, xlabel="Location", ylabel="Time", title="Infectious")

This notebook was generated using Literate.jl.