Modeling a SIR model on a LxL lattice

Hey guys,

I need to create a SIR model on an LxL (100x100) Lattice, where it starts out with 10 infected sites on the lattice, and then spreads to the neighboring sites with the β probability. I just don't know how to make the actual lattice. I can model the SIR model normally, but I am unsure of how to transfer this system to a lattice. Any help on creating the actual lattice would be appreciated.

%SIR_Model

%Models the spread of a disease

%Variables S, I and R

%Set parameters

I0 = 0.01; %Initial infected

beta = 0.8; %Infection coefficient in weeks^-1

gamma = 0.5; %Removal coefficient in weeks^-1

tmax = 52; %number of weeks

dt = 0.01; %Size of time steps in weeks

Imax = 1.1; %Max number of infected for graph

%initial vectors

t = 0:dt:tmax; %Time vector

Nt = length(t); %Number of time steps

I = zeros(1,Nt); %Infection vector

S = zeros(1,Nt); %susceptiple vector

R = zeros(1,Nt); %Removed vector

I(1) = I0; %Set initial infection rate

%Calculations

for it = 1:Nt-1

S(it) = 1-I(it)-R(it);

dI = beta*I(it) * S(it)-gamma*I(it); %rate of change pr. week

I(it+1) = I(it) + dI * dt; %Adding how many new infections we have to the ones we already have

dR = gamma*I(it); %How many are added into the removed category pr. unit of time

R(it+1) = R(it)+dR*dt; %Refreshing the Removed

end

S(Nt) = 1-I(Nt)-R(Nt);

Answer 1

[Aashray]

Hi @Raazi,

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:

1. I placed all the individuals on a 100×100 grid.
2. Each person (cell) is either: Susceptible (S) = 0, Infected (I) = 1, Recovered (R) = 2.
3. Initially, I choose 10 random cells (individuals) to infect.
4. Then, at every time step, each infected cell attempts to infect its 4 neighbours (up, down, left, right) with a probability β.
5. 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