15. Cellular automota#
Intro#
Cellular automata suppose that a set of hosts is arranged on a grid, called a lattice. For a cellular automata, we are asked to simulate the number of individuals in each of the possible supposed disease states over time. Most—if not all—cellular automata are stochastic. The two cellular automata models that we will look at (inspired again by Keeling and Rohani) are the contact process model introduced by Theodore Harris in 1974 and the “Forest Fire” model introduced by Bak et. al in 1990.
To be clear, cellular automata typically do not accurately represent real-life epidemic spread. Instead, we can extract from these models principles about how disease spreads, and attempt to find these (potentially unifying) principles in nature.
Contact process model#
Consider a single row of \(L\) squares, called a 1-dimensional lattice. Each square in the lattice represents on host who can be in one of two disease states: susceptible (S) and infected (I). An infector will transmit disease at the rate \(\beta\) infections per time unit, and the infectors will recover at a rate \(\gamma\) recoveries per time unit. Importantly, a susceptible host can only be infected if they are next to an infector—otherwise they cannot be infected. This makes clear that the (cellular) contact process model is a spatial model.
To begin this process, we select, at random, a single square to be infected.
How should we proceed for simulating this stochastic process?
Gillespie’s algorithm is a straightforward approach for simulating this process. We first need to define the events that can take place. In the previous models for which we used Gillespie’s algorithm, all the individuals were assumed to be mixed and in contact with everyone in the system. Because of this assumption we needed to define an event for infection and an event for recovery.
However, the (cellular) contact process model is spatial, and the events now depend on where they take place. At each time unit \(t\), we can collect the events that take place by looking at each host (each square in the model) and assessing whether that host can be infected or, if infected, then can recover.
Start with an empty list of events and begin iterating through hosts. If a host is next to an infector then they can be infected with rate \(\beta\), and we will append that event to the list If a host is infected already then they can recover with rate \(\gamma\), and we will append that event to the list.
An example list of events, \(E\), may look like \(E\) = [ (3,”recover”), (2,”infect”), (4,”infect”) ]. Here there are three events: host 3 is infected and can recover at rate \(\gamma\). Because host 3 is infected, host 2 is susceptible and can be infected at rate \(\beta\) and host 4 is susceptible and can be infected at rate \(\beta\).
This process can be implemented in Python
import numpy as np
import matplotlib.pyplot as plt
L = 50 #--L sites
beta = 2.0 #--transmission rate
gamma = 1.0 #--recovery rate
#--Pick how long to run the simulation.
T = 100
time = 0
sites = np.zeros(L) #--Long grid of zeros
#--pick at random one site to be infected.
#--Infections will be assigned the value one and susceptible the value zero.
pick_first_infection_site = np.random.choice(L)
sites[pick_first_infection_site] = 1.
#--We will record the number of infectors over time and time increments
I = [1]
times = [0]
while time < T:
# sum up rates for all possible events
events = []
event_rates = []
for site in range(L):
if sites[site]==0:
if sites[max(0,site-1)]==1 or sites[min(site+1,L-1)]==1:
events.append( (site,'infected') )
event_rates.append( beta )
elif sites[site] == 1:
events.append( (site,'recover') )
event_rates.append( gamma )
#--Extinction
if events == []:
break
total = sum(event_rates)
time_increment = np.random.exponential(1./total)
time+=time_increment
event_chosen = np.random.choice(len(event_rates), p = np.array(event_rates)/total )
site, event_type = events[event_chosen]
if event_type == "infected":
sites[site]=1
elif event_type=="recover":
sites[site]=0
times.append(time)
I.append(sum(sites))
fig, ax = plt.subplots()
ax.plot(times,I)
ax.set_xlabel("Time units")
ax.set_ylabel("Infectors")
fig.set_tight_layout(True)
plt.show()
Importance of spatial structure and visualization#
Though the number of infectors over time is often of interest, the (cellular) contact process model exhibits complex spatial structures for a simple—remember this model is 1-dimensional—model of an infection process.
To visualize infections in space and time we will record the number of infected across the \(L\) lattice sites for each time step \(t\). The end result will be a grid where the rows represent time and the columns represent the \(L\) hosts in the susceptible disease state (white squares) or the infected disease state (black squares).
def contact_process(L=50,beta=2,gamma=1, T = 100):
import numpy as np
time = 0
sites = np.zeros(L) #--Long grid of zeros
#--pick at random one site to be infected.
#--Infections will be assigned the value one and susceptible the value zero.
pick_first_infection_site = np.random.choice(L)
sites[pick_first_infection_site] = 1.
#--We will record the number of infectors over time and time increments
I = [1]
times = [0]
infectors = [ sites.copy() ]
while time < T:
# sum up rates for all possible events
events = []
event_rates = []
for site in range(L):
if sites[site]==0:
if sites[max(0,site-1)]==1 or sites[min(site+1,L-1)]==1:
events.append( (site,'infected') )
event_rates.append( beta )
elif sites[site] == 1:
events.append( (site,'recover') )
event_rates.append( gamma )
#--Extinction
if events == []:
break
total = sum(event_rates)
time_increment = np.random.exponential(1./total)
time+=time_increment
event_chosen = np.random.choice(len(event_rates), p = np.array(event_rates)/total )
site, event_type = events[event_chosen]
if event_type == "infected":
sites[site]=1
elif event_type=="recover":
sites[site]=0
times.append(time)
I.append(sum(sites))
infectors.append(sites.copy().astype(int))
return times,I,infectors
#--plot the spatial dynamics over time
plt.figure(figsize=(12, 12))
plt.imshow(infectors, cmap='gray_r', aspect='auto',interpolation='none')
plt.xticks([]) # remove ticks
plt.yticks([])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 53
51 #--plot the spatial dynamics over time
52 plt.figure(figsize=(12, 12))
---> 53 plt.imshow(infectors, cmap='gray_r', aspect='auto',interpolation='none')
54 plt.xticks([]) # remove ticks
55 plt.yticks([])
NameError: name 'infectors' is not defined
<Figure size 1200x1200 with 0 Axes>
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
def contact_process(L=50, beta=2, gamma=1, T=50):
sites = np.zeros(L, dtype=int)
sites[np.random.choice(L)] = 1
infectors = [sites.copy()]
time = 0
while time < T:
events, rates = [], []
for i in range(L):
if sites[i] == 1:
events.append((i, 'recover'))
rates.append(gamma)
else:
left = sites[i - 1] if i > 0 else 0
right = sites[i + 1] if i < L - 1 else 0
if left == 1 or right == 1:
events.append((i, 'infect'))
rates.append(beta)
if not events:
break
total = sum(rates)
dt = np.random.exponential(1 / total)
time += dt
chosen = np.random.choice(len(events), p=np.array(rates) / total)
i, action = events[chosen]
sites[i] = 1 if action == 'infect' else 0
infectors.append(sites.copy())
return infectors
infectors = contact_process(L=50, beta=2, gamma=1, T=30)
# Setup the figure and axis
fig, ax = plt.subplots(figsize=(8, 6))
data = np.zeros((1, len(infectors[0]))) # initialize empty frame
im = ax.imshow(data, cmap='gray_r', aspect='auto', vmin=0, vmax=1,interpolation='none')
ax.set_xlabel('Site')
ax.set_ylabel('Time step')
title = ax.set_title("Building Infection Timeline")
# Update function
def update(frame):
current = np.array(infectors[:frame + 1])
im.set_data(current)
im.set_extent([0, len(infectors[0]), frame + 0.5, -0.5]) # scrolls down
title.set_text(f"Time step {frame}")
return [im, title]
ani = animation.FuncAnimation(fig, update, frames=len(infectors), interval=100, blit=False)
# Render in Jupyter
HTML(ani.to_jshtml())
Forest-fire model: 2D Lattice simulation of an outbreak#
def from_empty_to_sus(hosts,p):
n_empty_slots = sum(sum(hosts==0))
sus = np.random.binomial(n_empty_slots, p)
empty_slots = [ (row,col) for row,persons in enumerate(hosts) for col,person in enumerate(persons) if person==0]
idxs = np.random.choice( len(empty_slots), size=sus )
for i in idxs:
row,col = empty_slots[i]
hosts[row,col]=1
return hosts
def spread_infection(hosts):
def determine_sus_in_contact(hosts):
sus_in_contact = []
infectors = [ (row,col) for row,persons in enumerate(hosts) for col,person in enumerate(persons) if person==2]
ROW,COL= hosts.shape
ROW-=1
COL-=1
for (row,col) in infectors:
if hosts[np.clip(row-1,0,ROW),np.clip(col-1,0,COL)] == 1:
sus_in_contact.append( ( np.clip(row-1,0,ROW) ,np.clip(col-1,0,COL)) )
if hosts[np.clip(row-0,0,ROW),np.clip(col-1,0,COL)] == 1:
sus_in_contact.append( (np.clip(row-0,0,ROW),np.clip(col-1,0,COL)) )
if hosts[np.clip(row+1,0,ROW),np.clip(col-1,0,COL)] == 1:
sus_in_contact.append( (np.clip(row+1,0,ROW),np.clip(col-1,0,COL)) )
if hosts[np.clip(row-1,0,ROW),np.clip(col,0,COL)] == 1:
sus_in_contact.append( (np.clip(row-1,0,ROW),np.clip(col,0,COL)) )
if hosts[np.clip(row+1,0,ROW),np.clip(col,0,COL)] == 1:
sus_in_contact.append( (np.clip(row+1,0,ROW),np.clip(col,0,COL)) )
if hosts[np.clip(row-1,0,ROW),np.clip(col+1,0,COL)] == 1:
sus_in_contact.append( (np.clip(row-1,0,ROW),np.clip(col+1,0,COL)) )
if hosts[np.clip(row-0,0,ROW),np.clip(col+1,0,COL)] == 1:
sus_in_contact.append( (np.clip(row-0,0,ROW),np.clip(col+1,0,COL)) )
if hosts[np.clip(row+1,0,ROW),np.clip(col+1,0,COL)] == 1:
sus_in_contact.append( (np.clip(row+1,0,ROW),np.clip(col+1,0,COL)) )
return sus_in_contact
sus_in_contact = determine_sus_in_contact(hosts)
while len(sus_in_contact)>0:
for (row,col) in sus_in_contact:
hosts[row,col]=2
sus_in_contact = determine_sus_in_contact(hosts)
return hosts
def from_infect_t_empty(hosts,current_infectors):
for (row,col) in current_infectors:
hosts[row,col]=0
return hosts
def spontaneous_infection(hosts,f):
sus_slots = [ (row,col) for row,persons in enumerate(hosts) for col,person in enumerate(persons) if person==1]
def neighbor(hosts, row,col):
ROW,COL= hosts.shape
ROW-=1
COL-=1
if hosts[np.clip(row-1,0,ROW),np.clip(col-1,0,COL)] == 2:
return True
if hosts[np.clip(row-0,0,ROW),np.clip(col-1,0,COL)] == 2:
return True
if hosts[np.clip(row+1,0,ROW),np.clip(col-1,0,COL)] == 2:
return True
if hosts[np.clip(row-1,0,ROW),np.clip(col,0,COL)] == 2:
return True
if hosts[np.clip(row+1,0,ROW),np.clip(col,0,COL)] == 2:
return True
if hosts[np.clip(row-1,0,ROW),np.clip(col+1,0,COL)] == 2:
return True
if hosts[np.clip(row-0,0,ROW),np.clip(col+1,0,COL)] == 2:
return True
if hosts[np.clip(row+1,0,ROW),np.clip(col+1,0,COL)] == 2:
return True
candidates = []
for (row,col) in sus_slots:
if neighbor( hosts, row,col ):
continue
candidates.append((row,col))
if np.random.random() < f:
idx = np.random.choice( len(candidates) )
row,col = candidates[idx]
hosts[row,col]=2
return hosts
#--initialize
L = 40
p = 0.01
f = 0.1
#--
# 0 means empty
# 1 means sus
# 2 means infected
hosts = np.zeros((L,L))
random_row = np.random.choice(L)
random_col = np.random.choice(L)
hosts[ random_row, random_col] = 2
infector_locs = [ (random_row,random_col) ]
T = 500
time = 0
times_record = [0]
infectors_record = [1]
frames = [hosts.copy()]
for t in range(T):
current_infectors = [ (row,col) for row,persons in enumerate(hosts) for col,person in enumerate(persons) if person==2]
hosts = spread_infection(hosts)
hosts = from_empty_to_sus(hosts,p)
hosts = from_infect_t_empty(hosts,current_infectors)
hosts = spontaneous_infection(hosts,f)
frames.append( hosts.copy() )
n_infectors = len([ (row,col) for row,persons in enumerate(hosts) for col,person in enumerate(persons) if person==2])
infectors_record.append(n_infectors)
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
from IPython.display import HTML
from matplotlib.colors import ListedColormap
cmap = ListedColormap(['white', 'red', 'blue'])
fig, ax = plt.subplots()
img = ax.imshow(frames[0], cmap=cmap, interpolation='none')
ax.axis('off')
plt.title("Infection Spread")
def update(frame_idx):
img.set_data(frames[frame_idx])
return [img]
ani = FuncAnimation(fig, update, frames=len(frames), interval=200, blit=True)
HTML(ani.to_jshtml())
Homework#
Code-up a stochastic SIS model. This model assumes that everyone comes in contact with everyone else at each time step. When an individual recovers they become susceptible again immediately. We assume that infections occur at a rate of \(\beta S I /N\). Run this model once for \(N=100\) individuals, transmission rate \(\beta=1.5\), recovery rate \(\gamma = 1\), and for 250 time steps. Begin this stochastic process with \(I_{0}=1\) infector in the system.
Plot the number of infectors over time (to ensure you are coding this up correctly).We will compare the frequency and size of outbreaks between the SIS model above (in 1) and the (cellular) contact process model. Run the SIS model 500 times and record the number of infectors at the final time step. Complete the same process for the (cellular) contact process model. Produce a complementary cumulative distribution function for the number of final infectors. Why are these curves different? What are the potential contributing factors towards this observed difference?
Read the paper On the Critical Behaviour of Simple Epidemics on course site
What is the relationship between the size-distribution exponent (what they call \(\tau\)) and infectious disease dynamics? Higher \(\tau\) implies what about transmission?
The authors discuss the role of dynamics versus random (or stochastic effects) that govern the number of observed cases of an infectious disease. Do main results from this paper support the idea that outbreaks could occur mostly from stochastic effects or could occur from dynamics?
The authors attempt to reproduce their observed data with a SEIR model but cannot. They claim that homogeneous mixing is to blame for the failed attempt to reproduce the data. Why might homogeneous mixing prevent the high frequency of small outbreaks interspersed with unlikely, but massive outbreaks that we observe in the forest fire model?
For the forest-fire model above, we will conduct a data analysis to measure the power law behavior from this simulation.
Please plot the number of infectors over time.
Please produce a complementary cumulative distributions plot of the number of infections over time.
Create a function that inputs a list (or array) of numbers
For every value \(x\) in the array compute the proportion of observing the value \(x\) or greater and call this px
Return a sorted list of \(x\) values and associated proportions (px)
Make two plots
plot x vs px
plot log(x+1) vs log(px)
Overlay a linear regression line on the (log(x+1),log(px)) “log-log” plot
Report the estimated slope and present the value \(\tau\) as was done in the paper above