Dynamical Modelling

better understanding through mathematics

Dynamics is the study of change: why it happens, where it takes us to. In contrast, the study of things that don’t change is called statics, and is popular among people who want to build things they hope won’t change with time, like aqueducts and arches.

Reasoning informally about dynamics is hard. People did it for thousands of years before Newton figured out how to do it properly, and he invented a formal language--differential calculus--to help him along. Believe it or not, mathematical reasoning about dynamics is easier than informal reasoning about dynamics. It just feels harder because the tools are less familiar.

My aim today is to familiarize you with some of the basic tools, and work through an example of how to apply them to pandemics by building a formal description of the dynamics: a dynamical model.

The basic tool of this formal language is the differential equation, which expresses a relationship between values and rates of change. It’s less scary than it sounds.

It’s long, so I’ve broken it up into bits.

Differential Equations for Beginners

Here's a differential equation, written in non-differential form: at a hundred kilometres per hour it'll take two and a half hours to drive from Kingston to Toronto.

If you can understand that statement, you can understand differential equations, because all they do is related a rate of change--moving at 100 km/hr, say--to a particular outcome, like: being in Toronto two and half hours from now if you're in Kingston now.

The last bit, "if you're in Kingston now" is called a "boundary value" or "initial value": it's the value of the thing you're trying to predict at the start or edge of the process. If you don't know where you're starting from, then knowing how fast things are changing is not going to tell you where you'll be sometime later on.

The range and diversity of differential equations is large, but boundary value problems, where the initial or edge conditions are specified, make up a large and interesting class. The goal is usually to integrate the equations to find out what's happening later (or elsewhere, for problems in space rather than time). I'm concerned here only with what are called "ordinary differential equations" (ODEs) which are ones where causes in one dimension result in effects in the same dimension, and nothing else: if I stretch out a spring along the x-axis, say, it pulls back only along the x-axis.

There is a vastly more complex and difficult set of problems that involve causes in one dimension causing effects in others, like when you squish a pillow (changing its height) and it fattens out on the sides (changing its width and/or length). Even though it's pretty simple for us to see what's going on, the mathematical description rapidly gets quite exciting, and involves things called "tensors", which describe the way a cause in one dimension creates effects in others. Fortunately we don't need to worry about any of that today.

The changes we're concerned with when dealing with ODEs are almost always relative to either time and/or position. Rate of change of position with respect to time is called speed or velocity, rate of change of velocity with respect to time is called acceleration. The rate of change of acceleration with time also has a name, although it's rarely used: it's called "jerk", which kind of makes sense if you think about giving a sharp yank on something. It goes from zero acceleration to a very high acceleration and back to zero acceleration again, now at a much higher velocity.

When the change is with respect to position, the first derivative is called slope--imagine the level of a floor changing as you walk across it... that's a slope--and the second derivative (the rate of change of slope) is called curvature: imagine an arch, say. There's probably some special name for the third space derivative but I don't know what it is.

A lot of the challenges people face with mathematical descriptions of the world is that we use multiple names for basically similar things. Here’s a table to relate them, where “derivative” just means “how many rates of change we take”, so the third derivative is a rate of change (1) of a rate of change (2) of a rate of change (3).

Thing        | Changing With | Name         | Derivative
--------------------------------------------------------
Value        | Time          | Velocity     | First
Velocity     | Time          | Acceleration | Second
Acceleration | Time          | Jerk         | Third
Value        | Position      | Slope        | First
Slope        | Position      | Curvature    | Second
Curvature    | Position      | ???          | Third

It is a truism of physics that God loves second order differential equations--ones that only involve first and second derivatives--because we can describe practically the whole world with them. These are equations that involve a quantity like position, it's rate of change (the first derivative), and the rate of change of that rate of change (the second derivative). What I talk about in detail below only deals with first derivatives, but I'll likely discuss second order equations at some later date.

The "differential" part of "differential equations" is just the notation that mathematicians use to write them down. Rather than write "speed" or "acceleration" they write dx/dt and d2x/dt2, but "speed" and "acceleration" is all that they mean. Physicists use conventions for co-ordinate names: x, y and z are positions and t is a time, so dx/dt means "the rate of change along the x-axis with time" (that is, the x-axis velocity) and d2x/dt2 means "the rate of change of the x-axis velocity with time".

I'm only going to be concerned with one-dimensional problems because once you understand those everything else is book-keeping, unless you're working on something genuinely hard, like fluid mechanics, which is described by the kind of differential equation that involves tensors (which for completeness I'll mention are called "partial differential equations" or PDEs).

A Pandemic Model

When we build a dynamical model of a system we want to relate the cause of change to the rate of change. When we're talking about pandemics, there are a bunch of rates of change we're concerned with because people go through different stages or states: they start out in the "Susceptible" state, jump to the "Exposed" (but not infectious) state when they first catch it, transition to the "Infectious" state when they become able to spread the disease themselves, and eventually end up in the "Recovered" state--where they are immune to infection... for a time. They may also end up in the "Dead" state, where they are immune to infection forever. Finally, because covid immunity wanes in the matter of months, everyone in the Recovered state winds up back in the Susceptible state. Dynamical models of diseases that use these stages are called SEIRS models, and that's what I'm building here.

What we want to know is the rates of change in the number of people in each of these various states. Then, knowing those rates of change and the initial conditions (the boundary values) we can sum up the changes day by day and see how the number of people in each state increases or decreases. This is called "integrating the differential equations".

First, consider the number of newly exposed people per unit time (generally the "unit time" is a day, so I'll just say "day" from here on in.) The cause of the change is the number of infectious people around. The more infectious people there are today, the more people will be newly exposed today. So we have a proportionality:

dE'/dt = K*I

where E' is the number of people newly exposed to the disease on a given day, I is the number of infectious people on that day, and K is some konstant.

Physicists have a convention that constants are named K... "c" is reserved for the speed of light, and stands for "celerity", meaning "swift" in Latin, although the convention was introduced by Weber and Kohlrausch, who are Germans. German was the language of physics for most of the 19th century, which is probably where the K for constant comes from.

While the number of newly exposed people is interesting, what we really want to know is the number of people in the Exposed state on any given day, as well as the number who are Infectious, Recovered, and Susceptible. Each stage of the disease has both an entry rate and an exit rate, and the difference between these rates is what determines the overall rate of change, like a sink that has water coming in from a tap and running out the drain. If the tap runs faster, the sink fills. If the drain is faster, the sink empties.

Modelling these transitions between states involves a bit of book-keeping, because people don't become infectious until a few days after being exposed, and people only recover a week or so after that, and people only become susceptible again a month or so further down the road. So our equations about what is happening today have to be written in terms of the state of the world days or weeks ago. In particular, the number of newly exposed people T days ago determines what happens today, where "T" has various possible values.

To keep the model simple, I'm going to assume that all people spend exactly Te days in the exposed state, Ti' days in the infectious state, and Tr' days in the recovered state before becoming susceptible again. This requirement could be relaxed, but then it would involve another set of weighted sums over probability distributions, which are tedious and don't add much to our conceptual understanding, which is what I'm trying to communicate.

The apostrophe on the variables Ti' and Tr' are pronounced "prime", so Ti' is read "Ti-prime". Primed values are usually meant to be special, which is the case here: although they are useful for thinking about the disease, they are not the values that determine the dynamics, as we'll see below.

On any given day, the number of people in the Exposed state changes according to:

dE/dt = (number of newly exposed people today) - (number of exposed people who transition to infectious today)

which in formal terms is:

dE/dt = K*I - E'[Te] (1)

where E'[Te] is the number of people who were newly exposed Te days ago: because people are in the Exposed state for exactly Te days, the cohort who got infected Te days ago will all transition to the Infectious state today, so they will disappear from the Exposed population. The prime in this case indicates “newly exposed” rather than “total exposed” at that time.

Likewise, the number of people in the Infectious state changes accord to:

dI/dt = E'[Te] - E'[Ti] (2)

where Ti = Te + Ti'

Remember, everyone spends Ti' days in the Infectious state, so the number of people who transition to Recovered today must be equal to the number of people who were newly exposed Te+Ti' days ago, which is E'[Ti]. And the number who are entering the Infectious state must be equal to the number of people leaving the Exposed state, at least so long as we are respecting the Law of Conservation of People: conservation laws are very important to building dynamical models, because they ensure that our book-keeping is exact.

Recovered works on the same principle:

dR/dt = E'[Ti] - E'[Tr] (3)

where Tr = Te + Ti' + Tr'

Again, it is the number of people who were newly exposed Tr days ago that determine the number of people whose immunity wanes today. Conservation of People ensures that everyone is shunted along the viral tracks like the cars of a train. There is nowhere else for them to go. Only if I included death in the model, which I've not done, would there be any need for a siding.

Finally, the number of susceptible people has to go up by the number of people whose immunity has waned, which is E'[Tr] and down by the number of newly exposed people, which is S*K*I/N, so:

dS/dt = E'[Tr] - S*K*I/N (4)

And there we have it: four equations that describe a pandemic, up to the meaning of our constant K, which will nail things down entirely. [Note: I managed to leave out the factor of S/N from this equation in the original version of this post, although it’s correct in the code I’ve been using. C’est la vie.]

Notice that there are so far no very free parameters in this model: the number of days someone spends Exposed, Infectious, or Recovered are all pretty well known. This is good. The fewer free parameters we have, and the easier the parameters are to measure and validate, the better. The uncertainty on Tr is high, but it is still known to be a month or so, not years or days.

But what is K?

It relates the newly exposed to the number of infectious people. For a simple model like this there are only two factors that come into it: the rate of sharing air between people, and the probability of infection given two people share air:

K = n*p

where n = number of people an infectious person shares air with in a day, and p = the probability of infecting someone you're sharing air with.

One of the big goals of dynamical modelling is to make it easier to reason about the reality the model describes. This model does that. The times are all measurable, and properties of the disease itself. They are determined by the viral progression through the body.

The constant p is a combination of how infectious the virus is, and the average conditions under which people share air. The value of p can be raised or lowered by ventilation, filtration, masking, and so on. For any given encounter, p is probably a complicated function of air changes per hour, HEPA filtration, far-UV exposure, quality of masks or respirators worn, and so on. We average over all of that and just assign it a single number, although we can also model mitgations by asking what effect cutting p by a factor of two or three would be.

The constant n--how many people breathe the same unfiltered air as an infectious person--is purely a function of population structure and population dynamics: how people live and how they move around. It is the only place in all of this that population dynamics comes into the picture at all, and as we'll see the role it plays is remarkably minor.

Conclusion

The process of developing a set of dynamical equations like this is usually described as "writing them down", although it will hopefully be clear there's a good deal of thinking that goes into doing so. But we don't "derive" them in the formal sense of mathematical deduction. We figure them out, and write them down. We're scientists, not mathematicians.

Having written down these equations that describe a pandemic in simple terms, the next trick is to solve them. I may say something about how we do that next week, or may jump straight in to what the solutions tell us about the state of the world.

Like any model, like any word, like any thought, this is all an approximation to reality. It does not follow from this that thinking or speaking or modelling is useless or pointless. It means we have to think or speak or model with care and understanding.

If we do that--if we interpret our models with care and understanding--it turns out we can learn a great deal about the world that they describe.