### Archive

Posts Tagged ‘Differential Equations’

## On Lanchester’s Differential Equations and WWII: Modeling the Iwo Jima Battle

August 18th, 2010 1 comment

Almost at the end of World War II, a battle to the death between US and Japanese forces took place off Japan in the island of Iwo Jima.  We can definitely apply a Lanchesterian analysis to this encounter; in fact, it has already been done by Martin Braun in Differential Equations and their Applications, 2nd ed. (New York: Springer and Verlag, 1975).  If $x$ represents the number of US soldiers and $y$ represents the number of Japanese troops, Braun estimates that $a = 0.05$ and $b = 0.01$ approximately.  Assuming then that the initial strength of the US forces was hovering around 54,000, and that of the Japanese was 21,500, we'd be interested in knowing the trajectory of the battle and see who wins.  I've used the method described in my previous post to obtain estimates for the proportion of  soldiers alive in $X = US$ and in $Y = Japan$, namely the Markovization or discretization (or Eulerization) method, and this is what we come up with.  We assume the battle is fought to the death and there are no reinforcements from either side:

$L = \left[ \begin{array}{ccc} 1 - 0.05 \cdot \frac{p_{alive,y}}{p_{alive,x}} & 0 & 0.05 \cdot \frac{p_{alive,y}}{p_{alive,x}} \\ 0 & 1 - 0.01 \cdot \frac{p_{alive,x}}{p_{alive,y}} & 0.01 \cdot \frac{p_{alive,x}}{p_{alive,y}} \\ 0 & 0 & 1 \end{array} \right]$

The attrition constants are within the range so that all entries are less than one as required by a Markov transition matrix provided the excess ratio of one army population to another is not too large (also, each row sums to 1).  The matrix is time-dependent, however, because each time the proportion of alives in any army changes according to the strength of the other army (proportional to both numerical superiority and technology, as can be seen in the matrix above, thus, from the Markov chain vantage point, the probability of death at each time step changes), as modeled by Lanchester's continuous differential equations.  Here are the proportion trajectories I came up with, computing at each time step the above matrix and pulling the original vector through:

Iwo Jima battle

It's pretty clear that, despite having the technology advantage, the Japanese lose because of numerical inferiority after about 59 time-steps.  The approximate proportion of US soldiers at the end of the battle is 0.35 (of the total soldier population), or about 26,276 soldiers, a HUGE decline from the initial 54,000 (roughly 0.72 of the total soldier population).  The Japanese army essentially obliterated half the soldier population of the US, a testament to their ferocity in battle.

My Calculus book (Deborah Hughes-Hallet, Andrew M. Gleason, et al., New York, John Wiley & Sons, Inc., 1994, p.551) suggests that the US in fact would have had 19,000 soldiers that could have stepped in as reinforcements: this does not matter really if the reinforcements come at the end of the battle, the battle has already been won with 54,000 soldiers, but it does alter the outcome significantly if they are added to the original soldiers at the beginning of battle (namely the US has 73,000 troops, or 0.78 of the total soldier population): the Japanese army would lose in roughly 34 time-steps (to the death), and the number of US soldiers left standing would be about 55,420, or 0.59 of the total soldier population: the Japanese obliterate only about 25% of the US troops.

In my next post, I want to discuss Dan's (from Giant Battling Robots) pdf document that he provided, essentially looking at a different stochastic model related to warfare and battle.  It's pretty neat, although it counts each potential outcome as a Markovian state (and each state is linked "downward" by a transition probability so that a given probability governs the death of one soldier in either army).  Thanks Dan!

## On Stochastic Processes

August 11th, 2010 1 comment

Why can't I shake the feeling that a stochastic process really... isn't? In previous posts we have been able to frame a deterministic system in terms of its Markov transform matrix, and going backward really doesn't seem like a problem. So what's going on, I can't help but wonder?  Are deterministic systems a subclass of random processes? Are they isomorphic spaces, or convergent at the limit as we take smaller pieces of time?  Was Einstein right, that it all converges to a determined state, and "God does not roll dice?"  What do you think?

## On Lanchester's Differential Equations and their Transform into a Markov Transition Matrix (or, the Markovization of the Lanchester Equations)

I have been very curious as to whether the Markovization of differential equations is applicable to all systems, linear or not.  In the case of the SIR differential equations, Markovization was possible even when the system is both non-linear and coupled.  I think that perhaps there is an algorithm to frame the differential equations thusly... I haven't fleshed it out completely (other than by using quite a bit of intuition, as in the SIR case).  Markovization, I think, is distinct from putting differential equations in state-space form; the vector pulled through the state-space matrix represents the trajectories of a function and its derivatives, not the trajectories of a set of functions related in that their sum is one: probabilities or proportions.  Who knows, though, a happy mix may be possible.

I was examining a Calculus book, specifically by Deborah Hughes-Hallet, Andrew M. Gleason, et al. (New York, John Wiley & Sons, Inc., 1994), and came across this very interesting problem that models how armies or soldiers fight it out... and tries to predict who wins.  I hope you don't judge: even after having read the most complicated books in mathematics, I honestly believe that revisiting more basic books can bring about insights that one didn't think of before.  Maybe even ground the complicated knowledge one acquired, or refine it.  Such has been my own experience... so I think it's a useful exercise, for all of those who have resisted "going back to the basics." If you are one, perhaps you are missing the opportunity to bridge the gaps you may (or not) still have!  :-\

The problem on page 550 says that, if $x$ represents the number of soldiers from army $X$ and $y$ represents the number of soldiers from army $Y$, the rate at which soldiers in one army are injured or die can be thought of as proportional to the number of "firepower" the other army can concentrate on them.  The model was first thought by F.W. Lanchester:

$\frac{dx}{dt} = -ay$

$\frac{dy}{dt} = -bx$ with $a, b > 0$.

In order to Markovize Lanchester's differential equations, we need to flesh out what the states are: clearly, you can be alive and belong to army $X$, or army $Y$, or be dead (from a different vantage point, we could have "Dead from $X$" and "Dead from $Y$" if we wanted to count each individually).  The Markov transition matrix should look something like this, then:

$L = \left[ \begin{array}{ccc} something & 0 & something \\ 0 & something & something \\ 0 & 0 & 1 \end{array} \right]$

where $L_{1,1}$ represents the probability of being alive in army $X$ (if you are from army $X$) after a time step (or the proportion of alive soldiers in $X$ after a time step), $L_{1,2}$ represents the probability of "defecting" to army $Y$ if you are in army $X$ (I assume such doesn't happen), $L_{1,3}$ represents the probability of dying if you are in army $X$.  Then $L_{2,1}$ represents the probability of defecting to army $X$ if you are in army $Y$ (I also assume this doesn't happen), $L_{2,2}$ represents the probability of being alive in army $Y$, if you are from army $Y$, and $L_{2,3}$ is the probability of dying if you are in army $Y$.  The third row of $L$ simply represents the fact that one cannot come back from the dead once dead:  there are no zombies (such may be in fact possible in videogames, for example, hence why I make a point of saying this).

The next order of business is to note that the total population does not change (alives and deads).  Thus we would have that:

$p_{alive,x} = x/N$

$p_{alive,y} = y/N$, and

$p_{alive,x} + p_{alive,y} + p_{dead} = 1$.  Rewriting Lanchester's equations in terms of proportions, we would have that:

$\frac{d p_{alive,x}}{dt} = -a p_{alive,y}$

$\frac{d p_{alive,y}}{dt} = -b p_{alive,x}$

To be complete, we must model what the rate of change of the dead is, which we can figure out by differentiating in time the equation $p_{alive,x} + p_{alive,y} + p_{dead} = 1$.  In effect, what I'm saying is that, Lanchester's equations are incomplete - they are missing the key information that:

$\frac{d p_{dead}}{dt} = a p_{alive,y} + b p_{alive,x}$

This equation basically states that the dead increases as the rate at which soldiers in army $X$ are put down and the rate at which soldiers in army $Y$ are put down.

Euler's method suggests that to write the approximation to these differential equations, we need think of the derivative not at the limit, but taking tiny time steps.  So in the following equation

$lim_{\Delta t \rightarrow \infty} \frac{\Delta p_{alive,x}}{\Delta t} = -a p_{alive,y}$

let's forget about the limit for a second, and solve thus:

$\Delta p_{alive,x} \approx -a p_{alive,y} \Delta t$

and

$p_{alive,x} (t+\Delta t) \approx p_{alive,x} (t) - a p_{alive,y} \Delta t$.  Writing all differential equations thus:

$p_{alive,y} (t+\Delta t) \approx p_{alive,y} (t) -b p_{alive,x} \Delta t$ and

$p_{dead} (t+\Delta t) \approx p_{dead} (t) + a p_{alive,y} \Delta t + b p_{alive,x} \Delta t$

Euler's method gives us just what we need for our transition matrix construct.  If the approximations are in fact equalities, from the point-of-view of Markov chain theory, the left-hand side is the next-step vector, and the right-hand side explains how this vector would be obtained by a transformation on the at-step vector.  In terms of Markov dynamics:

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

Playing a bit with the matrix $L$, it should look like:

$L = \left[ \begin{array}{ccc} 1 - a \cdot \frac{p_{alive,y}}{p_{alive,x}} \Delta t & 0 & a \cdot \frac{p_{alive,y}}{p_{alive,x}} \Delta t \\ 0 & 1 - b \cdot \frac{p_{alive,x}}{p_{alive,y}} \Delta t & b \cdot \frac{p_{alive,x}}{p_{alive,y}} \Delta t \\ 0 & 0 & 1 \end{array} \right]$

Letting $\Delta t$ be unity,

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

with

$L = \left[ \begin{array}{ccc} 1 - a \cdot \frac{p_{alive,y}}{p_{alive,x}} & 0 & a \cdot \frac{p_{alive,y}}{p_{alive,x}} \\ 0 & 1 - b \cdot \frac{p_{alive,x}}{p_{alive,y}} & b \cdot \frac{p_{alive,x}}{p_{alive,y}} \\ 0 & 0 & 1 \end{array} \right]$

and, I think, we have successfully Markovized the differential equations, in effect discretizing them at unit time steps.  We should be cautious though, because now both $a, b$ have additional restrictions in order that the entries of $L$ remain below or equal to 1 (and specifically, the sum of each row must be equal to one).

Euler's method in effect is the link here between the continuous differential equations and the Markov transition matrix.  In the SIR discretization of my previous post, I went backwards, showing how we can write continuous differential equations from the Markov transition matrix.

Several interesting things to note about the matrix $L$: the probability of death, if one is in army $X$ is proportional to the attrition $a$ and by how many soldiers army $Y$ exceeds $X$ at the time.  A similar argument works for the death probability of soldiers in $Y$, with attrition $b$ and excedent of $X$ to $Y$.

As with the SIR equations, since $L$ is really changing every time step, it is $L(t)$.  So, from initial conditions, the amount of soldiers alive in army $X$, army $Y$, and dead at time step $n$ can be calculated using the discretized matrix by pulling the original vector through, as:

$\textbf{p(t +} n \Delta t \textbf{)} = \textbf{p(t)} \cdot L(t) \cdot L(t + \Delta t) \cdot \ldots \cdot L(t + (n-1) \Delta t)$

In my next post, I want to model, as in the Calculus book I mentioned, the Iwo Jima encounter during World War II, the fight between Japan and the US, and solving the trajectory numerically using the discretization above.

An additional note, from examining our assumptions of matrix $L$, we realize now simple ways in which we could extend Lanchester's differential equations, by filling in "probabilities that a soldier will defect to the other army."  We have indeed now enriched the model in the sense that we can now take into account other complicated dynamics (the least of which is the "creation of zombies").  I think this is very useful in any model, because it allows us to contemplate transitions into states we wouldn't have considered otherwise, and, more powerfully, allows us to transform our mindset to the realm of differential equations, where we can have DE-related insights, and back again, where we can have Markovian insights.

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

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)$

$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.

## On the Swine Flu

So I've been MIA for the last month-and-a-half because of several exciting things that have happened to me, and I apologize for the hiatus... but it may be prolonged for a while still.  However, in the meantime, since a lot of my Mexican colleagues and students have been asking me about the swine flu, in somewhat of a panic, I thought I would try to put some numbers into some well-known SIR equations and see what insights can be obtained regarding the dynamics and proportions of the purported epidemic.  I believe this kind of analysis is routine in CDC in the US and in Secretaria de Salud in Mexico, although I'm not too sure how mathematically equipped the latter institution is in my home country. In much of what follows, I have to do some hand-waving because of the unavailability of accurate information in the web-o-sphere, or my inability to access true statistics.  Much of what I know I got from CNN and NYTIMES reports of today.  I find such fuzzy math unacceptable, personally, except to derive a notion of the magnitude of a problem, and so it would be a mistake to take the calculations as set fact or hard evidence.  Caveat in mind, I don't derive the model (it's up to the reader to find several excellent sites that might explain it), but instead delve to account for my assumptions and my results.

The SIR model (Susceptibles-Infecteds-Recovereds) uses coupled differential equations to analyze the progress of a disease in a closed population.  I focus on the population of Mexico City, assuming that the count of infected individuals nationally can be mostly found there.  Thus, out of Mexico City's 20 million people, most of the suspected 1600 currently-infecteds are in that city (all for simplicity).  Recovered individuals are the sum of those dead (149) plus those that did not die.  The coupled time-differential equations are:

$\frac{dS}{dt} = -a S I$

$\frac{dR}{dt} = b I$

$\frac{dI}{dt} = a S I - b I$

which basically says that the change in susceptible individuals is proportional to the amount infected and the amount currently susceptible, that the change in recovereds (removed from infection) is proportional to those who are already infected, and that the change in infecteds depends is the rate at which susceptibles get sick minus the rate at which infecteds get removed from the population.  The proportionality constants can be calculated with some cleverness (though not necessarily accuracy), like this:

$a = \frac{-\frac{dS}{dt}}{S I}$

and so we must guesstimate $\frac{dS}{dt}$ as well as $S$ and $I$.  If there are currently 1600 individuals that are infected, linearly, 400 have been infected per day since this became news four days ago.  So my guesstimate is that the current rate $\frac{dS}{dt} = -400$ individuals per day. The number of susceptibles is the population of Mexico City, so all 20 million, minus infecteds, about (I'm thinking naturally immune individuals are so few that the population number doesn't change much).  Finally, the number of current infecteds (Mexico-wide? focused in Mexico City) is 1600 according to the NYTIMES article I've been linking to.  This gives a value of $a = 1.25 \times 10^{-8}$.  Calculating $b$ is a bit trickier.  The datum says that 149 people have died from the swine flu, but I don't think all people infected with the swine flu die.  I'm going to guesstimate that approximately 20% of the infecteds either die or recover (since apparently 10% of them die) in a day.  There's no reason for this except my hunch. So $b \approx .2$.

By the chain rule, $\frac{dI}{dt} = \frac{dI}{dS} \cdot \frac{dS}{dt}$, and so $\frac{dI}{dS} = \frac{\frac{dI}{dt}}{\frac{dS}{dt}}$.  With $I$ not zero, this means

$\frac{dI}{dS} = \frac{1.25 \times 10^{-8} S I - 0.2 I}{-1.25 \times 10^{-8} S I} = -1 + \frac{16,000,000}{S}$.  The partial is zero at $S_* = 16,000,000$, and has initial conditions $S_0 = 19,998,400$ and $I_0 = 1600$ (twenty million total).

The value $S_*$ is called the threshold value, and it is less than the initial condition $S_0 \approx 20,000,000$.  This suggests an epidemic in fact occurs.

Luckily, $\frac{dI}{dS}$ can be solved in closed form, as

$I = -S + 1.6 \times 10^{7} ln(S) + C$.

The value of $C$ is of course determined by the initial conditions, as $C = I_0 + S_0 - 1.6 \times 10^7 ln(S_0) \approx -2.489 \times 10^8$. Then

$I = -S + 1.6 \times 10^7 ln (S) - 2.489 \times 10^8$.

With this in mind, the maximum number of infecteds at a time occurs at $S_* = 16,000,000$ and is

$I_* \approx 500,000$,

or about half a million people, equivalent to about 2.5% of Mexico City's population.

If there is enough interest, I may calculate the time dynamics (how long the epidemic lasts, etc.) with numerical methods (as by Euler's method), unfortunately by hand since access to fast computers and cool software is limited to me at present.

------

UPDATE May 5, 2009.

So it appears that the foundational numbers above were vastly overstated (since Mexico hadn't confirmed the particular strains of the alleged infecteds due to under-equipment): from the number of actual infecteds to the actual number of deaths related to the illness.  It now appears that the progression of the swine flu is a lot slower, and also that our derived coefficients are vastly different than originally thought.  Still, a happy exercise using the SIR equations.  I may yet post a new derivation that reflects reality more truly.

Categories: