Simulating state transitions
So far, we've covered the underlying mathematical principles that allow us to convert between continuous and discrete time, deterministic and stochastic dynamics, and how to interpret rates and probabilities in disease models. Now, let's see how we can simulate these processes in practice.
Linear approximation of transition probabilities
To perfectly simulate a linear time-invariant system, characterized by a transition rate \(\lambda\), we now know that each at-risk individual should have a transition probability over a step of length \(dt\) of
Calculating this probability requires an exponential function, which is computationally expensive. However, it is sometimes appropriate to use a more computationally efficient linear approximation.
To understand this, we can use the Taylor series expansion of the exponential function:
Thus, if \(\lambda dt\) is sufficiently small, so that (\(\lambda dt)^2\) and higher order terms are negligible, we can approximate the transition probability as:
In the next section, we will see how to choose the time step \(dt\) to ensure that this approximation is valid and that the transition probability never exceeds 1.0.
Choosing the time step
Above, we see that the transition probability \(p\) depends on the product of the rate \(\lambda\) and the time step \(dt\). This means that when we are implementing a discrete time model, we need to choose a time step that is sensitive to the rates.
If we select a time step that is too large, the transition probability \(p\) may be close to 1.0 (or exceed 1.0 if inappropriately using the linear approximation). Choosing a time step that is too large will lead to inaccurate results.
However, choosing a time step that is too small will lead to unnecessarily long computation times. So how do we choose an appropriate time step?
To avoid inaccuracies, start by identifying the largest (fastest) rate in your model, call it \(\lambda_{max}\). This rate can be identified mathematically for simple systems, but in practice you can think through the model piece by piece to identify the part that changes the fastest.
Intuitively, we will choose a time step that is inversely proportional to this maximum rate because the rate is multiplied by the time step to determine the transition probability. Thus we want to choose:
where \(\epsilon \ll 1\) is a small number that ensures that the transition probability remains sufficiently small. A reasonable rule of thumb that balances accuracy and computation time is to choose \(\epsilon = 0.1\). In other words, the time step should be an order of magnitude slower than the fastest rate in the system.
Notice that we're choosing the time step based on the fastest rate in the system. A system is called stiff if it has a wide range of rates. In this case, the fastest rate will determine the time step, which may be unnecessarily small for the slower rates. Stiff systems are challenging to simulate accurately and efficiently. If your system is stiff, consider updating different parts of the system at different frequencies or even using an event-based simulation approach.
Finally, there's no reason that the time step has to be constant. If properly adjusting rates into transition probabilities as a function of the time step, \(dt\), you can change the step duration throughout the simulation as needed to balance accuracy and computation time.