SIRS

An epidemic whose immunity wanes, which returns in damped waves to an endemic equilibrium.

• Name: sirs
• Field: biology/epidemiology
• Formalism: ordinary differential equations
• References: Kermack and McKendrick 1932; Argun et al. 2021

Factory

@register_model(
    field='biology/epidemiology',
    tags=('compartments', 'ode', 'jump', 'endemic', 'waning immunity'),
    references=('Kermack and McKendrick 1932', 'Argun et al. 2021'),
    defaults={'until': 300.0, 'dt': 1.0},
)
def sirs(
    susceptible: float = 990.0,
    infected: float = 10.0,
    recovered: float = 0.0,
    *,
    beta: float = 0.0003,
    gamma: float = 0.1,
    xi: float = 0.02,
    stochastic: bool = False,
) -> System:
    """Return an epidemic whose immunity wanes, which returns in damped waves to an endemic equilibrium.

    The equilibrium has `S = gamma / beta` susceptible and `I = xi (N - S) / (gamma + xi)` infected.

    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.

        xi:
            The rate at which immunity wanes.

        stochastic:
            Whether infections, recoveries, and losses of immunity are random events of single individuals.

    Returns:
        system:
            The model, with the compartments `S`, `I`, and `R`.

    Examples:
        >>> from skcomplex.models.biology.epidemiology import sirs
        >>> sirs(xi=0.1).interactions[-1].xi
        0.1
    """
    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),
            ImmunityLossInteraction('waning', recovered='R', susceptible='S', xi=xi, 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)

class ImmunityLossInteraction(ReactionInteraction):
    """Loss of immunity, the reaction `R -> S` at the rate `xi R` from recovered back to susceptible.

    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.

        recovered:
            The path of the recovered compartment.

        susceptible:
            The path of the susceptible compartment.

        xi:
            The rate at which immunity wanes.

        stochastic:
            Whether losses of immunity are random events.

    Examples:
        >>> from skcomplex.models.biology.epidemiology import ImmunityLossInteraction
        >>> ImmunityLossInteraction('waning', recovered='R', susceptible='S', xi=0.05).check().xi
        0.05
    """

    def __init__(
        self,
        label: str,
        *,
        recovered: Annotated[str, RoleSpec()],
        susceptible: Annotated[str, RoleSpec()],
        xi: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
        stochastic: bool = False,
    ) -> None:
        super().__init__(label)
        self.recovered = recovered
        self.susceptible = susceptible
        self.xi = xi
        self.stochastic = stochastic

    def _get_reactants(self) -> Stoichiometry:
        """Return the states `n` of the compartments the reaction consumes."""
        return {(self.recovered, 'n'): 1}

    def _get_products(self) -> Stoichiometry:
        """Return the states `n` of the compartments the reaction produces."""
        return {(self.susceptible, '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 `xi`.

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

Equations

S.n: diff(n_S(t), t) == -beta_infection*n_I(t)*n_S(t) + xi_waning*n_R(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) - xi_waning*n_R(t)