Networks

🧠20% 🤖80%

Core Starsim’s networks (random contacts, partnerships, mother and child, mixing pools, and contact matrices), and how to write a network of your own, are on Mixing. The Starsim library (import starsim.library as ssl; see Library) has more, for households and other clusters, for space, and for theoretical work:

Class	Based on	What it is	Main parameters
ssl.ClusterNet	ss.Network	Fully connected clusters, such as households	cluster_size, who, dynamic
ssl.HybridNet	ss.Network	Covasim’s households, schools, workplaces, and community	household_size, n_contacts, school_ages, work_ages, clustered
ssl.HouseholdNet	ss.Network	Households built from survey data in the format of the Demographic and Health Surveys (DHS), optionally changing over time	dhs_data, dynamic, prob_move_out
ssl.DiskNet	ss.Network	People who move in a unit square, in contact with everyone within a radius	r, v
ssl.StaticNet	ss.Network	A fixed network from a NetworkX graph, or a function that makes one	graph, p, n_contacts
ssl.GeoNet	ss.RandomNet	People in real places, with most contacts near home and the rest in places that are large and near	places, level, local, c, n_contacts
ssl.GridNet	ssl.StaticNet	People on a square grid, in contact with their neighbors	diagonal, periodic
ssl.SmallWorldNet	ssl.StaticNet	A ring with a few long-range contacts (Watts-Strogatz)	n_contacts, rewire
ssl.ScaleFreeNet	ssl.StaticNet	A few people with many contacts (Barabási-Albert)	n_contacts
ssl.NullNet	ss.Network	Each person in contact only with themselves, with no transmission, as a placeholder	none

Like every network, these carry agents, so a model with one runs under the agent-based solver; Mixing says how each would be approximated by a compartmental one. The tutorial T4 – Mixing introduces networks with a worked example.

Households and other clusters

Much transmission happens within households, schools, and workplaces: small groups whose members meet each other again and again. ssl.ClusterNet divides people (optionally only those who= gives) into clusters of sizes drawn from cluster_size (by default Poisson with a mean of 4), with everyone in a cluster in contact with everyone else. Its contacts are made once, at the start (or redrawn every step, with dynamic=True).

Here the same SIR disease spreads with about 6 contacts per person: all of them random, or 1 random and the rest in households of about 6, with each scenario run 3 times:

import numpy as np
import matplotlib.pyplot as plt
import starsim as ss
import starsim.library as ssl

sir = lambda: ss.SIR(beta=ss.perday(0.05), dur_inf=ss.days(10), init_prev=0.005)
kw = dict(n_agents=2000, start='2000-01-01', dur=ss.days(200), dt=ss.days(1))
scenarios = {
    'Random contacts only': ss.Sim(diseases=sir(), networks=ss.RandomNet(n_contacts=6), **kw),
    'Households and random contacts': ss.Sim(diseases=sir(), networks=[ssl.ClusterNet(cluster_size=ss.poisson(6)), ss.RandomNet(n_contacts=1)], **kw),
}
msim = ss.parallel(scenarios, n_runs=3, workers=1)
msim.plot(['sir.I', 'sir.R'])

MultiSim: 6 runs of 2 scenarios (Random contacts only, Households and random contacts)
  Done in 2.6 s

ssl.HybridNet is the population of Covasim: households (h, a ClusterNet with sizes drawn from a Poisson distribution with a mean of 2, as in Covasim, which gives about 2 contacts per person), schools for ages 6 to 22 (s), workplaces for ages 22 to 65 (w), and the community (c). Schools and workplaces are clusters too, and the community is a random network, with a Poisson number of contacts in each (n_contacts, by default ss.poisson(20), ss.poisson(16), and ss.poisson(20); a number is the same for everyone, as in ss.RandomNet); clustered=False makes schools and workplaces random networks, as in Covasim. When added to a sim, it becomes four networks, each reached by name. How much a contact in each transmits is the disease’s, with one β per network: Covasim’s relative weights are ssl.HybridNet.rel_betas, 3, 0.6, 0.6, and 0.3, so that a contact at home counts for 10 times one in the community. ssl.HybridNet.betas() multiplies them by a β per contact (by default Covasim’s, 0.016 a day), giving a dict to use as a disease’s β: for a β of 0.02 a day, ss.SIR(beta=ssl.HybridNet.betas(ss.perday(0.02))) is ss.SIR(beta=dict(h=ss.perday(0.06), s=ss.perday(0.012), w=ss.perday(0.012), c=ss.perday(0.006))).

The layers make questions about settings answerable. Here an influenza-like infection spreads with schools open, and with schools closed for the whole outbreak (a β of 0 on the school layer):

flu = lambda **beta: ss.SIR(beta=ssl.HybridNet.betas(ss.perday(0.02)) | beta, dur_inf=ss.days(5), init_prev=0.005) # Covasim's weights: 10 times as much at home as in the community
kw = dict(n_agents=2000, start='2000', dur=ss.days(120), dt=ss.days(1))
scenarios = {
    'Schools open':   ss.Sim(diseases=flu(), networks=ssl.HybridNet(), **kw),
    'Schools closed': ss.Sim(diseases=flu(s=ss.perday(0)), networks=ssl.HybridNet(), **kw),
}
msim = ss.parallel(scenarios, n_runs=5, workers=1)

sim = msim.sims[0]
for name, net in sim.networks.items():
    print(name, 'contacts per person:', round(2*len(net) / sim.people.alive.count(), 1))
msim.plot('sir.I')

MultiSim: 10 runs of 2 scenarios (Schools open, Schools closed)
  Done in 3.5 s

h contacts per person: 2.1
s contacts per person: 5.7
w contacts per person: 8.7
c contacts per person: 19.9

A closure for a while, rather than the whole outbreak, is a modifier on the school layer, ss.modify(sir.pars.beta, multiply=0, network='s') with the times it is closed (see Interventions).

Households from survey data

ssl.HouseholdNet builds households from DHS household survey data, setting each agent’s age and sex to match a randomly drawn real household, so that the ages of the people living together are realistic (children with their parents, grandparents with both). In a sim with ss.Pregnancy, households also change (dynamic=True, which needs it; without it, they are fixed): one woman is the head of each household, pregnant women other than the head may move out to form new households (with probability prob_move_out), and babies join their mother’s household. To change the households less often than every step, give the network a step of its own, e.g. dt=ss.years(1).

Real DHS data requires registration, so dhs_data='default' makes synthetic data with the same structure, for demonstrations and testing:

sim = ss.Sim(networks=ssl.HouseholdNet(dhs_data='default'), demographics=ss.Pregnancy(), start='2000', stop='2010', dt=ss.months(1), n_agents=2000)
sim.run()

net = sim.networks.householdnet # The sim's own copy of the network
hh_ids = net.household_ids[net.household_ids.notnan] # Unborn agents have no household
sizes = np.bincount(hh_ids.astype(int))
sizes = sizes[sizes > 0]
print('Households:', len(sizes))
print('Mean household size:', round(sizes.mean(), 1))
print('Largest household:', sizes.max())

Sim: pregnancy, householdnet, as 2,000 agents, 2000-01-01 to 2010-01-01
  Running as agents: n_agents was given
  Note: pregnancy: women pregnant at the start, conceived over the 9 steps before it (burnin=True)

  Done in 0.56 s

Households: 761
Mean household size: 3.3
Largest household: 9

To use real data, download a Household Recode (HR) dataset in Stata format, and load it with HouseholdNet.load_dhs(); any table with a household ID (hh_id) and the ages of its members as text (ages, such as '72, 17, 30'), and optionally their sexes (sexes, 1 for male and 2 for female), will also do:

dhs_data = ssl.HouseholdNet.load_dhs('XXHR7xDT/XXHR7xFL.DTA')
sim = ss.Sim(networks=ssl.HouseholdNet(dhs_data=dhs_data))

Warning

HouseholdNet sets the age and sex of every agent (other than unborn agents) when it starts. Take care when using it with other modules that set or depend on age and sex, such as ss.population('Kenya') (see People and demographics).

Space

ssl.DiskNet is a spatial network: each agent has a position (x, y) in a unit square and a heading, moves at speed v each step (by default 0.05 a day), bounces off the walls, and is in contact with every agent within radius r (by default 0.1). Because contacts depend on position, an epidemic spreads outwards from each case as a wave. The positions are states of the network, so they can be plotted at any point in a run; here the sim runs in stages, and plots who is susceptible, infected, and recovered at each:

sir = lambda: ss.SIR(beta=ss.perday(0.1), dur_inf=ss.days(10), init_n=3)
kw = dict(n_agents=2000, start='2000-01-01', dur=ss.days(150), dt=ss.days(1))
disk = ss.Sim(diseases=sir(), networks=ssl.DiskNet(r=0.04, v=ss.freqperday(0.001)), **kw)

states = dict(S=('Susceptible', 'lightgray'), I=('Infected', 'tab:red'), R=('Recovered', 'tab:blue'))
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
for ax, day in zip(axes, ['2000-01-31', '2000-03-01', '2000-04-15']):
    disk.run(until=day) # Run up to this day, then continue from here the next time
    net, sir_states = disk.networks.disknet, disk.diseases.sir
    for state, (label, color) in states.items():
        uids = getattr(sir_states, state).uids
        ax.scatter(net.x[uids], net.y[uids], s=4, color=color, label=label)
    ax.set_title(day)
    ax.set_aspect('equal')
    ax.set_xticks([])
    ax.set_yticks([])
axes[0].legend(loc='lower left')
plt.show()

Sim: sir, as 2,000 agents, 2000-01-01 to 2000-05-30
  Running as agents: n_agents was given

  ━━━━━━────────────────────────  19%  2000-01-29  0.8 s

  ━━━━━━━━━━━━──────────────────  38%  2000-02-27  0.7 s

  ━━━━━━━━━━━━━━━━━━━━──────────  68%  2000-04-12  1.5 s

With a radius of 0.04, each agent has about 10 contacts (π × 0.04² × 2000); a random network with as many gives a much faster, higher epidemic, since every contact can be anywhere.

For movement between distinct places (villages, districts, countries) rather than continuous space, see Places and data.

Theoretical networks

ssl.StaticNet makes a fixed network from a NetworkX graph: either a graph made beforehand, with one node for each agent, or a function that makes one, given with its arguments (n, the number of nodes, is filled in, and seed=True uses the sim’s random numbers). It is the way to use any of the many graphs that NetworkX can make: small-world, scale-free, regular, and so on. The most common have classes of their own: ssl.SmallWorldNet (Watts–Strogatz), ssl.ScaleFreeNet, and ssl.GridNet. Here a small-world graph (each person in contact with their 10 nearest neighbors on a ring, with 5% of contacts rewired) is compared with a random regular graph, made by StaticNet from NetworkX’s function, in which each person also has exactly 10 contacts but anywhere:

import networkx as nx

sir = lambda: ss.SIR(beta=ss.perday(0.05), dur_inf=ss.days(10), init_prev=0.005)
kw = dict(n_agents=2000, start='2000-01-01', dur=ss.days(200), dt=ss.days(1))
scenarios = {
    'Small world': ss.Sim(diseases=sir(), networks=ssl.SmallWorldNet(n_contacts=10, rewire=0.05), **kw),
    'Random regular': ss.Sim(diseases=sir(), networks=ssl.StaticNet(graph=nx.random_regular_graph, d=10), **kw),
}
msim = ss.parallel(scenarios, workers=1)
msim.plot('sir.I')

MultiSim: 2 runs of 2 scenarios (Small world, Random regular)
  Done in 0.64 s

With no graph, ssl.StaticNet is the classic random (Erdős–Rényi) graph: each possible pair of people is connected with probability p (or n_contacts over the population, by default 10 contacts each), so each person’s number of contacts is binomial; it is useful for comparison with analytic results. For random contacts drawn again on every step, ss.RandomNet is the one to use. ssl.NullNet connects each agent only to themselves, with a weight of 0, so it never transmits: a placeholder while developing a model, or for testing a disease’s natural history alone.

Advanced: contacts when people are added

ss.RandomNet(method='safe') is a random network in which each agent’s contacts depend only on that agent’s own random numbers (see Random numbers). Here the same network is made for 1,000 and 1,001 agents, counting how many of the contacts among the first 1,000 survive the extra person:

def contacts(net, n_agents):
    """ The contacts among the first 1,000 agents after one step, as a set of pairs """
    sim = ss.Sim(networks=net, n_agents=n_agents, start='2000-01-01', dur=ss.days(1), dt=ss.days(1), verbose=0)
    sim.run()
    edges = sim.networks[0].to_df()
    edges = edges[(edges.p1 < 1000) & (edges.p2 < 1000)]
    return set(zip(np.minimum(edges.p1, edges.p2), np.maximum(edges.p1, edges.p2)))

networks = {
    'ss.RandomNet': lambda: ss.RandomNet(n_contacts=10),
    'ss.RandomNet (exact)': lambda: ss.RandomNet(n_contacts=10, method='exact'),
    'ss.RandomNet (safe)': lambda: ss.RandomNet(n_contacts=10, method='safe'),
}
for label, net in networks.items():
    before = contacts(net(), 1000)
    after = contacts(net(), 1001)
    print(label, 'contacts kept (%):', round(100*len(before & after)/len(before)))

ss.RandomNet contacts kept (%): 51
ss.RandomNet (exact) contacts kept (%): 1
ss.RandomNet (safe) contacts kept (%): 100

For comparing scenarios that add or remove people (births, deaths, migration), the safe network keeps the comparison sharp; for everything else, ss.RandomNet is faster. The appendix says which of the library’s networks keep their random numbers in this way.

Further reading

• NetworkX: graph generators: the graphs that ssl.StaticNet can use
• EpiModel: network modeling for epidemics: epidemics on dynamic networks estimated from data (exponential random graph models)
• Covasim: tutorial on people and networks: the hybrid population that ssl.HybridNet comes from
• DHS Program: data: the household surveys that ssl.HouseholdNet reads
• Mesa: agent-based modeling in Python: agents in space and on networks, more generally

Appendix: every library network

Every network in ssl.networks, generated from the classes themselves: the class it builds on, its parameters and arguments with their defaults, and what its documentation says about common random numbers (“CRN safe”: whether adding or removing a person leaves everyone else’s contacts unchanged; see Random numbers). The networks built on ssl.StaticNet are fixed when the sim starts, so people born later have no contacts in them.

import re
import inspect
import pandas as pd

def crn_note(cls):
    """ The sentence of the docstring that says whether the network keeps common random numbers """
    sentences = re.split(r'(?<=\.)\s+', ' '.join(inspect.getdoc(cls).split()))
    found = [s for s in sentences if 'CRN safe' in s]
    return found[0].replace('Note: ', '').rstrip('.') if found else ''

def defaults(cls):
    """ The network's own parameters and arguments, with their defaults """
    net = cls(dhs_data='default') if cls is ssl.HouseholdNet else cls() # HouseholdNet needs data: here, synthetic
    pars = dict(net.pars)
    for key, arg in inspect.signature(cls).parameters.items():
        if arg.default is not inspect.Parameter.empty and key not in ['pars', 'dhs_data'] and key not in pars:
            pars[key] = arg.default
    return ', '.join(f'{k}={v}' for k,v in pars.items())

rows = []
for name, cls in vars(ssl.networks).items():
    if inspect.isclass(cls):
        rows.append(dict(network=f'ssl.{name}', based_on=cls.__bases__[0].__name__, parameters=defaults(cls), random_numbers=crn_note(cls)))

pd.set_option('display.max_colwidth', None) # Show each entry in full
pd.DataFrame(rows).set_index('network')

	based_on	parameters	random_numbers
network
ssl.HouseholdNet	Network	prob_move_out=ss.bernoulli(p=0.7), dynamic=None	not CRN safe unless `dynamic=False` (the partners who move out are drawn together)
ssl.ClusterNet	Network	cluster_size=ss.poisson(mean=4), who=None, dynamic=False	CRN safe, unless `dynamic=True` and people are born or die
ssl.HybridNet	Network	household_size=ss.poisson(mean=2.0), n_contacts={'s': ss.poisson(mean=20), 'w': ss.poisson(mean=16), 'c': ss.poisson(mean=20)}, school_ages=[6, 22], work_ages=[22, 65], dynamic=False, method=fast, clustered=True	CRN safe, unless `dynamic=True` and people are born or die
ssl.NullNet	Network	n_people=None	CRN safe
ssl.StaticNet	Network	seed=True, p=None, n_contacts=10, graph=None	CRN safe
ssl.SmallWorldNet	StaticNet	n_contacts=4, rewire=0.1	CRN safe
ssl.ScaleFreeNet	StaticNet	n_contacts=4	CRN safe
ssl.DiskNet	Network	r=0.1, v=0.05	CRN safe
ssl.GeoNet	RandomNet	memory=None, n_contacts=ss.constant(value=10), dur=0, who=None, dynamic=True, method=fast, remove_self_loops=False, places=Kenya, level=1, local=0.9, c=2.0	Like `ss.RandomNet`, this is not CRN safe when people are born or die
ssl.GridNet	StaticNet	diagonal=False, periodic=False	CRN safe