[stem-dev] Handling more than a single circulating disease

As a first step towards handling a true multi-serotype scenario in STEM, we need to correctly handle two or more isolated diseases circulating among the population in the same scenario. While it is possible to add more than one disease in a scenario today, there's some problems:

1. Since birth/death rates are defined inside each disease, they can be different so the models would disagree on how many people are around at a given time.
2. The disease death rate for one of the diseases would not be taken into account in the other disease, so it would tend to overestimate the number of people availabe.

Solving this is not straightforward at all, here's some thoughts.

Tackling the first problem, the diseases are not truly isolated since they should work off the same assumptions about birth and death rate, and one disease needs to know how many people the other disease killed off. So one idea would be to pull out the birth/death rate into a separate decorator that models the demographic changes for a scenario. This "DemographicModel" (I'm open to better name suggestions, it's just a decorator) simply determine how many individuals are born and die at a given time (for any species, humans, mosquitos etc.). This is not necessarily just a fixed birth/death rate. For instance, I can envision writing a special DemographicModel for mosquitos increasing birth rate during rainy conditions etc. The DemographicModel will also consult all disease models in the scenario at each step to see how many people were killed by diseases to make sure population count is correct.

In this sense, I think a DemographicModel is similar to an infector (or innoculator). We would be able to create one using the wizard, specify parameters (birth/death rates) and a region. A top level region would propagate down to lower level regions just like for innoculators today. If we can get global demographic data for various regions we can even add a DemographicModel to the STEM library for a bunch of regions.

As for the second problem, I propose we implement (change) the API's:

A DiseaseModelLabelValue would only have these variables:

1. S, I, R etc. depending on the type of the disease
2. DiseaseDeaths. Each disease model is responsble for determining how many people die from the disease
3. Incidence How many people become infected in this time period. We still need to keep this around since it's an important epi parameter

A DemographicLabelValue has these variables:

1. Population. Current total size of the population
2, OriginalPopulation. Needed so we can reset population back to the original value when a simulation is restarted
3. Births. How many people are born this time period
4. Deaths. How many people die (for any reason) this time period

A DiseaseModel has only two methods it needs to implement:

1. calculateDelta(STEMTime time, long timeDelta, EList<DiseaseModelLabel>list). When this method returns, the delta value of each of the passed in labels must have been set so that it can be applied and is ready to advance the state of the disease from the current state to the next. At this stage, this method includes lots of things like it does today like transportation/mixing etc. It does NOT handle births or deaths (from other than deaths by the disese).
2. applyBirthDeaths(EList<DemographicLabelValue>birthDeaths, EList<DiseaseModelLabelLabel>list). This is the method that will add the births and deaths from the passed in list in the first argument to the delta value for each label in the second argument. Only labels with overlapping nodes are changed. Important With the algorithm below the deaths will include disease deaths from this disease model so it will need to be substracted before applying.

A DemographicModel has two methods to implement:

1. calculateDelta(STEMTime, long timeDelta, EList<DemographicLabel>list). Similar to the disease model, but the delta value calculated is the birth and death given the rates associated with the demographic model.
2. applyDeaths(EList<DiseaseModelLabel>list, EList<DemographicLabelValue>birthDeaths) Additional deaths are applied to the delta value here from the passed in list of disease model labels. Only labels with overlapping nodes are changed. This is to adjust for deaths caused by diseases.

Assuming we have two disease models DM1 and DM2, and one demographic decorator D, the process of figuring out the change (delta) for DM1, DM2 and D for one time step would be something like:

DM1.calculateDelta(time, timeDelta, DM1List); // DM1List is every label updated by DM1
DM2.calculateDelta(time, timeDelta, DM2List);
// Now the deltas has been set for every label for both decorators, but the birth/deaths haven't been applied yet
D.calculateDelta(time, timeDelta, DList); // DList is every label updated by D
// Now we need to correct for for the disease deaths for each decorator
D.applyDeaths(DM1List, DList);
D.applyDeaths(DM2List, DList);
// And finally adjust the birth/deaths for each disease model
DM1.applyBirthDeaths(DList, DM1List);
DM2.applyBirthDeaths(DList, DM2List);

The error in the time step (for the integration solution) would now be considered for the delta calculated for each label by DM1, DM2 and D, and if any exceed the tolerance we'd reduce the step size etc.

Still trying to think if there's a more elegant solution to this, but I think the method above would work...

/ Stefan

Stefan Edlund
Public Health and Computer Science Research