Movement Basics

In the geography chapter, we discussed how spatially-explicit models represent a finite set of discrete locations.

An illustration of locations as discrete places containing groups of individuals.

Every individual in the simulation is at one of these locations at every point in time, and within a location all individuals are assumed to have a roughly equal chance of interacting with each other (they are well-mixed). This is a fundamental assumption of the mathematics of compartmental modeling, which originated from chemical kinetics simulations from chemistry.

The inevitable question, then, is how can a pathogen transmit between discrete locations? In epymorph, we assume this is due entirely to the movement of infected individuals in the model. For example: a commuter with the flu travels from their home location to their work location and has a chance to infect others there before returning home. Of course in the real world, movement patterns can be arbitrarily complex. Not all of this complexity, however, is useful in our effort to model and predict the spread of a pathogen, so we simplify in an attempt to be parsimonious.

epymorph represents movement mechanics as the Movement Model: a set of rules that govern how individuals move between locations in the simulation. With movement defined as a module, you can swap for more- or less-complex movement models to decide how much complexity is warranted for your use-case. Or even to test different hypotheses about the impact of movement on disease spread.

The Effect of Movement

Let’s repeat the Getting Started example (which models four US states) with no movement. We can simply replace the mm.Centroids() configuration with another model, in this case mm.No().

from epymorph.kit import *
from epymorph.adrio import us_tiger, acs5


rume = SingleStrataRUME.build(
    ipm=ipm.SIRS(),
    mm=mm.No(),  # <--- UPDATED!
    init=init.SingleLocation(location=0, seed_size=100),
    scope=StateScope.in_states(["AZ", "CO", "NM", "UT"], year=2020),
    time_frame=TimeFrame.rangex("2020-01-01", "2021-01-01"),
    params={
        "beta": 0.3,
        "gamma": 1/5,
        "xi": 1/90,
        "population": acs5.Population(),
        "centroid": us_tiger.InternalPoint(),
    },
)

sim = BasicSimulator(rume)
output = sim.run(rng_factory=default_rng(42))

output.plot.line(
    geo=output.rume.scope.select.all(),
    time=output.rume.time_frame.select.all(),
    quantity=output.rume.ipm.select.compartments("I"),
)

Infection spreads where it’s initialized, in Arizona, but the other states are completely isolated.

When we use a movement model that causes the populations to interact, we see that infection “spills over” to the other locations after some time because they are exchanging individuals.

rume = SingleStrataRUME.build(
    ...
    mm=mm.Centroids(),
    ...
)

# ... (unchanged lines omitted for brevity)

The Centroids movement model implements basic commuter movement where a fixed proportion of the population commutes every day. Commuters travel to another location for 1/3 of a day and then return home for the remaining 2/3 of the day. Exactly where commuters wind up is randomly determined based on a probability which is proportional to the distance between locations. This is not a strict model of real commuter behavior (people tend to have real jobs they go back to every day!) but it’s a reasonable approximation. Remember the goal is parsimony, which is not just about the time it takes to compute the simulation, but also about how much data you have to provide to run the simulation. (I don’t have a perfectly-accurate list of everyone’s jobs in my back pocket!)

The Centroids movement model takes a number of parameters which can modify the results you get out of it. For example, commuter_proportion sets the fraction of the population that commute daily. It has a default of 10%, but if we set it higher we’ll see the effects of movement happen more rapidly.

rume = SingleStrataRUME.build(
    ...
    mm=mm.Centroids(),
    ...
    params={
        ...
        "commuter_proportion": 0.5,
    },
)

In the prior simulation, Colorado’s infection peaked around day 260. With the increased movement proportion it peaks around day 210. Arizona’s trajectory is almost unchanged, though, as its infection is influenced far more by its initial infection than subsequent movement. (Also, by the latter half of the simulation a large number of the population have temporary immunity from prior infection; the “R” in “SIRS”.)

So we see movement plays a vital role in infectious disease simulation.

Movement’s Role in the Simulation

Traditional approaches to compartmental modeling often combine the effects of movement and the disease mechanics in the same set of differential equations. epymorph decomposes these into separate modules, gaining a number of benefits including ease of development and comparison, and so that experimenting with different theories of movement is as simple as swapping one model for another. How does this work?

As discussed in IPM Basics, the disease model is still represented as a system of differential equations. Computing the trajectory of a pathogen model requires integrating these differential equations forward in time from their initial values. epymorph uses the tau leaping algorithm. This method does not require solving the equations analytically, which can be computationally prohibitive. It computes integer-valued results so we do not need to interpret fractional numbers of individuals. And it provides a source of stochasticity: tau-leaping is based on poisson random variables. But for movement, the most important advantage is that tau leaping provides a natural “seam”, a stopping point, in which to interpose additional computation. It is in these seams that epymorph applies the effects of movement.

The base unit of time in epymorph is a single day. But the Movement Model divides that day up into parts. During each part of a day, the disease model is advanced by computing a tau leap. But before each part, we apply the effects of movement. It is a simplifying assumption that movement happens instantaneously, and there’s not any transitional space (like an airport) that movers exist in for any period of time.

So epymorph’s simulations follow a simple loop: first apply movement, potentially shuffling individuals between locations; then advance the disease model using tau leaping. Repeat for all tau steps to complete a single day. Finally, repeat all of the above for every day of the simulation time frame. We record the state of the simulation after each tau step to produce the simulation’s output.

Figure: An example simulation day using two tau steps of lengths one-third and two-thirds, respectively

Why are tau steps defined by the Movement Model? One of our target use-cases for epymorph is to model movement patterns like human commuting, which involves spending a portion of the day, the “work day”, at one’s work location and the rest of the day at home. We need the “seams” to line up with the work day so we can process movement at the appropriate times. More complex movement models might divide the day into more parts, with more rules governing the movement of individuals. This gives movement models in epymorph a lot of flexibility. Of course some use-cases for movement modeling might have no need to divide the day at all. (The “No” movement model uses the full day as its only tau step, for example.) The longest tau step possible in epymorph is one day. There is no fixed limit on how short tau steps can be, nor how many tau steps you can have in your model. However beware that more tau steps means more computation time and an extremely small tau step could cause other issues. (Though we’ve tested up to about a million tau steps per day with good results, so normal use-cases should be well covered!)

TipBut is this equivalent?

It’s natural to question if epymorph’s approach produces results which are comparable to the more traditional approaches mentioned. Our experimentation involved implementing epymorph-equivalent models to those in commonly-cited literature. We concluded that epymorph’s approach does in fact produce qualitatively similar results to pure ODE modeling involving movement.

So now we know when movement happens, but what does it mean to apply movement? Movement Models define a set of rules called clauses which govern movement. Clauses combine configuration (like when should this rule apply?) and a function that computes how many individuals should move between locations as a result of the clause. epymorph will check all of the clauses in the movement model and apply their effects one by one: randomly selecting individuals to move, placing them in a group (a cohort), and moving them to their destination. Cohorts keep track of whether or not they need to return home, and if so, when (their “travel itinerary”, which is determined by the clause configuration.) The final step is an implicit clause, to check for any cohorts which are due to return home and to move them back to their home location.

Let’s look at an example using a hypothetical movement clause:

epymorph cohort movement example

In this example, a movement clause requested 123 people to move from Arizona to Colorado. epymorph first randomly selects that number of individuals from Arizona’s resident population. In this case it selected 100 Susceptible, 20 Infectious, and 3 Recovered individuals. It subtracts the selected individuals from the resident cohort and creates a new cohort to contain the movers. This new cohort records its return itinerary and is placed in Colorado’s list of cohorts. This basic process is repeated for every pair of locations in the simulation, including folks moving from Colorado to Arizona.

As long as this cohort is in Colorado it participates in disease progression, potentially spreading infection or becoming infected itself. After two ticks, the cohort is due to return home. It will be moved back to Arizona and merged back into Arizona’s resident cohort.

All of this cohort book-keeping is handled automatically by epymorph. During the tau leap, when the disease model is advanced, all cohorts at a location are considered well-mixed as though they were part of the same population. However IPM transition events get distributed among the cohorts proportionally, so they participate in the disease progression as you would expect while maintaining their unique identity and movement scheduling. It’s some complicated work, but that’s why epymorph is here!

Movement Models and Clauses

So far we’ve talked about Movement Models in the abstract, but let’s look at them in a bit more detail. Movement Models are made up of one or more clauses. When evaluated, a clause’s function returns a matrix describing how many individuals should move from every location to every other location. It is a node-by-node (NxN) numpy array of integers. “Should” is a key distinction: your model is making a request for movement and the system determines if and how to execute that movement. If your model requests 1000 people to move from a location but there are only 100 people available to move, the request will be only partially fulfilled. (100 people would be moved and the remainder is ignored.)

The Centroids model used above has a single clause, that is, just one rule for how people move around. Here’s an example result from that clause during the simulation above:

array([[3585154,      28,    1379,     471],
       [     20, 2840726,     996,     721],
       [    434,     404, 1047657,      15],
       [    227,     388,      23, 1574981]])

In this array, the first axis (rows) is the “from” location and the second axis (columns) is “to”. The nodes are in geo scope order: AZ, CO, NM, UT (which, for US Census geo scopes, is their geo ID ordering, which is the same as their alphabetical ordering.) So on this tick, the model has requested 28 people move from Arizona to Colorado. Note that values on the diagonal of the matrix are much larger than the off-diagonals. These values represent self-movement (e.g., within Arizona) and are a quirk of how the Centroids model works. epymorph effectively ignores self-movement, so the movement request is interpreted more like this:

array([[   0,   28, 1379,  471],
       [  20,    0,  996,  721],
       [ 434,  404,    0,   15],
       [ 227,  388,   23,    0]])

So this is our clause’s request. epymorph then decides how to fulfill that request. It follows these basic steps for each origin location:

  • Compute the requested total number of individuals leaving this origin, \(m\).
  • If the requested total is larger than the number of available movers, reduce the requested values via random selection and cap the value of \(m\).
  • Randomly select \(m\) individuals from the pool of available movers at this origin, keeping track of their disease compartment.
  • Randomly distribute the selected individuals to destinations, interpreting the requested movement array as a set of destination weights.

The ultimate result is a node-by-node-by-compartment array (NxNxC). It tells us not only how many individuals travel from each node to each other node, but also which disease compartments they belong to.

An important consequence to keep in mind is that movement requests will be approximately fulfilled, and not exactly fulfilled. Even considering situations where movement is not capped due to a lack of available movers, while the total number of movers in your request is preserved, the number of movers going to each destination is preserved only in a probabilistic sense.

Movement fulfillment is another important source of stochasticity in epymorph.

NoteCohort Book-keeping

Although you will not need to interact with this directly, it can be helpful to understand how the movement system deals with cohorts from a book-keeping perspective. In epymorph’s internal model of the world, each location maintains a list of cohorts which are present at that location at the current time.

Cohorts that have a return itinerary are considered to be temporarily visiting another location and are not subject to further movement selection. Cohorts with no itinerary are considered to be “home” and this makes them generally available to be selected for movement. Cohorts that have different itineraries must be tracked separately, while cohorts with the same itineraries (which are at the same location) get merged.

For commuter-style movement these properties guarantee that your commuters don’t lose track of where they came from by getting moved twice. If you want to model migratory-style movement, you must make sure your movement clause defines no return rule, so that those cohorts will merge into the destination population and be subject to further movement in future ticks.

Incidentally, this kind of book keeping is how we keep track of who is where when computing “data by home” vs. “data by visit”, as discussed in the chapter on Basic Simulator’s Output Data Aggregation Modes. Both data modes count every person, but they are either counted by their origin (“by home”) or by whichever location they happen to be at, at each point in time (“by visit”).

Evaluating the clause’s function gives us the requested magnitude of movement, but the clause also defines configuration which determines when this movement occurs and when they should return home (if at all).

Here is a simplified definition of the Centroids movement model and its clause:

class Centroids(MovementModel):
    steps = (1 / 3, 2 / 3)
    clauses = (CentroidsClause(),)

class CentroidsClause(MovementClause):
    # ... (omitted: data requirements definitions)

    predicate = EveryDay()
    leaves = TickIndex(step=0)
    returns = TickDelta(step=1, days=0)

    # ... (omitted: the evaluate function)

The MovementModel implementation declares the tau steps for the model: one-third and two-thirds. It then lists the clauses it contains: just one in this case. The MovementClause implementation has three critical attributes:

  • predicate: Does this clause apply to the current time step? This is mostly used when movement differs by day of the week.
  • leaves: On which tau step does this clause apply? Tau steps are indexed from zero, so the zeroth tau step is the first one. In this case people move before the first tau step.
  • returns: On which tau step should these movers return home? In this case, with a delta of zero days and tau step index one, people moved by this clause will return before the second tau step of the same day.

Thus the Centroids model defines a system of movement where our movers leave their homes to be at another location (work) for the first third of the day, before returning home for the remaining two thirds of the day.

We could imagine adding a “vacation” clause in which people leave at the start of the day on a Saturday and return home at the start of the Saturday a week later. That might look like:

class VacationClause(MovementClause):
    # ... (omitted: data requirements definitions)

    predicate = DayIs(["Sa"])
    leaves = TickIndex(step=0)
    returns = TickDelta(step=0, days=7)

    # ... (omitted: the evaluate function)

We’ll discuss creating custom movement models in detail in the next chapter.

Immobility Due to Disease State

One additional feature blurs the line between epymorph’s compartment model and movement model. We may wish for certain disease states to make the individual immobile (for example, people in the hospital rarely move). IPMs can apply tags when defining their compartments. The special tag “immobile” indicates that individuals in this state should not be subjected to movement, and thus are ignored during movement calculations.

See the section on defining IPMs for more detail.