I’m running a simple agent-based disease spread model on a 2D lattice, where disease transmission in a population depends on each individual’s chance of avoiding infection at every time step, which decreases with the local density of the infected and transmission rate.

import numpy as np
import random
from math import *
from scipy.signal import fftconvolve
import matplotlib.pyplot as plt


''' Landscape parameters ''' # slightly unusual setup
L = 10. 
nx = 100
dx = L/nx 
hs = 0.5 * dx  # half-step 
ulist = np.linspace(hs, L-hs, nx)

''' Epidemiological parameters ''' 
alpha = .1  # spatial scale used for convolution when approx. local densities
delta = .5  # transmission rate

''' Infection kernel '''  # kernel function is a bivariate Gaussian with zero covariance
i, j = ulist.reshape(nx,1), ulist.reshape(1,nx) 
ii = np.minimum(i-ulist[0], L-i+ulist[0]) 
jj = np.minimum(j-ulist[0], L-j+ulist[0])

# bivariate Gaussian: [functional form][1] 
norm_constant = 1/(2*pi*alpha**2) 
xmu = (ii-0)/alpha; ymu = (jj-0)/alpha 
rhs = np.exp(-.5 * (xmu**2 + ymu**2)) 
z = norm_constant * rhs  

ker = np.roll(z, np.int(nx//2-1), 0) 
ker = np.roll(ker,np.int(nx//2-1), 1)

''' Infection function '''
def infect2(ido_, idx_):  # ido: "healthy" individuals; idx: "infected" individuals
    og = np.zeros((nx,nx))
    np.add.at(og, (idx_[:,0], idx_[:,1]), 1)  # tally the local number of infected

    dd = fftconvolve(og, ker, mode='same')  # 2D convolution

    dx_lo = dd[ido_[:,0], ido_[:,1]] # density of infected local to the "healthy" individuals 
    transmit = 1-np.exp(-delta  * dx_lo)
    mask = np.random.random(len(ido_)) > transmit
    ido_new = ido_[mask]
    idx_new = np.vstack((idx_, ido_[~mask]))
    return ido_new, idx_new, dd

''' Population distribution '''    
ido = np.random.randint(nx, size=(19999,2)).astype(np.int) 
idx = np.array([[np.floor(nx/2), nx-np.ceil(.5*nx)]]).astype(np.int)

''' Simulation ''' 
for t in range(15):
    print 't=', t
    ido, idx, dd = infect2(ido, idx)

    plt.plot(ulist[ido[:,1]], ulist[ido[:,0]], 'ko', alpha=.1)
    plt.plot(ulist[idx[:,1]], ulist[idx[:,0]], 'ro', alpha=.3)
    print len(idx)
    plt.xlim([0., L]); plt.ylim([0., L])

    if len(ido) == 0: break

Increasing lattice resolution by setting nx=200 appeared to also change the resulting dynamics. I tried to fix this by spatially scaling the transmission rate setting delta = .5 /dx**2, I'm not sure about the reasoning behind this scaling method but it did make the results closer to those under nx=100. How can I further improve on the similarity between systems under different spatial resolutions? And what would be the mathematical/numerical justifications for those fixes?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.