[A series of posts explaining a paper on the mathematical modeling of the spread of antiviral resistance. Links to other posts in the series by clicking tags, "Math model series" or "Antiviral model series" under Categories, left sidebar. Preliminary post here. Table of contents at end of this post.]
In the previous posts we have walked through the Introduction and Methods sections of the paper, "Antiviral resistance and the control of pandemic influenza," by Lipsitch et al., published in PLoS Medicine. The Methods section sets out the detailed model, which is summarized in the Figure in the upper left and the equations on the lower right of page 3 of the .pdf version of the freely downloadable paper. In this penultimate post on the Methods section, we'll clean up a few details and give an overview. Then, on to Results.
Reading the figure:
To summarize, here is how to read figure 1. Start on day zero with a population of susceptible people, all in box X. On day one seed the population with a single person, infected with a strain of the virus sensitive to Tamiflu. That person is initially untreated. Thus we start the model running with X-1 people in box X and one infected person in the box YSU, an untreated index case. On each succeeding day the rule is that you move a fraction of people from box X to each of the Y boxes and the Z box and a fraction from each of the Y boxes to the Z box by the cumulated factors on each of the arrows. Thus to get from X to YST you traverse two pathways, one along the first arrow out of the X box at 3 o'clock, βSTYST + βSUYSU, then either along the top branch fp(1-epcp) or the middle brach, which has two legs, giving (1-fp) times fT(1-cT). Each of these is a fraction of the number of people in the X box on that day, so the total number of people who move is X times βSTYST + βSUYSU times fp(1-ep-cp) + X times βSTYST + βSUYSU times fT(1-cT). Looking at the YST box you see that it also loses people to the Z box at the rate of vT per day (just the rate is on the arrow label) for a total of vT times YST people per day.
The figure and the equation are the same:
The quantity (βST YST + βSU YSU) is clumsy to write over and over again, so the authors abbreviate it in equation 2 with another Greek letter, λS. You can see that abbreviation as the second last line of the equation at the bottom right of page 3. Instead of an equal sign with two bars there is one with three bars. This indicates it is a definition. A lot of mathematics just looks hard and the use of Greek letters often makes it look even harder. But what we have done so far is just bookkeeping, using the Law of Mass Action to move things around. Even equation (2) is just bookkeeping. Using the abbreviation, the total number of people added to and lost from compartment YST on the next day is:
λS[fp(1-ep-cp) + (1-fp)fT(1-cT)]X - vSTYST.
I wrote all this out so you can see that these are exactly the terms in the third line of equation 2 to the right of the equal sign.
The left of the equal sign, dYST/dt, you can consider shorthand for the amount of change in the YSU box each day, the arrow going in, minus the arrow going out. The number of people each day riding the arrow going in is just λS[fp(1-ep-cp) + (1-fp)fT(1-cT)]X and the number of people per day riding the arrow going out is vSTYST (which has a minus sign n front of it because they are leaving the box YST).
Thus the figure and equation 2 say the same thing. It is a good exercise to follow the arrows and check that each line of equation 2 (except for the last two) is just another way to write the pathways between the boxes shown in figure 1, each of which was explained in detail in the last two posts.
A little detail: equation (1):
We need to make a special comment about equation 1 and the last line of equation 2. Here the authors introduce the Greek character ξ, used to prevent tiny fractions of a resistant case from affecting the findings initially.
The authors want to model what happens if the emergence of fit and resistant viruses does occur, but only very rarely. The way the model is set up, extremely small fractions of a case would still "transmit" infection, but in real life you need at least one full case for that to happen. .001 cases (a thousandth of a case) doesn't transmit infection. Equation 1, with its Greek letter ξ, is a standard way to prevent any transmission until the cumulated fractions of resistant cases amounts to at least one full case. That is also what the last line of equation 2 says: don't let resistant viruses be transmitted until there is at least one full case of resistance in the population.
How to account for everything else - R0:
The last detail to take care of is how to handle all the other things that are going on to manage a pandemic, like social distancing, handwashing hygiene, isolation, etc. We also don't know how good an emerging pandemic strain will be in reproducing itself in the population. Since the object of the model is to examine antiviral resistance, all these other factors are lumped into a single component, the basic reproductive rate of the virus, R0. R0 expresses (on average) how many new cases a newly introduced infectious case would produce in a totally susceptible population, that is, how many new cases you would expect on day 2. Many people think of this as some kind of intrinsic property of the virus, but it is really a property of the virus, the host it infects and the environment in which they live (see our entry at The Flu Wiki on this). Lipsitch et al. consider a range of possible R0s (1.2, 1.5, 2.0, 2.2) and point out that you can consider a virus with an R0=1.5 to be either one with an R0=1.5 without intervention, or one with an R0=2.2 with an intervention like social distancing. In other words, they are lumping all the non drug interventions into the size of R0. More on this later.
We are now ready to move on to Results in the next post.
Table of contents for posts in the series:
The Introduction. What's the paper about?
Sidebar: thinking mathematically
The rule book in equation form
Effects of treatment and prophylaxis on resistance
Effects of Tamiflu use and non drug interventions
Effects of fitness costs of resistance
- Log in to post comments
Good post, but the 'regular' text is kinda tiny.
Mike: Are you using Firefox? It should be "normal" text as the equation is a .jpg and not formatted as text.
Overnight Revere it has reverted to microtext onmy computer too. Cant post either
Small quibble (again - sorry): is it the norm in the field to use d/dt - continuous change - for expressions that actually are discrete-time? That was one of the things that tripped me up when trying to decipher the paper by myself.
Janne: Good question (again). In short, yes, it is the norm. The models are usually continuous, not discrete. Of course the actual implementation in the computer is discrete, but the mathematical analysis (e.g., phase plane analysis) is usually on the continuous model. This allows bringing in all the tools of qualitative analysis of nonllinear ODEs.
There are potentially tricky things that could happen when you move to difference equations (think of the difference equation version of the logistic equation that has peridic and chaotic regimes but the the ODE version is always stable and well behaved) but these problems are usually ignored. They can also happen with ODEs, so looking for bifurcations and strange behavior is warranted. Chaotic behavior in nature isn't very easy to demonstrate, but our tools are pretty blunt in that regard.
Your problem is you know too much.
Randy: That's weird. It still looks normal on mine using either Firefox or Safari. Anyone else having that trouble?
Randy, Mike: Found an unclosed subscript tag (the subscripts are small font) and that might have been the problem. Does it look all right now?
Might have reset something on the Seedpage without knowing it when you imported the stuff above in. You are starting to answer some of the questions that I posed yesterday. I meant what I said though, the modelers if they have their act together and full body suits and masks can use the media, personal observations etc. to create a whole new modeling equation that the world could use. It could be used for everything from BF, SARS, HIV, stocks, food, Katrina and Tsunami's to produce the fine line outcomes for the individual and the group if they can stay alive if it comes. Every facet of a major disaster could be unfolded and all they have to do is input the expected, the actuals and it squirts out a real live no kidding picture. Major opportunity here and I hope someone besides me sees it.
Very neat stuff and very impressive. Jeez,
Thanks. And yes, I know well how difference equations are messier and more difficult to solve than their continuous counterparts, with ghost solutions and all; and yet, I find difference equations easier to grasp. Go figure.
But, in this case, does it mean the model in the paper is actually continuous and you just present it in terms of day-by-day steps?
Janne: Yes the model is continuous, only discretized in the numerical solution (as it must be).
When you said you were going to post this series I thought I had better read the paper beforehand. Whilst I saw the logic in defining a term to prevent fractions of resistance remaining in the equation and being multiplied until they became significant, I was confused by the paragraph covering the use of this term once the basic model was expanded to include the age group interactions matrix. I had hoped this would become clearer once you reached that section in your posts (I have copied the relevant section below). My problem is this; my reading of the paragraph is that a de novo resistant person needs to occur in a group to be counted. I would expect this to be approximately 1/6 as all you need is emergence of one person in any of the groups to get +1 in the Yr box. Had they broken the ages into on million gradations then it would a million times more difficult to achieve de novo resistance emergence in any one of the age compartments.
Have I just misread this or am I missing something else?
In generating the numerical results, we divided the population into six age groups, i ¼ 1 to , each with its own initially susceptible population size Xi(0) and force of infection kSi and kRi. The restriction on transmission of resistant strains in this model is implemented separately for each age class; hence, once the expected number of resistant infections in class i exceeds 1, these infections can be transmitted to individuals in all six age classes. Before that time, resistant infections accumulate in an age class only by de novo acquired resistance and by transmission from other age classes j for which nj . 1.
JJackson: Ah, yes. A great question. So I put it to Prof. Lipsitch, who also thought it was a great question: Here is his reply to me, which he has given me permission to quote:
This answers your question but may or may not satisfy you. Science is a process of continually trying to get the world to give up its secrets and we try all sorts of tricks to do it. At any rate, you are clearly following this closely and using good intuition.
Maybe you should start a new career as a modeler?
Thank you for following up on my question (and if you communicate with Prof. Lipsitch again please thank him for me).
I had not considered the case where each of the groups is isolated and now see (qualitatively) how the value would fall between 1/n (n being the number of age groups) and 1.
Thanks again and I thinks your pessimism in predicting this series of posts would be a failed experiment was unjustified.
I kind of thought after reading the paper that the single set of ODEs represented essentially s fully mixed situation and the six independent sets a segregated. Doing both gives a range of results, probably including reality in the middle. If that's a hack it's one we use a lot in engineering.
Yup, that seems to have gotten it Revere. The size is back to normal.