Discrete models

This week, we convert last week's boxes into discrete dynamical equations that we can do maths and compute. We start with the SIR model, and show how we can make it more heterogenous by introducing, for instance, different degrees of susceptibility, distinct removal status (e.g., long lasting sickness, dead, recovered), or infection awareness (e.g., asymptomatic).

In this video, LHD walkthroughs the SIS model and introduces the mean-field approximation:

We are now ready to convert the mathematics into code:

As an example, LHD imagines a cloning factory as idealized through the diagram below:

where is the rate at which people clone themselves, and and are the rates at which people arrive and leave the factory, respectively. Events within each of these three processes occur at a fixed rate, hence independently from each other. These are Poisson processes, perhaps the simplest type of random process.

Back to example above, every now and then (e.g., everyday if we assume ) the factory (weirdly enough) thus ask you to follow this procedure:

# Julia code
for t=1:Δt:T
  if rand() < a*Δt
    N_t += 1         # 'arrivals' (people getting into the building)
  end
  for n=1:N_t 
    if rand() < l*Δt
      N_t -= 1     # 'departures'  (people leaving it)
    else
      if rand() < c*Δt
        N_t += 1       # 'cloning' (people cloning themselves)
      end
    end
  end
end

A couple of observations here. What if a transition rate, say , is too high so to make the associated transition probability larger than 1? This occurs because the duration of the time step you are considering is too large. For instance, if there are, on average, three people entering the building in a day and stands for one day, then . To solve the problem, take for example , so that now and .

Also, order matters. We are testing departure first, assuming that people can replicate only if they have not decided to leave the building. One could think of alternative situations where, during the same time step (which, remember, being finite allows for multiple processes to take place within it), people first try to clone themselves and then consider to leave; or, where recently arrived people are not allow to replicate or leave immediately after arrival (which would mean moving the `arrivals' code block to the end). In all cases, two events involving the same individual can never be simultaneous (in the factory, e.g., you can't leave while arriving or cloning yourself while leaving). Time discreteness forces us to hardwire the order of the potential events.

In the clip below, LHD shows how both these issues – transition probabilities larger than one and time ordering of competing events – are solved by considering small enough (eventually infinitesimal) time steps.

...

Key properties of Poisson processes

P1: The rate of a sum of Poisson processes of rates , , , is the sum of the rates,

P2: Given competing Poisson processes, the probability of Poisson process of rate occuring first is

Case study: discrete SIR

We review this week's key concepts by analyzing the discrete-time SIR model. For each contact between a susceptible and an infected individual, an infection event occurs with probability . Susceptible individuals then recover (more generally, the compartment can stand for 'removed', for individuals there are effectively outside the dynamics, either because they gain definitive immunity or die) with probability . We can resume these assumptions using a flow diagram:

We already saw in the second clip how to derive the probability that a susceptible individual gets the infection via at least one infectious contact. It reads , namely, [probability of not getting infected by any of the contacts]. (Say 4% transmission probability and 5 infected contacts, you get % of chance of not being infected.)

We thus get the following system of three equations

In code (taking ):

function runMAthDiscreteSIR(steps, N, β, γ) {
    // Initialize
    let I = 1, R = 0
    let S = N - I - R;
    
    // Pre-allocated arrays (faster)
    let St = new Array(steps);
    let It = new Array(steps);
    let Rt = new Array(steps);
    
    // Observe
    for (let step = 0; step < steps; step++) {
      St[step] = S;
      It[step] = I;
      Rt[step] = R;

      // Update
      let delta_I = S*(1 - Math.pow(1 - β, I))
      let delta_R = γ*I  

      S -= delta_I;
      I += delta_I - delta_R;
      R += delta_R;

    }
    return [St, It, Rt];
}

...

The model above is deterministic as it only tracks the expected values of the number of individuals in each compartment. Such formulation is exact only for a system of infinite size, as fluctuation around the mean becomes negligible. The stochasticity of the contagion process thus never disappears for a system of finite size (aka, real-world systems). We can see this by simulating the process using random variables. In this case, we rely on the Binomial distribution, which helps us determine how many successful events occur – such as how many susceptible individuals become infected – based on the probability of infection derived from the expected rate. The code remains largely the same, but we adjust how we update key quantities.

function runCompDiscreteSIR(steps, N, β, γ) {
    
    [...]
    
    for (let step = 0; step < steps; step++) {
    
      [...]
      
      // We use our previously defined math as probability of infection in the Binomial distribution
      let p_inf = 1 - Math.pow(1 - β, I); 
      let new_I = Binomial(S, p_inf); 
      let new_R = Binomial(I, γ);
      
      // Update
      S -= new_I;
      I += new_I - new_R;
      R += new_R;

      [...]
}

Click the button below to run a simulation and compare with the prediction made by the mathematical model:

Because of stochasticity, there are variation in draws. Notice, in particular, how from time to time the stochastic process goes extinct before an exponential growth can take place. With only the derministic model in our hand, when should we expect the process to be likely to take off? Imposing in Equation (1a) or in Equation (1c) requires . Any state configuration is thus a fixed point of the dynamics. This is not really informative. The revelant question is instead whether or not the system experiences a rise in infections before relaxing to an inactive equilibrium. To get an answer, we need to estimate the stability of those fixed points.

Linear stability analysis and criticality

This is the art of nudging a little our system around a fixed point to determine the (local) stability of that point. If following the perturbation the system tends to go back to that point, then the latter is said to be stable. It is unstable if the systems tends to move away from it, or neutrally stable if the perturbation neither grows or shrinks.

For our case study, we want to know the conditions under which the generic fixed point is stable or not. Nudging the system then means introducing an infinitesimal amount in the infected compartment and see whether such amount will grow to larger values or shrink towards zero. To do this, we first compute Equation (1b) in the limit. From , we get

It follows that – the perturbation shrinks – if and only if . The infection-free state is stable in this case. Otherwise, if , the perturbation expands and an outbreak takes place. We can restate these conditions in a more meaningful way by defining the basic reproduction number . Given that is the average infection period, and that is the average rate at which an infected individual produces new infections, represents the average number of new infections (secondary cases) produced by an infected individual (zero/primary case) in an otherwise infection-free population. The infection-free state is thus stable for and unstable for . The condition marks a critical point: outbreaks are expected above it but not below it. This is shown in the plot below, which reports results from the simulated SIR process.

We considered a population of size in a quasi-susceptible initial state (), and took . Imposing yields as the critical value for the infection rate. This is the value where a phase transition occurs, separating the inactive phase (left) from the active one (right).

We can get a more detailed view of this transition if we look at the distribution of outbreak sizes. That is, we choose a value for (notice the slider below) and run lots of simulations to record how many outbreaks of each possible size occurred. As we cross the critical point , a bulk of large outbreaks appears on the right. These are indeed the large outbreaks we expect to see above the critical point. The adjective large here has a specific meaning: it indicates outbreaks whose size is proportional to the size of the system – double the system's size and those large outbreaks will double their size too. The small outbreaks stored in the left chunk of the distribution, instead, do not scale and so their size becomes negligible as you consider bigger and bigger populations. Because of stochastic fluctuations, such small, localized infection chains are there even far above the critical point, although they become less and less frequent.

You may have noticed that the plot above has both its axis in logarithmic scale. The reason is that it allows to get a first, visual check of whether the system is or not at a critical point. How? Setting you should be able to see (either making an effort of imagination or running many more simulations to remove some noise) that the tail of the distribution is fitted by a straight line. A straight line in a log-log plot is a power-law in a linear plot, i.e., a curve of the form . At the critical point, the tail of the cluster size distribution thus follows a power-law. This is a general signature of criticality that you may encounter in any kind of system undergoing a phase transition. Also observe that moving just a little away from the critical point, either below or above it, makes the tail of the (small outbreaks) distribution decaying exponentially.

Bonus content

Lhd introducing Taylor Series and L'Hospital rule

Obligatory reference to 3b1b on L'Hôpital's rule

And Taylor Series