Prototyping Polio‐Inspired Features
The goal here was to show we could have (vaccine) interventions that are targeted in time and space, based on measured incidence or prevalence, in addition to or in contrast to pre-scheduled interventions.
This is relatively straightforward, assuming we are reporting incidence and prevalence by time and space, mid-simulation, and can adapt our SIA (or related) intervention to have access to the latest reports.
In my prototype, Incidence and Prevalence reporting were already being done and we already had an SIA component. Reporting can be a performance hit in LASER if done as a standalone component or step. But in my laser-polio prototype, I was already doing reporting of S,E,I(infectious),R,W and Incidence by node and timestep, in a performant way. These report channels are in the centralized model instance which is passed to components and so the necessary data is available to the SIA component. It's possible, depending on exactly how things are ordered, that the intervention component has access to the previous timestep and not the current timestep. I imagine either approach is usable, but the modeler probably wants to know exactly which they are using. That will be in their control.
In my prototype, I added code to the SIA component in which I checked every week (every 7 time steps) for incidence by node over the previous week, and distributed vaccine interventions to the nodes with over 1000 infections based on a (nominal and deterministic) ratio of of 1000:1 for infections to observed cases. These vaccines were distributed to 75% of susceptibles. None of that is very hard and JENNER can give numpy code to do that if one wants.
A quick note on age-targeting. There is no inherent need to have an age or even date_of_birth property in LASER, so we only add if we are going to need it and use it. For scheduled, age-targeted interventions, we can and may prefer to use timers. And even for age-based transmission/mixing, we may use timers. But for reactive interventions, if we want to target agents specifically by age, and not just target susceptibles as a proxy for age, we will almost certainly want an age or dob property. There is a discussion elsewhere on the tradeoffs involved in those choices. In short, dob requires less maintenance (routine updating) if you find it sufficiently intuitive to work with.
In the plots above, we see some nodes where outbreaks start to take off, and get "nipped in the bud", at least at a qualitative level. I did not do any quantitative analysis or testing.
I didn't measure exact performance impact but no noticeable slowdowns were perceived.
Each individual will get a personal transmission heterogeneity factor (multiplier), drawn from 0.0 to floatmax and which stays constant, and determines their actual shed amount when infected. It isn't correlated with any other property.
In my laser-polio prototype, infectious shedding is done in the Transmission component. Previously the "contagion" was the sum of infectious individuals by node, multiplied by an infectiousness value derived from beta or R-nought. This could be the output of a seasonality calculation.
I started with the intuitive code:
# Sum the tx_hetero_factor values for the infectious individuals, grouped by node
infectious_mask = model.population.itimer > 0
contagion = np.zeros(len(nodes))
np.add.at(contagion, model.population.nodeid[infectious_mask], model.population.tx_hetero_factor[infectious_mask])
But GPT recommended the faster:
contagious_indices = model.population.nodeid[infectious_mask]
values = model.population.tx_hetero_factor[infectious_mask]
contagion = np.bincount(contagious_indices, weights=values, minlength=len(nodes)).astype(np.float32)
The transmission component tends to be the overall performance bottleneck, which isn't surprising. It's where our agents are doing all the interaction with each other. And any time we are doing centralized counting or accumulating, we have to put in a little more effort for smart and performant parallelization. So it's easy for a line of code like the above to be surprisingly slower than expected. Hence the bin count approach.
This seems to work as expected. There's not a lot to it. Though I haven't thought of a good way to show that it's working.
Depending on how good our bincount solution performs, we may want to opt for a numba or ctypes implementation. E.g.,
C
#include <omp.h>
#include <stdint.h>
void accumulate_contagion(int *nodeid, double *tx_hetero_factor, uint32_t *contagion, size_t num_agents, size_t num_nodes) {
#pragma omp parallel for
for (size_t i = 0; i < num_nodes; i++) {
// For each node, accumulate the values for all the agents that are assigned to it
for (size_t j = 0; j < num_agents; j++) {
if (nodeid[j] == i) {
contagion[i] += tx_hetero_factor[j];
}
}
}
}
We could also check for nodes of 0 prevalence and skip those altogether.
We want to move from homogeneous mixing by age to age-based mixing (transmission/acquisition).
We can either try to do something very HINT-like (EMOD) but we need to check perf impact of doing that in LASER. We may want to resort to some age-based transmission and/or acquisition modifiers. Also, we may want to use timers to transition people between age buckets rather than age-bucketizing everyone every transmission step. (Looking up everyone's age every time they transmit, and every time they might acquire, is a non-zero cost.) It's also probably safe to say that, generally speaking, heterogeneity in acquisition is more expensive than heterogeneity in transmission.
Once again, the Transmission component tends to be our bottleneck and a lot of optimization was done in laser-polio to make this fast. It's easy to make it slower again. Age-based mixing by its nature is going to get right into the heart of the transmission code. Also, as discussed in the first section, we don't inherently have ages in our agents so the approach we use to agent ages is pretty wide open.
In the approach I took, I opted to prototype this in a standalone, rather than add to the existing laser-polio prototype. It's a runnable notebook here.
We have a normalized age-based mixing matrix, such as:
self.age_mixing_matrix = np.array([[0.5, 0.3, 0.2],
[0.3, 0.5, 0.2],
[0.2, 0.2, 0.6]])
Our age bins were defined as:
- Age bin 0: 0-5
- Age bin 1: 5-20
- Age bin 2: 20+
We opted to use two timer properties and 1 age_bin property rather than an age property. What that means is that everyone was initialized with the age_bin property according to their initial age (or date of birth). And then their age_bin_countdown_timer was initialized, depending on their current bin. And we added a AgeBinComponent which counts down everyone's (non-zero) timer -- which is independent per agent's and can be done infrequently. When a timer expires, the age bin gets incremented, and reinitialized if applicable. All very trivial math and logic. The mixing matrix value lookups are done based on an agent's age_bin and we don't have to compare current ages to 2 age boundaries every time.
We shed every infectious person's contagion into a bucket based on their age_bin. We expose people to each bin, using the multipliers appropriate to their age_bin, and combine. Normalization questions are left up to the implementing modeler. The biggest thing to draw attention to here is that I prefer to do transmission by calculating the number of new infections (by node) first, and then picking that number of infectees (hopefully at random!) from the current susceptibles. This tends to be faster than combining those steps, and easier to implement with fixed memory. The code to do all this ended up being relatively simple and easy to follow.
See the notebook. Here's a plot which shows age-based prevalence for 1 simulation.
It looks like age-based mixing, with 3 age buckets, can triple the time we spend in the transmission step, at least in an all-python standalone prototype. Even applying speedups.
Instead of modeling infections and vaccinations as separate things, both of which can induce immunity, we want to model vaccines as a type of infection, such that the vaccines can transmit (at a lower rate) and revert, as they circulate, at which point they become identical to infections.
Each individual will have a vax_delivered_timestep, set to the current timestep when someone gets vaxed (SIA). During transmission, the infectivity of an infection will be a binary function of whether their vax has reached reversion duration (low for vax, high for reverted). Upon acquisition, a person will get "vaxed" or "infected" depending on reversion of transmitted strain.
This was actually a little tricky to implement. Details TBD. Also, without an environmental reservoir, I had to hack some disease parameter values in order to demonstrate that the code was working as intended.
A particular design change, and this is reminiscent of EMOD, is that with "vaccine infections" and "real infections" co-circulating in the same node, we have to use the total infectivity to calculate the number of new infections -- though this is driven primarily by the "real infections" -- and then we have to probabilistically apportion new infectees between the two kinds of infections.
Another design point is that I chose to use the timestep of initial vaccine distribution as the "strain id" such that the strain id can also be used as the baseline for the reversion timer. This enabled me to avoid creating any "master table" of vaccines and reversion timers.
DeepNote notebook with code you can run and existing plots from prior runs here.
TBD. This definitely needs to be explored.