Metapopulation reaction-diffusion

SIR epidemics in subpopulations at the nodes of a network, linked by individuals who travel.

• Name: metapopulation
• Field: biology/epidemiology
• Formalism: ordinary differential equations
• References: Colizza, Pastor-Satorras, and Vespignani 2007; Pastor-Satorras et al. 2015

Factory

@register_model(
    field='biology/epidemiology',
    tags=('network', 'metapopulation', 'reaction-diffusion', 'ode', 'field'),
    references=('Colizza, Pastor-Satorras, and Vespignani 2007', 'Pastor-Satorras et al. 2015'),
    defaults={'until': 150.0, 'dt': 1.0},
)
def metapopulation(
    graph: object = None,
    *,
    population: float = 1000.0,
    infected: float = 10.0,
    beta: float = 0.3,
    gamma: float = 0.1,
    D: float = 0.01,
) -> System:
    """Return SIR epidemics in subpopulations at the nodes of a network, linked by individuals who travel.

    The infection starts in the first node and spreads to the others with the travellers. The total number of
    individuals stays constant.

    Args:
        graph:
            The NetworkX graph, or `None` for a Barabási-Albert graph of 30 nodes.

        population:
            The initial number of individuals of every node.

        infected:
            The initial number of infected individuals in the first node.

        beta:
            The infection rate.

        gamma:
            The recovery rate.

        D:
            The rate at which an individual hops along one edge.

    Returns:
        system:
            The model, with the compartments `cities.S`, `cities.I`, and `cities.R`, one value per node.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.biology.epidemiology import metapopulation
        >>> metapopulation(nx.path_graph(4)).elements[0].check().prepare().size_
        4
    """
    graph = nx.barabasi_albert_graph(30, 2, seed=0) if graph is None else graph
    size = nx.number_of_nodes(graph)
    sick = np.zeros(size)
    sick[0] = infected
    values = {
        'S': np.full(size, population) - sick,
        'I': sick,
        'R': np.zeros(size),
        'beta': beta,
        'gamma': gamma,
        'D': D,
    }
    return System(
        elements=[FieldElement('cities', element=SubpopulationElement, space=NetworkSpace(graph), values=values)],
    )

Components

class SubpopulationElement(Element):
    """Well-mixed SIR population at a node of a network, whose individuals hop along each edge at the rate `D`.

    Inside the node, `beta S I / N` individuals are infected and `gamma I` recover per unit time, with
    `N = S + I + R`. Hopping along every edge adds `D` times the graph Laplacian to every compartment.

    Args:
        label:
            The label of the subpopulation or the field.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        S:
            The initial number of susceptible individuals.

        I:
            The initial number of infected individuals.

        R:
            The initial number of recovered individuals.

        beta:
            The infection rate.

        gamma:
            The recovery rate.

        D:
            The rate at which an individual hops along one edge.

    Examples:
        >>> from skcomplex.models.biology.epidemiology import SubpopulationElement
        >>> SubpopulationElement('city', S=1000.0).beta
        0.3
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        S: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1000.0,
        I: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        R: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
        beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
        gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    ) -> None:
        super().__init__(label, visual=visual)
        self.S = S
        self.I = I
        self.R = R
        self.beta = beta
        self.gamma = gamma
        self.D = D

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the infections, recoveries, and hops of the node.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        infections = s.beta * s.S * s.I / (s.S + s.I + s.R)
        recoveries = s.gamma * s.I
        return [
            RateLaw(s.S, -infections + s.D * laplacian(s.S)),
            RateLaw(s.I, infections - recoveries + s.D * laplacian(s.I)),
            RateLaw(s.R, recoveries + s.D * laplacian(s.R)),
        ]

Equations

cities.S[i]: diff(S_cities(t), t) == -beta_cities*I_cities(t)*S_cities(t)/(I_cities(t) + R_cities(t) + S_cities(t)) + D_cities*laplacian(S_cities(t))
cities.I[i]: diff(I_cities(t), t) == -gamma_cities*I_cities(t) + beta_cities*I_cities(t)*S_cities(t)/(I_cities(t) + R_cities(t) + S_cities(t)) + D_cities*laplacian(I_cities(t))
cities.R[i]: diff(R_cities(t), t) == gamma_cities*I_cities(t) + D_cities*laplacian(R_cities(t))