Multistrata IPMs

In the previous chapter we looked at modeling with single-stratum compartmental models. These models express a hypothesis that the pathogen affects everyone in the same way. But in reality, infectious diseases often behave differently based on characteristics of the individuals exposed to them. In these cases, single-stratum models are a simplification — which can be a useful position to explore — while more complex models might disentangle the critical individual factors, and in doing so, gain expressive power. In these cases, classic epidemiological theory says that parts of the population belong to meaningfully different “risk classes” that should be accounted for explicitly.

For example, an influenza strain might be more likely to hospitalize older patients than younger patients. Locations can have significant differences in demographic makeup. A location with a population that skews older would expect to see a higher proportion of its population wind up in the hospital when compared with another location whose population skews younger, assuming both locations are equally impacted by infection. A model that represents a single, combined population would not be able to explain these differences. Modeling these subpopulations separately could improve the model’s fidelity to reality.

Multistrata models exist to address these situations. They divide populations, stratifying them by characteristics thought to be critical to the progression of the infectious pathogen. Characteristics may include demographic factors, like age, sex at birth, or racial or ethnic identity; economic and social factors, like wealth and poverty, educational attainment, or access to healthcare facilities; behavioral risk categorization, like occupation in a high-risk industry, diet and exercise activity, or sexual preferences; or any number of factors. Anything which might alter a group’s participation in the modeled system and which is independent from strictly geographic factors, like climate.

Multi-species models are another important use-case. Disease systems often involve interactions between more than one species, yet the way the disease affects each species differs in critical ways. West Nile Virus, for instance, involves pathogen transmission between mosquitoes, birds, and humans. We can use multistrata modeling in epymorph to represent these as well.

Another advantage of multistrata modeling in epymorph is that strata can have unique movement dynamics, in addition to unique transmission and disease dynamics. This modularization allows for enhanced flexibility in simulation design, allowing for sophisticated movement behavior.

Multistrata simulations in epymorph

epymorph provides first-class support for simulations involving multiple strata. Strictly speaking, it would be possible to use epymorph’s basic functionality to implement multistrata systems; it would just be extremely tedious and inflexible. Instead we included features to handle the common heavy lifting and extended epymorph’s modular approach so that multistrata models are nearly as easy to work with as simpler models. With epymorph you can quickly assemble experiments that have heterogenous disease models between strata, different rules for movement, and specialized initialization, and can easily provide the data required for these components.

The basic concept is this: a multistrata IPM merges multiple base IPMs (one for each strata), and adds meta edges and meta parameters to form one big combined IPM.

But assembling a multistrata model is more of a “whole RUME” task, not limited to just the compartmental model of your simulation. Adapting IPMs for a multistrata simulation is typically the most involved part, so while this chapter will focus primarily on details relevant to the IPM, it will by necessity discuss the other components as well.

When designing your simulation, identify which strata you are interested in studying and which factors you want to vary between strata. Do the individuals in each stratum progress through different infection phases (IPM compartments), or is it sufficient to use the same structure but vary the parameterization? How do individuals in one stratum influence the spread of a pathogen in other strata? Do individuals move between locations according to different rules? How do you want to initialize the starting populations of the strata?

Example

Let’s set up a simulation for a hypothetical pathogen, perhaps a virus transmitted through aerosolized droplets, where younger people are significantly less likely to become infected than older people. Assume the disease affects both groups in roughly the same way after infection, though. We will use an SIRS model for this. And to simplify things, we will assume there are no differences in movement between the groups. We will need population data subdivided by age, but epymorph provides a data source for this.

Defining a multistrata RUME builder

We start by implementing a subclass of MultiStrataRUMEBuilder and defining the strata as GPM objects:

from epymorph.kit import ipm, mm, init
from epymorph.rume import GPM, MultiStrataRUMEBuilder

class ExampleRUME(MultiStrataRUMEBuilder):
    strata = [
        GPM(
            name="youth",
            ipm=ipm.SIRS(),
            mm=mm.Centroids(),
            init=init.SingleLocation(location=0, seed_size=1_000),
        ),
        GPM(
            name="adult",
            ipm=ipm.SIRS(),
            mm=mm.Centroids(),
            init=init.SingleLocation(location=0, seed_size=1_000),
        ),
    ]

The GPM packages the modules (IPM, MM, Initializer) that apply to just a single stratum, along with a name for the stratum.

If we don’t add anything more to this setup it would be like running two entirely separate simulations for the two strata. How do people from the youth stratum interact with people from the adult stratum? To answer this, we layer on top of our base IPMs a “meta” layer of additional parameters and transitions. The base IPMs define the compartments and keep their original transitions to describe how individuals interact within that stratum, and we add transitions on top of those for interactions between different strata.

Add “meta edges” to model interactions between strata.

In this example, our base models allow susceptible youths to become infected due to contact with an infected youth, and for adults to become infected due to contact with an infected adult. (These are effects within strata.) But we also want to model that adults can become infected due to contact with an infected youth, and for youths to become infected due to contact with an infected adult. (Effects between strata.) That’s four total transition edges, two in the base model and two in the meta model. And each of these should have a unique infectivity parameter — \(\beta\) — so we can tune the infection transmission rate in each of these cases independently. This diagram lays out the beta parameters for our IPM with the new ones highlighted.

A matrix describing the four beta parameters for all possible interactions among two strata

It’s conventional to annotate multistrata models using subscripts for the strata. If we use “a” for adult and “y” for youth, beta_ay (or \(\beta_{ay}\) in mathematical notation) parameterizes the transition when a susceptible adult becomes infected due to contact with an infectious youth. (Back to school season can be particularly rough for parents!)

We can also think of \(\beta\) as taking the form of a transmission matrix: with width and height equal to the number of strata. Values on the diagonal represent transmission within a stratum, while off-diagonal values represent transmission between different strata. The values we choose for this matrix may or may not be symmetrical about the diagonal. Transmission from one stratum to another may occur at a higher rate than infection with the strata reversed. So take care to specify these values correctly, keeping in mind which parameters map to which strata: who is “becoming infected” and who is “passing on the infection”. (All of our examples use the subscript ordering convention that the recipient stratum is listed first and the source stratum is listed second.)

Continuing the class definition from above:

# additional imports....
from sympy import Max

from epymorph.attribute import AttributeDef
from epymorph.compartment_model import edge
from epymorph.data_shape import Shapes


class ExampleRUME(MultiStrataRUMEBuilder):
    strata = [...]  # strata as defined above

    # Add "meta requirements" for new parameters, just like adding
    # requirements to an IPM.
    meta_requirements = [
        AttributeDef("beta_ya", float, Shapes.TxN),
        AttributeDef("beta_ay", float, Shapes.TxN),
    ]

    # Add "meta edges" for new IPM transitions due to cross-strata
    # interactions, very similar to adding edges to a regular IPM.
    # Remember, our base IPM transitions will be automatically included
    # in the full model, we don't have to redefine them.
    def meta_edges(self, symbols):
        # extract compartment symbols by strata
        S_y, I_y, R_y = symbols.strata_compartments("youth")
        S_a, I_a, R_a = symbols.strata_compartments("adult")

        # extract meta attributes
        beta_ya, beta_ay = symbols.all_meta_requirements

        # compute totals by strata
        N_y = Max(1, S_y + I_y + R_y)
        N_a = Max(1, S_a + I_a + R_a)

        return [
            # youth infection due to contact with adult
            edge(S_y, I_y, rate=beta_ya * S_y * I_a / N_a), 
            # adult infection due to contact with youth
            edge(S_a, I_a, rate=beta_ay * S_a * I_y / N_y),
        ]

We added two edges: \(S_y \to I_y\) and \(S_a \to I_a\). These edges are “parallel” (i.e., the edges share the same source and destination) with the \(S \to I\) edges in the base models. Parallel edges are maintained as-is, because it can be useful to track these as separate events. Your results will be able to distinguish between the four infection cases, although because parallel edges use the same name you have to use positional notation to access them.

IPM rate denominator \(N\)

In the \(S_y \to I_y\) transition rate specified above:

\[\frac{\beta_{ya} S_y I_a}{N_a}\]

Why is the denominator \(N_a\) instead of \(N_y\)? Why is there a denominator at all?

You will often see epidemiology papers whose transmission rate equations do not divide by population totals. Some model fractional populations, where the total population is \(1\) and compartments divide the population into fractions. Some model infection as dependent upon population density rather than expected contact frequency. Others may calibrate parameter values to assume a fixed population size.

First, epymorph’s built-in IPMs generally use frequency-dependent transmission as a convention. Second, epymorph is designed to model discrete population sizes and aims to maximize modularity. When comparing parameters between different geographies, it would be inconvenient to have to account for population differences. Rates like \(\frac{\beta S I}{N}\), define \(\beta\) so as to be scaled for normalized populations. In multistrata models this is a bit more complicated: how should we define \(N\)?

It may help to think about the value of \(\beta\) as being a placeholder for several factors. How many close contacts does someone have in a day, on average? Social distancing policies were designed to reduce this factor. What mitigations are people using when they do interact? Mask policies were designed to reduce this factor. What is the biological mechanism for infection and how potent is it? The impact of pathogen biology is subtle and can shift over time. And so on. These factors, individually, are difficult to quantify and reason about, so we combine their effects into one number. (You don’t typically come up with the value of this number simply by reasoning about it; you measure it by comparing theoretical results against real world results. This is what parameter fitting is all about.) But if we could actually quantify all of these factors individually, we might hypothetically express the transition rate like this:

\[\frac{a b c S_y I_a}{N_a}\]

Which breaks down \(\beta_{ya}\) into separate “quantifiable factors” \(a\), \(b\), and \(c\).

\(a\) could represent the biology of the pathogen: how easily will an infectious adult transmit the pathogen to a susceptible youth? \(b\) could represent behavioral factors: how closely are people using good hygiene, isolating at home when sick, covering coughs, washing hands, etc. And finally \(c\) could represent how many adults will the average youth come into contact with each day? Youths might interact with a lot of their peers (e.g., in school) and a relatively smaller circle of adults (e.g., teachers and parents).

With these definitions, the quantity \(c S_y\) is the total number of interactions between a susceptible youth and any adult. \(\frac{I_a}{N_a}\) represents the probability that a randomly selected adult is infectious. Multiply those two quantities and we have the total number of interactions between a susceptible youth and an infectious adult. Finally multiply that by the quantity \(a b\) (the pathogen’s transmissibility and relevant behavioral factors) and we get the total expected number of new infections — exactly what we’re trying to compute!

So ultimately the answer is dependent upon what we intend a quantity like \(\beta_{ya}\) to represent. We, as the authors, get to choose the definitions of our terms. You may have reasons to use a different set of conventions in specifying model rates, which is just one more reason why epymorph gives you so much flexibility when it comes to defining IPMs.

Events for the combined IPM are maintained in a specific order:

First events from each base IPM (in their definition order)
- with each stratum listed in strata definition order,
then meta events as they were defined in the RUME builder.

So for our example the events are (zero indexed):

\(S_{youth} \to I_{youth}\)
\(I_{youth} \to R_{youth}\)
\(R_{youth} \to S_{youth}\)
\(S_{adult} \to I_{adult}\)
\(I_{adult} \to R_{adult}\)
\(R_{adult} \to S_{adult}\)
\(S_{youth} \to I_{youth}\) (meta edge)
\(S_{adult} \to I_{adult}\) (meta edge)

The work is not terribly different from defining a single-stratum IPM, there are simply more compartments and new parameters involved. You can also extract the parameters from the base models, if you need to use these as well in the new meta edges. However you cannot alter the edges from the base models; if you need to do that, you should swap or redefine the base IPM(s) directly.

Cross-strata edges

It is possible to define edges which cross between strata. That is, rather than the edge rate simply being influenced by another stratum, the source and destination compartment for the edge are in different strata. This could be used, for example, to model the population aging in very long term simulations. If your IPM uses this, care should be taken that all components (like your movement model) can handle this kind of population shift.

Otherwise, and especially with multispecies simulations where cross-strata edges would perhaps not make sense, we recommend checking that the source and destination strata match for all meta edges.

Instantiating the RUME

Having defined our RUME builder, we should now use it to build our RUME. A MultiStrataRUMEBuilder is just what it says: a class which builds a RUME. It is not a RUME in and of itself.

To build a RUME, call the build() method and supply a scope, time frame, and parameter values. These are the remaining parts of a RUME which were not already defined above.

Scope and time frame are straightforward, but it may not be obvious which parameters you need to provide. You can query a RUME for this information, but you do have to construct it first. We can temporarily leave the parameters dictionary empty while setting the scope and time frame (counties in Arizona and the year 2020).

from epymorph.geography.us_census import CountyScope
from epymorph.time import TimeFrame

rume = ExampleRUME().build(
    scope=CountyScope.in_states(["AZ"], year=2020),
    time_frame=TimeFrame.rangex("2020-01-01", "2020-03-01"),
    params={}, # <-- empty for now...
)

print(rume.params_description())

gpm:youth::ipm::beta (type: float, shape: TxN)
    infectivity

gpm:youth::ipm::gamma (type: float, shape: TxN)
    progression from infected to recovered

gpm:youth::ipm::xi (type: float, shape: TxN)
    progression from recovered to susceptible

gpm:adult::ipm::beta (type: float, shape: TxN)
    infectivity

gpm:adult::ipm::gamma (type: float, shape: TxN)
    progression from infected to recovered

gpm:adult::ipm::xi (type: float, shape: TxN)
    progression from recovered to susceptible

meta::ipm::beta_ya (type: float, shape: TxN)

meta::ipm::beta_ay (type: float, shape: TxN)

gpm:youth::mm::centroid (type: [(longitude, float), (latitude, float)], shape: N)
    The centroids for each node as (longitude, latitude) tuples.

gpm:youth::mm::phi (type: float, shape: S, default: 40.0)
    Influences the distance that movers tend to travel.

gpm:youth::mm::commuter_proportion (type: float, shape: S, default: 0.1)
    The proportion of the total population which commutes.

gpm:youth::init::population (type: int, shape: N)
    The population at each geo node.

gpm:adult::mm::centroid (type: [(longitude, float), (latitude, float)], shape: N)
    The centroids for each node as (longitude, latitude) tuples.

gpm:adult::mm::phi (type: float, shape: S, default: 40.0)
    Influences the distance that movers tend to travel.

gpm:adult::mm::commuter_proportion (type: float, shape: S, default: 0.1)
    The proportion of the total population which commutes.

gpm:adult::init::population (type: int, shape: N)
    The population at each geo node.

This parameters list is pretty long, but it’s basically just the same set of parameters we would expect for a single strata RUME multiplied by the number of strata, plus our meta requirements. When specifying values for these, we can use asterisks as wild cards to cover multiple strata, modules, requirements, or any combination of those, so we don’t have to repeat ourselves.

Fully-specified parameter names are in the format:

"<stratum name>::<module name>::<requirement name>"

That’s three parts — stratum, module, requirement — separated by double-colons. We use gpm:youth to identify the youth stratum and gpm:adult to identify the adult stratum. There is also the meta stratum for our new meta parameters. There are parameters for modules ipm (the pathogen model), mm (the movement model), and init (the initializer).

Now we’ll reconstruct the RUME with the parameters filled in.

# additional imports....
from epymorph.adrio import acs5, us_tiger
from epymorph.log.messaging import sim_messaging

rume = ExampleRUME().build(
    scope=CountyScope.in_states(["AZ"], year=2020),
    time_frame=TimeFrame.rangex("2020-01-01", "2021-01-01"),
    params={
        # SIRS parameters for the youth base IPM
        "gpm:youth::ipm::beta": 0.2,
        "gpm:youth::ipm::gamma": 1 / 4,
        "gpm:youth::ipm::xi": 1 / 45,
        # SIRS parameters for the adult base IPM
        "gpm:adult::ipm::beta": 0.5,
        "gpm:adult::ipm::gamma": 1 / 4,
        "gpm:adult::ipm::xi": 1 / 45,
        # Parameters for the "meta edges" of our IPM
        "meta::ipm::beta_ya": 0.2,
        "meta::ipm::beta_ay": 0.5,
        # Movement and Initialization parameters
        "*::*::centroid": us_tiger.InternalPoint(),
        "*::*::population_by_age_table": acs5.PopulationByAgeTable(),
        "gpm:youth::*::population": acs5.PopulationByAge(0, 19),
        "gpm:adult::*::population": acs5.PopulationByAge(20, None),
    },
)

I’ve used a beta value of 0.2 when it concerns youth infection, and 0.5 when it concerns adult infection so as to model the age differences: a higher beta value corresponds to a higher transmission rate. Note that we have the freedom to use different beta values for, e.g., youth-adult versus youth-youth infection. For this example, we’ll keep them the same.

Next, for centroids, it would be very strange if the geographic centroids of our counties were different between strata or modules, so we use a double wild card for the centroid parameter. It’s the same values no matter where it’s used.

Lastly, the population parameters deserve a closer look. These parameters are pretty unique: they’re used both by the movement model and by the initializer, but the values used should be the same in each case. Hence I’ve used a module wild card. Populations do vary between age groups, though, so I’ve been explicit about which stratum each applies to. As with all parameters, I could have provided an array of values myself, but it’s easier to use our built-in ADRIOs for fetching population estimates from the US Census Bureau’s ACS 5-year data. Each instance of a PopulationByAge ADRIO accepts arguments specifying the desired age range. These data will be fetched and combined so as to produce the population total for that age bracket, if possible. My chosen age brackets have to align with the brackets provided in the ACS data. (See the ADRIO guide for more info.)

Using ACS Population by Age data

Beware a common pitfall when using PopulationByAge ADRIOs. You might think you could use these ADRIOs on their own, but this will result in an error when you run a simulation. For example:

SimValidationError: RUME attribute requirements were not met. See errors:
- Cannot evaluate parameters, there are missing values:
gpm:youth::mm::population_by_age_table
gpm:youth::init::population_by_age_table
gpm:adult::mm::population_by_age_table
gpm:adult::init::population_by_age_table

PopulationByAge is a rare ADRIO that, aside from being used to fulfill RUME requirements, also has requirements of its own. It’s expecting to find a parameter called "population_by_age_table", and this value is intended to be shared between all strata. We need to add this parameter to the RUME with the value being a PopulationByAgeTable ADRIO (or an equivalent data source). Then the PopulationByAge ADRIOs can do their work without loading the same dataset twice.

If you look at the example above you’ll see we have, in fact, provided the necessary “table” parameter.

This is, incidentally, a quirk that we are hoping to smooth out in future epymorph versions.

Running simulations

From here, running a simulation is similar to the single stratum case. We will visualize the results by aggregating the state of Arizona and showing the Infectious compartments from both strata.

from epymorph.kit import BasicSimulator, default_rng, sim_messaging

with sim_messaging():
    out = BasicSimulator(rume).run(
        rng_factory=default_rng(seed=42),
    )

out.plot.line(
    geo=rume.scope.select.all().group("state").sum(),
    time=rume.time_frame.select.all(),
    quantity=rume.ipm.select.compartments("I_youth", "I_adult"),
)

Loading gpm:youth::init::population_by_age_table (epymorph.adrio.acs5.PopulationByAgeTable):
  |####################| 100%  (1.496s)
Running simulation (BasicSimulator):
• 2020-01-01 to 2020-12-31 (366 days)
• 15 geo nodes
  |####################| 100% 
Runtime: 2.910s

Notice that a multistrata RUME automatically appends the strata names to compartments as subscripts. The equivalent thing happens to event names. When making selections you can also use wild cards: "I_*" selects all compartments with “I” as a base name.

Our results show there much less infection in the youth population. This may not be surprising: the youth population is much smaller. But if we transform the chart to “per 100k population”, the relationship remains:

# Calculate strata populations by doing a slice-and-sum on the initials matrix.
# Our populations totals are stable throughout the simulation, so this is valid.
# Both youth and adult strata have 3 compartments each.
# (If that changed you would have to change these indices.)
youth_pop = out.initial[:, :3].sum()
adult_pop = out.initial[:, 3:].sum()


def per_100k(df):
    qty = df["quantity"].iloc[0]
    p = (youth_pop if qty == "I_youth" else adult_pop) / 100_000
    return df.assign(value=df["value"] / p)


out.plot.line(
    geo=rume.scope.select.all().group("state").sum(),
    time=rume.time_frame.select.all(),
    quantity=rume.ipm.select.compartments("I_*"),
    transform=per_100k,
)

The combined IPM

As mentioned, when you create a multistrata RUME you are, in effect, creating one large combined IPM. You can get access to the resulting IPM for querying the compartments and events, or drawing diagrams, just as you would a single-stratum IPM.

print([str(c.name) for c in rume.ipm.compartments])

['S_youth', 'I_youth', 'R_youth', 'S_adult', 'I_adult', 'R_adult']

rume.ipm.diagram()

Movement and time

The movement model you specify for a stratum applies only to the individuals of that stratum. However there is one detail worth noting: a RUME’s movement model defines the tau steps to use, and it’s possible to have different movement models in your multistrata RUME that use different tau steps. In this case, the multistrata RUME computes a new set of tau steps that act as the combination of the tau steps from the movement models, and reworks the movement models to use these modified tau steps.

For example, if you are building a model with two strata and assign a movement model to each, where:

Stratum One’s movement model uses tau steps: (1/3, 2/3);
Stratum Two’s movement model uses tau steps: (2/3, 1/3);

Then the combined movement model will use tau steps: (1/3, 1/3, 1/3).

A diagram of the combination of different tau step schemes.

This combination preserves all of the dividing lines of the underlying models, so movement clauses can still fire when they were supposed to. For example, if Stratum Two’s original movement model has a clause which fires on the second tau step (when we are two-thirds through the day), in the combined movement model it will fire on the third tau step (still two-thirds through the day.)

The combination and necessary adjustments are done automatically, but you should be aware of this because the new tau steps will be reflected in the structure of the simulation output. If you’re not sure what the combined tau step structure will be, you can query your RUME to find out with the tau_step_lengths and num_tau_steps properties.

Naturally, if all of the movement models in your RUME use the same set of tau steps there are no changes necessary: it will just use those tau steps as-is.

Initialization

The initializer you specify for a stratum applies only to the individuals of that stratum. There’s nothing special to note here; it’s just epymorph’s modular design used to maximum advantage. All initializers will work the same way for a single-stratum RUME as for a multistrata RUME.