Numeric Models

Provides an API to define epidemiological models.

class epipack.numeric_epi_models.ConstantBirthRate(rate)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters: rate (float) -- Constant rate value

class epipack.numeric_epi_models.ConstantLinearRate(rate, comp0)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters

rate (float) -- Constant rate value
comp0 (int) -- Index of the corresponding reacting component. The incidence of this component will be multiplied with the value of rate.

class epipack.numeric_epi_models.ConstantQuadraticRate(rate, comp0, comp1)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters

rate (float) -- Constant rate value
comp0 (int) -- Index of one of the reacting components. The incidence of this component will be multiplied with the value of rate.
comp1 (int) -- Index of the other reacting component. The incidence of this component will be multiplied with the value of rate.

class epipack.numeric_epi_models.DynamicBirthRate(rate)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters: rate (function) -- Function of time t and state y

class epipack.numeric_epi_models.DynamicLinearRate(rate, comp0)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters

rate (function) -- Function of time t and state y
comp0 (int) -- Index of the corresponding reacting component. The incidence of this component will be multiplied with the value of rate.

class epipack.numeric_epi_models.DynamicQuadraticRate(rate, comp0, comp1)[source]

Bases: object

Will be used as a function of time t and state y, returning a rate value.

Parameters

rate (function) -- Function of time t and state y
comp0 (int) -- Index of one of the reacting components. The incidence of this component will be multiplied with the value of rate.
comp1 (int) -- Index of the other reacting component. The incidence of this component will be multiplied with the value of rate.

class epipack.numeric_epi_models.EpiModel(compartments, initial_population_size=1, correct_for_dynamical_population_size=False, integral_solver='solve_ivp')[source]

Bases: IntegrationMixin

A general class to define a standard mean-field compartmental epidemiological model, based on reaction events.

Parameters

compartments (list of string) -- A list containing compartment strings.
initial_population_size (float, default = 1.0) -- The population size at \(t = 0\).
correct_for_dynamical_population_size (bool, default = False) -- If True, the quadratic coupling terms will be divided by the population size.
integral_solver (str, default = 'solve_ivp') -- Whether or not to use the initial-value solver solve_ivp. to determine a time leap for time-varying rates. If not 'solve_ivp', a Newton-Raphson method will be used on the upper bound of a quad-integrator.

compartments

A list containing strings or hashable types that describe each compartment, (e.g. "S", "I", etc.).

Type: list

compartment_ids

Maps compartments to their indices in self.compartments.

Type: dict

N_comp

Number of compartments (including population number)

Type: int

initial_population_size

The population size at \(t = 0\).

Type: float

correct_for_dynamical_population_size

If True, the quadratic coupling terms will be divided by the sum of all compartments, otherwise they will be divided by the initial population size.

Type: bool

birth_rate_functions

A list of functions that return rate values based on time t and state vector y. Each entry corresponds to an event update in self.birth_event_updates.

Type: list of ConstantBirthRate or DynamicBirthRate

birth_event_updates

A list of vectors. Each entry corresponds to a rate in birth_rate_functions and quantifies the change in individual counts in the compartments.

Type: list of numpy.ndarray

linear_rate_functions

A list of functions that return rate values based on time t and state vector y. Each entry corresponds to an event update in self.linear_event_updates.

Type: list of ConstantLinearRate or DynamicLinearRate

linear_event_updates

A list of vectors. Each entry corresponds to a rate in linear_rate_functions and quantifies the change in individual counts in the compartments.

Type: list of numpy.ndarray

quadratic_rate_functions

A list of functions that return rate values based on time t and state vector y. Each entry corresponds to an event update in self.quadratic_event_updates.

Type: list of ConstantQuadraticRate or DynamicQuadraticRate

quadratic_event_updates

A list of vectors. Each entry corresponds to a rate in quadratic_rate_functions and quantifies the change in individual counts in the compartments.

Type: list of numpy.ndarray

y0

The initial conditions.

Type: numpy.ndarray

rates_have_explicit_time_dependence

Internal switch that's flipped when any functional rate is passed to the model.

Type: bool

rates_have_explicit_time_dependence

Internal switch that's flipped when a non-constant rate is passed to the model.

Type: bool

use_ivp_solver

Whether or not to use the initial-value solver to determine a time leap for time-varying rates. If False, a Newton-Raphson method will be used on the upper bound of a quad-integrator.

Type: bool

Example

>>> epi = EpiModel(["S","I","R"])
>>> print(epi.compartments)
[ "S", "I", "R" ]

add_fission_processes(process_list)[source]

Define linear fission processes between compartments.

Parameters

process_list (list of tuple) --

A list of tuples that contains fission rates in the following format:

[
    ("source_compartment", rate, "target_compartment_0", "target_compartment_1" ),
    ...
]

Example

For pure exponential growth of compartment B.

epi.add_fission_processes([
    ("B", growth_event, "B", "B" ),
])

add_fusion_processes(process_list)[source]

Define fusion processes between compartments.

Parameters

process_list (list of tuple) --

A list of tuples that contains fission rates in the following format:

[
    ("coupling_compartment_0", "coupling_compartment_1", rate, "target_compartment_0" ),
    ...
]

Example

Fusion of reactants "A", and "B" to form "C".

epi.add_fusion_processes([
    ("A", "B", reaction_rate, "C" ),
])

add_linear_events(event_list, allow_nonzero_column_sums=False)[source]: Add linear events without resetting the existing event terms. See epipack.numeric_epi_models.EpiModel.set_linear_events() for docstring.

add_quadratic_events(event_list, allow_nonzero_column_sums=False)[source]: Add quadratic events without resetting the existing event terms. See epipack.numeric_epi_models.EpiModel.set_quadratic_events() for docstring.

add_transition_processes(process_list)[source]

Define the linear transition processes between compartments.

Parameters

process_list (list of tuple) --

A list of tuples that contains transitions events in the following format:

[
    ( source_compartment, rate, target_compartment ),
    ...
]

Example

For an SEIR model.

epi.add_transition_processes([
    ("E", symptomatic_rate, "I" ),
    ("I", recovery_rate, "R" ),
])

add_transmission_processes(process_list)[source]

A wrapper to define quadratic process rates through transmission reaction equations. Note that in stochastic network/agent simulations, the transmission rate is equal to a rate per link. For the mean-field ODEs, the rates provided to this function will just be equal to the prefactor of the respective quadratic terms.

on a network of mean degree \(k_0\), a basic reproduction number \(R_0\), and a recovery rate \(\mu\), you would define the single link transmission process as

("I", "S", R_0/k_0 * mu, "I", "I")

For the mean-field system here, the corresponding reaction equation would read

("I", "S", R_0 * mu, "I", "I")

Parameters

process_list (list of tuple) --

A list of tuples that contains transitions rates in the following format:

[
    ("source_compartment",
     "target_compartment_initial",
     rate,
     "source_compartment",
     "target_compartment_final",
     ),
    ...
]

Example

For an SEIR model.

epi.add_transmission_processes([
    ("I", "S", +1, "I", "E" ),
])

dydt(t, y)[source]

Compute the current momenta of the epidemiological model.

Parameters

t (float) -- Current time
y (numpy.ndarray) -- The entries correspond to the compartment frequencies (or counts, depending on population size).

get_compartment(iC)[source]: Get the compartment, given an integer ID iC

get_compartment_changes(rates)[source]

Sample a state change vector with probability proportional to its rate in rates.

Needed for stochastic simulations.

Parameters: rates (numpy.ndarray) -- A non-zero list of rates. Expects rates to be sorted according to self.birth_event_updates + self.linear_event_updates + self.quadratic_event_updates.
Returns: dy -- A state change vector.
Return type: numpy.ndarray

get_compartment_id(C)[source]: Get the integer ID of a compartment C

get_event_rates(t, y)[source]

Get a list of rate values corresponding to the previously set events.

Parameters

t (float) -- Current time
y (numpy.ndarray) -- Current state vector

Returns

rates -- A list of rate values corresponding to rates. Ordered as birth_rate_functions + linear_rate_functions + quadratic_rate_functions.

Return type

list

get_numerical_dydt()[source]: Return a function that obtains t and y as an input and returns dydt of this system

get_numerical_event_and_rate_functions()[source]

This function is needed to generalize stochastic simulations for child classes.

Returns

get_event_rates (func) -- A function that takes the current time t and state vector y and returns numerical event rate lists.
get_compartment_changes (funx) -- A function that takes a numerical list of event rates and returns a random event state change vector with probability proportional to its entry in rates.

get_time_leap_and_proposed_compartment_changes(t, current_event_rates=None, get_event_rates=None, get_compartment_changes=None, use_ivp_solver=None, rates_have_explicit_time_dependence=None)[source]

For the current event rates, obtain a proposed time leap and concurrent state change vector.

This method is needed for stochastic simulations.

Parameters

t (float) -- current time
current_event_rates (list, default = None) -- A list of constant rate values. Will be ignored if self.rates_have_explicit_time_dependence is True, which is why None is a valid value.
get_event_rates (function, default = None) -- A function that takes time t and current state y as input and computes the rates of all possible events. If None, will attempt to set this to self.get_event_rates().
get_compartment_changes (function, default = None) -- A function that takes computed event rates and returns a random state change with probability proportional to its rate. If None, will attempt to set this to self.get_compartment_changes().
use_ivp_solver (bool, default = None) -- Whether or not to use an initial value problem solver to obtain a time leap in explicitly time-dependent problems. If None, will use the value of the class attribute self.use_ivp_solver.
rates_have_explicit_time_dependence (bool, default = None) -- Whether or not the problem is explicitly time-dependent. If None, will use the value of the class attribute self.rates_have_explicit_time_dependence.

Returns

tau (float) -- A time leap.
dy (numpy.ndarray) -- A state change vector.

set_linear_events(event_list, allow_nonzero_column_sums=False, reset_events=True)[source]

Define the linear transition events between compartments.

Parameters

event_list (list of tuple) --

A list of tuples that contains transition events in the following format:

[
    (
        ("affected_compartment_0",),
        rate,
        [
            ("affected_compartment_0", dN0),
            ("affected_compartment_1", dN1),
            ...
        ],
     ),
    ...
]

allow_nonzero_column_sums (bool, default : False) -- Traditionally, epidemiological models preserve the total population size. If that's not the case, switch off testing for this.
reset_events (bool, default : True) -- Whether to reset all linear events to zero before converting those.

Example

For an SEIR model with infectious period tau and incubation period theta.

epi.set_linear_events([
    ( ("E",),
      1/theta,
      [ ("E", -1), ("I", +1) ]
    ),
    ( ("I",),
      1/tau,
      [ ("I", -1), ("R", +1) ]
    ),
])

Read as "compartment E reacts with rate \(1/\theta\) which leads to the decay of one E particle to one I particle."

set_processes(process_list, allow_nonzero_column_sums=False, reset_events=True, ignore_rate_position_checks=False)[source]

Converts a list of reaction process tuples to event tuples and sets the rates for this model.

Parameters

process_list (list of tuple) --

A list containing reaction processes in terms of tuples.

[
    # transition process
    ( source_compartment, rate, target_compartment),

    # transmission process
    ( coupling_compartment_0, coupling_compartment_1, rate, target_compartment_0, target_ccompartment_1),

    # fission process
    ( source_compartment, rate, target_compartment_0, target_ccompartment_1),

    # fusion process
    ( source_compartment_0, source_compartment_1, rate, target_compartment),

    # death process
    ( source_compartment, rate, None),

    # birth process
    ( None, rate, target_compartment),
]

allow_nonzero_column_sums (bool, default : False) -- Traditionally, epidemiological models preserve the total population size. If that's not the case, switch off testing for this.
reset_events (bool, default : True) -- If this is True, reset all events to zero before setting the new ones.
ignore_rate_position_checks (bool, default = False) -- This function usually checks whether the rate of a reaction is positioned correctly. You can turn this behavior off for transition, birth, death, and transmission processes. (Useful if you want to define symbolic transmission processes that are compartment-dependent).

set_quadratic_events(event_list, allow_nonzero_column_sums=False, reset_events=True, initial_time_for_column_sum_test=0)[source]

Define quadratic transition events between compartments.

Parameters

event_list (list of tuple) --

A list of tuples that contains transmission events in the following format:

[
    (
        ("coupling_compartment_0", "coupling_compartment_1"),
        rate,
        [
            ("affected_compartment_0", dN0),
            ("affected_compartment_1", dN1),
            ...
        ],
     ),
    ...
]

allow_nonzero_column_sums (bool, default : False) -- Traditionally, epidemiological models preserve the total population size. If that's not the case, switch off testing for this.
reset_events (bool, default : True) -- Whether to reset all linear events to zero before converting those.

Example

For an SEIR model with infection rate eta.

epi.set_quadratic_events([
    ( ("S", "I"),
      eta,
      [ ("S", -1), ("E", +1) ]
    ),
])

Read as

"Coupling of S and I leads to the decay of one S particle to one E particle with rate \(\eta\).".

simulate(tmax, return_compartments=None, sampling_dt=None, sampling_callback=None, adopt_final_state=False, use_ivp_solver=None, rates_have_explicit_time_dependence=None, ignore_warnings=False)[source]

Returns values of the given compartments at the demanded time points (as a numpy.ndarray of shape (return_compartments), len(time_points).

If return_compartments is None, all compartments will be returned.

Parameters

tmax (float) -- maximum length of the simulation
return_compartments (list of compartments, default = None:) -- The compartments for which to return time series. If None, all compartments will be returned.
sampling_dt (float, default = None) -- Temporal distance between samples of the compartment counts. If None, every change will be returned.
sampling_callback (funtion, default = None) -- A function that's called when a sample is taken
use_ivp_solver (bool, default = None) -- Wether or not to use an initial value problem solver to obtain a time leap in explicitly time-dependent problems. If None, will use the value of the class attribute self.use_ivp_solver.
rates_have_explicit_time_dependence (bool, default = None) -- Wether or not the problem is explicitly time-dependent. If None, will use the value of the class attribute self.rates_have_explicit_time_dependence.
ignore_warnings (bool, default = False) -- wether or not to raise warnings about unset explicit time.

Returns

t (numpy.ndarray) -- times at which compartment counts have been sampled
result (dict) -- Dictionary mapping a compartment to a time series of its count.

class epipack.numeric_epi_models.SEIRModel(infection_rate, recovery_rate, symptomatic_rate, initial_population_size=1.0)[source]

Bases: EpiModel

An SEIR model derived from epipack.numeric_epi_models.EpiModel.

class epipack.numeric_epi_models.SIModel(infection_rate, initial_population_size=1.0)[source]

Bases: EpiModel

An SI model derived from epipack.numeric_epi_models.EpiModel.

class epipack.numeric_epi_models.SIRModel(infection_rate, recovery_rate, initial_population_size=1.0)[source]

Bases: EpiModel

An SIR model derived from epipack.numeric_epi_models.EpiModel.

class epipack.numeric_epi_models.SIRSModel(infection_rate, recovery_rate, waning_immunity_rate, initial_population_size=1.0)[source]

Bases: EpiModel

An SIRS model derived from epipack.numeric_epi_models.EpiModel.

class epipack.numeric_epi_models.SISModel(infection_rate, recovery_rate, initial_population_size=1.0)[source]

Bases: EpiModel

An SIS model derived from epipack.numeric_epi_models.EpiModel.

Parameters

R0 (float) -- The basic reproduction number. From this, the infection rate is computed as infection_rate = R0 * recovery_rate
recovery_rate (float) -- Recovery rate
population_size (float, default = 1.0) -- Number of people in the population.

epipack.numeric_epi_models.custom_choice(p)[source]: Return an index of normalized probability vector p with probability equal to the corresponding entry in p.