Problem 1: Exploring the SIR mean-field model

In this problem, we will examine how the network structure affects how well the Erdos-Renyi SIR mean field model captures the dynamics. Feel free to adapt the code we used in class.

Figure 1. Example small-world network generated using the Watts-Strogatz algorithm (image from Wikipedia).

For this problem, we will use both the Erdos-Renyi graph, as well as the Watts-Strogatz small world network. As a reminder, the way the Watts-Strogatz method works is that we start with a ring of nodes, connected to their \(k\) nearest neighbors. We then rewire one end of each edge to a new, randomly chosen neighbor with probability \(\beta\). This results in a network that has mostly local connections, but some long range connections (depending on \(\beta\)), leading to the small world phenomenon of high clustering but low path lengths between nodes.

Another useful feature of this model is that as you adjust \(\beta\) between \(0\) and \(1\), you can move between a very localized ring network when \(\beta = 0\) (matching what we used in cellular automata), to a more random network when \(\beta = 1\) (similar to an Erdos-Renyi network). Note that while the random network we achieve when \(\beta = 1\) (i.e. we rewire all edges) is very similar to an Erdos-Renyi network, it is not quite the same. This is because we only rewire one end of each edge (e.g. we only rewire edges going to clockwise nodes). This rewiring rule means that every node must have at least \(k/2\) edges (unlike in an Erdos-Renyi network, where the degree can potentially range from zero to all other nodes).

Problem 1A) Finding an ‘equivalent’ small-world network. Let’s see if we can find a sort of ‘equivalent’ small world network to match the Erdos-Renyi network we looked at in class. We can do this by setting the average degree of the small world network, \(k\), equal to that of the Erdos-Renyi network. Calculate the average degree of an Erdos-Renyi network with \(N\) nodes and edge probability \(p_e\). Be sure to explain your reasoning as well as giving a formula for the average degree.

Problem 1B) Implementing the simulation on a small-world network. Now that we know what average degree to use, adapt the SIR model on an Erdos-Renyi network that we used in class to run on a Watts-Strogatz network. You can do this as follows:

Implement the model and run it a couple of times for a few different \(\beta\) values ranging from 0 to 1 (e.g. at least include \(\beta\) values that are low, medium, and high, and you may want more values of \(\beta\) to fully see how it works!). How well does the Erdos-Renyi mean field model perform as you change \(\beta\) from 0 to 1? In your write-up, be sure to include a couple of plots of the prevalence (fraction of the population that is infected) over time, comparing the Erdos-Renyi mean field prediction to the Watts-Strogatz network simulation.

Problem 1C) Exploring the Erdos-Renyi mean-field model performance. Running one or two simulations is often not enough to really explore the model behavior, particularly for a stochastic model like this one. Adjust your code so that for a given parameter set, you run at least \(100\) simulations for \(100\) timesteps each. For this you will probably want to skip the visualization at each time step (and the pycxsimulator function if you’re using it), and instead run your model in a for loop without observing it until the simulation has finished. Some example code for this might be something like:

timesteps = 100  # number of timesteps per simulation
numsims = 100    # number of simulations to run

# this array will store the time series of fraction infected for each simulation
prevarray = np.zeros([timesteps,numsims])  

for i in range(0,numsims):        # loop over all simulations
    initialize()                  # initialize each simulation
    for j in range(1,timesteps):  # loop over the timesteps for that simulation
        update()                  # update
    prevarray[:,i] = prev         # store the resulting simulation in prevarray

Note however, that if you do it this way you will also need to move the prev = [], prev_mf = [], and susc_mf = [] lines into the initialize() function! Otherwise the code as written will keep appending new prevalences from one simulation to the next, so that when we try to store each run into our prevarray, the prev variable will be the wrong length.

Try running your code for at least five different values of \(\beta\), and be sure to include 0 and 1 (i.e. no rewiring and 100% rewiring). For each choice of \(\beta\) make a plot of all of your simulations, and plot the mean-field model prediction on top of this as a line (or use a distinct color, so we can tell the difference between the mean field and the network simulations). Include your plots in your writeup. You should see something like this (plotted for \(\beta = 0.5\), with the black line as the mean field):