Reaction networks for disease dynamics

Some simulations of the effect of an asymptomatic phase

Last week, I changed a bit of my modelling lecture to discuss the notion ot $R_0$, and how can estimate it in simple models of epidemics. And then the question of “what happens if there are no symptoms” came up, which is an interesting pretense for a short simulation. In this entry, I will look at using discrete stochastic differential equations to see how our understanding of the disease dynamics changes over time. This is as good a place as any to state that what follows is purely a theoretical exercise and does not have any consequences or implications, practical imagined or otherwise, on anything.

A good-enough model for this situation is as follows. We have a susceptible population S, an asymptomatic population A, a symptomatic population Y and a recovered population R. Both A and Y can transmit the disease to S. A eventually progresses to Y, and Y eventually recovers. We assume no demography, and no mortality due to infection, so N = S+A+Y+R remains constant. This differs from a SEIR model in that both sub-component of the infectious population can transmit the disease.

We can write this model as a series of ODEs, in this case Ṡ = - β(A+Y)S, Ȧ = βAS - ρA, Ẏ = ρA - ρY, and Ṙ = ρY. Note that for the sake of simplification, we assume that both A and Y have the same transmission rate β, and the duration for the symptomatic and asymptomatic steps are the same (ρ⁻¹).

We could write this as a series of ODEs, but we won’t. Instead, we will represent this as a reaction network, using the DiffEqBiological package for Julia, which offers a really intuitive domain specific language for this type of problems, using the rate, inputs → products syntax (and more!). So let us start by re-expressin our model as as series of reactions.

First, when individuals from A and S meet, they react as a rate β to change the individual from S into an individual from A, so, the net result of this interaction is 2S. We can therefore represent this first reaction as β, S+A → 2A. Just about the same thing happens when individuals from S and Y meet, except that we keep the individual from Y, and the individual from S becomes A, which yields the rule β, S+Y → A+Y. Finally, at rate ρ, the asymptomatic cases start showing symptoms, so ρ, A → Y, and the symptomatic cases recover, so ρ, Y → R. We have defined the reaction network for our model.

Let’s write the code:

sayr = @reaction_network begin
  β, S+A → 2A
  β, S+Y → A+Y
  ρ, A → Y
  ρ,Y → R
end β ρ
(::Main.WeaveSandBox1.reaction_network) (generic function with 2 methods)

I love how expressive the syntax is! We can take this network, and convert it into a problem for DifferentialEquations to solve. As always, this requires u₀, t, and the parameters:

p = (0.02, 1.0/3.5)
u0 = [999,1,0,0]
t = (0., 40.)
timepoints = t[1]:1:t[2]

We can now transform our network into a JumpProblem (it’s all in the documentation!), and solve it:

hyperparameters = DiscreteProblem(u0, t, p)
problem = JumpProblem(hyperparameters, Direct(), sayr)
timeseries = solve(problem, SSAStepper(), saveat=timepoints)

Our model has time units of days, and we get the results for every day between 0 and 40. What does it looks like?

plot(timeseries.t, timeseries[1,:], c=:black, label="Susceptible")
plot!(timeseries.t, timeseries[2,:], c=:orange, ls=:dash, label="Asymptomatic")
plot!(timeseries.t, timeseries[3,:], c=:orange, label="Symptomatic")
plot!(timeseries.t, timeseries[4,:], c=:teal, label="Recovered")

This looks more or less like what we expected - the proportion of the population that is recoved increases steadily, and eventually the transmission stops. Now, the discussion that motivated this was, is our understanding of the epidemics changing if we miss the information about the asymptomatic transmission?

To do so, let’s look at measured prevalence over time (Y/N), and the actual prevalence over time ((A+Y)/N):

plot(timeseries.t, timeseries[3,:]./sum(u0), c=:grey, lab="Measured on symptomatic")
plot!(timeseries.t, (timeseries[2,:].+timeseries[3,:])./sum(u0), c=:black, lw=2, lab="Actual")

As we might expect, the prevalence measured on the symptomatic cases only would be much lower, and would peak much later. How much later? Well, this is a good question, and we might think about it for a minute and get an approximation, or we can run the model thousands of times and measure how many days after the actual peak we dected the peak in symptomatic cases:

peak_measures = zeros(Float64, 500)
peak_actual = similar(peak_measures)
for i in eachindex(peak_measures)
  ti = solve(problem, SSAStepper())
  prev_mes = ti[3,:]./sum(u0)
  prev_act = (ti[2,:].+ti[3,:])./sum(u0)
  peak_measures[i] = ti.t[last(findmax(prev_mes))]
  peak_actual[i] = ti.t[last(findmax(prev_act))]
Δ = peak_measures .- peak_actual

Let’s look at the average error we would make:

import Statistics

Now, we could have gotten a very reasonable estimate by thinking about the time it takes to start showing symptoms, which is ρ⁻¹ days, and using our parameters this turns out to be about 3.5 days, a result to which our result based on simulations is very close (by about 8 hours).

So, what have we learned? Not a lot on epidemics that we would not have intuited. The main point of this exercise was to play with the reaction networks for differential equations, and I must say that I am now fully convinced. For some problems, this is a very natural way of representing changes, in addition to being easy to write compared to the full differential equations.