Skip to content

Tutorial

Jump to bottom
rincedd edited this page Sep 22, 2012 · 12 revisions

Building an SIS-model simulation

The model

In this short tutorial, you will create a simulation of a simple epidemiological model, the susceptible-infected-susceptible model of disease propagation on a random network. The model consists of a set of agents, represented by the nodes in a network, which can be in either of two states: susceptible to a communicable disease, S, or infected by the disease, I. The links in the network represent contacts among the agents. The disease can be transmitted through these contacts, so that for each link connecting an infected node to a susceptible node, the latter becomes infected with a fixed probability p per unit time. In addition to that, infected nodes can recover with a probability r per unit time.

The model dynamics is thus prescribed by two rules:

  • infection: S-I -> I-I at rate p
  • recovery: I -> S at rate r

In the following, these two rules will be implemented in a stochastic simulation on a random network.

Preliminaries

We will write the SIS model simulation in C++ using the largenet2 library. To keep things simple, we will put everything in a single file sis.cpp. This is of course not recommended for larger projects.

Create a file sis.cpp with the following content:

#include <iostream>
#include <largenet.h>

using namespace std;
using namespace largenet;

int main()
{
  return 0;
}

So far, this does nothing. It should compile with the following command on the command line (see also using largenet2 as shared library):

g++ -o sis sis.cpp -I <path to largenet2 headers> -L <path to liblargenet2> -llargenet2-0.3

where you have to supply the paths where the largenet2 include files (headers) and the library itself have been installed (e.g. $HOME/local/include and $HOME/local/lib).

Setting up the network

First, we will set up a random network consisting of susceptible and infected nodes. Node and edge states are represented by integers in the largenet library. For our convenience, we define symbolic names for the states:

enum NodeState { S, I };
enum EdgeState { SS, SI, II };

Let's define some constants for our simulation,

const node_size_t numberOfNodes = 1000;
const double averageDegree = 4.0;
const edge_size_t numberOfEdges = averageDegree * numberOfNodes / 2.0;

const double infectionRate = 0.03;
const double recoveryRate = 0.02;

const double initiallyInfectedFraction = 0.1;
const double maximumSimulationTime = 1000;

Now, we create a graph object with two possible node states and three possible edge states:

Graph graph(2, 3);

Graph listeners

To monitor changes in the graph we can use graph listeners. These are custom objects derived from the GraphListener class. Graph listeners can be attached to a graph object and will be notified of any event in the graph, such as a change in a node's or edge's state, the addition or removal of an edge and so on. In the following we will use a special predefined graph listener to ensure consistency between node and edge states in our model.

In our model, an edge's state is determined by the states of the nodes it connects, i.e. an edge connecting two S-nodes is in state SS and so on. In largenet, edge and node states are independent of each other. To ensure that the edge states are always consistent with the node states, we can use a StateConsistencyListener that automatically updates the states of adjacent edges whenever a node changes state:

class EdgeStateCalc
{
public:
  edge_state_t operator()(node_state_t source, node_state_t target) const
  {
    if (source == target)
    {
      if (source == S)
        return SS;
      else
        return II;
    }
    else
      return SI;
  }
};

StateConsistencyListener<EdgeStateCalc> consistencyListener(auto_ptr<EdgeStateCalc>(new EdgeStateCalc));
graph.addGraphListener(&consistencyListener);

The StateConsistencyListener uses the user-defined EdgeStateCalc object, which is a simple functor computing an edge's state from the states of the nodes it connects. (This object is passed to the consistencyListener through a shared auto_ptr for convenience.) Then, the consistencyListener is attached to the graph and will subsequently ensure the consistency of edge and node states for us.

Creating nodes and edges

The graph object does not yet contain any nodes or edges. We will create a random graph, i.e. an Erdős-Rényi network, using one of largenet's network generators. To that end, we need a random number generator. In principle, any random number generator can be used with largenet. Here, we will use the implementation of MyRNG that comes with the largenet examples in examples/lib. In fact, we recommend MyRNG as a powerful random number generator for most scientific purposes. We use the WELL1024a generator to create a random graph:

typedef myrng::RandomVariates<myrng::WELLEngine> RandomNumberGenerator;
RandomNumberGenerator randomNumberGenerator;

generators::randomGnm(graph, numberOfNodes, numberOfEdges, randomNumberGenerator, false);

Iterating over nodes

We now have a Graph object in which each node is in state 0, i.e. in the state S. We have to set the initial condition containing a number of infected nodes:

node_size_t numberOfInitiallyInfectedNodes = initiallyInfectedFraction * graph.numberOfNodes();
for (Graph::NodeIterator n = graph.nodes().first; n != graph.nodes().second; ++n)
  graph.setNodeState(n.id(), I);

Here, we used a Graph::NodeIterator object to iterate over all nodes in the network. Largenet provides a number of different such iterators for nodes and edges, similar to the STL iterators.

The simulation loop

We simulate the processes of infection and recovery using Gillespie's algorithm. To that end, we compute the instantaneous rates of both processes and draw the time after which the next event occurs from an exponential distribution. Then, infection occurs with the relative probability given by the ratio of the infection rate and the total rate of anything happening. Accordingly, recovery occurs with the complementary probability.

double t = 0;
printOutputLine(t, graph);
while (t < maximumSimulationTime)
{
  if (graph.numberOfNodes(I) == 0)
    break;
  double instantaneousInfectionRate = infectionRate * graph.numberOfEdges(SI);
  double instantaneousRecoveryRate = recoveryRate * graph.numberOfNodes(I);
  double totalRate = instantaneousInfectionRate + instantaneousRecoveryRate;
  if (totalRate == 0.0)
    break;
  t += randomNumberGenerator.Exponential(1.0 / totalRate);
  if (randomNumberGenerator.Chance(instantaneousInfectionRate / totalRate))
    infectOnRandomSIEdge(graph, randomNumberGenerator);
  else
    recoverRandomINode(graph, randomNumberGenerator);
  printOutputLine(t, graph);
}

Infection and recovery

In the case of an infection event, we select a random SI-edge from the graph and set the state of the S-node to I:

void infectOnRandomSIEdge(Graph& graph, RandomNumberGenerator& randomNumberGenerator)
{
  Graph::EdgeStateIterator edge = myrng::util::random_from(graph.edges(SI), randomNumberGenerator);
  node_id_t susceptibleNode = 0;
  if (graph.nodeState(edge->source()->id()) == S)
    susceptibleNode = edge->source()->id();
  else
    susceptibleNode = edge->target()->id();
  graph.setNodeState(susceptibleNode, I);
}

Here, we used the utility function random_from that comes with the MyRNG random number generator to select a random edge from the graph.

In the case of a recovery event, we select a random I-node from the graph and set its state to S:

void recoverRandomINode(Graph& graph, RandomNumberGenerator& randomNumberGenerator)
{
  Graph::NodeStateIterator node = myrng::util::random_from(graph.nodes(I), randomNumberGenerator);
  graph.setNodeState(node->id(), S);
}

Again, we used the random_from utility function.

Output

We simply write the numbers of infected and susceptible nodes as well as the numbers of edges in each state to the standard output stream:

void printOutputLine(double t, const Graph& graph)
{
  cout << t << "\t" << graph.numberOfNodes(S) << "\t"
      << graph.numberOfNodes(I) << "\t" << graph.numberOfEdges(SS) << "\t"
      << graph.numberOfEdges(SI) << "\t" << graph.numberOfEdges(II) << "\n";
}

To make things work, we have to include the following headers at the top of sis.cpp in addition to the ones already there:

#include <largenet2/StateConsistencyListener.h>
#include <largenet2/generators/generators.h>
#include <lib/RandomVariates.h>
#include <lib/WELLEngine.h>
#include <lib/util.h>

Additionally, we have to tell the compiler where to find the MyRNG headers, so we compile sis.cpp using

g++ -o sis sis.cpp -I <path to largenet2 headers> -I <path to largenet2>/examples -L <path to liblargenet2> -llargenet2-0.3

The full example sis.cpp can be found as examples/simple-sis/simple-sis.cpp in the examples directory. In the examples directory, you can find two more complete examples dealing with adaptive networks in which also the network structure changes and edges are rewired dynamically.

Clone this wiki locally