Stochastic simulations

Stochastic simulations with compartment models and age structure. Has Gillespie and tau-leaping implemented.

A list of methods for stochastic simulations of age-structured compartment models along with link to notebook examples is given below.

Model

class pyross.stochastic.Model

Generic user-defined epidemic model.

…

Parameters:
  • model_spec (dict) – A dictionary specifying the model. See Examples.
  • parameters (dict) – A dictionary containing the model parameters. All parameters can be float if not age-dependent, and np.array(M,) if age-dependent. An optional element with the key ‘seed’ may be supplied as a seed for the pseudo-random number generator.
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(3*M, )) – Initial number in each compartment and class
  • time_dep_param_mapping (python function, optional) – A user-defined function that takes a dictionary of time-independent parameters and time as an argument, and returns a dictionary of the parameters of model_spec. Default: Identical mapping of the dictionary at all times.

Examples

An example of model_spec and parameters for SIR class with a constant influx

>>> model_spec = {
        "classes" : ["S", "I"],
        "S" : {
            "constant"  : [ ["k"] ],
            "infection" : [ ["I", "S", "-beta"] ]
        },
        "I" : {
            "linear"    : [ ["I", "-gamma"] ],
            "infection" : [ ["I", "S", "beta"] ]
        }
    }
>>> parameters = {
        'beta': 0.1,
        'gamma': 0.1,
        'seed': 1234,
        'k': 1,
    }
model_class_data()
Parameters:data (dict) – The object returned by simulate.
Returns:
Return type:The population of class model_class_key as a time series
simulate()

Performs the Stochastic Simulation Algorithm (SSA)

Parameters:
  • x0 (np.array) – Initial condition.
  • contactMatrix (python function(t)) – The social contact matrix C_{ij} denotes the average number of contacts made per day by an individual in class i with an individual in class j
  • Tf (float) – Final time of integrator
  • Nf (Int) – Number of time points to evaluate.
  • method (str, optional) – SSA to use, either ‘gillespie’ or ‘tau_leaping’. The default is ‘gillespie’.
  • nc (TYPE, optional) –
  • epsilon (float, optional) –

    The acceptable relative change of the rates during each tau-leaping step, as defined in Cao et al:

    The default is 0.03

  • tau_update_frequency (TYPE, optional) –
Returns:

X: output path from integrator, t : time points evaluated at, ‘event_occured’ , ‘param’: input param to integrator.
Return type:

dict

Spp

class pyross.stochastic.Spp

This is a slightly more specific version of the class Model.

Spp is still supported for backward compatibility.

Model class is recommended over Spp for new users.

The Spp class works like Model but infection terms use a single class S

…

Parameters:
  • model_spec (dict) – A dictionary specifying the model. See Examples.
  • parameters (dict) – A dictionary containing the model parameters. All parameters can be float if not age-dependent, and np.array(M,) if age-dependent. An optional element with the key ‘seed’ may be supplied as a seed for the pseudo-random number generator.
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(3*M, )) – Initial number in each compartment and class
  • time_dep_param_mapping (python function, optional) – A user-defined function that takes a dictionary of time-independent parameters and time as an argument, and returns a dictionary of the parameters of model_spec. Default: Identical mapping of the dictionary at all times.

Examples

An example of model_spec and parameters for SIR class with a constant influx

>>> model_spec = {
        "classes" : ["S", "I"],
        "S" : {
            "constant"  : [ ["k"] ],
            "infection" : [ ["I", "-beta"] ]
        },
        "I" : {
            "linear"    : [ ["I", "-gamma"] ],
            "infection" : [ ["I", "beta"] ]
        }
    }
>>> parameters = {
        'beta': 0.1,
        'gamma': 0.1,
        'seed': 1234,
        'k': 1,
    }

SppQ

class pyross.stochastic.SppQ

User-defined epidemic model with quarantine stage.

This is a slightly more specific version of the class Model.

SppQ is still supported for backward compatibility.

Model class is recommended over SppQ for new users.

To initialise the SppQ model,

…

Parameters:
  • model_spec (dict) – A dictionary specifying the model. See Examples.
  • parameters (dict) – A dictionary containing the model parameters. All parameters can be float if not age-dependent, and np.array(M,) if age-dependent. An optional element with the key ‘seed’ may be supplied as a seed for the pseudo-random number generator.
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(3*M, )) – Initial number in each compartment and class
  • time_dep_param_mapping (python function, optional) – A user-defined function that takes a dictionary of time-independent parameters and time as an argument, and returns a dictionary of the parameters of model_spec. Default: Identical mapping of the dictionary at all times.

Examples

An example of model_spec and parameters for SIR class with random testing (without false positives/negatives) and quarantine

>>> model_spec = {
        "classes" : ["S", "I"],
        "S" : {
            "infection" : [ ["I", "-beta"] ]
        },
        "I" : {
            "linear"    : [ ["I", "-gamma"] ],
            "infection" : [ ["I", "beta"] ]
        },
        "test_pos"  : [ "p_falsepos", "p_truepos", "p_falsepos"] ,
        "test_freq" : [ "tf", "tf", "tf"]
    }
>>> parameters = {
        'beta': 0.1,
        'gamma': 0.1,
        'seed': 1234,
        'p_falsepos': 0,
        'p_truepos': 1,
        'tf': 1
    }
model_class_data()
Parameters:data (dict) – The object returned by simulate.
Returns:
Return type:The population of class model_class_key as a time series
simulate()

Performs the Stochastic Simulation Algorithm (SSA)

Parameters:
  • x0 (np.array) – Initial condition.
  • contactMatrix (python function(t)) – The social contact matrix C_{ij} denotes the average number of contacts made per day by an individual in class i with an individual in class j
  • testRate (python function(t)) – The total number of PCR tests performed per day
  • Tf (float) – Final time of integrator
  • Nf (Int) – Number of time points to evaluate.
  • method (str, optional) – SSA to use, either ‘gillespie’ or ‘tau_leaping’. The default is ‘gillespie’.
  • nc (TYPE, optional) –
  • epsilon (float, optional) –

    The acceptable relative change of the rates during each tau-leaping step, as defined in Cao et al:

    The default is 0.03

  • tau_update_frequency (TYPE, optional) –
Returns:

X: output path from integrator, t : time points evaluated at, ‘event_occured’ , ‘param’: input param to integrator.
Return type:

dict

SIR

class pyross.stochastic.SIR

Susceptible, Infected, Removed (SIR)

  • Ia: asymptomatic
  • Is: symptomatic

…

Parameters:
  • parameters (dict) –

    Contains the following keys:

    alpha: float, np.array (M,)
    fraction of infected who are asymptomatic.
    beta: float
    rate of spread of infection.
    gIa: float
    rate of removal from asymptomatic individuals.
    gIs: float
    rate of removal from symptomatic individuals.
    fsa: float
    Fraction by which symptomatic individuals do not self-isolate.
    seed: long
    seed for pseudo-random number generator (optional).
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(3*M, )) – Initial number in each compartment and class
rate_matrix:

Calculates the rate constant for each reaction channel.

simulate:

Performs stochastic numerical integration.

simulate()

Performs the Stochastic Simulation Algorithm (SSA)

Parameters:
  • S0 (np.array) – Initial number of susceptables.
  • Ia0 (np.array) – Initial number of asymptomatic infectives.
  • Is0 (np.array) – Initial number of symptomatic infectives.
  • contactMatrix (python function(t)) – The social contact matrix C_{ij} denotes the average number of contacts made per day by an individual in class i with an individual in class j
  • Tf (float) – Final time of integrator
  • Nf (Int) – Number of time points to evaluate.
  • method (str, optional) – SSA to use, either ‘gillespie’ or ‘tau_leaping’. The default is ‘gillespie’.
  • nc (TYPE, optional) –
  • epsilon (float, optional) –

    The acceptable relative change of the rates during each tau-leaping step, as defined in Cao et al:

    The default is 0.03

  • tau_update_frequency (TYPE, optional) –
Returns:

X: output path from integrator, t : time points evaluated at, ‘event_occured’ , ‘param’: input param to integrator.
Return type:

dict

SIkR

class pyross.stochastic.SIkR

Susceptible, Infected, Removed (SIkR). Method of k-stages of I

…

Parameters:
  • parameters (dict) –

    Contains the following keys:

    beta: float
    rate of spread of infection.
    gI: float
    rate of removal from infectives.
    fsa: float
    Fraction by which symptomatic individuals do not self-isolate.
    kI: int
    number of stages of infection.
    seed: long
    seed for pseudo-random number generator (optional).
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array((kI + 1)*M, )) – Initial number in each compartment and class

SEIR

class pyross.stochastic.SEIR

Susceptible, Exposed, Infected, Removed (SEIR)

  • Ia: asymptomatic
  • Is: symptomatic
  • E: exposed

…

Parameters:
  • parameters (dict) –

    Contains the following keys:

    alpha: float, np.array (M,)
    fraction of infected who are asymptomatic.
    beta: float
    rate of spread of infection.
    gIa: float
    rate of removal from asymptomatic individuals.
    gIs: float
    rate of removal from symptomatic individuals.
    fsa: float
    Fraction by which symptomatic individuals do not self-isolate.
    gE: float
    rate of removal from exposed individuals.
    seed: long
    seed for pseudo-random number generator (optional).
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(4*M, )) – Initial number in each compartment and class

SEAIRQ

class pyross.stochastic.SEAIRQ

Susceptible, Exposed, Asymptomatic and infected, Infected, Removed, Quarantined (SEAIRQ)

  • Ia: asymptomatic
  • Is: symptomatic
  • E: exposed
  • A: asymptomatic and infectious
  • Q: quarantined

…

Parameters:
  • parameters (dict) –

    Contains the following keys:

    alpha: float, np.array(M,)
    fraction of infected who are asymptomatic.
    beta: float
    rate of spread of infection.
    gIa: float
    rate of removal from asymptomatic individuals.
    gIs: float
    rate of removal from symptomatic individuals.
    gE: float
    rate of removal from exposed individuals.
    gA: float
    rate of removal from activated individuals.
    fsa: float
    Fraction by which symptomatic individuals do not self-isolate.
    tE : float
    testing rate and contact tracing of exposeds
    tA : float
    testing rate and contact tracing of activateds
    tIa: float
    testing rate and contact tracing of asymptomatics
    tIs: float
    testing rate and contact tracing of symptomatics
    seed: long
    seed for pseudo-random number generator (optional).
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(6*M, )) – Initial number in each compartment and class

SEAIRQ_testing

class pyross.stochastic.SEAIRQ_testing

Susceptible, Exposed, Asymptomatic and infected, Infected, Removed, Quarantined (SEAIRQ)

  • E: exposed
  • A: Asymptomatic and infectious
  • Ia: asymptomatic
  • Is: symptomatic
  • Q: quarantined

…

Parameters:
  • parameters (dict) –

    Contains the following keys:

    alpha: float, np.array(M,)
    fraction of infected who are asymptomatic.
    beta: float
    rate of spread of infection.
    gIa: float
    rate of removal from asymptomatic individuals.
    gIs: float
    rate of removal from symptomatic individuals.
    gE: float
    rate of removal from exposed individuals.
    gA: float
    rate of removal from activated individuals.
    fsa: float
    Fraction by which symptomatic individuals do not self-isolate.
    ars: float
    fraction of population admissible for random and symptomatic tests
    kapE: float
    fraction of positive tests for exposed individuals
    seed: long
    seed for pseudo-random number generator (optional).
  • M (int) – Number of compartments of individual for each class. I.e len(contactMatrix)
  • Ni (np.array(6*M, )) – Initial number in each compartment and class
  • testRate (python function(t)) – number of tests per day and age group

stochastic_integration

class pyross.stochastic.stochastic_integration
Integrators used by stochastic models:
Gillespie and tau-leaping
A()
Parameters:data (Data dict) –
Returns:A
Return type:Activated population time series
E()
Parameters:data (Data dict) –
Returns:E
Return type:Exposed population time series
I()
Parameters:data (Data dict) –
Returns:Ia
Return type:Asymptomatics population time series
Ia()
Parameters:data (Data dict) –
Returns:Ia
Return type:Asymptomatics population time series
Is()
Parameters:data (Data dict) –
Returns:Is
Return type:symptomatics population time series
R()
Parameters:data (Data dict) –
Returns:R
Return type:Removed population time series
S()
Parameters:data (Data dict) –
Returns:S
Return type:Susceptible population time series
Sx()
Parameters:data (Data dict) –
Returns:
Return type:Generic compartment Sx
check_for_event()
simulate_gillespie()

Performs the stochastic simulation using the Gillespie algorithm.

  1. Rates for each reaction channel r_i calculated from current state.
  2. The timestep tau is chosen randomly from an exponential distribution P ~ e^(-W tau).
  3. A single reaction occurs with probablity proportional to its fractional rate constant r_i/W.
  4. The state is updated to reflect this reaction occuring and time is propagated forward by tau

Stops if population becomes too small.

Parameters:
  • contactMatrix (python function(t)) – The social contact matrix C_{ij} denotes the average number of contacts made per day by an individual in class i with an individual in class j
  • Tf (float) – Final time of integrator
  • Nf (Int) – Number of time points to evaluate.
Returns:

  • t_arr (np.array(Nf,)) – Array of time points at which the integrator was evaluated.
  • out_arr (np.array) – Output path from integrator.
simulate_tau_leaping()

Tau leaping algorithm for producing stochastically correct trajectories Based on Cao et al (2006): https://doi.org/10.1063/1.2159468 This method can run much faster than the Gillespie algorithm

  1. Rates for each reaction channel r_i calculated from current state.
  2. Timestep tau chosen such that Delta r_i < epsilon Sum r_i
  3. Number of reactions that occur in channel i ~Poisson(r_i tau)
  4. Update state by this amount
Parameters:
  • contactMatrix (python function(t)) – The social contact matrix C_{ij} denotes the average number of contacts made per day by an individual in class i with an individual in class j
  • Tf (float) – Final time of integrator
  • Nf (Int) – Number of time points to evaluate.
  • nc (optional) – The default is 30
  • epsilon (float, optional) – The acceptable relative change of the rates during each tau-leaping step, as defined in Cao et al. The default is 0.03
  • tau_update_frequency (optional) –
Returns:

  • t_arr (np.array(Nf,)) – Array of time points at which the integrator was evaluated.
  • out_arr (np.array) – Output path from integrator.