Boundary conditions (?) in MethodOfLines.jl PDE

[sdwfrost]

Hi everyone,

I’m trying out the new capability to deal with integrals in MethodOfLines.jl (in the new main, v0.7.6). I’m trying to simply expand an ODE example into a PDE. For those of you interested, the ODE is a simple SIR compartmental model, and the PDE is one where individuals are structured by age. It’s very close to working, except the sum of the states over time should remain constant, yet they decline over time (see the last line). I thought that setting the derivative wrt a to 0 at the upper boundary would prevent mass from leaving the system. Can anyone spot where I’m going wrong?

using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using Plots

β = 0.05
c = 10.0
γ = 0.25
S₀ = 990.0
I₀ = 10.0
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 26.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
N(t) = Ia(S(a,t)) + Ia(I(a,t)) + Ia(R(a,t))

eqs = [Dt(S(a,t)) + Da(I(a,t)) ~ - β * c * S(a,t) * Ia(I(a,t))/N(t),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * c * S(a,t) * Ia(I(a,t))/N(t) - γ*I(a,t),
       Dt(R(a,t)) + Da(R(a,t)) ~ γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        R(a,0.0) ~ 0.0,
        R(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0,
        Da(R(amax,t)) ~ 0.0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t), R(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, Tsit5(), saveat = dt)

Smat = sol.u[S(a, t)]
Imat = sol.u[I(a, t)]
Rmat = sol.u[R(a, t)]

sum(Smat .+ Imat .+ Rmat, dims=1)

[RobertGregg]

I haven’t checked any of these solutions but you could try:

1. Lowering the tolerance of the ODE solver (abstol and reltol)
2. Changing the solver to one better suited for conservation laws or PDEs ( SSPRK22, ROCK4)
3. Using a Manifold Projection Callback

[sdwfrost]

Thanks! Switching the solver didn’t help. Changing the parameter values does affect the conservation of mass, but I worry that I haven’t specified the model correctly, so using a ManifoldProjection may just be a sticking plaster.

[Sevi]

Is there a typo in your equations? The advection term for S(a,t) is Da(I(a,t)), but should it not be Da(S(a,t)) ? I.e. people in all three categories S, I, R just age independently, all at the same rate.

Also, I would argue that the equations plus boundary/initial conditions are simply not mass-conserving. If you consider the temporal change of the mass N(a,t) = S(a,t) + I(a,t) + R(a,t) in the system, i.e. Dt(N(a,t)) = -Da(N(a,t)) and integrate it over the whole domain (to get the change of total mass), you find that the result at any time is given by N(amin,t) - N(amax,t). Since from the very beginning there is some mass at the right boundary, and the advection term will transport the mass there anyway, the loss of mass should never be zero.

From this calculation, to get zero loss of mass, the right boundary conditions would be that N(amax,t) ~ 0 although that does seem a bit counter intuitive to me as well – a steady influx and outflux with some “birth rate”/“death rate” would probably make more sense physically. Or localizing the initial condition well below amax and stopping the simulation when the first people reach that age … then either choice of b.c. doesn’t really matter and mass conservation should hold until that time is reached.

Otherwise the equations seem reasonable as far as I can tell…

[xtalax]

This may be numerical dispersion from the upwind scheme, used for the advection terms - an unfortunate numerical artefact. You might try advection_scheme=WENOScheme() which is better, though is often unstable when used without a periodic condition, or at least 2 bcs on each boundary - we need better boundary handling for this.

[sdwfrost]

Hi @Sevi

Thanks - you were right about the typo - the equations should be

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * c * S(a,t) * Ia(I(a,t))/N(t),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * c * S(a,t) * Ia(I(a,t))/N(t) - γ*I(a,t),
       Dt(R(a,t)) + Da(R(a,t)) ~ γ*I(a,t)]

as you say. However, this doesn’t fix the problem!

You’re also right about the loss of mass from the right boundary - however, I thought that I had dealt with this by setting the boundary conditions as follows:

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        R(a,0.0) ~ 0.0,
        R(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0,
        Da(R(amax,t)) ~ 0.0
        ]

The last 3 BCs are intended to prevent further flow out of the state, so over time, the mass at N(amax,t) should increase. The above also results in da(N(amax,t)) ~ 0 (or at least, that’s the intention).

[sdwfrost]

Thanks for the tip about changing the advection scheme - still no luck however. Current code below:

using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using Plots

β = 0.05
c = 10.0
γ = 0.25
S₀ = 990.0
I₀ = 10.0
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
N(t) = Ia(S(a,t)) + Ia(I(a,t)) + Ia(R(a,t))

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * c * S(a,t) * Ia(I(a,t))/N(t),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * c * S(a,t) * Ia(I(a,t))/N(t) - γ*I(a,t),
       Dt(R(a,t)) + Da(R(a,t)) ~ γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        R(a,0.0) ~ 0.0,
        R(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0,
        Da(R(amax,t)) ~ 0.0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t), R(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob = MethodOfLines.discretize(pde_system, discretization, advection_scheme=WENOScheme())
sol = solve(prob, Tsit5(), saveat = dt)

Smat = sol.u[S(a, t)]
Imat = sol.u[I(a, t)]
Rmat = sol.u[R(a, t)]

plot(sol.t, sum(Smat,dims=1)',label="S",xlabel="Time",ylabel="Number")
plot!(sol.t, sum(Imat,dims=1)',label="I")
plot!(sol.t, sum(Rmat,dims=1)',label="R")

sum(Smat .+ Imat .+ Rmat, dims=1)

[Sevi]

I see. But my point was that Da(N) ~ 0 is not the right boundary condition to begin with. Consider the change in N(t), so

\frac{\mathrm{d}}{\mathrm{d}t} N(t) = \int_{a_{\text{min}}}^{a_{\text{max}}} \mathrm{d}a \; \partial_t N(a,t) = - \int_{a_{\text{min}}}^{a_{\text{max}}} \mathrm{d}a \; \partial_a N(a,t) = N(a_{\text{min}},t) - N(a_{\text{max}},t)

(We just used the definition of N(t) and the equations above, where all the reaction terms cancel on the right-hand side).
If we want this to be zero, then not \partial_a N(a,t) has to vanish on the boundaries, but N(a,t) itself.

However, you will still see that mass will leave the system with these (“absorbing”) boundary conditions, N(a_{\text{max}},t) = 0, but that’s not a problem. The missing mass is just in the point a_{\text{max}}. We could also specify N(a_{\text{max}},t) = N_{0} - \int_{a_{\text{min}}}^{a_{\text{max}}} \mathrm{d}a \; N(a,t) to get the correct value, but that just makes the numerical solution more tricky I believe.

Not sure this is the only way of interpreting your equations, but it would make sense to me intuitively. This also comes up when studying first-passage-time problems, where the equations are somewhat similar. The picture there is that once an individual has reached a_{\text{max}} it cannot come back, so it should be “absorbed” at that boundary.

Also physically, this is what should happen. Your population doesn’t suddenly stop ageing at a certain point – instead they will age past the range that you can describe with the given model. But we can still just (manually) summarize all of these individuals and define their mass to be whatever is missing from the simulation.

[sdwfrost]

I agree with what you’re saying, but my issue is how to implement these absorbing boundary conditions in the right way within MethodOfLines.jl. The textbook definitions of this model (see here) don’t mention the upper boundary condition (as they rarely give details of the numerical simulation of the model).

I tried a simpler example based on the linked paper (alongside the equivalent ODE) that treats everything as proportions, and so skips computation of , and I’m running into other problems now. I’m sure it’s just a silly mistake, but I can’t see it…

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots

β = 0.5
γ = 0.25
S₀ = 0.99
I₀ = 0.01
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * S(a,t) * Ia(I(a,t)),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * S(a,t) * Ia(I(a,t)) - γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob_pde = MethodOfLines.discretize(pde_system, discretization, advection_scheme=WENOScheme())
sol_pde = solve(prob, Tsit5(), saveat = dt)

Smat = sol_pde.u[S(a, t)]
Imat = sol_pde.u[I(a, t)]

plot(sol_pde.t, sum(Smat,dims=1)',label="S",xlabel="Time",ylabel="Number")
plot!(sol_pde.t, sum(Imat,dims=1)',label="I")

function sir_ode(u,p,t)
    (S, I) = u
    (β, γ) = p
    dS = -β*S*I
    dI = β*S*I - γ*I
    [dS, dI]
end

prob_ode = ODEProblem(sir_ode, [S₀,I₀], (tmin,tmax), [β,γ])
sol_ode = solve(prob_ode, Tsit5(), saveat=dt)

[xtalax]

@sdwfrost

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots

β = 0.5
γ = 0.25
S₀ = 0.99
I₀ = 0.01
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * S(a,t) * Ia(I(a,t)),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * S(a,t) * Ia(I(a,t)) - γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob_pde = MethodOfLines.discretize(pde_system, discretization, advection_scheme=WENOScheme())
sol_pde = solve(prob, Tsit5(), saveat = dt)

Smat = sol_pde.u[S(a, t)]
Imat = sol_pde.u[I(a, t)]

plot(sol_pde.t, sum(Smat,dims=1)',label="S",xlabel="Time",ylabel="Number")
plot!(sol_pde.t, sum(Imat,dims=1)',label="I")

function sir_ode(u,p,t)
    (S, I) = u
    (β, γ) = p
    dS = -β*S*I
    dI = β*S*I - γ*I
    [dS, dI]
end

prob_ode = ODEProblem(sir_ode, [S₀,I₀], (tmin,tmax), [β,γ])
sol_ode = solve(prob_ode, Tsit5(), saveat=dt)

Are you getting errors or what are the problems you’re getting?

By the way, try updating - whole domain integrals now correctly reduce dimensionality and I fixed the problems you found with using them in bcs

[BernhardAhrens]

In your new example, the wrong problem is solved. Instead of

prob_pde = MethodOfLines.discretize(pde_system, discretization, advection_scheme=WENOScheme())
sol_pde = solve(prob, Tsit5(), saveat = dt)

it should read

prob_pde = MethodOfLines.discretize(pde_system, discretization, advection_scheme=WENOScheme())
sol_pde = solve(prob_pde, Tsit5(), saveat = dt)

[xtalax]

Also I’m not sure the advection scheme works as a kwarg for discretize, it has to be in your MOLFiniteDifference

[sdwfrost]

Here’s one with the typo fixed.

# https://courses-archive.maths.ox.ac.uk/node/view_material/52367
using DifferentialEquations
using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots

β = 0.5
γ = 0.25
S₀ = 0.99
I₀ = 0.01
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * S(a,t) * Ia(I(a,t)),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * S(a,t) * Ia(I(a,t)) - γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob_pde = MethodOfLines.discretize(pde_system, discretization);
sol_pde = solve(prob_pde, saveat = dt)

Smat = sol_pde.u[S(a, t)]
Imat = sol_pde.u[I(a, t)]

plot(sol_pde.t, sum(Smat,dims=1)',label="S",xlabel="Time",ylabel="Number")
plot!(sol_pde.t, sum(Imat,dims=1)',label="I")

function sir_ode(u,p,t)
    (S, I) = u
    (β, γ) = p
    dS = -β*S*I
    dI = β*S*I - γ*I
    [dS, dI]
end

prob_ode = ODEProblem(sir_ode, [S₀,I₀], (tmin,tmax), [β,γ])
sol_ode = solve(prob_ode, Tsit5(), saveat=dt)

The PDE gives this trajectory:

The trajectory should look closer to this:

My original 3-state model (from which the above was derived) is close to what it should be:

using ModelingToolkit
using MethodOfLines
using DomainSets
using OrdinaryDiffEq
using Plots

β = 0.05
c = 10.0
γ = 0.25
S₀ = 990.0
I₀ = 10.0
dt = 0.1
tmin = 0.0
tmax = 40.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
N(t) = Ia(S(a,t)) + Ia(I(a,t)) + Ia(R(a,t))

eqs = [Dt(S(a,t)) + Da(S(a,t)) ~ - β * c * S(a,t) * Ia(I(a,t))/N(t),
       Dt(I(a,t)) + Da(I(a,t)) ~ β * c * S(a,t) * Ia(I(a,t))/N(t) - γ*I(a,t),
       Dt(R(a,t)) + Da(R(a,t)) ~ γ*I(a,t)]

bcs = [
        S(a,0.0) ~ S₀*da/(amax-amin),
        S(0.0,t) ~ 0.0,
        I(a,0.0) ~ I₀*da/(amax-amin),
        I(0.0,t) ~ 0.0,
        R(a,0.0) ~ 0.0,
        R(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0,
        Da(I(amax,t)) ~ 0,
        Da(R(amax,t)) ~ 0
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t), R(a, t)])
discretization = MOLFiniteDifference([a=>da],t)

prob = MethodOfLines.discretize(pde_system, discretization);
sol = solve(prob, Tsit5(), saveat = dt);

Smat = sol.u[S(a, t)]
Imat = sol.u[I(a, t)]
Rmat = sol.u[R(a, t)]

plot(sol.t, sum(Smat,dims=1)',label="S",xlabel="Time",ylabel="Number")
plot!(sol.t, sum(Imat,dims=1)',label="I")
plot!(sol.t, sum(Rmat,dims=1)',label="R")

sum(Smat .+ Imat .+ Rmat, dims=1)

The difference is probably because my system is leaking mass from the upper boundary at amax, and I have to figure out how to model that in ModelingToolkit.

[Sevi]

Thanks for the reference, that was quite helpful in understanding the problem better. Here are some thoughts

• the reference doesn’t include upper boundary conditions, because they are not needed for the analytical solution – because the boundary is at infinity, we can get away with specifying no boundary condition
• of course we can’t solve an infinite system like this on the computer, but if it should resemble the results from the reference, then the upper boundary should be very large compared to the maximum age at which individuals die – this way, there is no individual reaching the boundary and the boundary condition doesn’t matter. In the current setup, the initial condition spans the whole space, hence there are always weird effects from the boundaries

@sdwfrost

The PDE gives this trajectory:

@sdwfrost

The trajectory should look closer to this:

It just seems like an issue of the parameters. You might try decreasing the infection rate and death rate.

@sdwfrost

The difference is probably because my system is leaking mass from the upper boundary at amax, and I have to figure out how to model that in ModelingToolkit.

Is this from the ODE or PDE ? I’m guessing ODE, because of the initial conditions you chose.

[Sevi]

I made some edits to your code (mainly changing the initial conditions – see below).

With that, we get from the full PDE solution basically the same behavior that you showed here:

@sdwfrost

My original 3-state model (from which the above was derived) is close to what it should be:

Until the time when the first individuals hit the upper boundary, the whole system evolves like the ODE would predict, then we “leak” some mass. Of course you can choose any values, the only important thing is that the system boundary is well above the region where something actually happens. In the model from your reference, that is the case as the boundary is at infinity. This is the integrated mass:

And the individual masses (“space” / age vs current time). The initial condition is just a rectangular function, which gets advected in time and performs the SIR model individually. Everything behaves as expected until the wave hits the upper boundary. There is just some “smearing out” of the wave wich is probably the numerical dispersion that @xtalax mentioned:

And here the code to reproduce it

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots

β = 0.05
γ = 0.005
S₀ = 1.0
I₀ = 0.001
dt = 0.1
tmin = 0.0
tmax = 4000.0
da = 52.0
amin = 0.0
amax = 52.0*75.0

min_age = 52.0 * 5.0
max_age = 52.0 * 40.0

theta(a) = (a > 0) ? 1 : 0
@register_symbolic theta(a)

@parameters t a
@variables S(..) I(..) R(..)
Dt = Differential(t)
Da = Differential(a)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

eqs = [
    Dt(S(a,t)) + Da(S(a,t)) ~ - β * S(a,t) * Ia(I(a,t)),
    Dt(I(a,t)) + Da(I(a,t)) ~   β * S(a,t) * Ia(I(a,t)) - γ*I(a,t),
    Dt(R(a,t)) + Da(R(a,t)) ~ γ * I(a,t),
    ]

bcs = [
        S(a,0.0) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age),
        I(a,0.0) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age),
        R(a,0.0) ~ 0.0,
        S(0.0,t) ~ 0.0,
        I(0.0,t) ~ 0.0,
        R(0.0,t) ~ 0.0,
        Da(S(amax,t)) ~ 0.0,
        Da(I(amax,t)) ~ 0.0,
        Da(R(amax,t)) ~ 0.0,
        ]

domains = [t ∈ (tmin,tmax), a ∈ (amin,amax)]

@named pde_system = PDESystem(eqs, bcs, domains, [a, t], [S(a, t), I(a, t), R(a,t)])
discretization = MOLFiniteDifference([a=>da],t)

prob_pde = discretize(pde_system, discretization)
sol_pde = solve(prob_pde, Tsit5(), saveat = dt)

integrate(component) = da .* sum(sol_pde[component], dims=1) |> transpose

plot(sol_pde[t], integrate(S(a,t)), label="S", xlabel="Time", ylabel="Number")
plot!(sol_pde[t], integrate(I(a,t)), label="I")
plot!(sol_pde[t], integrate(R(a,t)), label="R")
plot!(sol_pde[t], integrate(S(a,t)) .+ integrate(I(a,t)) .+ integrate(R(a,t)), label="I + S + R")

plot(
    heatmap(sol_pde[t], sol_pde[a], sol_pde[S(a,t)], title="S(a,t)"),
    heatmap(sol_pde[t], sol_pde[a], sol_pde[I(a,t)], title="I(a,t)"),
    heatmap(sol_pde[t], sol_pde[a], sol_pde[R(a,t)], title="R(a,t)"),
    xlabel="time passed (weeks)",
    ylabel="age (weeks)",
    layout=(3, 1),
    size=(900,1000),
)

[Sevi]

All in all, the numerical solution makes sense to me (I didn’t try the WENO scheme yet, maybe it reduces the numerical dispersion?), but I’m not sure what you want to simulate in terms of physics/biology in the end.

I hope the comments above are helpful anyway, but feel free to tell me if not

[sdwfrost]

Dear @sevimora

Thanks! This looks great - I thought all that was needed was to (in effect) increase the boundary conditions to prevent leaking of mass over the (tmin,tmax) period, but your code on how to implement that in the BCs is super useful.

Now, back to trying to understand why the reduced system doesn’t work…

[mkretzschmar]

Dear all above,
I have taken the above code for the age structured SIR system (thanks for publishing it!) and tried to extend it. I added waning immunity and an age-dependent death rate, which worked well (is set to a constant at the moment). I also added inflow into the system at the a=0 boundary and got some results that I did not understand. I am just starting with Julia, so not sure whether I coded the boundary conditions correctly, but I could not really find any documentation on it. So I thought I could ask here.
I thought I provided a constant inflow at the a=0 boundary and an initial age distribution that is concentrated in the interval [20,30]. The latter seems to work, but then in the heatmap for the total population I get the following:

I don’t understand why there is only inflow during a limited time interval, even though my boundary condition is simply a constant.

Here is the code:

using OrdinaryDiffEq
using ModelingToolkit
using MethodOfLines
using DomainSets
using Plots
using Distributions

tmin = 0.0 
tmax = 100.0 
amin = 0.0 
amax = 100.0

#Model parameters
β = 3.0    #transmission rate
γ = 1.0      # recovery rate
alpha = 0.05    #waning rate
N = 100.0     #population size
k = 10.0
lambda = 90.0

function mu(age) #mortality function
	0.01
 #   pdf(Weibull(k,lambda),age)/(1-cdf(Weibull(k, lambda),age))
end

@register_symbolic mu(age)

# Parameters, variables, and derivatives
@parameters t a 
@variables S(..) I(..) R(..) 
Dt = Differential(t) 
Da = Differential(a) 
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

#Heaviside function 
theta(a) = (a > 0) ? 1 : 0 
@register_symbolic theta(a)


# PDE and boundary conditions
eqs = [ 
Dt(S(t,a)) + Da(S(t,a)) ~ - β * S(t,a) * Ia(I(t,a))/N + alpha * R(t,a) - mu(a)*S(t,a), 
Dt(I(t,a)) + Da(I(t,a)) ~ β * S(t,a) * Ia(I(t,a))/N - γ*I(t,a) - mu(a)*I(t,a), 
Dt(R(t,a)) + Da(R(t,a)) ~ γ * I(t,a) - alpha * R(t,a) - mu(a)*R(t,a)] 

S₀ = N * 0.99 #initial number susceptible
I₀ = N * 0.01 #initial number infected
R₀ = 0.0 	 #initial number recovered
min_age = 20.0 #minimum age
max_age = 30.0 #maximum age 
bcs = [ 
	S(0.0,a) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age), 
	I(0.0,a) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age), 
	R(0.0,a) ~ 0.0, 
	S(t,0.0) ~ N/(max_age - min_age),  
	I(t,0.0) ~ 0.0, 
	R(t,0.0) ~ 0.0 
]

# Space and age domains
domains = [t ∈ Interval(tmin,tmax), a ∈ Interval(amin,amax)]	

# PDE system
@named SIRmodel = PDESystem(eqs, bcs, domains, [t, a], [S(t, a), I(t,a), R(t,a)]) 

# Method of lines discretization
da = 1.0 
discretization = MOLFiniteDifference([a=>da],t)

# Convert the PDE problem into an ODE problem
prob = discretize(SIRmodel, discretization) 

# Solve ODE problem
dt = tmax/100.0
sol_pde = solve(prob, Tsit5(), saveat = dt)  

discrete_a = sol_pde[a]
discrete_t = sol_pde[t]
sols = sol_pde[S(t, a)]
soli = sol_pde[I(t, a)]
solr = sol_pde[R(t, a)]
soln = sol_pde[S(t, a)]+sol_pde[I(t,a)]+sol_pde[R(t,a)]

#Plotting results
x = range(0, 100, length=100)
y = mu.(x)
pltmort=plot(x,y, xaxis = "Age (a)", yaxis = "Mortality rate")
display(pltmort)

n=10

plt4 = plot()
for i in 1:11
	plot!(discrete_a, soln[1+n*(i-1), :],label="t=$(discrete_t[1+n*(i-1)])", title = "Population age distribution",
	xaxis = "Age (a)", yaxis = "Number")
end

display(plt4)
savefig(plt4, "AgeDistribution.pdf")


integrate(component) = da .* sum(sol_pde[component], dims=1) |> transpose
cumulativeplot=plot(sol_pde[t], integrate(S(t,a)), label="S", xlabel="Time", ylabel="Number") 
plot!(sol_pde[t], integrate(I(t,a)), label="I") 
plot!(sol_pde[t], integrate(R(t,a)), label="R") 
plot!(sol_pde[t], integrate(S(t,a)) .+ integrate(I(t,a)) .+ integrate(R(t,a)), label="S + I + R",legend=:outerbottom,legendcolumns=4)

display(cumulativeplot)
savefig(cumulativeplot, "cumulativeplot.pdf")

plotall=plot( 
heatmap(sol_pde[t], sol_pde[a], sol_pde[S(t,a)], title="S(t,a)"), 
heatmap(sol_pde[t], sol_pde[a], sol_pde[I(t,a)], title="I(t,a)"), 
heatmap(sol_pde[t], sol_pde[a], sol_pde[R(t,a)], title="R(t,a)"), 
xlabel="Time", 
ylabel="Age", 
layout=(3, 1), 
size=(900,1000),
tspan=(0.0,50.0) 
) 
display(plotall)
savefig(plotall, "heatmaps.pdf")

plotpopdist=plot( 
heatmap(sol_pde[t], sol_pde[a], sol_pde[S(t,a)]+sol_pde[I(t,a)]+sol_pde[R(t,a)], title="N(t,a)"), 
xlabel="Time", 
ylabel="Age",  
size=(900,300)
) 
display(plotpopdist)
savefig(plotpopdist, "heatmapN.png")

[Sevi]

Welcome to the forum @mkretzschmar

@mkretzschmar

I don’t understand why there is only inflow during a limited time interval, even though my boundary condition is simply a constant.

Are the axes perhaps labelled the wrong way? If we assume that the x-axis represents age and the y-axis time, then I would say we see exactly what you describe:

• a plateau at age 20-30 at t=0 which is advected
• a constant inflow from the “left” (age 0) by fixing the value there

The confusion might have been that in the above code examples I used the arguments in the order S(a,t), but you used S(t,a), etc.

[mkretzschmar]

Thanks for your reply, Sevi. I have also thought about that, and it is true that I switched t and a. But in the heatmap routine first t is called and then a, so that should plot t on the x-axis and a on the y axis. To check that, I now switched the order of t and a in the plotting routine to
plotpopdist=plot(
heatmap(sol_pde[a], sol_pde[t], sol_pde[S(t,a)]+sol_pde[I(t,a)]+sol_pde[R(t,a)], title=“N(t,a)”),
xlabel=“Time”,
ylabel=“Age”,
size=(1000,300)
)
display(plotpopdist)

and strangely, I get exactly the same plot.