uofuepibio/epiworld

A general framework for quick epidemiological ABM models

epiworld

This C++ library provides a general framework for epidemiologic simulation. The core principle of epiworld is fast epidemiological prototyping for building complex models quickly. Here are some of its main features:

• It only depends on the standard library (C++11 required.)
• It is a template library.
• It is header-only (single file).
• Models can have an arbitrary set of states.
• Viruses and tools (e.g., vaccines, mask-wearing) can be designed to have arbitrary features.
• Multiple tools and viruses can live in the same simulation.
• It is FAST: About 30 Million person/day simulations per second (see example below).

Various examples can be found in the examples folder.

Hello world

Here is a simple SIR model implemented with epiworld. The source code can be found here, and you can compile the code as follows:

g++ -std=c++17 -O2 readme.cpp -o readme.o

As you can see in helloworld.cpp, to use epiworld you only need to incorporate the single header file epiworld.hpp:

#include "epiworld.hpp"

using namespace epiworld;

int main()
{

    // epiworld already comes with a couple
    // of models, like the SIR
    epimodels::ModelSIR<> hello(
        "COVID-19", // Name of the virus
        0.01,        // Initial prevalence
        0.9,        // Transmission probability
        0.3         // Recovery probability
        );

    // We can simulate agents using a smallworld network
    // with 100,000 individuals, in this case
    hello.agents_smallworld(100000, 4L, false, .01);

    // Running the model and printing the results
    // Setting the number of days (100) and seed (122)
    hello.run(100, 122);
    hello.print();

    return 0;

}

Compiling (with make helloworld.o) and running the problem yields the following result:

_________________________________________________________________________
Running the model...
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
 done.
________________________________________________________________________________
________________________________________________________________________________
SIMULATION STUDY

Name of the model   : Susceptible-Infected-Recovered (SIR)
Population size     : 100000
Agents' data        : (none)
Number of entities  : 0
Days (duration)     : 100 (of 100)
Number of viruses   : 1
Last run elapsed t  : 103.00ms
Last run speed      : 96.34 million agents x day / second
Rewiring            : off

Global events:
 (none)

Virus(es):
 - COVID-19 (baseline prevalence: 1.00%)

Tool(s):
 (none)

Model parameters:
 - Recovery rate     : 0.3000
 - Transmission rate : 0.9000

Distribution of the population at time 100:
  - (0) Susceptible :  99000 -> 2565
  - (1) Infected    :   1000 -> 366
  - (2) Recovered   :      0 -> 97069

Transition Probabilities:
 - Susceptible  0.96  0.04  0.00
 - Infected     0.00  0.70  0.30
 - Recovered    0.00  0.00  1.00

Building from scratch

One of the best things about epiworld is the capability to build models from scratch. Here is one example (readme.cpp):

#include "epiworld.hpp"

using namespace epiworld;

int main()
{

    // Creating a model with three statuses:
    // - Susceptible: Status 0
    // - Infected: Status 1
    // - Recovered: Status 2
    Model<> model;
    model.add_status("Susceptible", default_update_susceptible<>);
    model.add_status("Infected", default_update_exposed<>);
    model.add_status("Recovered");

    // Desgining a virus: This virus will:
    // - Have a 90% transmission rate
    // - Have a 30% recovery rate
    // - Infected individuals become "Infected" (status 1)
    // - Recovered individuals become "Recovered" (status 2)
    // 100 individuals will have the virus from the beginning.
    Virus<> virus("covid 19");

    virus.set_prob_infecting(.90);
    virus.set_prob_recovery(.30);
    
    virus.set_status(1, 2);

    model.default_add_virus<TSeq>n(virus, 1000);
    
    // Generating a random pop from a smallworld network
    model.agents_smallworld(100000, 4L, false, .01);

    // Initializing setting days and seed
    model.init(100, 122);

    // Running the model
    model.run();
    model.print();
  
}

Which should print something like the following:

Running the model...
_________________________________________________________________________
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.

________________________________________________________________________________
SIMULATION STUDY

Name of the model   : (none)
Population size     : 100000
Number of entitites : 0
Days (duration)     : 100 (of 100)
Number of variants  : 1
Last run elapsed t  : 209.00ms
Last run speed      : 47.64 million agents x day / second
Rewiring            : off

Virus(es):
 - covid 19 (baseline prevalence: 1000 seeds)

Tool(s):
 (none)

Model parameters:
 (none)

Distribution of the population at time 100:
 - (0) Susceptible :  99000 -> 2565
 - (1) Infected    :   1000 -> 366
 - (2) Recovered   :      0 -> 97069

Transition Probabilities:
 - Susceptible  0.96  0.04  0.00
 - Infected     0.00  0.70  0.30
 - Recovered    0.00  0.00  1.00

Which took about 0.209 seconds (~ 47 million ppl x day / second).

Simulation Steps

The core logic of the model relies on user-defined statuses and their corresponding update functions. In particular, the model does not have a predefined set of statuses, e.g., susceptible, infected, recovered; the user establishes them. This provides flexibility as models in epiworld can have an arbitrary set of statuses.

Like most other ABMs, epiworld simulates the evolution of a system in discrete steps. Each step represents a day in the system, and changes are reflected at the beginning of the following day. Therefore, agents can recover and transmit a virus on the same day. A single step of epiworld features the following procedures:

Status update: Agents are updated according to their status.

2. (optional) Execute Global events: A call of user-defined functions affecting the system. These can make any type of change in the system.

3. (optional) Apply rewiring algorithm: When specified, the network is rewired according to a user-defined function.

4. Lock the results: The current date is incremented in one unit, and the changes (exposition, new infections, recoveries, etc.) are recorded in the database.

5. (optional) Mutate Variants: When defined, variants can mutate, with the new variants appearing the next day.

To speed up computations, epiworld uses a queuing system that decides which agents will be active during each step and which will not. Agents are active when either they or at least one of their neighbors has a virus active. Agents' updates are triggered only for those in the queue, accelerating the completion of the current step.

Agents

Agents carry two sets of important information: viruses and tools. Each agent can have multiple instances of them, meaning multiple viruses and tools can coexist in a model. At each step of the simulation, an agent can face the following changes:

• Acquire a virus (add_virus()): Become exposed to a particular virus+host.

• Lose a virus (rm_virus()): Removing a virus from the agent. Losing a virus triggers a call to the virus's postrecovery() function, which can, for example, result in gaining immunity to that variant.

• Acquire a tool (add_tool()): For example, mask-wearing, vaccines, etc.

• Lose a tool (rm_tool()): For example, stop wearing masks, lose immunity, etc.

• Change status (change_status()): An arbitrary change in the status of the agent. Examples of this are moving from "exposed" to "infected," from "infected" to "ICU," etc.

• Become removed (rm_agent_by_virus()): An agent becomes inactive after its condition worsens. In such a case, all viruses attached to the agent are also removed.

Any action in the model can trigger a change in its queuing system. By default, becoming exposed makes the agent (and its neighbors) active in the queuing system. Likewise, losing all viruses could make the agent and its neighbors inactive.

Contagion

Susceptible individuals can acquire a virus from any of their infected connections. The probability that susceptible individual i gets the virus v from individual j depends on how three things:

1. The transmissibility of the virus, ,
2. The contagion reduction factor of i, , and
3. The host's transmission reduction factor, .

The last two are computed from and 's tools. Ultimately, the probability of getting virus $v$ from equals:

Nonetheless, the default behavior of the simulation model is to assume that individuals can acquire only one disease at a time, if any. This way, the actual probability is:

The latter is calculated using Bayes' rule

Where

This way, viruses with higher transmissibility will be more likely to be acquired when competing with other variants.

epiworld's Issues

Implement missing model.entity methods

The following lack of implementation:

• Model.rm_entity
• Model.add_entity_n
• Model.add_entity_f

Recompute connected model transmission

The rationale to compute transmission should be different.

1. Using contact_rate, sample neighbors from the population.
2. Follow the same steps as in the ModelSIR class, i.e., regular transmission prob.

This alternative should be used whenever any of the follow conditions is true:

• Transmission probability in the virus is a function.
• There are tools involved.

Add WebAssembly Demo

• Research feasibility of WebAssembly v.s. Shiny (Python/Rlang)
• Implement

Always reset the state of the agents

Currently, when reset()=true and there's a backup population, the agents' states are not reset. This should not be the case, as the backup should only preserve network structure (and data), not states, tools, and viruses.

Design and implement print method for entities and entity

Both lack of a print method.

Global action changing status

A contradiction can occur when a global action changes an agent to a status that depends on a virus, e.g.:

1. Update status moves agent to recovered, and the virus is to be removed.
2. The global action sees the agent still in the model, and moves the agent to infectious.
3. The virus will be removed, but the status is set as the last one, infectious.
4. The agent will then be incorporated into the model as infectious, but with no virus, causing an error.

Pair status update with status

The Person::update_status function uses the IN template function which right now is causing a bit of a bottleneck. Instead of checking set inclusion, we could leverage that statuses are unsigned ints, and thus an updating function could be reached by simply accessing by index, i.e., update_status_funs[p->status_current]();.