SIR
An epidemic that moves individuals from susceptible to infected to recovered.
- Name:
sir - Field:
biology/epidemiology - Formalism: ordinary differential equations
- References: Kermack and McKendrick 1927
Factory
@register_model(
field='biology/epidemiology',
tags=('compartments', 'ode', 'jump'),
references=('Kermack and McKendrick 1927',),
defaults={'until': 160.0, 'dt': 1.0},
)
def sir(
susceptible: float = 990.0,
infected: float = 10.0,
recovered: float = 0.0,
*,
beta: float = 0.0003,
gamma: float = 0.1,
stochastic: bool = False,
) -> System:
"""Return an epidemic that moves individuals from susceptible to infected to recovered.
Args:
susceptible:
The initial number of susceptible individuals.
infected:
The initial number of infected individuals.
recovered:
The initial number of recovered individuals.
beta:
The infection rate per susceptible-infected pair.
gamma:
The recovery rate.
stochastic:
Whether infections and recoveries are random events of single individuals.
Returns:
system:
The model, with the compartments `S`, `I`, and `R`.
Examples:
>>> from skcomplex.models.biology.epidemiology import sir
>>> [element.label for element in sir().elements]
['S', 'I', 'R']
"""
return System(
elements=[
CompartmentElement('S', n=susceptible),
CompartmentElement('I', n=infected),
CompartmentElement('R', n=recovered),
],
interactions=[
InfectionInteraction('infection', susceptible='S', infected='I', beta=beta, stochastic=stochastic),
RecoveryInteraction('recovery', infected='I', recovered='R', gamma=gamma, stochastic=stochastic),
],
)
Components
class CompartmentElement(Element):
"""Compartment, the amount of units in one named state, as a count or a fraction.
The label of the element names the compartment, such as `'A'`, `'B'`, or `'source'`. A compartment has no laws
of its own, and reactions and flows between compartments move its units.
Args:
label:
The label of the compartment, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
n:
The initial amount, a nonnegative count or fraction.
Attributes:
default_visual_:
The visual of the element when it sets none.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> system = System(elements=[CompartmentElement('A', n=990), CompartmentElement('B', n=10)])
>>> Simulation(system, dt=1.0).check().prepare().run(until=2).trajectory_.values('A.n').tolist()
[990.0, 990.0, 990.0]
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.n = n
class InfectionInteraction(ReactionInteraction):
"""Mass-action infection, the reaction `S + I -> 2 I` at the rate `beta S I` from susceptible to infected.
It is a `ReactionInteraction` between the states `n` of the two compartments its roles bind, a rate equation, or
in stochastic mode a random event that moves one individual.
Args:
label:
The label of the interaction.
susceptible:
The path of the susceptible compartment.
infected:
The path of the infected compartment.
beta:
The transmission rate.
stochastic:
Whether infections are random events.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.epidemiology import InfectionInteraction
>>> system = System(
... elements=[CompartmentElement('S', n=99.0), CompartmentElement('I', n=1.0)],
... interactions=[InfectionInteraction('infection', susceptible='S', infected='I', beta=0.01)],
... )
>>> trajectory = Simulation(system, dt=0.1).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('S.n')[-1] + trajectory.values('I.n')[-1]), 6)
100.0
"""
def __init__(
self,
label: str,
*,
susceptible: Annotated[str, RoleSpec()],
infected: Annotated[str, RoleSpec()],
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.001,
stochastic: bool = False,
) -> None:
super().__init__(label)
self.susceptible = susceptible
self.infected = infected
self.beta = beta
self.stochastic = stochastic
def _get_reactants(self) -> Stoichiometry:
"""Return the states `n` of the compartments the reaction consumes."""
return {(self.susceptible, 'n'): 1, (self.infected, 'n'): 1}
def _get_products(self) -> Stoichiometry:
"""Return the states `n` of the compartments the reaction produces."""
return {(self.infected, 'n'): 2}
def _get_mode(self) -> str:
"""Return the mode of the reaction, random events when stochastic and a rate equation otherwise."""
return MODES[self.stochastic]
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the mass action of the reaction at the rate constant `beta`.
Args:
namespace:
The namespace of the symbols the rate can use, with the compartments.
Returns:
rate:
The propensity of the reaction.
"""
return self._get_mass_action(namespace, namespace.beta)
class RecoveryInteraction(ReactionInteraction):
"""Recovery, the reaction `I -> R` at the rate `gamma I` from infected to recovered.
It is a `ReactionInteraction` between the states `n` of the two compartments its roles bind, a rate equation, or
in stochastic mode a random event that moves one individual.
Args:
label:
The label of the interaction.
infected:
The path of the infected compartment.
recovered:
The path of the recovered compartment.
gamma:
The recovery rate.
stochastic:
Whether recoveries are random events.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.epidemiology import RecoveryInteraction
>>> system = System(
... elements=[CompartmentElement('I', n=10.0), CompartmentElement('R')],
... interactions=[RecoveryInteraction('recovery', infected='I', recovered='R', gamma=0.5)],
... )
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('I.n')[-1]), 3)
3.679
"""
def __init__(
self,
label: str,
*,
infected: Annotated[str, RoleSpec()],
recovered: Annotated[str, RoleSpec()],
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
stochastic: bool = False,
) -> None:
super().__init__(label)
self.infected = infected
self.recovered = recovered
self.gamma = gamma
self.stochastic = stochastic
def _get_reactants(self) -> Stoichiometry:
"""Return the states `n` of the compartments the reaction consumes."""
return {(self.infected, 'n'): 1}
def _get_products(self) -> Stoichiometry:
"""Return the states `n` of the compartments the reaction produces."""
return {(self.recovered, 'n'): 1}
def _get_mode(self) -> str:
"""Return the mode of the reaction, random events when stochastic and a rate equation otherwise."""
return MODES[self.stochastic]
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the mass action of the reaction at the rate constant `gamma`.
Args:
namespace:
The namespace of the symbols the rate can use, with the compartments.
Returns:
rate:
The propensity of the reaction.
"""
return self._get_mass_action(namespace, namespace.gamma)
Equations
S.n: diff(n_S(t), t) == -beta_infection*n_I(t)*n_S(t)
I.n: diff(n_I(t), t) == beta_infection*n_I(t)*n_S(t) - gamma_recovery*n_I(t)
R.n: diff(n_R(t), t) == gamma_recovery*n_I(t)