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:

- In your
`initialize()`

function, change the network from the`nx.erdos_renyi_graph()`

function, to use the`nx.watts_strogatz_graph()`

function. This function takes three arguments: \(N\), the number of nodes, \(k\), the average degree (i.e.Â number of ring neighbors to connect to), and \(\beta\), the rewiring probability. - For the number of neighbors to connect to, calculate the average
degree using the formula you derived in 1A, using an edge probability of
\(p_e = 0.06\).
**Note**that you will need to round your answer to the nearest integer, since you can only connect to an integer number of neighbors. - For the other parameters, take \(p_i = p_r = 0.1\) and \(N = 100\). We will try out a range of different \(\beta\) values, so set it to whatever you like for now.

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):