Overview

The subpop_setup section of the configuration file is where users can input the information required to define a population structure on which to simulate the model. The options allow the user to determine the population size of each subpopulation that makes up the overall population, and to specify the amount of mixing that occurs between each pair of subpopulations.

An example configuration file with the global header and the spatial_setup section is below:

name: test_simulation
data_path: data
model_output_dirname: model_output
start_date: 2020-01-01
end_date: 2020-12-31
nslots: 100

subpop_setup:
  geodata: model_input/geodata.csv
  mobility: model_input/mobility.csv

Items and options

geodata file

• geodata is a .csv with column headers, with at least two columns: subpop and population.

• nodenames is the name of a column in geodata that specifies unique geographical identification strings for each subpopulation.

• selected is the list of selected locations in geodata to be modeled

Example geodata file format

subpop,population
10001,1000
20002,2000

mobility file

The mobility file is a .csv file (it has to contain .csv as extension) with long form comma separated values. Columns have to be named ori, dest, amount, with amount being the average number individuals moving from the origin subpopulation ori to destination subpopulation dest on any given day. Details on the mathematics of this model of contact are explained in the Model Description section. Unassigned relations are assumed to be zero. The location entries in the ori and dest columns should match exactly the subpop column in geodata.csv

Example mobility file format

ori, dest, amount
10001, 20002, 3
20002, 10001, 3

It is also possible, but not recommended to specify the mobility file as a .txt with space-separated values in the shape of a matrix. This matrix is symmetric and of size K x K, with K being the number of rows in geodata. The above example corresponds to

0 3
3 0

Examples

Example 1

To simulate a simple population structure with two subpopulations, a large province with 10,000 individuals and a small province with only 1,000 individuals, where every day 100 residents of the large province travel to the small province and interact with residents there, and 50 residents of the small province visit the large province

subpop_setup:
  geodata: model_input/geodata.csv
  mobility: model_input/mobility.csv

geodata.csv contains the population structure (with columns subpop and population)

subpop,          population
large_province, 10000
small_province, 1000

mobility.csv contains

ori,            dest,           amount
large_province, small_province, 100
small_province, large_province, 50

Within flepiMoP, gempyor is an open-source Python package that constructs and simulates compartmental infectious disease dynamics models. gempyor is meant to be used within flepiMoP, where it integrates with parameter inference and data processing scripts, but can also be run standalone with a command-line interface, generating simulations of disease incidence under different scenario assumptions.

To simulate an infectious disease dynamics problem, the following building blocks needs to be defined:

• The population structure over which the disease is transmitted

• The transmission model, defining the compartments and the transitions between compartments

• An observation model, defining different observable outcomes (serology, hospitalization, deaths, cases) from the transmission model

• The parameters and modifiers that apply to them

Generalized compartmental infection model

At the core of our pipeline is a dynamic mathematical model that categorizes individuals in the population into a discrete set of states ('compartments') and describes the rates at which transitions between states can occur. Our modeling approach was developed to describe classic infectious disease transmission models, like the SIR model, but is much more general. It can encode any compartmental model in which transitions between states are of the form

where $X$, $Y$, and $Z$ are time-dependent variables describing the number of individuals in each state, $b$ is a rate parameter (units of time$^{-1}$) and $a$ is a scaling parameter (unitless). $Z$ may be $X$, $Y$, a different variable, or 1, and the rate may also be the sum of terms of this form. Rates that involve non-linear functions or more than two variables are currently not possible. For simplicity, we omitted the time dependencies on parameters (e.g $X$ is in fact $X(t)$ and $a$,$b$ are $a(t)$,$b(t)$).

The model can be simulated as a continuous-time, deterministic process (i.e., a set of ordinary differential equations), which in this example would be in the form

Details on the numerical integration procedure for simulating such an equation is given in the Advanced section.

SEIR model

and as a stochastic process, it is

A common COVID-19 model is a variation of this SEIR model that incorporates:

1. multiple identical stages of the infectious period, which allows us to model gamma-distributed durations of infectiousness, and

A three-stage infectious period model is given by

The flepiMoP model structure is specifically designed to make it simple to encode the type of more complex "stratified'' models that often arise in infectious disease dynamics. The following are some examples of possible stratifications.

Age groups

To describe an SEIR-type disease that spreads and progresses differently among children versus adults, one may want to repeat each compartment of the model for each of the two age groups (C – Children, A – Adults), creating an age-stratified model

Vaccination status

Vaccination status could influence disease progression and infectiousness, and could also change over time as individuals choose to get the vaccine (V – vaccinated, U – unvaccinated)

Pathogen strain

Another common stratification would be pathogen strain, such as COVID-19 variants. Individuals may be infected with one of several variants, strains, or serotypes. Our framework can easily create multistrain models, for example

All combinations of these situations can be quickly specified in flepiMoP. Details on how to encode these models is provided in the Model Implementation section, with examples given in the Tutorials section.

Clinical outcomes and observations model

The pipeline allows for an additional type of dynamic state variable beyond those included in the mathematical model. We refer to these extra variables as "Outcomes" or "Observations". Outcome variables can be functions of model variables, but do not feed back into the model by influencing other state variables. Typically, we use outcome variables to describe the process through which some subset of individuals in a compartment are "observed'' and become part of the data to which models are compared and attempt to predict. For example, in the context of a model for an infectious disease like COVID-19, outcome variables include reported cases, hospitalizations, and deaths.

In addition to single values (drawn from a distribution), the duration and delay can be inputted as distributions, producing a convolution of the output.

Formally, for a deterministic, continuous-time model

For a discrete-time, stochastic model

Outcomes can also be constructed as functions of other outcomes. For example, a fraction of hospitalized patients may end up in the intensive care unit (ICU).

There are several benefits to separating outcome variables from the mathematical model. Firstly, these variables can be calculated after the model is run, and only at the timepoints of interest, which can dramatically reduce the memory needed during model simulation. Secondly, outcome variables can be fully stochastic even when the mathematical model is simulated deterministically. This becomes useful when an infection might be at high enough prevalence that a deterministic simulation is appropriate, but when there is a rare and therefore quite stochastic outcome reported in the data (e.g., severe cases) that the model is tasked with predicting. Thirdly, outcome variables can have arbitrary delay distributions, to take into account the complexities of health reporting practices, whereas our mathematical modeling framework is designed mainly for exponentially distributed delays and only easily permits extensions to gamma-distributed delays. Finally, this separation keeps the pipeline modular and allow for easy editing of one component of the model without disrupting the other.

Details on how to specify these outcomes in the model configuration files is provided in the Model Implementation section, with examples given in the Tutorials section.

Population structure

The pipeline was designed specifically to simulate infection dynamics in a set of connected subpopulations. These subpopulations could represent geographic divisions, like countries, states, provinces, or neighborhoods, or demographic groups, or potentially even different host species. The equations and parameters of the transmission and outcomes models are repeated for each subpopulation, but the values of the parameters can differ by location. Within each subpopulation, infection is equally likely to spread between any pair of susceptible/infected individuals after accounting for their infection class, whereas between subpopulations there may be varying levels of mixing.

Formally, this type of population structure is often referred to as a “metapopulation”, and each subpopulation may be called a “deme”.

The following properties may be different between subpopulations:

• the population size

• the parameters of the transmission model (see LINK)

• the parameters of the outcomes model (see LINK)

• the amount of transmission that occurs within this subpopulation versus from any other subpopulation (see LINK)

• the timing and extent of any interventions that modify these parameters (see LINK)

• the initial timing and number of external introductions of infections into the population (see LINK)

• the ground truth timeseries data used to compare to model output and infer model parameters (see LINK)

Currently, the following properties must be the same across all subpopulations:

• the compartmental model structure

• the form of the likelihood function used to estimate parameters by fitting the model to data (LINK)

• ...

Mixing between subpopulations

The generalized compartmental model allows for second order “interaction” terms that describe transitions between model states that depend on interactions between pairs of individuals. For example, in the context of a classical SIR model, the rate of new infections depends on interactions between susceptible and infectious individuals and the transmission rate

For a model with multiple subpopulations, each of these interactions can occur either between individuals in the same or different subpopulations, with specific rate parameters for each combination of individual locations

The transmission matrix is then

The list of all pairwise mobility values and the interaction scaling factor are model input parameters. Details on how to specify them are given in the Model Implementation section.

If an alternative compartmental disease model is created that has other interactions (second order terms), then the same mobility values are used to determine the degree of interaction between each pair of subpopulations.

Initial conditions

It might also be necessary to model instantaneous changes in values of model variables at any time during a simulation. We call this 'seeding'. For example, individuals may import infection from other external populations, or instantaneous mutations may occur, leading to new variants of the pathogen. These processes can be modeled with seeding, allowing individuals to change state at specified times independently of model equations.

We also note that seeding can also be used as a convenient way to specify initial conditions, particularly ealy in an outbreak where the outbreak is triggered by a few 'seedings'.

Time-dependent interventions

Parameters in the disease transmission model or the observation model may change over time. These changes could be, for example: environmental drivers of disease seasonality; “non-pharmaceutical interventions” like social distancing, isolation policies, or wearing of personal protective equipment; “pharmaceutical interventions” like vaccination, prophylaxis, or therapeutics; changes in healthcare seeking behavior like testing and diagnosis; changes in case reporting, etc.

The model allows for any parameter of the disease transmission model or the observation model to change to a new value for a time interval specified by start and end times (or multiple start and end times, for interventions that are recurring). Each change may be subpopulation-specific or apply to the entire population. Changes may be overlapping in time.

The magnitude of these changes are themselves model parameters, and thus may inferred along with other parameters when the model is fit to data. Currently, the start and end times of interventions must be fixed and cannot be varied or inferred.

• StackedModifier - TBA

Overview

flepiMoP allows users to specify instantaneous changes in values of model variables, at any time during the simulation. We call this "seeding". For example, some individuals in the population may travel or otherwise acquire infection from outside the population throughout the epidemic, and this importation of infection could be specified with the seeding option. As another example, new genetic variants of the pathogen may arise due to mutation and selection that occurs within infected individuals, and this generation of new strains can also be modeled with seeding. Seeding allows individuals to change state at specified times in ways that do not depend on the model equations. In the first example, the individuals would be "seeded" into the infected compartment from the susceptible compartment, and in the second example, individuals would be seeded into the "infected with new variant" compartment from the "infected with wild type" compartment.

The seeding option can also be used as a convenient alternative way to specify initial conditions. By default, flepiMoP initiates models by putting the entire population size (specified in the geodata file) in the first model compartment. If the desired initial condition is only slightly different than the default state, it may be more convenient to specify it with a few "seedings" that occur on the first day of the simulation. For example, for a simple SIR model where the desired initial condition is just a small number of infected individuals, this could be specified by a single seeding into the infected compartment from the susceptible compartment at time zero, instead of specifying the initial values of three separate compartments. For larger models, the difference becomes more relevant.

Specifying model seeding

The configuration items in the seeding section of the config file are

seeding:method Must be either "NoSeeding", "FromFile", "PoissonDistributed", "NegativeBinomialDistributed", or "FolderDraw".

seeding::seeding_file Only required for method: “FromFile”. Path to a .csv file containing the list of seeding events

seeding::lambda_file Only required for methods "PoissonDistributed" or "NegativeBinomialDistributed". Path to a .csv file containing the list of the events from which the actual seeding will be randomly drawn.

seeding::seeding_file_type Only required for method "FolderDraw". Either seir or seed

Details on implementing each seeding method and the options that go along with it are below.

seeding::method

NoSeeding

If there is no seeding, then the amount of individuals in each compartment will be initiated using the values specified in theinitial_conditions section and will only be changed at later times based on the equations defined in the seir section. No other arguments are needed in the seeding section in this case

Example

seeding:
    method: “NoSeeding”

FromFile

This seeding method reads in a user-defined file with a list of seeding events (instantaneous transitions of individuals between compartments) including the time of the event and subpopulation where it occurs, and the source and destination compartment of the individuals. For example, for the simple two-subpopulation SIR model where the outbreak starts with 5 individuals in the small province being infected from a source outside the population, the seeding section of the config could be specified as

seeding:
  method: "FromFile"
  seeding_file: seeding_2pop.csv

Where seeding.csv contains

subpop, date, amount, source_infection_stage, destination_infection_stage
small_province, 2020-02-01, 5, S, E

seeding::seeding_file must contain the following columns:

• subpop – the name of the subpopulation in which the seeding event takes place. Seeding cannot move individuals between different subpopulations.

• date – the date the seeding event occurs, in YYYY-MM-DD format

• amount – an integer value for the amount of individuals who transition between states in the seeding event

• source_* and destination_* – For each compartment group (i.e., infection stage, vaccination stage, age group), a different column describes the status of individuals before and after the transition described by the seeding event. For example, for a model where individuals are stratified by age and vaccination status, and a 1-day vaccination campaign for young children and the elderly moves a large number of individuals into a vaccinated state, this file could be something like

subpop, date, amount, source_infection_stage, source_vaccine_doses, source_age_group, destination_infection_stage, destination_vaccine_doses, destination_age_group
anytown, 1950-03-15, 452, S, 0dose, under5years, S, 1dose, under5years
anytown, 1950-03-16, 527, S, 0dose, 5_10years, S, 1dose, 5_10years
anytown, 1950-03-17, 1153, S, 0dose, over65years, S, 1dose, over65years

PoissonDistributed or NegativeBinomialDistributed

These methods are very similar to FromFile, except the seeding value used in the simulation is randomly drawn from the seeding value specified in the file, with an average value equal to the file value. These methods can be useful when the true seeded value is unknown, and only an observed value is available which is assumed to be observed with some uncertainty. The input requirements are the same for both distributions

seeding:
  method: "PoissonDistributed"
  lambda_file: seeding.csv

or

seeding:
  method: "NegativeBinomialDistributed"
  lambda_file: seeding.csv

and the lambda_file has the same format requirements as the seeding_file for the FromFile method described above.

For method::PoissonDistributed, the seeding value for each seeding event is drawn from a Poisson distribution with mean and variance equal to the value in the amount column. Formethod::NegativeBinomialDistributed, seeding is drawn from a negative binomial distribution with mean amount and variance amount+5 (so identical to "PoissonDistributed" for large values of amount but has higher variance for small values).

FolderDraw

TB ;

Thinking about outcomes variables

Our pipeline allows users to encode state variables describing the infection status of individuals in the population in two different ways. The first way is via the state variables and transitions of the compartmental model of disease transmission, which are specified in the compartments and seir sections of the config. This model should include all variables that influence the natural course of the epidemic (i.e., all variables that feed back into the model by influencing the rate of change of other variables). For example, the number of infected individuals influences the rate at which new infections occur, and the number of immune individuals influences the number of individuals at risk of acquiring infection.

However, these intrinsic model variables may be difficult to observe in the real world and so directly comparing model predictions about the values of these variables to data might not make sense. Instead, the observable outcomes of infection may include only a subset of individuals in any state, and may only be observed with a time delay. Thus, we allow users to define new outcome variables that are functions of the underlying model variables. Commonly used examples include detected cases or hospitalizations ;

Variables should not be included as outcomes if they influence the infection trajectory. The choice of what variables to include in the compartmental disease model vs. the outcomes section may be very model specific. For example, hospitalizations due to infection could be encoded as an outcome variable that is some fraction of infections, but if we believe hospitalized individuals are isolated from the population and don't contribute to onward infection, or that the number of hospitalizations feeds back into the population's perception of risk of infection and influences everyone's contact behavior, this would not be the best choice. Similarly, we could include deaths due to infection as an outcome variable that is also some fraction of infections, but unless death is a very rare outcome of infection and we aren't worried about actually removing deceased individuals from the modeled populations, deaths should be in the compartmental model instead.

The outcomes section is not required in the config. However, there are benefits to including it, even if the only outcome variable is set to be equivalent to one of the infection model variables. If the compartmental model is complicated but you only want to visualize a few output variables, the outcomes output file will be much easier to work with. Outcome variables always occur with some fixed delay from their source infection model variable, which can be more convenient than the exponential distribution underlying the infection model. Outcome variables can be created to automatically sum over multiple compartments of the infection model, removing the need for post-processing code to do this. If the model is being fit to data, then the outcomes section is required, as only outcome variables can be compared to data.

As an example, imagine we are simulating an SIR-style model and want to compare it to real epidemic data in which cases of infection and death from infection are reported. Our model doesn't explicitly include death, but suppose we know that 1% of all infections eventually lead to hospitalization, and that hospitalization occurs on average 1 week after infection. We know that not all infections are reported as cases, and assume that only 50% are detected and are reported 2 days after infection begins. The model and outcomes section of the config for these outcomes, which we call incidC (daily incidence of cases) and incidH (daily incidence of hospital admission) would be

compartments:
  infection_stage: ["S", "I", "R"]
  
seir:
  transitions:
    # infection
    - source: [S]
      destination: [I]
      proportional_to: [[S], [I]]
      rate: [beta]
      proportion_exponent: 1
    # recovery
    - source: [I]
      destination: [R]
      proportional_to: [[I]]
      rate: [gamma]
      proportion_exponent: 1
  parameters:
    beta: 
      value: 0.2
    gamma: 
      value: 0.1

outcomes:
  settings:
    method: delayframe
  outcomes:
    incidC:
      source:
        incidence:
          infection_stage: "I"
      probability: 
        value: 0.5
      delay: 
        value: 2
    incidH:
      source:
        incidence:
          infection_stage: "I"
      probability: 
        value: 0.01
      delay: 
        value: 21 

in the following sections we describe in more detail how this specification works

Specifying outcomes in the configuration file

The outcomes config section consists of a list of defined outcome variables (observables), which are defined by a user-created name (e.g., "incidH"). For each of these outcome variables, the user defines the source compartment(s) in the infectious disease model that they draw from and whether they draw from the incidence (new individuals entering into that compartment) or prevalence (total current individuals in that compartment). Each new outcome variable is always associated with two mandatory parameters ;

• probability of being counted in this outcome variable if in the source compartment

• ;delay between when an individual enters the source compartment and when they are counted in the outcome variable

and one optional parameter

• duration after entering that an individual is counted as part of the outcome variable

The value of the probability, delay, and duration parameters can be a single value or come from distribution ;

Outcome model parameters probability, delay, and distribution can have an additional attribute beyond value called modifier_key. This value is explained in the section on coding time-dependent parameter modifications (also known as "modifiers") as it provides a way to have the same modifier act on multiple different outcomes ;

Example

// Some code

`source ;

Required, unless sum option is used instead. This sub-section describes the compartment(s) in the infectious disease model from which this outcome variable is drawn. Outcome variables can be drawn from the incidence of a variable - meaning that some fraction of new individuals entering the infection model state each day are chosen to contribute to the outcome variable - or from the prevalence, meaning that each day some fraction of individuals currently in the infection state are chosen to contribute to the outcome variable. Note that whatever the source type, the named outcome variable itself is always a measure of incidence ;

To specify which compartment(s) contribute the user must specify the state(s) within each model stratification. For stratifications not mentioned, the outcome will sum over that states in all strata ;

For example, consider a configuration in which the compartmental model was constructed to track infection status stratified by vaccination status and age group. The following code would be used to create an outcome called incidH_child (incidence of hospitalization for children) and incidH_adult (incidence of hospitalization for adults) where some fraction of infected individuals would become hospitalized and we wanted to track separately track pediatric vs adult hospitalizations, but did not care about tracking the vaccination status of hospitalized individuals as in reality it was not tracked by the hospitals ;

 compartments:
   infection_state: ["S", "I", "R"]
   age_group: ["child", "adult"]
   vaccination_status: ["unvaxxed", "vaxxed"]
   
outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    ...
  incidH_adult:
    source:
      incidence:
        infection_state: "I"
        age_group: "adult"
    ...
  incidH_all:
    source:
      incidence:
        infection_state: "I"
    ...

to instead create an outcome variable for cases where on each day of infection there is some probability of testing positive (for example, for the situation of an asymptomatic infection where testing is administered totally randomly), the following code would be used

 compartments:
   infection_state: ["S", "I", "R"]
   age_group: ["child", "adult"]
   vaccination_status: ["unvaxxed", "vaxxed"]
   
outcomes:
  incidC:
    source:
      prevalence:
        infection_state: "I"
    ...

The source of an outcome variable can also be a previous defined outcome variable. For example, t to create a new variable for the number of individuals recruited to be part of a contact tracing program (incidT), which is just some fraction of diagnosed cases ;

outcomes:
  incidC:
    source:
      prevalence:
        infection_state: "I"
    ...
  incidT:
    source: incidC
    ...

`probability ;

Required, unless sum option is used instead. Probability is the fraction of individuals in the source compartment who are counted as part of this outcome variable (if the source is incidence; if the source is prevalence it is the fraction of individuals per day). It must be between 0 and 1 ;

Specifying the probability creates a parameter called outcome_name::probability that can be referred to in the outcome_modifiers section of the config. The value of this parameter can be changed using the probability::intervention_param_name option ;

For example, to track the incidence of hospitalization when 5% of children but only 1% of adults infected require hospitalization, and to create a modifier_key such that both of these rates could be modified by the same amount during some time period using the outcomes_modifier section:

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 
      value: 0.05
      modifier_key: hosp_rate
  incidH_adult:
    source:
      incidence:
        infection_state: "I"
        age_group: "adult"
    probability: 
      value: 0.01
      modifier_key: hosp_rate

To track the incidence of diagnosed cases iterating over uncertainty in the case detection rate (ranging 20% to 30%), and naming this parameter "case_detect_rate"

outcomes:
  incidC:
    source:
      prevalence:
        infection_state: "I"
    probability:
      value:
        distribution: uniform
        low: 
          value: 0.2
        high: 
          value: 0.3
      intervention_param_name: "case_detect_rate"

Each time the model is run a new random value for the probability of case detection will be chosen ;

Delay

Required, unless sum option is used instead. delay is the time delay between when individuals are chosen from the source compartment and when they are counted as part of this outcome variable ;

For example, to track the incidence of hospitalization when 5% of children are hospitalized and hospitalization occurs 7 days after infection:

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 
      value: 0.05
    delay: 
      value: 7

To iterate over uncertainty in the exact delay time, we could include some variation between simulations in the delay time using a normal distribution with standard deviation of 2 (truncating to make sure the delay does not become negative). Note that a delay distribution here does not mean that the delay time varies between individuals - it is identical) ;

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 
      value: 0.05
    delay: 
      value: 
        distribution: truncnorm
        mean: 7
        sd: 2
        a: 0
        b: Inf

Duration

By default, all outcome variables describe incidence (new individuals entering each day). However, they can also track an associated "prevalence" if the user specifies how long individuals will stay classified as the outcome state the outcome variable describes. This is the duration parameter ;

When the duration parameter is set, a new outcome variable is automatically created and named with the name of the original outcome variable + "_curr". This name can be changed using the duration::name option ;

For example, to track the incidence and prevalence of hospitalization when 5% of children are hospitalized, hospitalization occurs 7 days after infection, and the duration of hospitalization is 3 days:

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 
      value: 0.05
    delay: 
      value: 7
    duration: 
      value: 3

which creates the variable "incidH_child_curr" to track all currently hospitalized children. Since it doesn't make sense to call this new outcome variable an incidence, as it is a prevalence, we could instead rename it:

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 
      value: 0.05
    delay: 
      value: 7
    duration: 
      value: 3
      name: "hosp_child_curr"

Sum

Optional. sum is used to create new outcome variables that are sums over other previously defined outcome variables ;

If sum is included, source, probability, delay, and duration will be ignored ;

For example, to track new hospital admissions and current hospitalizations separately for children and adults, as well as for all ages combined

outcomes:
  incidH_child:
    source:
      incidence:
        infection_state: "I"
        age_group: "child"
    probability: 0.05
    delay: 6
    duration: 
      value: 14
      name: "hosp_child_curr"
  incidH_adult:
    source:
      incidence:
        infection_state: "I"
        age_group: "adult"
    probability: 0.01
    delay: 8
    duration:
      value: 7
      name: "hosp_adult_curr"
  incidH_total: 
    sum: ["incidH_child","incidH_adult"]
  hosp_curr_total:   
    sum: ["hosp_child_curr","hosp_adult_curr"]

outcomes::settings

There are other required and optional configuration items for the outcomes section which can be specified under outcomes::settings:

method: delayframe.This is the mathematical method used to create the outcomes variable values from the transmission model variables. Currently, the only model supported is delayframe, which .. ;

param_from_file: Optional, TRUE or FALSE. It is possible to allow any of the outcomes variables to have values that vary across the subpopulations. For example, disease severity rates or diagnosis rates may differ by demographic group. In this case, all the outcome parameter values defined in outcomes::outcomes will represent baseline values, and then you can define a relative change from this baseline for any particular subpopulation using the paths section. If params_from_file: TRUE is specified, then these relative values will be read from the params_subpop_file. Otherwise, if params_from_file: FALSE or is not listed at all, all subpopulations will have the same values for the outcome parameters, defined below ;

param_subpop_file: Required if params_from_file: TRUE. The path to a .csv or .parquet file that contains the relative amount by which a given outcome variable is shifted relative to baseline in each subpopulation. File must contain the following columns:

• subpop: The subpopulation for which the parameter change applies. Must be a subpopulation defined in the geodata file. For example, small_province

• parameter: The outcomes parameter which will be altered for this subpopulation. For example, incidH_child: probability

• value: The amount by which the baseline value will be multiplied, for example, 0.75 or 1.1

Examples

Consider a disease described by an SIR model in a population that is divided into two age groups, adults and children, which experience the disease separately. We are interested in comparing the predictions of the model to real world data, but we know we cannot observe every infected individual. Instead, we have two types of outcomes that are observed.

First, via syndromic surveillance, we have a database that records how many individuals in the population are experiencing symptoms from the disease at any given time. Suppose careful cohort studies have shown that 50% of infected adults and 80% of infected children will develop symptoms, and that symptoms occur in both age groups around 3 days after infection (following a log-normal distribution with log mean X and log standard deviation of Y). The duration that symptoms persist is also a variable, following a ...

Secondly, via laboratory surveillance we have a database of every positive test result for the infection. We assume the test is 100% sensitive and specific. Only individuals with symptoms are tested, and they are always tested exactly 1 day after their symptom onset. We are unsure what portion of symptomatic individuals are seeking out testing, but are interested in considering two extreme scenarios: 95% of symptomatic individuals are tested, or only 75% of individuals are tested.

The configuration file we could use to model this situation includes

// Some code

About configuration files

flepiMop is set up so that all parameters and other options for running the pipeline can be specified in a single "configuration" file (aka "config"). Users do not need to edit any other code files, or even be aware of their contents, to create and run complex model scenarios. Configuration files also provide a convenient record of model options and promote reproducibility of model results.

We use the YAML language syntax to write config files, which are typically named something like config.yml. The file has simple plain text contents and follows a tabbed outline structure. When config files are read by the model code, a data structure encoding the model options is created.

Comments can be added to the config file by starting with the hash key (#) then a space. Comments can start anywhere on a line and continue until the end, but if they run over to a new line, a new # must be used at the start of the new line.

Example

(give a simple configuration for a toy model with two subpopulations, SEIR, single "cases" outcome, single seeded infection, single NPI that starts after some time? this page is currently under development, please see our example repo _for some simple configurations) ;

When referring to config items (individual parameters), we use their full position in the outline. For example, in the sample config file above, we denote

subpop_setup:
  ...
  geodata: minimal

as subpop_setup::geodata having a value of minimal

Notation

Parameters and other options specified in the configuration files can take on a variety of types of values, using the following notations:

• dates are specified as [year]-[month]-[day]. (e.g., 2020-01-31)

• boolean values are either "TRUE" or "FALSE"

• files names are strings

• probability is a float between 0 and 1

• distribution is a probability distribution from which a random value for the parameter is drawn each time a new simulation is run (or chain, if doing inference). See here for the require schema.

Configuration files sections

Global header

These global configuration options typically sit at the top of the configuration file.

For example, for a configuration file to simulate the spread of COVID-19 in the US during 2020 and compare to data from March 1 onwards, with 1000 independent simulations, the header of the config might read:

name: USA_covid19_2020
model_output_dirname: model_output
start_date: 2020-01-01
end_date: 2020-12-31
start_date_groundtruth: 2020-03-01
end_date_groundtruth: 2020-12-31
nslots: 1000

subpop_setup section

This section specifies the population structure on which the model will be simulated, including the names and sizes of each subpopulation and the connectivity between them. More details here.

compartments section

This section is where users can specify the variables (infection states) that will be tracked in the infectious disease transmission model. More details can be found here. The other details of the model are specified in the seir section, including transitions between these compartments (seir::transitions), the names of the parameters governing the transitions (seir::parameters), and the numerical method used to simulate the equations over time (seir::integration). The initial conditions of the model can be specified in the initial_conditions section, and any other inputs into the model from external populations or instantaneous transitions between states that occur at later times can be specified in the seeding section. ;

seir section

This section is where users can specify the details of the infectious disease transmission model they wish to simulate (e.g., SEIR). This model describes the allowed transitions (seir::transitions) between the compartments that were specified in the compartments section, the values of the parameters involved in these transitions (seir::parameters), and the numerical method used to simulate the equations over time (seir::integration). More details here. The initial conditions of the model can be specified in the separate initial_conditions section, and any other inputs into the model from external populations or instantaneous transitions between states that occur at later times can be specified in the seeding section. ;

initial_conditions section

This section is used to specify the initial conditions of the model, which define how individuals are distributed between the model compartments at the time the model simulation begins. Importantly, the initial conditions specify the time and location where infection is first introduced. If this section is omitted, default values are used. If users want to add infections to the population at later times, or add or remove individuals from compartments separately from the model rules, they can do so via the related seeding section. More details here ;

seeding section

This section is used to specify how individuals are instantaneously "seeded" from one compartment to another, where they then continue to be governed by the model equations. For example, this seeding could be used to represent importations of infected individuals from an outside population, mutation events that create new strains, or vaccinations that alter disease susceptibility. Seeding events can occur at any time in the simulation. The seeding section specifies the numeric values added to or removed from any compartment of the model. More details here ;

outcomes section

This section is where users can define new variables representing the observed quantities and how they are related to the underlying state variables in the model (e.g., the fraction of infections that are detected as cases). More details here ;

interventions section

This section is where users can specify time-varying changes to parameters governing either the infectious disease model or the observational model. More details here ;

inference section

This section is where users can specify the details of how the model is fit to data, including what data streams they will be included and which outcome variables they represent and the likelihood functions describing the probability of the data given the model. More details here. ;

We want to allow users to work with a wide variety of infectious diseases or, one infectious disease under a wide variety of modeling assumptions. To facilitate this, we allow the user to specify their compartmental model of disease dynamics via the configuration file.

We originally considered asking users to specify each compartment and transition manually. However, we quickly found that this created long, confusing configuration files, and so we created a shorthand to more succinctly specify both compartments and transitions between them. This works especially well for models where individuals are stratified by other properties (like age, vaccination status, etc.) in addition to their infection status.

The model is specified in two separate sections of the configuration file. In the compartments section, users define the possible states individuals can be categorized into. Then in the seir section, users define the possible transitions between states, the values of parameters that govern the rates of these transitions, and the numerical method used to simulate the model.

An example section of a configuration file defining a simple SIR model is below.

compartments:
  infection_stage: ["S", "I", "R"]
  
seir:
  transitions:
    # infection
    - source: [S]
      destination: [I]
      proportional_to: [[S], [I]]
      rate: [beta]
      proportion_exponent: 1
    # recovery
    - source: [I]
      destination: [R]
      proportional_to: [[I]]
      rate: [gamma]
      proportion_exponent: 1
  parameters:
    beta: 0.1
    gamma: 0.2
  integration:
     method: rk4
     dt: 1.00

Specifying model compartments (compartments)

The first stage of specifying the model is to define the infection states (variables) that the model will track. These "compartments" are defined first in the compartments section of the config file, before describing the processes that lead to transitions between them. The compartments are defined separately from the rest of the model because they are also used by the seeding section that defines initial conditions and importations.

For simple disease models, the compartments can simply be listed with whatever notation the user chooses. For example, for a simple SIR model, the compartments could be ["S", "I", "R"]. The config also requires that there be a variable name for the property of the individual that these compartments describe, which for example in this case could be infection_stage

compartments:
  infection_stage: ["S", "I", "R"]

Our syntax allows for more complex models to be specified without much additional notation. For example, consider a model of a disease that followed SIR dynamics but for which individuals could receive vaccination, which might change how they experience infection.

In this case we can specify compartments as the cross product of multiple states of interest. For example:

 compartments:
   infection_stage: ["S", "I", "R"]
   vaccination_status: ["unvaccinated", "vaccinated"]

Corresponds to 6 compartments, which the code internally converts to this data frame

infection_stage, vaccination_status, compartment_name
S,               unvaccinated,       S_unvaccinated
I,               unvaccinated,       I_unvaccinated
R,               unvaccinated,       R_unvaccinated
S,               vaccinated,         S_vaccinated
I,               vaccinated,         I_vaccinated
R,               vaccinated,         R_vaccinated

In order to more easily describe transitions, we want to be able to refer to a compartment by its components, but then use it by its compartment name.

If the user wants to specify a model in which some compartments are repeated across states but others are not, there will be pros and cons of how the model is specified. Specifying it using the cross product notation is simpler, less error prone, and makes config files easier to read, and there is no issue with having compartments that have zero individuals in them throughout the model. However, for very large models, extra compartments increase the memory required to conduct the simulation, and so having unnecessary compartments tracked may not be desired.

For example, consider a model of a disease that follows SI dynamics in two separate age groups (children and adults), but for which only adults receive vaccination, with one or two doses of vaccine. With the simplified notation, this model could be specified as:

 compartments:
   infection_stage: ["S", "I"]
   age_group: ["child", "adult"]
   vaccination_status: ["unvaccinated", "1dose", "2dose"]

corresponding to 12 compartments, 4 of which are unnecessary to the model

infection_stage, age_group, vaccination_status, compartment_name
S,		 child,	    unvaccinated,	S_child_unvaccinated	
I,		 child,	    unvaccinated,	I_child_unvaccinated
S,		 adult,	    unvaccinated,	S_adult_unvaccinated
I,		 adult,	    unvaccinated,	I_adult_unvaccinated
S,		 child,	    1dose,		S_child_1dose
I,		 child,	    1dose,		I_child_1dose
S,		 adult,     1dose,		S_adult_1dose
I,		 adult,     1dose,		I_adult_1dose
S,		 child,     2dose,		S_child_2dose	
I,		 child,     2dose,		I_child_2dose
S,		 adult,	    2dose,		S_adult_2dose
I,		 adult,	    2dose,		I_adult_2dose

Or, it could be specified with the less concise notation

compartments:
   overall_state: ["S_child", "I_child", "S_adult_unvaccinated", "I_adult_unvaccinated", "S_adult_1dose", "I_adult_1dose", "S_adult_2dose", "I_adult_2dose"]

which does not result in any unnecessary compartments being included.

These compartments are referenced in multiple different subsequent sections of the config. In the seeding (LINK TBA) section the user can specify how the initial (or later imported) infections are distributed across compartments; in the seir section the user can specify the form and rate of the transitions between these compartments encoded by the model; in the outcomes section the user can specify how the observed variables are generated from the underlying model states.

Notation must be consistent between these sections.

Specifying compartmental model transitions (seir::transitions)

The way we specify transitions between compartments in the model is a bit more complicated than how the compartments themselves are specified, but allows users to specify complex stratified infectious disease models with minimal code. This makes checking, sharing, and updating models more efficient and less error-prone.

We specify one or more transition globs, each of which corresponds to one or more transitions. Since transition globs are shorthand for collections of transitions, we will first explain how to specify a single transition before discussing transition globs.

A transition has 5 pieces of associated information that a user can specify:

• source

• destination

• rate

• proportional_to

• proportion_exponent

For more details on the mathematical forms possible for transitions in our models, read the Model Description section.

We first consider a simple example of an SI model where individuals may either be vaccinated (v) or unvaccinated (u), but the vaccine does not change the susceptibility to infection nor the infectiousness of infected individuals.

We will focus on describing the first transition of this model, the rate at which unvaccinated individuals move from the susceptible to infected state.

Specifying a single transition

Source

The compartment the transition moves individuals out of (e.g., the source compartment) is an array. For example, to describe a transition that moves unvaccinated susceptible individuals to another state, we would write

[S,unvaccinated]

which corresponds to the compartment S_unvaccinated

Destination

The compartment the transition moves individuals into (e.g. the destination compartment) is an array. For example, to describe a transition that moves individuals into the unvaccinated but infected state, we would write

[I,unvaccinated]

which corresponds to the compartment I_unvaccinated

Rate

The rate constant specifies the probability per time that an individual in the source compartment changes state and moves to the destination compartment. For example, to describe a transition that occurs with rate 5/time, we would write:

5

instead, we could describe the rate using a parameter beta, which can be given a numeric value later:

beta

The interpretation and unit of the rate constant depend on the model details, as the rate may potentially also be per number (or proportion) of individuals in other compartments (see below).

Proportional to

A vector of groups of compartments (each of which is an array) that modify the overall rate of transition between the source and destination compartment. Each separate group of compartments in the vector are first summed, and then all entries of the vector are multiplied to get the rate modifier. For example, to specify that the transition rate depends on the product of the number of unvaccinated susceptible individuals and the total infected individuals (vaccinated and unvaccinated), we would write:

[[[S,unvaccinated]], [[I,unvaccinated], [I, vaccinated]]]

To understand this term, consider the compartments written out as strings

[[S_unvaccinated], [I_unvaccinated, I_vaccinated]]

and then sum the terms in each group

[S_unvaccinated, I_unvaccinated + I_vaccinated]

From here, we can say that the transition we are describing is proportional to S_unvaccinated and I_unvaccinated + I_vaccinated, i.e., the rate depends on the product S_unvaccinated * (I_unvaccinated + I_vaccinated).

For transitions that occur at a constant per-capita rate (ie, E -> I at rate $\gamma$ in an SEIR model), it is possible to simply write proportional_to: ["source"].

Proportion exponent

This is an exponent modifying each group of compartments that contribute to the rate. It is equivalent to the "order" term in chemical kinetics. For example, if the reaction rate for the model above depends linearly on the number of unvaccinated susceptible individuals but on the total infected individuals sub-linearly, for example to a power 0.9, we would write:

[1, 0.9]

or a power parameter alpha, which can be given a numeric value later:

[1, alpha]

The (top level) length of the proportion_exponent vector must be the same as the (top level) length of the proportional_to vector, even if the desire of the user is to have the same exponent for all terms being multiplied together to get the rate.

Summary

Putting it all together, the model transition is specified as

source: [S, unvaccinated]
destination: [I, unvaccinated]
proportional_to: [[[S,unvaccinated]], [[I,unvaccinated], [I,vaccinated]]]
rate: [5]
proportion_exponent: [1, 0.9]

would correspond to the following model if expressed as an ordinary differential equation

with parameter and parameter (we will describe how to use parameter symbols in the transitions and specify their numeric values separately in the section Specifying compartmental model parameters).

Transition globs

We now explain a shorthand we have developed for specifying multiple transitions that have similar forms all at once, via transition globs. The basic idea is that for each component of the single transitions described above where a term corresponded to a single model compartment, we can instead specify one or more compartment. Similarly, multiple rate values can be specified at once, for each involved compartment. From one transition glob, multiple individual transitions are created, by broadcasting across the specified compartments.

For transition globs, any time you could specify multiple arguments as a list, you may instead specify one argument as a non-list, which will be used for every broadcast. So [1,1,1] is equivalent to 1 if the dimension of that broadcast is 3.

We continue with the same SI model example, where individuals are stratified by vaccination status, but expand it to allow infection to occur at different rates in vaccinated and unvaccinated individuals:

Source

We allow one or more arguments to be specified for each compartment. So to specify the transitions out of both susceptible compartments (S_unvaccinated and S_unvaccinated), we would use

[[S], [unvaccinated,vaccinated]]

Destination

The destination variable should be the same shape as the source, and in the same relative order. So to specify a transition from S_unvaccinated to I_unvaccinated and S_vaccinated to I_vaccinated, we would write the destination as:

[[I], [unvaccinated,vaccinated]]

If instead we wrote:

[[I], [vaccinated,unvaccinated]]

we would have a transition from S_unvaccinated to I_vaccinated and S_vaccinated to I_unvaccinated.

Rate

The rate vector allows users to specify the rate constant for all the source -> destination transitions that are defined in a shorthand way, by instead specifying how the rate is altered depending on the compartment type. For example, the rate of transmission between a susceptible (S) and an infected (I) individual may vary depending on whether the susceptible individual is vaccinated or not and whether the infected individual is vaccinated or not. The overall rate constant is constructed by multiplying together or "broadcasting" all the compartment type-specific terms that are relevant to a given compartment.

For example,

rate: [[3], [0.6,0.5]]

This would mean our transition from S_unvaccinated to I_unvaccinated would have a rate of 3 * 0.6 while our transition from S_vaccinated to I_vaccinated would have a rate of 3 * 0.5.

The rate vector should be the same shape as source and destination and in the same relative order.

Note that if the desire is to make a model where the difference in the rate constants varies in a more complicated than multiplicative way between different compartment types, it would be better to specify separate transitions for each compartment type instead of using this shorthand.

Proportional to

The broadcasting here is a bit more complicated. In other cases, each broadcast is over a single component. However, in this case, we have a broadcast over a group of components. We allow a different group to be chosen for each broadcast.

[
  [[S,unvaccinated], [S,vaccinated]],
  [[I,unvaccinated],[I, vaccinated]], [[I,unvaccinated],[I, vaccinated]]
]

Again, let's unpack what it says. Since the broadcast is over groups, let's split the config back up

into those groups

[
  [S,unvaccinated],
  [[I,unvaccinated],[I, vaccinated]]
]
[
  [S,vaccinated],
  [[I,unvaccinated],[I, vaccinated]]
]

From here, we can say that we are describing two transitions. Both occur proportionally to the same compartments: S_unvaccinated and the total number of infections (I_unvaccinated+I_vaccinated).

If, for example, we want to model a situation where vaccinated susceptibles cannot be infected by unvaccinated individuals, we would instead write:

[
  [[S,unvaccinated], [S,vaccinated]],
  [[I,unvaccinated],[I, vaccinated]], [[I, vaccinated]]
]

Proportion exponent

Similarly to rate and proportional_to, we provide an exponent for each component and every group across the broadcast. So we could for example use:

[[1,1], [0.9,0.8]]

The (top level) length of the proportion_exponent vector must be the same as the (top level) length of the proportional_to vector, even if the desire of the user is to have the same exponent for all terms being multiplied together to get the rate. Within each vector entry, the arrays must have the same length as the source and destination vectors.

Summary

Putting it all together, the transition glob

seir:
  transitions:
    source: [[S],[unvaccinated,vaccinated]]
    destination: [[I],[unvaccinated,vaccinated]]
    proportional_to: [
                       [[S,unvaccinated], [S,vaccinated]],
                       [[I,unvaccinated],[I, vaccinated]], [[I, vaccinated]]
                     ]
    rate: [[3], [0.6,0.5]]
    proportion_exponent: [[1,1], [0.9,0.8]]

is equivalent to the following transitions

seir:
  transitions:
    - source: [S,unvaccinated]
      destination: [I,unvaccinated]
      proportional_to: [[[S,unvaccinated]], [[I,unvaccinated],[I, vaccinated]]]
      proportion_exponent: [1 * 0.9]
      rate: [3*0.6]
    - source: [S,vaccinated]
      destination: [I,vaccinated]
      proportional_to: [[[S,vaccinated]], [[I, vaccinated]]]
      proportion_exponent: [1 * 0.8]
      rate: [3*0.5]

Warning

We warn the user that with this shorthand, it is possible to specify large models with few lines of code in the configuration file. The more compartments and transitions you specify, the longer the model will take to run, and the more memory it will require.

Specifying compartmental model parameters (seir::parameters)

When the transitions of the compartmental model are specified as described above, they can either be entered as numeric values (e.g., 0.1) or as strings which can be assigned numeric values later (e.g., beta). We recommend the latter method for all but the simplest models, since parameters may recur in multiple transitions and so that parameter values may be edited without risk of editing the model structure itself. It also improves readability of the configuration files.

Parameters can take on three types of values:

• Fixed values

• Value drawn from distributions

• Values read from timeseries specified in a data file

Specifying fixed parameter values

Parameters can be assigned values by using the value argument after their name and then simply stating their numeric argument. For example, in a config describing a simple SIR model with transmission rate $\beta$ (beta) = 0.1/day and recovery rate $\gamma$ (gamma) = 0.2/day. This could be specified as

seir:
  parameters:
    beta: 
      value: 0.1
    gamma: 
      value: 0.2

The full model section of the config could then read

compartments:
  infection_state: ["S", "I", "R"]
  
seir:
  transitions:
    # infection
    - source: [S]
      destination: [I]
      proportional_to: [[S], [I]]
      rate: [beta]
      proportion_exponent: 1
    # recovery
    - source: [I]
      destination: [R]
      proportional_to: [[I]]
      rate: [gamma]
      proportion_exponent: [1,1]
  parameters:
    beta: 
      value: 0.1
    gamma: 
      value: 0.2

For the stratified SI model described above, this portion of the config would read

compartments:
  infection_stage: ["S", "I", "R"]
  vaccination_status: ["unvaccinated", "vaccinated"]
  
seir:
  transitions:
    source: [[S],[unvaccinated,vaccinated]]
    destination: [[I],[unvaccinated,vaccinated]]
    proportional_to: [
                       [[S,unvaccinated], [S,vaccinated]],
                       [[I,unvaccinated],[I, vaccinated]], [[I, vaccinated]]
                     ]
    rate: [[beta], [theta_u,theta_v]]
    proportion_exponent: [[1,1], [alpha_u,alpha_v]]
  parameters:
    beta: 
      value: 0.1
    theta_u: 
      value: 0.6
    theta_v: 
      value: 0.5
    alpha_u: 
      value: 0.9
    alpha_v: 
      value: 0.8

If there are no parameter values that need to be specified (all rates given numeric values when defining model transitions), the seir::parameters section of the config can be left blank or omitted.

Specifying parameters values from distributions

Parameter values can also be specified as random values drawn from a distribution, as a way of including uncertainty in parameters in the model output. In this case, every time the model is run independently, a new random value of the parameter is drawn. For example, to choose the same value of beta = 0.1 each time the model is run but to choose a random values of gamma with mean on a log scale of $e^{-1.6} = 0.2$ and standard deviation on a log scale of $e^{0.2} = 1.2$ (e.g., 1.2-fold variation):

seir:
  parameters:
    beta: 
      value:
        distribution: fixed
        value: 0.1
    gamma: 
      value:
        distribution: lognorm
        logmean: -1.6
        logsd: 0.2

Details on the possible distributions that are currently available, and how to specify their parameters, is provided in the Distributions section.

Note that understanding when a new parameter values from this distribution is drawn becomes more complicated when the model is run in Inference mode. In Inference mode, we distinguish model runs as occurring in different "slots" – i.e., completely independent model instances that could be run on different processing cores in a parallel computing environment – and different "iterations" of the model that occur sequentially when the model is being fit to data and update fitted parameters each time based on the fit quality found in the previous iteration. A new parameter values is only drawn from the above distribution once per slot. Within a slot, at each iteration during an inference run, the parameter is only changed if it is being fit and the inference algorithm decides to perturb it to test a possible improved fit. Otherwise, it would maintain the same value no matter how many times the model was run within a slot.

Specifying parameter values as timeseries from data files

Sometimes, we want to be able to specify model parameters that have different values at different timepoints. For example, the relative transmissibility may vary throughout the year based on the weather conditions, or the rate at which individuals are vaccinated may vary as vaccine programs are rolled out. One way to do this is to instead specify the parameter values as a timeseries.

This can be done by providing a data file in .csv or .parquet format that has a list of values of the parameter for a corresponding timepoint and subpopulation name. One column should be date, which should have an entry for every calendar day of the simulation, with the first and last date corresponding to the start_date and end_date for the simulation specified in the header of the config. There should be another column for each subpopulation, where the column name is the subpop name used in other files and the values are the desired parameter values for that subpopulation for the corresponding day. If any day or subpopulation is missing, an error will occur. However, if you want all subpopulations to have the same parameter value for every day, then only a single column in addition to date is needed, which can have any name, and will be applied to every subpop ;

For example, for an SIR model with a simple two-province population structure where the relative transmissibility peaks on January 1 then decreases linearly to a minimal value on June 1 then increases linearly again, but varies more in the small province than the large province, the theta parameter could be constructed from the file seasonal_transmission_2pop.csv with contents including

date,        small_province,    large_province
2022-01-01,  1.5,               1.3
.....
2022-05-01,  0.5,               0.7 
....
2022-12-31,  1.5,               1.3

as a part of a configuration file with the model sections:

compartments:
  infection_stage: ["S", "I", "R"]

seir:
  transitions:
    # infection
    - source: [S]
      destination: [I]
      proportional_to: [[S], [I]]
      rate: [beta*theta]
      proportion_exponent: 1
    # recovery
    - source: [I]
      destination: [R]
      proportional_to: [[I]]
      rate: [gamma]
      proportion_exponent: 1
  parameters:
    beta: 
      value: 0.1
    gamma: 
      value: 0.2
    theta:
       timeseries: data/seasonal_transmission.csv

Note that there is an alternative way to specify time dependence in parameter values that is described in the Specifying time-varying parameter modifications section. That method allows the user to define intervention parameters that apply specific additive or multiplicative shifts to other parameter values for a defined time interval. Interventions are useful if the parameter doesn't vary frequently and if the values of the shift is unknown and it is desired to either sample over uncertainty in it or try to estimate its value by fitting the model to data. If the parameter varies frequently and its value or relative value over time is known, specifying it as a timeseries is more efficient.

Compartmental model parameters can have an additional attribute beyond value or timeseries, which is called stacked_modifier_method. This value is explained in the section on coding time-dependent parameter modifications (also known as "modifiers") as it determines what happens when two different modifiers act on the same parameter at the same time (are they combined additively or multiplicatively?) ;

Specifying model simulation method (seir::integration)

A compartmental model defined using the notation in the previous sections describes rules for classifying individuals in the population based on infection state dynamically, but does not uniquely specify the mathematical framework that should be used to simulate the model.

Our framework allows for two major methods for implementing compartmental models of disease transmission:

• ordinary differential equations, which are completely deterministic, operate in continuous time (consider infinitesimally small timesteps), and allow for arbitrary fractions of the population (i.e., not just discrete individuals) to move between model compartments

• discrete-time stochastic process, which tracks discrete individuals and produces random variation in the number of individuals transitioning between states for any given rate, and which allows transitions between states only to occur at discrete time intervals

The mathematics behind each implementation is described in the Model Description section

For example, to simulate a model deterministically using the 4th order Runge-Kutta algorithm for numerical integration with a timestep of 1 day:

seir:
  integration:
     method: rk4
     dt: 1.00

Alternatively, to simulate a model stochastically with a timestep of 0.1 days

seir:
  integration:
     method: stochastic
     dt: 0.1

For any method, the results of the model will be more accurate when the timestep is smaller (i.e., output will more precisely match the mathematics of the model description and be invariant to the choice of timestep). However, the computing time required to simulate the model for a certain time range of interest increases with the number of timesteps required (i.e., with smaller timesteps). In our experience, the 4th order Runge-Kutta algorithm (for details see Advanced section) is a very accurate method of numerically integrating such models and can handle timesteps as large as roughly a day for models with the maximum per capita transition rates in this same order of magnitude. However, the discrete time stochastic model or the legacy method for integrating the model in deterministic mode require smaller timesteps to be accurate (around 0.1 for COVID-19-like dynamics in our experience.

inference settings

iterations_per_slot

do_inference

gt_data_path

With inference model runs, the number of simulations nsimulations refers to the number of final model simulations that will be produced. The filtering$simulations_per_slot setting refers to the number of iterative simulations that will be run in order to produce a single final simulation (i.e., number of simulations in a single MCMC chain).

f

inference::statistics options

required options

name

aggregator

period

sim_var

data_var

likelihood

The statistics specified here are used to calibrate the model to empirical data. If multiple statistics are specified, this inference is performed jointly and they are weighted in the likelihood according to the number of data points and the variance of the proposal distribution.

f

optional options ?

remove_na

add_one

gt_start_date

gt_end_date

Optional sections

inference::hierarchical_stats_geo

The hierarchical settings specified here are used to group the inference of certain parameters together (similar to inference in "hierarchical" or "fixed/group effects" models). For example, users may desire to group all counties in a given state because they are geograhically proximate and impacted by the same statewide policies. The effect should be to make these inferred parameters follow a normal distribution and to observe shrinkage among the variance in these grouped estimates.

inference::priors

It is now possible to specify prior values for inferred parameters. This will have the effect of speeding up model convergence.

Ground truth data

name

module

geo_group_col

transform

inference:::priors

inference::

Welcome to flepiMoP documentation!

The “FLexible EPIdemic MOdeling Pipeline” (flepiMoP; formerly known as the COVID Scenario Modeling Pipeline or CSP) is an open-source software suite designed by researchers in the Johns Hopkins Infectious Disease Dynamics Group and at UNC Chapel Hill to simulate a wide range of compartmental models of infectious disease transmission. The disease transmission and observation models are defined by a no-code configuration file, which allows models of varying complexity to be specified quickly and consistently, from simple problems described by SIR-style models in a single population to more complicated models of multiple pathogen strains transmitting between thousands of connected spatial divisions and age groups.

It was initially designed in early 2020 and was routinely used to provide projections of the emerging COVID-19 epidemic to health authorities worldwide. Currently, flepiMoP provides COVID-19 projections to the US CDC-funded model aggregation sites, the COVID-19 Forecast Hub and the COVID-19 Scenario Modeling Hub, influenza projections to FluSight and to the Flu Scenario Modeling Hub, and RSV projections to the RSV Scenario Modeling Hub.

However, the pipeline is much more general and can be used to simulate the dynamics of any infection that can be expressed as a compartmental epidemic model. These include applications in chemical reaction kinetics, pharmacokinetics, within-host disease dynamics, or applications in the social sciences.

In addition to producing forward simulations given a specified model and parameter values, the pipeline can also attempt to optimize unknown parameters (e.g., transmission rate, case detection rate, intervention efficacy) to fit the model to datasets the user provides (e.g., hospitalizations due to severe disease) using a Bayesian inference framework. This feature allows the pipeline to be utilized for short-term forecasting or longer-term scenario projections for ongoing epidemics, since it can simultaneously be fit to data for dates in the past and then use best-fit parameters to make projections into the future.

General description of flepiMoP

The main features of flepiMoP are:

• Open-source (GPL v3.0) infectious dynamics modeling software, written in R and Python

• Versatile, no-code design applicable for most compartmental models and outcome observation models, allowing for quick iteration in reaction to epidemic events (e.g., emergence of new variants, vaccines, non-pharmaceutical interventions (NPIs))

• Powerful, just-in-time compiled disease transmission model and distributed inference engine ready for large scale simulations on high-performance computing clusters or cloud workflows

• Adapted to small- and large-scale problems, from a simple SIR model to a complex model structure with hundreds of compartments on thousands of connected populations

• Strong emphasis on mechanistic processes, with a design aimed at leveraging domain knowledge in conjunction with statistical inference

• Portable for Windows WSL, MacOS, and Linux with the provided Docker image and an Anaconda environment

The mathematical model within the pipeline is a compartmental epidemic model embedded within a well-mixed metapopulation. A compartmental epidemic model is a model that divides all individuals in a population into a discrete set of states (e.g. “infected”, “recovered”) and tracks – over time – the number of individuals in each state and the rates at which individuals transition between these states. The well-known SIR model is a classic example of such a model, and much more complex versions of this model type have been simulated with this framework (for example, an SEIR-style model in which individuals are further subdivided into multiple age groups and vaccination statuses).

The structure of the desired model, as well as the parameter values and initial conditions, can be specified flexibly by the user in a no-code fashion. The pipeline allows for parameter values to change over time at discrete intervals, which can be used to specify time-dependent aspects of disease transmission and control (such as seasonality or vaccination campaigns).

The model is embedded within a meta-population structure, which consists of a series of distinct subpopulations (e.g. states, provinces, or other communities) in which the model structure is repeated, albeit with potentially different parameter values. The subpopulations can interact, either through the movement of individuals or the influence of individuals in one subpopulation on the transition rate of individuals in another ;

Within each subpopulation, the population is assumed to be well-mixed, meaning that interactions are assumed to be equally likely between any pair of individuals (since unique identities of individuals are not explicitly tracked). The same model structure can be simulated in a continuous-time deterministic or discrete-time stochastic manner ;

In addition to the variables described by the compartmental model, the model can track other observable variables (“outcomes”) that are functions of the basic model variables but do not themselves influence the dynamics (i.e., some portion of infections are reported as cases, depending on a testing rate). The model can be run iteratively to tune the values of certain parameters so that these outcome variables best match timeseries data provided by the user for a certain time period ;

Fitting is done using a Bayesian-like framework, where the user can specify the likelihood of observed outcomes in data given modeled outcomes, and the priors on any parameters to be fit. Multiple data streams (e.g., cases and deaths) can be fit simultaneously. A custom Markov Chain Monte Carlo method is used to sequentially propose and accept or reject parameter values based on the model fit to data, in a way that balances fit quality within each individual subpopulation with that of the total aggregate population, and that takes advantage of parallel computing environments.

The code is written in a combination of R and Python, and the vast majority of users only need to interact with the pipeline via the components written in R. It is structured in a modular fashion, such that individual components – such as the epidemic model, the observable variables, the population structure, or the parameters – can be edited or completely replaced without any handling of other parts of the code ;

When model simulation is combined with fitting to data, the code is designed to run most efficiently on a supercomputing cluster with many cores. We most commonly run the code on Amazon Web Services or on high-performance computers using SLURM. However, even relatively large models can be run efficiently on most personal computers. Typically, the memory of the machine will limit the number of compartments (i.e., variables) that can be included in the epidemic model, while the machine’s CPU will determine the speed at which each model run is completed and the number of iterations of the model that can be run during parameter searches when fitting the model to data. While the pipeline can be installed on any computer, it is sometime easier to use an Anaconda environment or the provided Docker container, where all the software dependencies (e.g., standardized R and Python versions along with required packages) are included, independent of the user’s local machine. All the code is maintained on our GitHub and shared with the GNU General Public License v3.0 license. It is build on top of a fully open-source software stack.

This documentation is organized as follows. The Model Description section describes the mathematical framework for the compartmental epidemic models that can be simulated forward in time by the pipeline. The Model Inference section describes the statistical framework for fitting the model to data. The Data and Parameter section describes the inputs the user must provide to the pipeline, in terms of the model structure and parameters, the population characteristics, the initial conditions, time-varying interventions, data to be fit, and more. The How to Run section provides concrete guidance on setting up and running the model and analyzing the output. The Quick Start Guide provides a simple example model setup. The Advanced section goes into more detail on specific features of the model and the code that are likely to only be of interest to users who want to run more complex models or data fitting routines or substantially edit the code. It includes a subsection describing each file and package used in the pipeline and their interactions during a model run.

Users who wish to jump to running the model themselves can see Quick Start Guide.

For questions about the pipeline or to report a bug, please use the “Issues” or "Discussions" feature on our GitHub.

Acknowledgments

flepiMoP is actively developed by its current contributors, including Joseph C Lemaitre, Sara L Loo, Emily Przykucki, Clifton McKee, Claire Smith, Sung-mok Jung, Koji Sato, Pengcheng Fang, Erica Carcelen, Alison Hill, Justin Lessler, and Shaun Truelove, affiliated with the ;

• Department of Epidemiology, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA for (JCL, JL)

• Johns Hopkins University International Vaccine Access Center, Department of International Health, Baltimore, MD, USA for (SLL, KJ, EC, ST)

• Department of Epidemiology, Johns Hopkins Bloomberg School of Public Health, Baltimore, Maryland, USA for (CM, CS, JL, ST)

• Carolina Population Center, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA for (S-m.J, JL)

• Institute for Computational Medicine, Johns Hopkins University, Baltimore, MD, USA for (AH).

The development of this model was supported by from funds the National Science Foundation (2127976; ST, CPS, JK, ECL, AH), Centers for Disease Control and Prevention (200-2016-91781; ST, CPS, JK, AH, JL, JCL, SL, CM, EC, KS, S-m.J), US Department of Health and Human Services / Department of Homeland Security (ST, CPS, JK, ECL, AH, JL), California Department of Public Health (ST, CPS, JK, ECL, JL), Johns Hopkins University (ST, CPS, JK, ECL, JL), Amazon Web Services (ST, CPS, JK, ECL, AH, JL, JCL), National Institutes of Health (R01GM140564; JL, 5R01AI102939; JCL), and the Swiss National Science Foundation (200021-172578; JCL)

We need to also acknowledge past contributions to the development of the COVID Scenario Pipeline, which evolved into flepiMoP. These include contributions by Heramb Gupta, Kyra H. Grantz, Hannah R. Meredith, Stephen A. Lauer, Lindsay T. Keegan, Sam Shah, Josh Wills, Kathryn Kaminsky, Javier Perez-Saez, Joshua Kaminsky, and Elizabeth C. Lee.

Config changes w.r.t classical inference

You need, under inference, to add method: emcee and modify the statistics: as shown in the diff below (basically: all resampling goes to one subsection, with some minor changes to names).

To see which llik options and regularization (e.g do you want to weigh more the last weeks for forecasts, or do you want to add the sum of all subpop) see files statistics.py.

Test on your computer

Install gempyor from branch emcee_batch . Test your config by running:

on your laptop. If it works, it should produce:

• plots of simulation directly from your config

• plots after the fits with the fits and the parameter chains

• and h5 file with all the chains

• and in model_output, the final hosp/snpi/seir/... files in the flepiMoP structure.

It will output something like

```

Here, it says the config fits 92 parameters, we'll keep that in mind and choose a number of walkers greater than (ideally 2 times) this number of parameters.

Run on cluster

Install gempyor on the cluster. test it with the above line, then modify this script:

so you need to have:

• -c (number of core) equal to roughly half the number of walkers (slots/parallel chains)

• mem to be around two times the number of walkers. Look at the computes nodes you have access to and make something that can be prioritized fast enough.

• nsamples is the number of final results you want, but it's fine not to care about it, I rerun the sampling from my computer.

• To resume from an existing run, add the previous line --resume and it 'll start from the last parameter values in the h5 files.

Postprocess EMCEE

To analyze run postprocessing/emcee_postprocess.ipynb First, this plots the chains and then it runs nsamples (you can choose it) projection with the end of the chains and does the plot of the fit, with and without projections

Updates to other files

\LLIK (inference runs only)

During inference runs, an additional file type, llik, is created, which is described in the section ;

These files contain the log-likelihoods of the model simulation for each subpopulation, as well as some diagnostics on acceptance.

The meanings of the columns are:

ll - These values are the log-likelihoods of data given the model and parameter values for a single subpopulation (in subpop column) ;

filename - ...

subpop - The values of this column are the names of the nodes from the geodata file.

accept - Either 0 or 1, depending on whether the parameters during this iteration of the simulation were accepted (1) or rejected (0) in that subpopulation ;

accept_avg ;

accept_prob ;

For inference runs, ... flepiMoP produces one file per parallel slot, for both global and chimeric outputs...

Some information to consider if you'd like your script to be run automatically after an inference run ;

• Most R/python packages are installed already installed. Try to run your script on the conda environment defined on the (or easier if you are not set up on MARCC, ask me)

• There will be some variables set in the environment. These variables are:

  • $CONFIG_PATH the path to the configuration fil ;

  • $FLEPI_RUN_INDEX the run id for this run (e.g `CH_R3_highVE_pesImm_2022_Jan29`

  • $JOB_NAME this job name (e.g USA-20230130T163847_inference_med)

  • $FS_RESULTS_PATH the path where lies the model results. It's a folder that contains the model_ouput/ as a subfolder

  • $FLEPI_PATH path of the flepiMoP repository.

  • $DATA_PATH path of the Data directory (e.g Flu_USA or COVID19_USA).

  • Anything you ask can theoretically be provided here.

• The script must run without any user intervention.

• The script is run from $DATA_PATH.

• Your script lies in the flepiMoP directory (preferably) or it's ok if it is in a data directory if it makes sense ;

• It is run on a 64Gb of RAM multicore machine. All scripts combined must complete under 4 hours, and you can use multiprocessing (48 cores)

• Outputs (pdf, csv, html, txt, png ...) must be saved in a directory named pplot/ (you can assume that it exists) in order to be sent to slack by FlepiBot 🤖 after the run.

• an example postprocessing script (in python) is .

• You can test your script on MARCC on a run that is already saved in /data/struelo1/flepimop-runs or I can do it for you.

• Once your script works, add (or ask to add) the command line to run in file batch/postprocessing_scripts.sh between the START and END lines, with a little comment about what your script does.

filtering section

The filtering section configures the settings for the inference algorithm. The below example shows the settings for some typical default settings, where the model is calibrated to the weekly incident deaths and weekly incident confirmed cases for each subpop. Statistics, hierarchical_stats_geo, and priors each have scenario names (e.g., sum_deaths, local_var_hierarchy, and local_var_prior, respectively).

filtering settings

With inference model runs, the number of simulations nsimulations refers to the number of final model simulations that will be produced. The filtering$simulations_per_slot setting refers to the number of iterative simulations that will be run in order to produce a single final simulation (i.e., number of simulations in a single MCMC chain).

filtering::statistics

The statistics specified here are used to calibrate the model to empirical data. If multiple statistics are specified, this inference is performed jointly and they are weighted in the likelihood according to the number of data points and the variance of the proposal distribution.

filtering::hierarchical_stats_geo

The hierarchical settings specified here are used to group the inference of certain parameters together (similar to inference in "hierarchical" or "fixed/group effects" models). For example, users may desire to group all counties in a given state because they are geograhically proximate and impacted by the same statewide policies. The effect should be to make these inferred parameters follow a normal distribution and to observe shrinkage among the variance in these grouped estimates.

filtering::priors

It is now possible to specify prior values for inferred parameters. This will have the effect of speeding up model convergence.

Ground truth data

Likelihood function

Fitting parameters

Ground truth data