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 , defining the compartments and the transitions between compartments

• An , 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 , , and are time-dependent variables describing the number of individuals in each state, is a rate parameter (units of time) and is a scaling parameter (unitless). may be , , 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 is in fact and , are ,).

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 section.

Alternatively, the model can be simulated as a discrete-time stochastic process, where the number of individuals transitioning between states and at time is a binomial random variable

where the second term is the expected fraction of individuals in the state at time who would transition to by time if there were no other changes to in this time, and time step is a chosen parameter that must be small for equivalence between continuous- and discrete-time versions of the model.

SEIR model

For example, an SEIR model – which describes an infection for which susceptible individuals () who are infected first pass through a latent or exposed () phase before becoming infectious () and that confers perfect lifelong immunity after recovery () – could be encoded as

where is the transmission rate (rate of infectious contact per infectious individual), is the rate of progression ( is the average latent/incubation period), and is the recovery rate ( is the average duration of the infectious period), and is the total population size (). In differential equation form, this model is

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

2. an infection rate modified by a 'mixing coefficient', , which is a rough heuristic for the slowdown in disease spread that occurs in realistically heterogeneous populations where more well-connected individuals are infected first.

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

where is the transmission rate between age X and Y, and we have assumed individuals do not age on the timescale relevant to the 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)

where is the vaccination rate (we assume that individuals do not receive the vaccine while they are exposed or infectious) and is the vaccine efficacy against infection. Similar structures could be used for other sources of prior immunity or other dynamic risk groups.

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

where is the immune cross-protection conferred from infection with strain A to subsequent infection with strain B. Co-infection is ignored. All individuals are assumed to be initially equally susceptible to both infections and are just categorized as (vs ) for convenience.

All combinations of these situations can be quickly specified in flepiMoP. Details on how to encode these models is provided in the section, with examples given in the 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.

An outcome variable can be generated from a state variable of the mathematical model using the following properties:

• The proportion of all individuals in who will be observed as ,

• The delay between when an individual enters state and when they are observed as , which can follow a class of probability distributions where is the parameters of the distribution (e.g., the mean and standard deviation of a normal distribution)

• (optional) the duration spent in observable , in which case the output will also contain the prevalence (number of individuals currently in in addition to the incidence into

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

The number of individuals in at time who become part of the outcome variable is a random variable, and individuals who are observed in at time could have entered at different times in the past.

Formally, for a deterministic, continuous-time model

For a discrete-time, stochastic model

Note that outcomes constructed in this way always represent incidence values; meaning they describe the number of individuals newly entering this state at time . If the model state is also an incidence, then is a unitless probability, whereas if is a prevalence (number of individuals currently in state at time ), then is instead a probability per time unit.

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 section, with examples given in the 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)

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

where is the per-contact per-time rate of disease transmission between an infected individual residing in subpopulation and a susceptible individual from subpopulation .

In general for infection models in connected subpopulations, the transmission rates can take on arbitrary values. In this pipeline, however, we impose an additional structure on these terms. We assume that interactions between subpopulations occur when individuals temporarily relocate to another subpopulation, where they interact with locals. We call this movement “mobility”, and it could be due to regular commuting, special travel, etc. There is a transmission rate () associated with each subpopulation , and individuals physically in that subpopulation – permanently or temporarily – are exposed and infected with this local rate whenever they encounter local susceptible individuals.

The transmission matrix is then

where is the onward transmission rate from infected individuals in subpopulation , is the number of individuals in subpopulation i who are interacting with individuals in subpopulation at any given time (for example, fraction who commute each day), and is a fractional scaling factor for the strength of inter-population contacts (for example, representing the fraction of hours in a day commuting individuals spend outside vs. inside their subpopulation).

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 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

Initial conditions can be specified by setting the values of the compartments in the disease transmission model at time zero, or the start of the simulation. For example, we might assume that for day zero of an outbreak the whole population is susceptible except for one single infected individual, i.e. and . Alternatively, we might assume that a certain proportion of the population has prior immunity from previous infection or vaccination.

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 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.

For example, the rate of transmission in subpopulation , , may be reduced by an intervention that acts between times and , and another intervention that acts between times and

In this case, and are both considered simple SinglePeriodModifier interventions. There are four possible types of interventions that can be included in the model

• SinglePeriodModifier - an intervention that leads to a fractional reduction in a parameter value in subpopulation (i.e., ) between two timepoints

• MultiPeriodModifier - an intervention that leads to a fractional reduction in a parameter value in subpopulation (i.e., ) value between multiple sets of timepoints

• ModifierModifier- an intervention that leads to a fractional reduction in the value of another intervention between two timepoints

• StackedModifier - TBA

Overview

In order for the models specified previously to be dynamically simulated, the user must provide initial conditions, in addition to the model structure and parameter values. Initial conditions describe the value of each variable in the model at the time point that the simulation is to start. For example, on day zero of an outbreak, we may assume that the entire population is susceptible except for one single infected individual. Alternatively, we could assume that some portion of the population already has prior immunity due to vaccination or previous infection. Different initial conditions lead to different model trajectories.

The initial_conditions section of the configuration file is detailed below. Note that in some cases, can replace or complement the initial condition, the table below provides a quick comparison of these sections.

Specifying model initial conditions

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

initial_conditions:method Must be either "Default", "SetInitialConditions", or "FromFile".

initial_conditions:initial_conditions_fileRequired for methods “SetInitialConditions” and “FromFile” . Path to a .csv or .parquet file containing the list of initial conditions for each compartment.

initial_conditions:initial_file_type Only required for method: “FolderDraw”. Description TBA

initial_conditions::allow_missing_subpops Optional for all methods, determines what will happen if initial_conditions_file is missing values for some subpopulations. If FALSE, the default behavior, or unspecified, an error will occur if subpopulations are missing. If TRUE, then for subpopulations missing from the initial_conditions file, it will be assumed that all individuals begin in the first compartment (the “first” compartment depends on how the model was specified, and will be the compartment that contains the first named category in each compartment group), unless another compartment is designated to hold the rest of the individuals ;

initial_conditions::allow_missing_compartments Optional for all methods. If FALSE, the default behavior, or unspecified, an error will occur if any compartments are missing for any subpopulation. If TRUE, then it will be assumed there are zero individuals in compartments missing from the initial_conditions file.

initial_conditions::proportional If TRUE, assume that the user has specified all input initial conditions as fractions of the population, instead of numbers of individuals (the default behavior, or if set to FALSE). Code will check that initial values in all compartments sum to 1.0 and throw an error if not, and then will multiply all values by the total population size for that subpopulation ;

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

initial_conditions::method

Default

The default initial conditions are that the initial value of all compartments for each subpopulation will be zero, except for the first compartment, whose value will be the population size. The “first” compartment depends on how the model was specified, and will be the compartment that contains the first named category in each compartment group.

For example, a model with the following compartments

with the accompanying geodata file

will be started with 1000 individuals in the S_child_unvaxxed in the "small province" and 10,000 in that compartment in the "large province".

SetInitialConditions

With this method users can specify arbitrary initial conditions in a convenient formatted input .csv or .parquet file.

For example, for a model with the following compartments and initial_conditions sections

with the accompanying geodata file

where initial_conditions.csv contains

the model will be started with half of the population of both subpopulations, consisting of children and the other half of adults, everyone unvaccinated, and 5 infections (in exposed-but-not-yet-infectious class) among the unvaccinated adults in the large province, with the remaining individuals susceptible (4995). All other compartments will contain zero individuals initially ;

initial_conditions::initial_conditions_file must contain the following columns:

• subpop – the name of the subpopulation for which the initial condition is being specified. By default, all subpopulations must be listed in this file, unless the allow_missing_subpops option is set to TRUE.

• mc_name – the concatenated name of the compartment for which an initial condition is being specified. The order of the compartment groups in the name must be the same as the order in which these groups are defined in the config for the model, e.g., you cannot say unvaccinated_S.

For each subpopulation, if there are compartments that are not listed in SetInitialConditions, an error will be thrown unless allow_missing_compartments is set to TRUE, in which case it will be assumed there are zero individuals in them. If the sum of the values of the initial conditions in all compartments in a location does not add up to the total population of that location (specified in the geodata file), an error will be thrown. To allocate all remaining individuals in a subpopulation (the difference between the total population size and those allocated by defined initial conditions) to a single pre-specified compartment, include this compartment in the initial_conditions_file but instead of a number in the amount column, put the word "rest" ;

If allow_missing_subpops is FALSE or unspecified, an error will occur if initial conditions for some subpopulations are missing. If TRUE, then for subpopulations missing from the initial_conditions file, it will be assumed that all individuals begin in the first compartment. (The “first” compartment depends on how the model was specified, and will be the compartment that contains the first named category in each compartment group.)

FromFile

Similar to "SetInitialConditions", with this method users can specify arbitrary initial conditions in a formatted .csv or .parquet input file. However, the format of the input file is different. The required file format is consistent with the from the compartmental model, so the user could take output from one simulation and use it as input into another simulation with the same model structure ;

For example, for an input configuration file containing

with the accompanying geodata file

where initial_conditions_from_previous.csv contains

The simulation would be initiated on 2021-06-01 with these values in each compartment (no children vaccinated, only adults in the small province vaccinated, some past and current infection in both compartments but ).

initial_conditions::initial_conditions_file must contain the following columns:

• mc_value_type – in model output files, this is either prevalence or incidence. Prevalence values only are selected to be used as initial conditions, since compartmental models described the prevalence (number of individuals at any given time) in each compartment. Prevalence is taken to be the value measured instantaneously at the start of the day

• mc_name – The name of the compartment for which the value is reported, which is a concatenation of the compartment status in each state type, e.g. "S_adult_unvaxxed" and must be in the same order as these groups are defined in the config for the model, e.g., you cannot say unvaxxed_S_adult.

SetInitialConditionsFolderDraw, FromFileFolderDraw

The way that initial conditions is specified with SetInitialConditions and FromFile results in a single value for each compartment and does not easily allow the user to instead specify a distribution (like is possible for compartmental or outcome model parameters). If a user wants to use different possible initial condition values each time the model is run, the way to do this is to instead specify a folder containing a set of file with initial condition values for each simulation that will be run. The user can do this using files with the format described in initial_conditions::method::SetInitialConditions using instead method::SetInitialConditionsFolder draw. Similarly, to provide a folder of initial condition files with the format described in initial_conditions::method:FromFile using instead method::FromFileFolderDraw ;

Each file in the folder needs to be named according to the same naming conventions as the model output files: run_number.runID.file_type.[csv or parquet] where ....[DESCRIBE] as it is now taking the place of the seeding files the model would normally outpu ;

Only one additional config argument is needed to use a FolderDraw method for initial conditions:

initial_file_type: either seir or seed

When using FolderDraw methods, initial_conditions_file should now be the path to the directory that contains the folder with all the initial conditions files. For example, if you are using output from another model run and so the files are in an seir folder within a model_output folder which is in within your project directory, you would use initial_conditions_file: model_outpu ;

flepiMoP can be used to conduct forward simulations of a model with user-defined parameter values, or, it can be used to iteratively run a model with some unknown parameters, compare the model output to ground truth data, and find parameter values that optimize the fit of the model to data (i.e., conduct model "inference"). We have developed a custom model inference method that is based on standard Markov Chain Monte Carlo (MCMC)-based approaches to Bayesian inference for dynamic models, but is adapted to deal with some of the particular challenges of large-scale epidemic models, including i) long times and high computational resources required to simulate single model runs, ii) multiple subpopulations with location-specific parameters but inter-location transmission, iii) a high-dimensional parameter space, iv) the need to produce real-time epidemic projections, and v) the availability of parallel computing resources.

Notation

• – A set of unknown model parameters to be estimated by fitting the model output to data. For a model with subpopulations each with their own parameters, this set includes all location-specific parameters .

• – The timeseries output of one or more of the state variables of the model under parameters For simplicity, we will often just use the notation . The value at a timepoint is . For a model with subpopulations for which there are different state variables, this becomes . (Note that for the general case when the dynamics in one location can effect the dynamics in another, the model state in one location depends on the full set of parameters, not just the location-specific parameters.)

• – The timeseries for the observed data (also referred to as "ground truth") that the model attempts to recreate. For a model with subpopulations each with their own observed data for variable , this becomes .

• – The likelihood of the observed data being produced by the model for an input parameter set . This is a probability density function over all possible values of the data being produced by the model, conditional on a fixed model parameter value ;

• – The prior probability distribution, which in Bayesian inference encodes beliefs about the possible values of the unknown parameter before any observed data is formally compared to the model.

• – The posterior probability distribution, which in Bayesian inference describes the updated probability of the parameters conditional on the observed data .

• – The proposal density, used in Metropolis-Hastings algorithms for Markov Chain Monte Carlo (MCMC) techniques for sampling the posterior distribution, describes the probability of proposing a new parameter set from a current accepted parameter set .

Background

This section can be skipped by those familiar with Markov Chain Monte Carlo approaches to Bayesian inference.

Bayesian inference

Our model fitting framework is based on the principles of Bayesian inference. Instead of estimating a single "best-fit" value of the unknown model parameters, our goal is to evaluate the consistency of every possible parameter value with the observed data, or in other words, to construct a distribution that describes the probability that a parameter has a certain value given the observations. This output is referred to as the posterior probability. This framework assumes that the model structure accurately describes the underlying generative process which created the data, but that the underlying parameters are unknown and that there can be some error in the observation of the data.

Bayes' Rule states that the posterior probability of a set of model parameters given observed data can be expressed as a function of the likelihood of observing the data under the model with those parameters () and the prior probability ascribed to those parameters before any data was observed ()

where the denominator is a constant factor – independent of – that only serves to normalize the posterior and thus can be ignored ;

The likelihood function can be defined for a model/data combination based on an understanding of both a) the distribution of model outcomes for a given set of input parameters (if output is stochastic), and b) the nature of the measurement error in observing the data (if relevant) ;

For complex models with many parameters like those used to simulate epidemic spread, it is generally impossible to construct the full posterior distribution either analytically or numerically. Instead, we rely on a class of methods called "Markov Chain Monte Carlo" (MCMC) that allows us to draw a random sample of parameters from the posterior distribution. Ideally, the statistics of the parameters drawn from this sample should be an unbiased estimate of those from the complete posterior.

Markov Chain Monte Carlo methods

In many Bayesian inference problems that arise in scientific model fitting, it is impossible to directly evaluate the full posterior distribution, since there are many parameters to be inferred (high dimensionality) and it is computationally costly to evaluate the model at any individual parameter set. Instead, it is common to employ Markov Chain Monte Carlo (MCMC) methods, which provide a way to iteratively construct a sequence of values that when taken together represent a sample from a desired probability distribution. In the limit of infinitely long sequences ("chains") of values, these methods are mathematically proven to converge to an unbiased sample from the distribution. There are many different MCMC algorithms, but each of them relies on some type of rule for generating a new "sampled" parameter set from an existing one. Our parameter inference method is based on the popular Metropolis-Hastings algorithm. Briefly, at every step of this iterative algorithm, a new set of parameters is jointly proposed, the model is evaluated at that proposed set, the value of the posterior (e.g., likelihood and prior) is evaluated at the proposed set, and if the posterior is improved compared to the previous step, the proposed parameters are "accepted" and become the next entry in the sequence, whereas if the value of the posterior is decreased, the proposed parameters are only accepted with some probability and otherwise rejected (in which case the next entry in the sequences becomes a repeat of the previous parameter set).

The full algorithm for Metropolis-Hastings Markov Chain Monte Carlo is:

• Generate initial set of parameters

• Evaluate the likelihood () and prior () at this parameter set

• For where is the length of the MCMC chain, add to the sequence of parameter values :

Inference algorithm

Likelihood

In our algorithm, model fitting involves comparing timeseries of variables produced by the model (either transmission model state variables or observable outcomes constructed from those variables) to timeseries of observed "ground truth" data with the same time points. For timeseries data that arises from a deterministic, dynamic model, then the overall likelihood can be calculated as the product of the likelihood of the model output at each timepoint (since we assume the data at each timepoint was measured independently). If there are multiple observed datastreams corresponding to multiple model outputs (e.g., cases and deaths) ;

For each subpopulation in the model, the likelihood of observing the "ground truth" data given the model parameters is

where describes the process by which the data is assumed to be observed/measured from the underlying true vales. For example, observations may be assumed to be normally distributed around the truth with a known variance, or, count data may be assumed to be generated by a Poisson process ;

And the overall likelihood taking into account all subpopulations, is the product of the individual likelihoods

Note that the likelihood for each subpopulation depends not only on the parameter values that act within that subpopulation, but on the entire parameter set , since in general the infection dynamics in one subpopulation are also affected by those in each other region. Also note that we assume that the parameters only impact the likelihood through the single model output timeseries . While this is exactly true for a deterministic model, we make the simplifying assumption that it is also true for stochastic models, instead of attempting to calculate the full distribution of possible trajectories for a given parameter set and include that in the likelihood as well.

Fitting algorithm

The method we use for estimating model parameters is based on the Metropolis-Hastings algorithm, which is a class of Markov Chain Monte Carlo (MCMC) methods for obtaining samples from a posterior probability distribution. We developed a custom version of this algorithm to deal with some of the particular mathematical properties and computational challenges of fitting large disease transmission models ;

There are to major unique features of our adapted algorithm:

• Parallelization – Generally MCMC methods starting from a single initial parameter set and generating an extremely long sequence of parameter samples such that the Markov process is acceptably close to a stationary state where it represents an unbiased sample from the posterior. Instead, we simulate multiple shorter chains in parallel, starting from different initial conditions, and pool the results. Due to the computational time required to simulate the epidemic model, and the timescale on which forecasts of epidemic trajectories are often needed (~weeks), it is not possible to sequentially simulate the model millions of times. However, modern supercomputers allow massively parallel computation. The hope of this algorithm is that the parallel chains sample different subspaces of the posterior distribution, and together represent a reasonable sample from the full posterior. To maximize the chance of at least local stationarity of these subsamples, we pool only the final values of each of the parallel chains.

• Multi-level – Our pipeline, and the fitting algorithm in particular, were designed to be able to simulate disease dynamics in a collection of linked subpopulations. This population structure creates challenges for model fitting. We want the model to be able to recreate the dynamics in each subpopulation, not just the overall summed dynamics. Each subpopulation has unique parameters, but due to the coupling between them (), the model outcomes in one subpopulation also depend on the parameter values in other subpopulations. For some subpopulations (), this coupling may effectively be weak and have little impact on dynamics, but for others (), spillover from another closely connected subpopulation may be the primary driver of the local dynamics. Thus, the model cannot be separately fit to each subpopulation, but must consider the combined likelihood. However, such an algorithm may be very slow to find parameters that optimize fits in all locations simultaneously, and may be predominantly drawn to fitting to the largest/most connected subpopulations. The avoid these issues, we

Note that while the traditional Metropolis-Hastings algorithm for MCMC will provably converge to a stationary distribution where the sequence of parameters represents a sample from the posterior distribution, no such claim has been mathematically proven for our method.

• For , where is the number of parallel MCMC chains (also known as slots)

  • Generate initial state

    • Generate an initial set of parameters , and copy this to both the global (

We consider the sequence to represent a sample from the posterior probability distribution, and use it to calculate statistics about the inferred parameter values and the epidemic trajectories resulting from them (e.g., mean, median, 95% credible intervals).

There are a few important notes/limitations about our method currently:

• All parameters to be fit must be location-specific. There is currently no way to fit a parameter that has the identical value across all locations

• The pipeline currently does not allow for fitting of the basic parameters of the compartmental epidemic model. Instead, these must be fixed, and the value of location-specific "interventions" acting to increase/reduce these parameters can be fit. All parameters related to the observational/outcomes model can be fit, as well as "interventions" acting to increase or reduce them ;

• At no point is the parameter fitting optimizing the fit of the summed total population data to total population model predictions. The "overall" likelihood function used to make "global" parameter acceptance decisions is the product of the individual subpopulations likelihoods (which are based on comparing location-specific data to location-specific model output), which is not equivalent to likelihood for the total population. For example, if overestimates of the model in some subpopulations were exactly balanced by underestimates in others, the total population estimate could be very accurate and the total population likelihood high, but the overall likelihood we use here would still be low.

Hierarchical parameters

The baseline likelihood function used in the fitting algorithm described above allows for parameter values to differ arbitrarily between different subpopulations. However, it may be desired to instead impose constraints on the best-fit parameters, such that subpopulations that are similar in some way, or belong to some pre-defined group, have parameters that are close to one another. Formally, this is typically done with group-level or hierarchical models that fit meta-parameters from which individual subpopulation parameters are assumed to draw. Here, we instead impose this group-level structure by adding an additional term to the likelihood that describes the probability that the set of parameters proposed for a group of subpopulations comes from a normal distribution. This term of the likelihood will be larger when the variance of this parameter set is smaller. Formally

where is a group of subpopulations, is one of the parameters in the set , is the probability density function of the normal distribution, and are the mean and standard deviation of all values of the parameter in the group . There is also the option to use a logit-normal distribution instead of a standard normal, which may be more appropriate if the parameter is a proportion bounded in [0,1].

flepiMoP is flexible pipeline for modeling epidemics. It has functionality for simulating epidemics as well as doing inference for simulation parameters and post-processing of simulation/inference outputs. It is written in a combination of python and R and uses anaconda to manage installations which allows flepiMoP to enforce version constraints across both languages.

Prerequisites

flepiMoP requires the following:

• , and

• .

If you do not have git installed you can go to to find the appropriate installation for your system. It's also recommended, but not required, to have a account. If you're totally new to git and GitHub, GitHub has a very nice introduction to the that is worth reading before continuing.

If you do not have conda installed you can go to to find the appropriate installation for your system. We would recommend selecting the Anaconda Distribution installer of conda.

Installing flepiMoP

Navigate to the parent location of where you would like to install flepiMoP, a subdirectory called flepiMoP will be created there. For example, if you navigate to ~/Desktop then flepiMoP will be installed to ~/Desktop/flepiMoP.

This installations script will guide you through a series of prompts to determine how and where to install flepiMoP. Loosely this script:

1. Determines what directory flepiMoP is being installed into,

2. Optionally gets a clone of flepiMoP if it is not present at the install location,

3. Creates a conda environment to house the installation,

4. Installs

For more help on how to use the installation script you can do ./flepimop-install -h to get help information. For first time users accepting the default prompt will be the best choice (as shown below):

Once the prompts are done the installer will output information about the installations that it is doing. After the installation has completed you should see an installation summary similar to:

This summary gives a brief overview of the R/python/package versions installed. If you encounter any issues with your installation please include this information with your issue report.

Activating A flepiMoP Installation

To activate flepiMoP you need to activate the conda environment that it is installed to with:

Or replacing flepimop-env with the appropriate conda environment if you decided on a non-default conda environment. Once you do this you should have the flepimop CLI available to you with:

Defining Environment Variables (Optional)

flepiMoP frequently uses two environment variables to refer to specific directories both as a default for CLI arguments and throughout the documentation:

1. FLEPI_PATH: Refers to the directory where flepiMoP is installed to, and

2. PROJECT_PATH: Refers to the directory where flepiMoP is being ran from.

Furthermore, you'll likely be navigating between these directories frequently in production usage so having these environment variables set can save some typing.

Continuing with the same paths from the installation example flepiMoP was installed to /Users/example/Desktop/flepiMoP. On Linux/MacOS or in linux shells on windows setting an environment variable can be done by:

Where /your/path/to is the directory containing flepiMoP. If you have already navigated to your flepiMoP directory you can just do:

You can check that the variables have been set by either typing env to see all defined environment variables, or typing echo $FLEPI_PATH and echo $PROJECT_PATH to see the values of FLEPI_PATH and PROJECT_PATH.

However, if you're on Windows:

Where /your/path/to is the directory containing flepiMoP. If you have already navigated to your flepiMoP directory you can just do:

You can check that the variables have been set by either typing set to see all defined environment variables, or typing echo $FLEPI_PATH$ and echo $PROJECT_PATH$ to see the values of FLEPI_PATH and PROJECT_PATH.

For more information on the usage of environment variables with flepiMoP please refer to the documentation.

Run flepiMoP

Now that flepiMoP has been successfully installed on your system you will be able to use the tool to model epidemics.

First, navigate to the PROJECT_PATH folder and make sure to delete any old model output files that are there:

Non-Inference Run

Stay in the PROJECT_PATH folder, and run a simulation directly from forward-simulation Python package gempyor. Call flepimop simulate providing the name of the configuration file you want to run. For example here, we use config_sample_2pop.yml.

This will produce a model_output folder, which you can look at using provided post-processing functions and scripts.

We recommend using model_output_notebook.Rmd as a starting point to interact with your model outputs. First, modify the YAML preamble in the notebook (make sure the configuration file listed matches the one used in simulation), then knit this markdown. This will produce plots of the prevalence of infection states over time. You can edit this markdown to produce any figures you'd like to explore your model output.

For your first flepiMoP run, it's better to run each command individually as described above to be sure each exits successfully. However, eventually you can put all these steps together in a script, seen below:

Note that you only have to re-run the installation steps once each time you update any of the files in the flepimop repository (either by pulling changes made by the developers and stored on Github, or by changing them yourself). If you're just running the same or different configuration file, just repeat the final steps:

Inference Run

An inference run requires a configuration file that has the inference section. Stay in the $PROJECT_PATH folder, and run the inference script, providing the name of the configuration file you want to run. For example here, we use config_sample_2pop_inference.yml.

This will run the model and create a lot of output files in $PROJECT_PATH/model_output/.

The last few lines visible on the command prompt should be:

[[1]]

[[1]][[1]]

[[1]][[1]][[1]]

NULL

If you want to quickly do runs with options different from those encoded in the configuration file, you can do that from the command line, for example

where:

• n is the number of parallel inference slots,

• j is the number of CPU cores to use on your machine (if j > n, only n cores will actually be used. If j <n, some cores will run multiple slots in sequence)

Again, it is helpful to run the model output notebook (model_output_notebook.Rmd to explore your model outputs. Knitting this file for an inference run will also provide an analysis of your fits: the acceptance probabilities, likelihoods overtime, and the fits against the provided ground truth.

For your first flepiMoP inference run, it's better to run each command individually as described above to be sure each exits successfully. However, eventually you can put all these steps together in a script, seen below:

Note that you only have to re-run the installation steps once each time you update any of the files in the flepimop repository (either by pulling changes made by the developers and stored on Github, or by changing them yourself). If you're just running the same or different configuration file, just repeat the final steps

Examining Model Output

If your run is successful, you should see your output files in the model_output folder. The structure of the files in this folder is described in the section. By default, all the output files are .parquet format (a compressed format which can be imported as dataframes using R's arrow package arrow::read_parquet or using the free desktop application for quick viewing). However, you can add the option --write-csv to the end of the commands to run the code (e.g., flepimop simulate --write-csv config.yml) to have everything saved as .csv files instead ;

Updating flepiMoP

You can use the flepimop-install script provided by the flepiMoP repository to update your install of flepiMoP with:

Or to reinstall flepiMoP from scratch (say if your conda environment is very out of date or in a bad state) you can do so with:

Next Steps

These configs and notebooks should be a good starting point for getting started with flepiMoP. To explore other running options, see .

These details cover how to install and initialize flepiMoP on an HPC environment and submit a job with slurm.

For getting access to one of the supported HPC environments please refer to the following documentation before continuing:

• for UNC users, or

• for JHU users.

External users will need to consult with their PI contact at the respective institution.

Installing flepiMoP

This task needs to be ran once to do the initial install of flepiMoP.

Download and run the the appropriate installation script with the following command:

Substituting <cluster-name> with either rockfish or longleaf. This script will install flepiMoP to the correct locations on the cluster. Once the installation is done the conda environment can be activated and the script can be removed with:

Updating flepiMoP

Updating flepiMoP is designed to work just the same as installing flepiMoP. First change directory to your flepiMoP installation and then make sure that your clone of the flepiMoP repository is set to the branch you are working with (if doing development or operations work) and then run the flepimop-install-<cluster-name> script, substituting <cluster-name> with either rockfish or longleaf.

Initialize The Created flepiMoP Environment

These steps to initialize the environment need to run on a per run or as needed basis.

Change directory to where a full clone of the flepiMoP repository was placed (it will state the location in the output of the script above). And then run the hpc_init script, substituting <cluster-name> with either rockfish or longleaf. This script will assume the same defaults as the script before for where the flepiMoP clone is and the name of the conda environment. This script will also ask about the path to your flepiMoP installation and project directory. It will also ask if you would like to set a default configuration file, if you plan to use the flepimop batch-calibrate command below we recommend pressing enter to skip setting this environment variable. If this is your first time initializing flepiMoP it might be helpful to use configs out of flepiMoP/examples/tutorials directory as a test.

Upon completing this script it will output a sample set of commands to run to quickly test if the installation/initialization has gone okay.

Submitting A Batch Inference Job To Slurm

The main entry point for submitting batch inference jobs is the flepimop batch-calibrate action. This CLI tool will let you submit a job to slurm once logged into a cluster. For details on the available options please refer to flepimop batch-calibrate --help. As a quick example let's submit an R inference and EMCEE inference job. For the R inference run execute the following once logged into either longleaf or rockfish:

This command will produce a large amount of output, due to -vvv. If you want to try the command without actually submitting the job you can pass the --dry-run option. This command will submit a job to calibrate the sample 2 population configuration which uses R inference. The R inference supports array jobs so each chain will be run on an individual node with 1 CPU and 1GB of memory a piece. Additionally the extra option allows you to provide additional info to the batch system, in this case what partition to submit the jobs to but email is also supported with slurm for notifications. After running this command you should notice the following output:

• config_sample_2pop-YYYYMMDDTHHMMSS.yml: This file contains the compiled config that is actually submitted for inference,

• manifest.json: This file contains a description of the submitted job with the command used, the job name, and flepiMoP and project git commit hashes,

• slurm-*_*.out: These files contain output from slurm for each of the array jobs submitted,

For operational runs these files should be committed to the checked out branch for archival/reproducibility reasons. Since this is just a test you can safely remove these files after inspecting them.

Now, let's submit an EMCEE inference job with the same tool. Importantly, the options we'll use won't change much because flepimop batch-calibrate is designed to provide a unified implementation independent interface.

One notable difference is, unlike R inference, EMCEE inference only supports running on 1 node so resources for this command are adjusted accordingly:

• Swapping 4 nodes with 1 cpu each to 1 node with 4 cpus, and

• Doubling the memory usage from 4 nodes with 1GB each for 4GB total to 1 node with 8GB for 8GB total.

The extra increase in memory is to run a configuration that is slightly more resource intense than the previous example. This command will also produce a similar set of record keeping files like before that you can safely remove after inspecting.

Estimating Required Resources For A Batch Inference Job

When inspecting the output of flepimop batch-calibrate --help you may have noticed several options named --estimate-*. While not required for the smaller jobs above this tool has the ability to estimate the required resources to run a larger batch estimation job. The tool does this by running smaller jobs and then projecting the required resources for a large job from those smaller jobs. To use this feature provide the --estimate flag, a job size of the targeted job, resources for test jobs, and the following estimation settings:

• --estimate-runs: The number of smaller jobs to run to estimate the required resources from,

• --estimate-interval: The size of the prediction interval to use for estimating the resource/time limit upper bounds,

• --estimate-vary: The job size elements to vary when generating smaller jobs,

Effectively using these options requires some knowledge of the underlying inference method. Sticking with the simple usa state level example above try submitting the following command (after cleaning up the output from the previous example):

In short, this command will submit 6 test jobs that will vary simulations and measure time and memory. The number of simulations will be used to project the required resources. The test jobs will range from 1/5 to 1/10 of the target job size. This command will take a bit to run because it needs to wait on these test jobs to finish running before it can do the analysis, so you can check on the progress by checking the output of the simple_usa_statelevel_estimation.log file.

Once this command finishes running you should notice a file called USA_influpaint_resources.json. This JSON file contains the estimated resources required to run the target job. You can submit the target job with the estimated resources by using the same command as before without the --estimate-* options and using the --from-estimate option to pull the information from the outputted file:

Saving Model Outputs On Batch Inference Job Finish

For production runs it is particularly helpful to save the calibration results after a successful run to long term storage for safe keeping. To accomplish this flepimop batch-calibrate can chain a call to flepimop sync after a successful run via the --sync-protocol option. For more details on the flepimop sync command in general please refer to the guide.

For a quick demonstration of how to use this option start with the config_sample_2pop_inference.yml configuration file and add the following section:

Where /path/to/an/example-folder and s3://my-bucket/and-sub-bucket are place holders for paths to your desired location. Importantly, note that there is no trailing slash on the model_output directory name. This will cause flepimop sync to sync the model_output directory itself and not just it's contents. You can also apply additional filters to the sync protocols here, say to limit the backed up model outputs to certain folders or exclude llik outputs, but the --sync-protocol option will add filters to limit the synced directories to those corresponding to the run submitted. Note that users do not need to specify run/job ids or configuration file names in the sync protocol. The flepimop batch-calibrate CLI will take advantage of flepimop sync's options to set paths appropriately to accommodate for run/job ids.

Modifying the first flepimop batch-calibrate command from before:

This command will submit an array job just like before, but will also add a dependent job with the same name prefixed with 'sync_'. This should looks like:

After those jobs finish the results can be found in a subdirectory named after the job and whose contents will look like:

Note that this contains the model_output directory but only limited to the batch run named 'sample_2pop-20250521T190823_Ro_all_test_limits' as well as a file called manifest.json which can be used to reproduce the run from scratch if needed.

Saving Model Outputs To S3 For Hopkins Users

For Hopkins affiliated users there is a configuration file patch included with flepiMoP that can be used to add S3 syncing for model outputs to s3://idd-inference-runs. Taking the example before of running the config_sample_2pop_inference.yml configuration we can slightly modify the command to:

This will take advantage of the patching abilities of the flepimop batch-calibrate to add a sync protocol named s3-idd-inference-runs that will save the results to the s3://idd-inference-runs bucket.

Access model files

See the Before any run section to ensure you have access to the correct files needed to run. On your local machine, determine the file paths to:

• the directory containing the flepimop code (likely the folder you cloned from Github), which we'll call <FLEPI_PATH>

• the directory containing your project code including input configuration file and population structure, which we'll call <PROJECT_PATH>

🧱 Set up Docker

is a software platform that allows you to build, test, and deploy applications quickly. Docker packages software into standardized units called containers that have everything the software needs to run including libraries, system tools, code, and runtime. This means you can run and install software without installing the dependencies in the local operating system.

A Docker container is an environment which is isolated from the rest of the operating system i.e. you can create files, programs, delete and everything but that will not affect your OS. It is a local virtual OS within your OS ;

For flepiMoP, we have a Docker container that will help you get running quickly ;

Make sure you have the Docker software installed, and then open your command prompt or terminal application ;

Run the Docker image

First, make sure you have the latest version of the flepimop Docker (hopkinsidd/flepimop) downloaded on your machine by opening your terminal application and entering:

Next, run the Docker image by entering the following, replace <FLEPI_PATH> and <PROJECT_PATH> with the path names for your machine (no quotes or brackets, just the path text):

In this command, we run the Docker container, creating a volume and mounting (-v) your code and project directories into the container. Creating a volume and mounting it to a container basically allocates space in Docker for it to mirror - and have read and write access - to files on your local machine ;

The folder with the flepiMoP code <PROJECT_PATH> will be on the path flepimop within the Docker environment, while the project folder will be at the path `drp. ;

{% hint style="success" %} You now have a local Docker container installed, which includes the R and Python versions required to run flepiMop with all the required packagers already installed ; {% endhint %}

{% hint style="info" %} You don't need to re-run the above steps every time you want to run the model. When you're done using Docker for the day, you can simply "detach" from the container and pause it, without deleting it from your machine. Then you can re-attach to it when you next want to run the model ; {% endhint %}

Define environment variables

Create environmental variables for the paths to the flepimop code folder and the project folder:

`

``bash export FLEPI_PATH=/home/app/flepimop/ export PROJECT_PATH=/home/app/drp/

Each installation step may take a few minutes to run.

Run the code

Everything is now ready 🎉 The next step depends on what sort of simulation you want to run: One that includes inference (fitting model to data) or only a forward simulation (non-inference). Inference is run from R, while forward-only simulations are run directly from the Python package gempyor.

In either case, navigate to the project folder and make sure to delete any old model output files that are there

Inference run

An inference run requires a configuration file that has the inference section. Stay in the $PROJECT_PATH folder, and run the inference script, providing the name of the configuration file you want to run (ex. config.yml ;

This will run the model and create a lot of output files in $PROJECT_PATH/model_output/ ;

The last few lines visible on the command prompt should be:

[[1]]

[[1]][[1]]

[[1]][[1]][[1]]

NULL

If you want to quickly do runs with options different from those encoded in the configuration file, you can do that from the command line, for example

where:

• n is the number of parallel inference slots,

• j is the number of CPU cores to use on your machine (if j > n, only n cores will actually be used. If j <n, some cores will run multiple slots in sequence)

You can put all of this together into a single script that can be run all at once ;

Non-inference run

Stay in the $PROJECT_PATH folder, and run a simulation directly from forward-simulation Python package gempyor,call flepimop simulate providing the name of the configuration file you want to run (ex. config.yml):

You can put all of this together into a single script that can be run all at once ;

Finishing up

You can avoid repeating all the above steps every time you want to run the code. When the docker run command creates an container, it is stored locally on your computer with all the installed packages/variables/etc you created. You can leave this container and come back to it whenever you want, without having to redo all this set up ;

When you're in the Docker container, figure out the name Docker has given to the container you created by typing

the output will be something silly like

write this down for later reference. You can also see the container name in the Docker Desktop app's Containers tab ;

To "detach" from the Docker container and stop it, type CTLR + c

The command prompt for your terminal application is now just running locally, not in the Docker container ;

Next time you want to re-start and "attach" the container, type

at the command line or hit the play button ▶️ beside the container's name in the Docker app. Replace container_name with the name for your old container ;

Then "attach" to the container by typing

The reason that stopping/starting a container is separate from detaching/attaching is that technically you can leave a container (and any processes within it) running in the background and exit it. In case you want to do that, detach and leave it running by typing CTRL + p then quickly CTRL + q. Then when you want to attach to it again, you don't need to do the part about starting the container ;

If you the core model code within the flepimop repository (flepimop/flepimop/gempyor_pkg/ or flepimop/flepimop/R_packages) has been edited since you created the contained, or if the R or Python package requirements have changed, then you'll have to re-run the steps to install the packages, but otherwise, you can just start running model code!

🖥 Start and access AWS submission box

Spin up an Ubuntu submission box if not already running. To do this, log onto AWS Console and start the EC2 instance.

Update IP address in .ssh/config file. To do this, open a terminal and type the command below. This will open your config file where you can change the IP to the IP4 assigned to the AWS EC2 instance (see AWS Console for this):

SSH into the box. In the terminal, SSH into your box. Typically we name these instances "staging", so usually the command is:

🧱 Setup

Now you should be logged onto the AWS submission box. If you haven't yet, set up your directory structure.

🗂 Create the directory structure (ONCE PER USER)

Type the following commands:

Have your Github ssh key passphrase handy so you can paste it when prompted (possibly multiple times) with the git pull command. Alternatively, you can add your github key to your batch box so you don't have to enter your token 6 times per day.

🚀 Run inference using AWS (do everytime)

🛳 Initiate the Docker

Start up and log into the docker container, and run setup scripts to setup the environment. This setup code links the docker directories to the existing directories on your box. As this is the case, you should not run job submission simultaneously using this setup, as one job submission might modify the data for another job submission.

Setup environment

To set up the environment for your run, run the following commands. These are specific to your run, i.e., change VALIDATION_DATE, FLEPI_RUN_INDEX and RESUME_LOCATION as required. If submitting multiple jobs, it is recommended to split jobs between 2 queues: Compartment-JQ-1588569569 and Compartment-JQ-1588569574.

NOTE: If you are not running a resume run, DO NOT export the environmental variable RESUME_LOCATION.

Additionally, if you want to profile how the model is using your memory resources during the run, run the following commands

Then prepare the pipeline directory (if you have already done that and the pipeline hasn't been updated (git pull says it's up to date). You need to set $PROJECT_PATH to your data folder. For a COVID-19 run, do:

for Flu do:

Now for any type of run:

For now, just in case: update the arrow package from 8.0.0 in the docker to 11.0.3 ;

Now flepiMoP is ready 🎉 ;

Do some clean-up before your run. The fast way is to restore the $PROJECT_PATH git repository to its blank states (⚠️ removes everything that does not come from git):

Then run the preparatory data building scripts and you are good

Now you may want to test that it works :

If this fails, you may want to investigate this error. In case this succeeds, then you can proceed by first deleting the model_output:

Launch your inference batch job on AWS

Assuming that the initial test simulation finishes successfully, you will now enter credentials and submit your job onto AWS batch. Enter the following command into the terminal:

You will be prompted to enter the following items. These can be found in a file you received from Shaun called new_user_credentials.csv.

• Access key ID when prompted

• Secret access key when prompted

• Default region name: us-west-2

• Default output: Leave blank when this is prompted and press enter (The Access Key ID and Secret Access Key will be given to you once in a file)

Now you're fully set to go 🎉

To launch the whole inference batch job, type the following command:

This command infers everything from you environment variables, if there is a resume or not, what is the run_id, etc., and the default is to carry seeding if it is a resume (see below for alternative options).

If you'd like to have more control, you can specify the arguments manually:

We allow for a number of different jobs, with different setups, e.g., you may not want to carry seeding. Some examples of appropriate setups are given below. No modification of these code chunks should be required ;

Document the submission

After the job is successfully submitted, you will now be in a new branch of the data repo. Commit the ground truth data files to the branch on github and then return to the main branch:

Send the submission information to slack so we can identify the job later. Example output: