Skip to content

SEIRSNetworkModel Class

Ryan Seamus McGee edited this page Aug 10, 2020 · 4 revisions

Contents:

Initializing the model

All model parameter values and properties, including the contact network are set in the call to the SEIRSNetworkModel constructor. The interaction network G and the basic rate parameters beta, sigma, and gamma are the only required arguments. All other arguments represent parameters for optional model features (e.g., handling quarantined individuals differently); these optional parameters take default values that bypass their corresponding features when not provided in the constructor.

All constructor parameters are listed and described below, followed by instructions for specifying contact networks and example instantiations of various elaborations of the model.

Specifying parameters

Individualized Parameter Values: Parameter values are assigned to members of the population on an individual basis. Parameter values can be specified by providing a single numerical value, which assigns that common value to all individuals, or by providing a list of values that specifies the N values to assign to each individual. The population may be either homogeneous or heterogeneous for a given parameter at the user's discretion. Some parameters may also be specified in a pairwise manner by providing a matrix of values.

Some parameters in the table are linked to additional details that appear below.

Required parameters

Constructor Argument Parameter Description Data Type Default Value
G graph specifying the contact network networkx Graph or numpy 2d array REQUIRED
beta transmissibility float or list-like REQUIRED
sigma rate of progression to infectiousness (inverse of latent period) float or list-like REQUIRED
gamma rate of recovery (inverse of infectious period) float or list-like REQUIRED

Additional parameters

Constructor Argument Parameter Description Data Type Default Value
p probability of global interactions (network locality) float or list-like 0.0
alpha susceptibility to infection float 1.0
beta_local transmissibility in local interactions (between close contacts)
See Transmissibility Parameters for more info
float or list-like None
(set equal to beta)
beta_pairwise_mode specifies how pairwise transmissibilities for local interactions are to be calculated from individual transmissibilities
See Transmissibility Parameters for more info
string "infected"
delta individuals' connectivity correction factor (for use in calculating propensity for exosure from local interactions)
See Connectivity Correction Factors for more info
float or list-like None
(use a log-degree term)
delta_pairwise_mode specifies how pairwise connectivity correction factors are to be calculated from individual connectivity correction factors
See Connectivity Correction Factors for more info
float or list-like None
(no correction factor used)
f probability of death for infected individuals (case fatality rate) float or list-like 0.0
mu_I rate of death for hospitalized individuals (inverse of admission-to-death period) float or list-like 0.0
xi rate of re-susceptibility (inverse of temporary immunity period; 0 if permanent immunity) float 0.0
mu_0 rate of baseline mortality float or list-like 0.0
nu rate of baseline birth float or list-like 0.0
G_Q graph specifying the quarantine contact network float or list-like None
q intensity of global interactions for quarantined individuals float or list-like 0.0
alpha_Q susceptibility to infection for quarantined individuals float or list-like None
(set equal to alpha)
beta_Q transmissibility of quarantined individuals float or list-like None
(set equal to beta)
beta_Q_local transmissibility of quarantined individuals in local interactions (between quarantine contacts)
See Transmissibility Parameters for more info
float or list-like None
(use beta_Q values)
delta_Q quarantined individuals' connectivity correction factor (for use in calculating propensity for exosure from local interactions)
See Connectivity Correction Factors for more info
float or list-like None
(use a log-degree term)
sigma_Q rate of progression to infectiousness for quarantined individuals (inverse of latent period) float or list-like None
(set equal to sigma)
gamma_Q rate of recovery for quarantined individuals (inverse of infectious period) float or list-like None
(set equal to gamma)
theta_E rate of testing for exposed individuals float or list-like 0.0
theta_I rate of testing for infectious individuals float or list-like 0.0
phi_E rate of contact tracing testing for exposed individuals float or list-like 0.0
phi_I rate of contact tracing testing for infectious individuals float or list-like 0.0
psi_E probability of positive test (sensitivity) for exposed individuals float or list-like 1.0
psi_I probability of positive test (sensitivity) for infectious individuals float or list-like 1.0
isolation_time amount of time to elapse before isolated individuals transition out of quarantine states float 14
initE initial number of exposed individuals int 0
initI initial number of infectious individuals int 0
initR initial number of recovered individuals int 0
initF initial number of dead individuals int 0
initQ_E initial number of quarantined exposed individuals int 0
initQ_I initial number of quarantined infectious individuals int 0
transition_mode method of calculating propensities of state transitions
See Transition Mode for more info
string "exponential_rates"
node_groups dictionary of lists of node indices specifying groups of nodes to record data for on a group-basis dict None
store_Xseries flag specifying if a timeseries of all node states should be stored bool False
seed ` seed for numpy random number generation int None

Transmissibility parameters

This model considers two modes of disease transmission: a well-mixed mode of global transmission and a contact-network-based mode of local transmission. The propensity for a given individual to become exposed due to global transmission depends on the mean transmissibility of all infectious individuals in the population; the propensity for a given individual to become exposed due to local transmission depends on the pairwise transmissibilities between the focal node and its infectious contacts in the network (see Transmission and Model Equations for more information about these calculations).

Individual Transmissibility Values

All individuals are assigned an Individual Transmissibility Value, which are stored in the beta attribute of the model object (e.g., model.beta[i] gives the Individual Transmissibility Value βi for the ith individual).

The means of the Individual Transmissibility Values for infectious subpopulations are used to calculate the global transmission terms. Individual Transmissibility Values may also be used to generate the Pairwise Transmissibility Values used to calculate the local transmission terms when no specific local/pairwise transmissibility values are provided by the user.

Individual Transmissibility Values can be specified to the model in two ways:

  • Providing a single numerical value to the beta constructor argument assigns the given value to all individuals in the population (homogeneous transmissibility).
    model = SEIRSNetworkModel(..., beta=0.2, ...)
  • Providing a list or array with N values assigns the values in the list to the corresponding individuals (allowing for heterogeneous transmissibility).
    model = SEIRSNetworkModel(..., beta=[0.2, 0.21, 0.2, ..., 0.19], ...)

Different Individual Transmissibility Values can be specified for individuals in quarantine states in the same fashion using the beta_Q constructor argument. When beta_Q values are not specified, the model defaults to using the values specified for beta for individuals in these states.

Pairwise Transmissibility Values

Local transmission from an infected individual $j$ to a susceptible individual i depends on a pairwise transmissibility βji that is assigned to this particular pairwise interaction. Internally, the model class generates a matrix of that gives the Pairwise Transmissibility Values for all interactions defined by the contact network.

This Pairwise Transmissibility Values matrix can be generated/specified in multiple ways using the beta_local constructor argument:

  • By default, when nothing is passed to the beta_local constructor argument, the matrix is generated automatically using the Individual Transmissibility Values given by beta (or beta_Q as applicable) using one of the rules specified by the beta_pairwise_mode constructor argument:
       model = SEIRSNetworkModel(..., beta_local=None, ...)
    • If beta_pairwise_mode = "infected" (default), each Pairwise Transmissibility Value is defined to be the Individual Transmissibility Value of the infected individual (i.e., βji = βj)
    • If beta_pairwise_mode = "infectee", each Pairwise Transmissibility Value is defined to be the Individual Transmissibility Value of the susceptible individual (i.e., βji = βi)
    • If beta_pairwise_mode = "min", each Pairwise Transmissibility Value is defined to be the minimum of the Individual Transmissibility Values of the interacting individuals (i.e., βji = min(βj, βi))
    • If beta_pairwise_mode = "max", each Pairwise Transmissibility Value is defined to be the maximum of the Individual Transmissibility Values of the interacting individuals (i.e., βji = max(βj, βi))
    • If beta_pairwise_mode = "mean", each Pairwise Transmissibility Value is defined to be the mean of the Individual Transmissibility Values of the interacting individuals (i.e., βji = mean(βj, βi))
  • Individual transmissibility values for use in local interactions can be provided to the beta_local constructor argument. In this case the matrix is generated automatically according to the beta_pairwise_mode rule as outlined above, but using the provided beta_local values instead of the beta values.
       model = SEIRSNetworkModel(..., beta_local=[0.24, 0.25, 0.23, ..., 0.25], ...)
  • A fully specified matrix of Pairwise Transmissibility Values can be provided. In this case the provided matrix is used directly.
       model = SEIRSNetworkModel(..., beta_local=numpy.array([[0, 0, 0.22, ..., 0],[0.24, 0, 0, ..., 0.25],...] ...)

Different Pairwise Transmissibility Values can be specified for individuals in quarantine states in the same fashion using the beta_Q_local constructor argument. When beta_Q_local values are not specified, the model defaults to using the matrix associated with beta_local for individuals in these states.

Connectivity correction factors

Another weighting factor 𝛿ji appears in the calculation of propensity for exposure due to local transmission. This pairwise factor can be used to weight the transmissibility of interactions according to the connectivity of the interacting individuals (or other arbitrary criteria).

Internally, the model class generates a matrix that gives the pairwise 𝛿ji values for all interactions defined by the contact network. This matrix can be generated/specified in multiple ways using the delta constructor argument. The SEIRSNetworkModel class implements several convenience options for generating delta matrices that weight the propensity of according to degree, described below (note: Dj and Di are the degrees of nodes j and i, respectively, and is the mean degree of the network):

  • By default, when nothing is passed to the delta constructor argument, the matrix is generated automatically using the using one of the rules specified by the delta_pairwise_mode constructor argument:
       model = SEIRSNetworkModel(..., delta=None, ...)
    • If delta_pairwise_mode = None (default), this factor is not used (as if 𝛿ji = 1 for all i,j)
    • If delta_pairwise_mode = "infected", this factor is defined as
    • If delta_pairwise_mode = "infectee", this factor is defined as
    • If delta_pairwise_mode = "min", this factor is defined as
    • If delta_pairwise_mode = "max", this factor is defined as
    • If delta_pairwise_mode = "mean" (mentioned here), this factor is defined as
  • Individual delta values can be provided to the delta constructor argument. In this case the matrix is generated automatically according to the delta_pairwise_mode rule as outlined above, but using the provided delta values (and their min/max/mean) instead of the log-degree terms shown in the equations above.
       model = SEIRSNetworkModel(..., delta=[1.0, 1.1, ..., 0.9], ...)
  • A fully specified matrix of pairwise *𝛿ji values can be provided. In this case the provided matrix is used directly.
       model = SEIRSNetworkModel(..., delta=numpy.array([[0, 0, 1.1, ..., 0],[0.9, 0, 0, ..., 1.04],...] ...)

Different 𝛿ji values can be specified for use with individuals in quarantine states using the delta_Q constructor argument. When delta_Q is not specified, the model defaults to using the matrix associated with delta for individuals in quarantine.

Transition mode

The Gillepsie algorithm that runs the stochastic dynamics of this model can be run in one of two modes, as specified by the choice of the following values for the transition_mode constructor argument:

  • "exponential_rates" (default): Propensities are calculated as described in Model Equations for all state transitions. This results in state transition times that are exponentially distributed with expected transition times given by the inverse of the specified rate parameter values (e.g., sigma, lamda, gamma, etc.). While stochastic, the Gillepsie algorithm method is rigorously analogous to the deterministic dynamics of standard mean-field compartment model differential equations, and the Gillepsie algorithm run using this mode will tend to converge on the same trajectories as deterministic models for large networks.
    Caveat: Since transition times are exponentially distributed (and given the shape of this distribution), there is some tendency for some individuals to transition between states significantly sooner than the expected transition time designated by the rate parameter. This is a known caveat for standard compartment models in epidemiology, which is worth considering in some but not all applications.
  • "time_in_state": This simulation mode is motivated by the caveat described above for the more standard "exponential_rates" mode. In this mode, the propensities for Susceptible → Exposed transitions are calculated according to the Pr(S→E) equation in the Model Equations. However, the propensities of all other transitions (i.e., disease progression transitions) are calculated as follows:
    • Every individual has a time-in-state counter that tracks how long they've been in their current state.
    • The propensity of a disease progression transition for a given individual is some very large number when that individual's time-in-state is greater than their expected transition time as designated by the corresponding rate parameter value; this propensity is 0 when their time-in-state is less than this expected time. This results in individuals transitioning between states very nearly exactly at the expected time they were assigned for that transition (note that individuals can be assigned heterogeneous expected transition times, such as from a gamma distribution).
      • Therefore, when an individual is due for a disease progression state transition they are essentially guaranteed to be the next transition in the population with a very small dt. When no individual is due for such a transition, the next transition will be a S→E transmission event with a dt calculated as usual. This amounts to trasmission events moving the simulation time forward as usual and other transitions occurring immediately after transmission events when they are due according to individuals' time-in-states. When the time between transmission events is small (as is the case in moderate and large networks), this is a good approximation of all desired event times (when the population is very small or there are very few cases, the time between transmission events can be relatively large which biases the time of all other transitions upwards as well).

Specifying contact networks

This class implements a model of SEIRS dynamics for populations with a structured contact network (as opposed to standard mean-field compartment models, which assume uniform mixing of the population). When using this network model, a graph specifying the contact network must be specified, where each node represents an individual in the population and edges connect individuals who have regular interactions.

This model also supports scenarios where individuals enter quarantine states in which their parameters and interactions may be different from baseline. Given that this model models transmission on contact networks, a separate graph defining the interactions for individuals in quarantine can be specified (i.e., the quarantine contact network).

  • The baseline contact network is specified by providing a networkx Graph object or a numpy 2d array representing the adjacency matrix to the G constructor argument.
  • The quarantine contact network is specified by providing a networkx Graph object or a numpy 2d array representing the adjacency matrix to the G_Q constructor argument.

Contact networks can be defined and generated by any method that is appropriate for representing the user's population and scenario of interest. Network generation is not a focus of the SEIRS+ package itself, but a few network tools have been included (see Networks for more information).

Epidemic scenarios of interest often involve interaction networks that change in time. Multiple interaction networks can be defined and used at different times in the model simulation using the checkpoints feature described below.

Note: Simulation time increases with network size. Small networks simulate quickly, but have more stochastic volatility. Networks with ~10,000 are large enough to produce per-capita population dynamics that are generally consistent with those of larger networks, but small enough to simulate quickly. We recommend using networks with ~10,000 nodes for prototyping parameters and scenarios, which can then be run on larger networks if more precision is required.

Example model instantiations

Here are a few examples of instantiating and initializing the SEIRS network model. This is far from an exhaustive list of use cases, and not all parameters and features are shown being used (the parameter values that are shown are just illustrative placeholders).

For an example of this model being used to simulate a sophisticated TTI scenario, see the Community TTI and Workplace TTI examples.

Minimal SEIR dynamics on a network
model = SEIRSNetworkModel(G=network, beta=0.45, sigma=1/5.1, gamma=1/7, initE=100)
Minimal SEIRS dynamics on a network
model = SEIRSNetworkModel(G=network, beta=0.45, sigma=1/5.1, gamma=1/7, xi=0.001, initE=100)
SEIRS dynamics on a network with rate-based testing and isolation (theta and psi testing params > 0, quarantine network Q provided)
model = SEIRSNetworkModel(G=network, p=0.2, beta=0.45, sigma=1/5.1 gamma=1/7,
                           Q=quarantine_network, q=0.05,
                           theta_E=0.02, theta_I=0.02, psi_E=0.1, psi_I=0.85, initE=100)
SEIRS dynamics on a network with rate-based testing, rate-based tracing , and isolation (theta and psi testing params > 0, phi tracing params > 0, quarantine network Q provided)
model = SEIRSNetworkModel(G=network, p=0.2, beta=0.45, sigma=1/5.1 gamma=1/7,
                           Q=quarantine_network, q=0.05,, 
                           theta_E=0.02, theta_I=0.02, psi_E=0.1, psi_I=0.85, phi_E=0.2, phi_I=0.2, 
                           initE=100)
SEIRS dynamics on a network with rate-based TTI, where transmissibility and susceptibility are reduced for quarantined individuals (beta_Q and alpha_Q specified)
model = SEIRSNetworkModel(G=network, p=0.2, beta=0.45, sigma=1/5.1 gamma=1/7,
                           Q=quarantine_network, q=0.05, beta_Q=0.25, alpha_Q=0.5, 
                           theta_E=0.02, theta_I=0.02, psi_E=0.1, psi_I=0.85, phi_E=0.2, phi_I=0.2, 
                           initE=100)

Running the model

Running a simulation

run() function

Once a model is initialized, a simulation can be run with a call to the following function:

model.run(T=300)

The run() function has the following arguments

Argument Description Data Type Default Value
T total runtime of simulation numeric REQUIRED
checkpoints dictionary of checkpoint lists (see section below) dictionary None
print_interval (network model only) time interval to print sim status to console numeric 10
verbose if True, print count in each state at print intervals, else just the time bool False

run_iteration() function

The run_iteration() function advances the model by a single update (this is called repeatedly within run()). The run_iteration() function can be used to develop custom simulation loops that execute additional simulation logic in between updates.

T = 300
running = True
while(running and model.t < T):
   # ... simulation logic ...
   running = model.run_iteration()
   # ... simulation logic ...

For an example of a sophisticated custom simulation loop that uses this feature, see TTI Simulation Loop

The run_iteration() function has no arguments.

Accessing & changing parameters during a simulation

All parameters are stored as attributes in the model object, and these attributes can be accessed and modified directly by your script after or during simulations.

For example:

# Get the list of susceptibility values for all individuals:
susceptibility_vals = model.alpha
# ...
# Change the transmissibility of a given individual (node ID 34)
model.beta[34] = 0.6

Checkpoints

Model parameters can be easily changed during a simulation run using "checkpoints". A dictionary holds a list of checkpoint times (checkpoints['t']) and lists of new values to assign to various model parameters at each checkpoint time.

Example of running a simulation with checkpoints:

checkpoints = {'t':       [20, 100], 
               'G':       [G_distancing, G_normal], 
               'p':       [0.1, 0.5], 
               'theta_E': [0.02, 0.02], 
               'theta_I': [0.02, 0.02], 
               'phi_E':   [0.2, 0.2], 
               'phi_I':   [0.2, 0.2]}

model.run(T=300, checkpoints=checkpoints)

Any model parameter listed in the model constructor can be updated in this way. Only model parameters that are included in the checkpoints dictionary have their values updated at the checkpoint times, all other parameters keep their pre-existing values.

Use cases of this feature include:

  • Changing parameters during a simulation, such as changing transition rates or testing parameters every day, week, on a specific sequence of dates, etc.
  • Starting and stopping interventions, such as social distancing (changing interaction network), testing and contact tracing (setting relevant parameters to non-zero or zero values), etc.

Manually updating parameters

One approach is to run the same model object multiple times. Each time the run() function of a given model object is called, it starts a simulation from the state it left off in any previous simulations. For example:*

model.run(T=100)    # simulate the model for 100 time units
# ... 
# do other things, such as processing simulation data or changing parameters 
model.beta[17] = 0.3    # example of changing the transmissibility of a particular individual (node 17)
# ...
model.run(T=200)    # simulate the model for an additional 200 time units, picking up where the first sim left off

Another option is to write a custom simulation loop using the run_iteration() function that defines logic for manually modifying model parameters at various updates.

Accessing & changing node states during a simulation

The run() and run_iteration() functions direct state updates according to the model dynamics. This section describes options for manually changing states outside of these normal dynamics.

The compartment states of individuals are stored in the model.X attribute of the model object, which is an array of each node's state. States are encoded as integers within the model object (these numerical state IDs are defined as object attributes, so you don't need to store or memorize them):

model.S       = 1
model.E       = 2
model.I       = 3
model.R       = 4
model.F       = 5
model.Q_E     = 6
model.Q_I     = 7

Therefore the state of a given node i can be manually accessed or changed by referencing the ith index of the X array, for example:

# Check the state of a given individual (node ID 23)(e.g., for custom testing logic):
if(model.X[23] == model.E or model.X[23] == model.I):
   # do something 

# Change the state of a given individual (node ID 117) to exposed:
model.X[23] = model.E

set_isolation() function

The set_isolate() convenience function makes it easy to manually move an individual into or out of isolation. This function implements the common logic for checking the current disease state of a node, moving it into or out of the corresponding quarantine state as directed, and resetting the node's isolation time counter.

The set_isolate() function has the following arguments:

Argument Description Data Type Default Value
node index of the node being whose quarantine status is being updated int REQUIRED
isolation True moves the given node into the quarantine state corresponding to its current disease state, False moves it out of the quarantine state bool REQUIRED

Accessing simulation results

Model parameter values and the state variable time series generated by the simulation are stored in the attributes of the model object and can be accessed directly as follows:

S = model.numS            # time series of S counts
E = model.numE            # time series of E counts
I = model.numI            # time series of I counts
R = model.numR            # time series of R counts
F = model.numF            # time series of F counts
Q_E = model.numQ_E        # time series of Q_E counts
Q_I = model.numQ_I        # time series of Q_I counts


t = model.tseries         # time values corresponding to the above time series

G_normal     = model.G    # interaction network graph
G_quarantine = model.Q    # quarantine interaction network graph

beta_vals = model.beta    # value of beta parameters
# similar for other parameters

Note: convenience methods for plotting these time series are included in the package. See below.

Node States Time Series

If the model object is initialized with the constructor argument store_Xseries = True, then the time series of all node states at all time points will be stored in the model.Xseries attribute of the model object. This data can be sizable for large networks and/or sims with many time points, so this data is not saved by default (the time series of aggregate number of nodes in each state are saved in any case, and can be accessed as shown in the code block above).

Node Groups

It may be of interest to track the trajectories for a particular subset of nodes (such as a subpopulation that you identify as being high risk or having some other property). A model can be set up to store time series of the counts of nodes in each state for user defined sets of nodes using the node_groups constructor argument. This argument takes a dict with key-value pairs corresponding to the name and node indices for the group(s), for example:

node_groups = { 'high_risk': [34, 117, 42, ...], 'low_risk': [23, 38, 10, ...] } 

Then time series data for these groups are stored in the nodeGroupData attribute of the model object, which is a dict that is accessed as follows:

high_risk_Eseries = model.nodeGroupData['high_risk']['numE']

This node groups feature is a convenience for recording time series for subpopulations of interest. When the user does not provide a dictionary of node group definitions nothing is stored in model.nodeGroupData. Of course, time series of the state counts are always stored for the full population.

Visualization

Visualizing the results

The SEIRSNetworkModel class has a plot() convenience function for plotting simulation results on a matplotlib axis. This function generates a line plot of the frequency of each model state in the population by default, but there are many optional arguments that can be used to customize the plot.

There are also convenience functions for generating a full figure out of model simulation results (optionally, arguments can be provided to customize the plots generated by these functions, see below).

  • figure_basic() calls the plot() function with default parameters to generate a line plot of the frequency of each state in the population.
  • figure_infections() calls the plot() function with default parameters to generate a stacked area plot of the frequency of only the infection states (E, I, DE, DI) in the population.

Parameters that can be passed to any of the above functions include:

Argument Description
plot_S 'line', shaded, 'stacked', or False
plot_E 'line', shaded, 'stacked', or False
plot_I 'line', shaded, 'stacked', or False
plot_R 'line', shaded, 'stacked', or False
plot_F 'line', shaded, 'stacked', or False
plot_D_E 'line', shaded, 'stacked', or False
plot_D_I 'line', shaded, 'stacked', or False
combine_D True or False
color_S matplotlib color of line or stacked area
color_E matplotlib color of line or stacked area
color_I matplotlib color of line or stacked area
color_R matplotlib color of line or stacked area
color_F matplotlib color of line or stacked area
color_D_E matplotlib color of line or stacked area
color_D_I matplotlib color of line or stacked area
color_reference matplotlib color of line or stacked area
dashed_reference_results seirsplus model object containing results to be plotted as a dashed-line reference curve
dashed_reference_label string for labeling the reference curve in the legend
shaded_reference_results seirsplus model object containing results to be plotted as a dashed-line reference curve
shaded_reference_label string for labeling the reference curve in the legend
vlines list of x positions at which to plot vertical lines
vline_colors list of matplotlib colors corresponding to the vertical lines
vline_styles list of matplotlib linestyle strings corresponding to the vertical lines
vline_labels list of string labels corresponding to the vertical lines
ylim max y-axis value
xlim max x-axis value
legend display legend, True or False
title string plot title
side_title string plot title along y-axis
plot_percentages if True plot percentage of population in each state, else plot absolute counts
figsize tuple specifying figure x and y dimensions
use_seaborn if True import seaborn and use seaborn styles