I like the fact that you are moving from traditional SIR model for a more spatially realistic simulation of disease spread. In the standard SIR model (which you have shared), the population is assumed to be well-mixed but in real life scenarios, people only infect their neighbours, and that is exactly what a lattice-based SIR model helps capture.
How I approached simulating the system is as follows:
- I placed all the individuals on a 100×100 grid.
- Each person (cell) is either: Susceptible (S) = 0, Infected (I) = 1, Recovered (R) = 2.
- Initially, I choose 10 random cells (individuals) to infect.
- Then, at every time step, each infected cell attempts to infect its 4 neighbours (up, down, left, right) with a probability β.
- Each infected cell has a probability γ of recovering.
I have also attached the code for simulation. I hope it helps you out!
% Defining all parameters initially
L = 100; % Grid size (L x L)
beta = 0.8; % Infection probability
gamma = 0.5; % Recovery probability
nInfectedStart = 10; % Initial number of infected sites
nSteps = 100; % Number of time steps to simulate
% numerating all the possible states
SUSCEPTIBLE = 0;
INFECTED = 1;
RECOVERED = 2;
% Initialize lattice
grid = zeros(L); % Everyone is initially susceptible (0)
% Infecting 10 random unique cells
idx = randperm(L*L, nInfectedStart);
grid(idx) = INFECTED;
% Displaying initial state
imagesc(grid);
colormap([1 1 1; 1 0 0; 0.52 0.74 1]); % white = S, red = I, blue = R
title('Step 0');
axis equal tight;
pause(0.5);
% Main Simulation loop
for step = 1:nSteps
newGrid = grid; % Copy of the current state
for i = 1:L
for j = 1:L
if grid(i,j) == INFECTED
% Try infecting 4 neighbors: up, down, left, right
% considering infect probability (beta).
neighbors = [i-1 j; i+1 j; i j-1; i j+1];
for k = 1:size(neighbors,1)
x = neighbors(k,1);
y = neighbors(k,2);
% Check boundaries
if x >= 1 && x <= L && y >= 1 && y <= L
if grid(x,y) == SUSCEPTIBLE
if rand < beta
newGrid(x,y) = INFECTED;
end
end
end
end
% Trying to recover the infected using recover probability
% (gamma)
if rand < gamma
newGrid(i,j) = RECOVERED;
end
end
end
end
grid = newGrid;
% Visualizating the grid
imagesc(grid);
colormap([1 1 1; 1 0 0; 0.52 0.74 1]); % white=S, red=I, light blue=R
title(['Step ' num2str(step)]);
axis equal tight;
pause(0.1);
end
