How to code SIR system to see state after time of each individual?

Let we have N individiviaduals in a system, out of which a random number of infividual are infected (say N/10) and the rest are susceptible. each suscptible can become infected with probability p_i and each infected can recover with probablity q_i. How can we code it to see the state after time t.

e.g let initially S=[0 1 1 1 1 1 1 1 1 1], I=[1 0 0 0 0 0 0 0 0 0] and R=[ 0 0 0 0 0 0 0 0 0 0].

I try to use Gillespie algorithm but its not working.

Answer 1

[Balavignesh]

Hi Wajid,

It is my understanding that you would like to simulate the spread of an infection in a system of individuals. The Gillespie algorithm is a stochastic simulation algorithm often used for such systems, especially in biochemical kinetics, but it might be overcomplicated for a basic simulation unless you're specifically interested in simulating the exact timings of events based on reaction rates. Instead of that, you could use a simple model based on probabilities for each event (infection and recovery).

Here is a rough implementation of the idea discussed:

% Example usage:
N = 10; % Total number of individuals
t = 5; % Number of time steps to simulate
p_i = 0.05; % Probability of infection
q_i = 0.1; % Probability of recovery

[S, I, R] = simulateInfection(N, t, p_i, q_i);
disp('Susceptible:'), disp(S)

Susceptible:
     1     1     1     1     1     1     0     1     1     0

disp('Infected:'), disp(I)

Infected:
     0     0     0     0     0     0     0     0     0     0

disp('Recovered:'), disp(R)

Recovered:
     0     0     0     0     0     0     1     0     0     1

function [S, I, R] = simulateInfection(N, t, p_i, q_i)
    % Initialize the state vectors for Susceptible (S), Infected (I), and Recovered (R)
    S = ones(1, N); % Initially, all individuals are susceptible
    I = zeros(1, N); % Initially, no one is infected
    R = zeros(1, N); % Initially, no one is recovered
    
    % Infect a random N/10 individuals
    infectedIndices = randperm(N, N/10);
    S(infectedIndices) = 0;
    I(infectedIndices) = 1;
    
    % Simulation loop for t time steps
    for step = 1:t
        for i = 1:N
            if S(i) == 1 && rand() < p_i
                % Susceptible individual becomes infected
                S(i) = 0;
                I(i) = 1;
            elseif I(i) == 1 && rand() < q_i
                % Infected individual recovers
                I(i) = 0;
                R(i) = 1;
            end
        end
    end
end

This code defines a function 'simulateInfection' that simulates the spread of an infection and recovery over a given number of time steps. It initializes the population as susceptible, randomly infects 'N/10' individuals, and then simulates the process of infection and recovery based on the provided probabilities 'p_i' and 'q_i'.

Note that this simulation assumes each time step allows for each individual to potentially change state based on the probabilities. The actual dynamics of your system could be more complex, and you might need to adjust the model accordingly (e.g., considering different probabilities for different individuals, incorporating a delay between infection and recovery, or using a more sophisticated model for the spread of infection).

The following documentation links might come handy:

• 'randperm' function : https://www.mathworks.com/help/matlab/ref/randperm.html
• 'ones' function: https://www.mathworks.com/help/matlab/ref/ones.html

Hope that helps!

Balavignesh