### Archive

Archive for August 6th, 2010

## On the SIR Equations and Their Transform into a Markov Transition Matrix (or, the Markovization of the SIR Model)

August 6th, 2010 No comments

It's been a while since I've had the time to write here, mostly for personal reasons.  However this period has not been without very exciting "discoveries," at least in my mathematics understanding.

Ever since last year in April with the Swine flu scare I attempted, at the very breakout, to model the progression of the disease directly using the differential equations and data provided by online articles and publications.  I'm not too sure if the raw data at the beginning of the scare was refined enough or accurate enough, but it was definitely enough to estimate some numbers and obtain a rough idea of the magnitude of the problem.  I suggested that, if anyone was interested, we could attempt to use Euler's method to complete a time evolution of the disease based on the coupled differential equations the model suggests.

Having read the excellent SIR explanation from Duke University, which differed slightly from another Calculus book in my possession in that it suggested using the proportions of susceptibles, infecteds, and recovereds as the basis for the calculations (than actual raw numbers), I began to think that, indeed, the SIR coupled-differential equation formulation was suspiciously very much like a Markovian formulation, and that, perhaps, there was an interesting link between the two.  If we think about it this way, a susceptible can (only) transition to being infected or stay as a susceptible, an infected can only transition to recovery or stay as an infected, and, once recovered, one can stay such.  In other words, S, I, and R are in effect Markovian "states."

Imagine for a moment that indeed we could somehow create an equivalence between a state-transition matrix and the actual differential equations, which, apparently we can, as I will explain shortly.  I couldn't at this point help but remember the Schrodinger vs. Heisenberg debate, in which Schrodinger's differential-equation formulation of quantum mechanics (which I know nothing about, by the way) was shown equivalent to Heisenberg's discrete-matrix approach (matrix mechanics).  At this point I began to wonder: can a system described by continuous differential equations be (always) transformed into matrix-form?  I'm not sure the answer is "yes" always, but I think I have succeeded transforming the SIR equations into matrix form, through Markovian reasoning.  Not that this is the way that matrix mechanics and wave mechanics were shown to be equivalent, that is definitely not my claim, nor do I know how they were shown equivalent either.

I started out reasoning that perhaps I could model the dynamics of susceptibles, infecteds, and recovereds by building a transition matrix as follows:

$A = \left[ \begin{array}{ccc} (1-a) & a & 0 \\ 0 & (1-b) & b \\ 0 & 0 & 1 \end{array} \right]$

where $a$ represents the probability (or, from another POV, the proportion) of susceptibles becoming infecteds, and $b$ represents the probability (proportion) of infecteds becoming recovereds.  From the Markovian point-of-view, the entry $A_{1,1}$ is the probability/proportion of susceptibles who stay susceptibles (at the next step), $A_{1,2}$ the proportion of susceptibles who become infecteds (at the next step), and, naturally, $A_{1,3}$ the proportion of susceptibles who, at the next step, become recovereds: zero.  Then $A_{2,1}$ is the proportion of infecteds who become susceptibles (zero), $A_{2,2}$ is the proportion of infecteds that stay infected in the next time period, and $A_{2,3}$ is the proportion of  infecteds recovering.  Finally, $A_{3,1}$ and $A_{3,2}$ represent the proportion of recovereds that become susceptible or infected (zero) respectively, and $A_{3,3}$ represents the proportion of recovereds who stay recovered.  In other words, the transition matrix represents a birth process of susceptibles into infecteds and infecteds into recovereds.

Following the typical Markov chain understanding, let the vector

$\textbf{p(0)} = \left[ \begin{array}{ccc} p_s(0) & p_i(0) & p_r(0) \end{array} \right]$

represent the initial proportions of susceptibles, infecteds, and recovereds.  Then the one-step forward proportion vector, $\textbf{p(1)}$, is given by

$\textbf{p(1)} = \textbf{p(0)} \cdot A$

or, following the argument through,

$\textbf{p(t+1)} = \textbf{p(t)} \cdot A$

where we necessarily assume that $t \in \mathbb{Z_+} \cup \{0\}$, or that time is discrete, which means we'd need or want to know or can calculate the proportion vector at discrete time intervals. Such is not necessarily the case (we may want to know the proportion vector at times on the real continuum), but we can resolve this objection (presently), by taking a hint from Queueing Theory.

At any rate, by writing the equations explicitly, we have:

$p_s(t+1) = p_s(t) \cdot (1-a) = p_s(t) - a \cdot p_s(t)$

$p_i(t+1) = a \cdot p_s(t) + (1-b) \cdot p_i(t) = a \cdot p_s(t) + p_i(t) - b \cdot p_i(t)$

$p_r(t+1) = b \cdot p_i(t) + p_r(t)$

Here's the amazing thing, we can write these as difference equations that are VERY much like the SIR differential equations (a discretized version I guess).  This would have been a first step to obtain the steady state of the transition matrix (assuming the transition matrix had such: regularity, and other conditions).

$p_s(t+1) - p_s(t) = - a \cdot p_s(t)$

$p_i(t+1) - p_i(t) = a \cdot p_s(t) - b \cdot p_i(t)$

$p_r(t+1) - p_r(t) = b \cdot p_i(t)$

Notice then the similarity (from Duke's website look that the SIR differential equations are formulated in terms of proportions, and from my previous post take notice of notation so that the analogy is clearer):

$\frac{ds}{dt} = -a s(t) i(t)$

$\frac{di}{dt} = a s(t) i(t) - b i(t)$

$\frac{dr}{dt} = b i(t)$

They are not exactly the same, of course: one is continuous, the other is not, and then the rate of change in the SIR equations of susceptibles also depends on the proportion of infecteds at the time, which in turn affects the rate of change of the infecteds.

Before taking on these objections, setting the difference equations equal to zero, in the usual Markov-chain context, gives the steady-state probability or proportion vector.  Thus, in the very long-run, it must be true that

$0 = - a \cdot p_s(\infty)$ and, with $a$ nonzero, $p_s(\infty) = 0$

$0 = a \cdot p_s(\infty) - b \cdot p_i(\infty) \Rightarrow p_i(\infty) = 0$ and

$0 = b \cdot p_i(t) \Rightarrow p_i(\infty) = 0$. This last bit is redundant, but the proportion $p_r(\infty)$ must be one (because the sum of all proportions must be one).

Now let's take on the objections: the SIR differential equations suggests that how susceptibles change in time is proportional to both the susceptibles and the infecteds at the time; this furthermore affects the rate of change of the infecteds.  Taking a clue, we can manipulate the "transition matrix" to incorporate this term.  Thus:

$B = \left[ \begin{array}{ccc} (1 - a \cdot p_i(t)) & a \cdot p_i(t) & 0 \\ 0 & (1-b) & b \\ 0 & 0 & 1 \end{array} \right]$

Notice that every row sums to 1 (provided appropriate conditions on $a, b$), as it should if this is a Markov chain transition matrix! This strongly suggests that we can follow a Markov chain treatment, and everything Markovian applies (particularly steady-state proportions). And so, the next-step equation is given by the usual Markov context as:

$\textbf{p(t+1)} = \textbf{p(t)} \cdot B$

or, explicitly, as

$p_s(t+1) = p_s(t) \cdot (1 - a \cdot p_i(t)) = p_s(t) - a \cdot p_i(t) \cdot p_s(t)$

$p_i(t+1) = a \cdot p_s(t) \cdot p_i(t) + (1-b) \cdot p_i(t) = a \cdot p_s(t) \cdot p_i(t) + p_i(t) - b \cdot p_i(t)$

$p_r(t+1) = b \cdot p_i(t) + p_r(t)$

with difference equations (now analogous to SIR differential equations)

$p_s(t+1) - p_s(t) = - a \cdot p_i(t) \cdot p_s(t)$

$p_i(t+1) - p_i(t) = a \cdot p_s(t) \cdot p_i(t) - b \cdot p_i(t)$

$p_r(t+1) - p_r(t) = b \cdot p_i(t)$

and steady state proportions

$0 = - a \cdot p_i(\infty) \cdot p_s(\infty)$

$0 = a \cdot p_s(\infty) \cdot p_i(\infty) - b \cdot p_i(\infty)$

$0 = b \cdot p_i(\infty)$

From the third equation it follows $p_i(\infty) = 0$, which now implies that both $p_s(\infty)$ and $p_r(\infty)$ can be anything at steady-state (as long, of course, that $p_s(\infty) + p_r(\infty) = 1$).

Examining the Duke University link to a SIR progression graph, this vector tendency can be seen pretty clearly at large $t$. Notice, for example, how $i(t)$ goes to zero while both $s(t), r(t)$ are anything.  We have in effect, proved that under the SIR model, the proportion of infecteds is asymptotic to 0. (Neat-o).

On to the continuity question: using the difference equations we are in effect assuming that the state vector $\textbf{p(t)}$ is known at discrete times.  In Queueing Theory when discussing birth-death processes, it is apparently conventional (in the standard treatment) to think that the probability in a given unit time interval tends to accumulate uniformly. Thus, we can rewrite the transition matrix as:

$C = \left[ \begin{array}{ccc} (1 - a \cdot p_i(t) \cdot \Delta t) & a \cdot p_i(t) \cdot \Delta t & 0 \\ 0 & (1-b \cdot \Delta t) & b \cdot \Delta t \\ 0 & 0 & 1 \end{array} \right]$

We want to write out the explicit equations of

$\textbf{p(t+} \Delta t \textbf{)} = \textbf{p(t)} \cdot C$

as

$p_s(t + \Delta t) = p_s(t) - a \cdot p_i(t) \cdot p_s(t) \cdot \Delta t$

$p_i(t + \Delta t) = a \cdot p_s(t) \cdot p_i(t) \cdot \Delta t + p_i(t) - b \cdot p_i(t) \cdot \Delta t$

$p_r(t + \Delta t) = b \cdot p_i(t) \cdot \Delta t + p_r(t)$

By further algebraic manipulation,

$\frac{p_s(t + \Delta t) - p_s(t)}{\Delta t} = - a \cdot p_i(t) \cdot p_s(t)$

$\frac{p_i(t + \Delta t) - p_i(t)}{\Delta t} = a \cdot p_s(t) \cdot p_i(t) - b \cdot p_i(t)$

$\frac{p_r(t + \Delta t) - p_r(t)}{\Delta t} = b \cdot p_i(t)$

Taking the limit as $\Delta t \rightarrow 0$, the difference quotient describes the very definition of the derivative, and we obtain the continuous SIR differential equations.  We have assumed, in effect, the vector $p(t)$ is known at all times $t \in \mathbb{R}^+ \cup \{0\}$.

Perhaps an obvious thing to note, yet worth mentioning, now that we have shown the equivalence between the SIR model and the Markov transition matrix, is that, in the SIR model the proporitonality constant $a$ is now shown to be related to the probability of transitioning from the susceptible state to the infected state, and the proportionality constant $b$ is now shown to be the probability of transitioning from the infected state to the recovered state.  This, in my opinion, is an amazing insight!

Before, since the SIR differential equations could not be solved explicitly, we would rely on Euler's method to get a numerical approximation of the solution curves $s(t), i(t), r(t)$.  The present analysis suggests we need only calculate powers of matrix $C$ to obtain a numerical sample of the all three solution curves concurrently (spaced as $\Delta t$), since, in the Markov chain context, the power of the matrix represents the $n \cdot \Delta t$ -step forward evolution of the initial (or current) state vector (in the context of a time-dependent matrix $C(t)$, I actually mean by "power" consecutive matrix multiplications of $C$, as by

$p(t + n \Delta t) = p(t) \cdot C(t) \cdot C(t + \Delta t) \cdot C(t + 2 \Delta t) \cdot \ldots \cdot C(t + (n-1) \Delta t)$,

not actual powers of $C(t)$).  Higher-rate samplings (higher-resolution solution curves) can be achieved by reducing $\Delta t$.  It may be worth exploring the link between Euler's method and this sampling-method (see my next post), in terms of computational requirements, or in terms of a mathematical equivalence.

Needless to say, this method can probably be extended to "Markovizate" (Markovize?) other differential equations, such as SIRD or Lotka-Volterra.  I can only speculate as to the limits of this process, but I am certainly cognizant of two things: several, if not all, regular Markov chains can be formulated as continuous differential equations, and, several, if not all continuous differential equations can be formulated as Markov transition matrices.  I'd love to explore more such, and will probably in the future.