Demo for SI model of disease on network with rewiring
Consider a network of N individuals, initially randomly connected with average node degree avgdegree. The paper [Gross et al, 2006] studies a disease spreading on the network governed by the following rules: Each node is either infected (I) or susceptible (S). At each time step nodes change their state and their connections according to the following rules:
- an infected node recovers, becoming susceptible with probability
(SfromI in the script)
- every link between between an infected and a susceptible node transmits the disease with probability
(IfromS in the script), turning the susceptible node into an infected node. Links transmit the disease independent of each other.
- every link between an infected and a susceptible node gets rewired with probability
(rewire in the script). If rewiring occurs for a link connecting a node S and a node I, then the link swaps its end point I, replacing it with a randomly chosen susceptible node (excluding self-connections and double-connections). This models the effect that individuals replace a connection with another infected individual by connecting with another healthy individual with a certain probability.
Gross et al (2006) also formulated a system of three differential equations for the macroscopic quantities i (fraction of infected), ss (density of S-S links) and ii (density of I-I links). They found bistability between an endemic state (a large fraction of nodes is infected) and the disease-free state (the diesase dies out). Note that the disease-free state is absorbing such that, strictly speaking, the endemic state is always only a long-lived transient, even if for the macroscopic model the disease-free state is unstable.
Contents
Provided functions
- network= GenerateNetwork(N,avgdegree) generates a random network with N nodes returning an average degree avgdegree.
- [network,nodes,ode]=SIstep(network,nodes,rates) performs a single time step, changing nodes and links folllowing the above rules with rates given in the structure rates. It updates network and nodes. ode is a structure containing the macroscopic quantities ode.I (fraction of infected), ode.lii (density of links from I to I) and ode.lss (density of links ofrom S to S).
- AssertNetwork(network) checks if the array network is consistent (see below): if node A has a link to B, B must have a link to A, no double links, and no self-links.
Data structures
The network is stored in the form of two arrays.
- network: an array of size N x maxdegree, lists in row i the indices of all nodes, which node i has a link to. If the node has fewer than maxdegree links, the remaining entries are not in (1:N). The column dimension is larger (not necessarily equal) to the maximum degree. (The function SIstep never reduces the column dimension of network.)
- nodes: a vector of length N, contains the states (S or I) of all nodes.
The demo below performs a parameter sweep from down to
, showing an abrupt transition from the endemic state to the disease-free state.
Definition of basic quantities of for the network
The rates for node- and link-changing probabilities are kept in the structure rates. Our initial condition is a random netowrk where a randomly chosen 80% of the nodes are infected.
clear addpath('../utilities/'); % some auxiliary functions used in the script % globaldefinitions; % S=0, I=1 N=1000; % number of nodes avgdegree=10; % average degree rates=struct('SfromI',0.002,'IfromS',0.004,'rewire',0.04); network=GenerateNetwork(N,avgdegree); nodes=zeros(N,1)+S; nI_ini=round(N*0.8); nodes(randperm(N,nI_ini))=I;
Setup of range for simulations
p=0.01:-0.001:0; % parameter range to be scanned nsteps=1000; % number of time steps to be simulated for each parameter fracI=NaN(nsteps,length(p)); % pre-allocated array for fraction of infected
Parameter sweep
clf for k=1:length(p) rates.IfromS=p(k); for i=1:nsteps [network,nodes,ode]=SIstep(network,nodes,rates,'assert',true); fracI(i,k)=sum(nodes==I)/N; end fprintf('k=%d of %d,p=%g, I/N=%g\n',k,length(p),p(k),mean(fracI(:,k))); plot(reshape(fracI(:,1:k),1,k*nsteps)); xlabel('steps'); ylabel('fraction of infected nodes'); title(sprintf(['network of N=%d nodes with step-wise ',... 'decreasing p from p=%g to p=%g'],N,p(1),p(end))); drawnow end
k=1 of 11,p=0.01, I/N=0.971416 k=2 of 11,p=0.009, I/N=0.971945 k=3 of 11,p=0.008, I/N=0.966145 k=4 of 11,p=0.007, I/N=0.960685 k=5 of 11,p=0.006, I/N=0.954369 k=6 of 11,p=0.005, I/N=0.94582 k=7 of 11,p=0.004, I/N=0.92678 k=8 of 11,p=0.003, I/N=0.900278 k=9 of 11,p=0.002, I/N=0.748936 k=10 of 11,p=0.001, I/N=0.371379 k=11 of 11,p=0, I/N=0.064215

Time profile of another sweep with longer wait for the transients
(nsteps=4000).
Corresponding pseudo bifurcation diagram from the parameter sweep
(nsteps=400). Black dots correspond to the mean of the observed fraction of infected. The colour shading represents the histogram (darker= value visited more frequently).