Numerical Integration of the Master Equation in Some Models of Stochastic Epidemiology

Garrett Jenkinson,John Goutsias,Stephen J Cornell

doi:10.1371/journal.pone.0036160

Abstract

The processes by which disease spreads in a population of individuals are inherently stochastic. The master equation has proven to be a useful tool for modeling such processes. Unfortunately, solving the master equation analytically is possible only in limited cases (e.g., when the model is linear), and thus numerical procedures or approximation methods must be employed. Available approximation methods, such as the system size expansion method of van Kampen, may fail to provide reliable solutions, whereas current numerical approaches can induce appreciable computational cost. In this paper, we propose a new numerical technique for solving the master equation. Our method is based on a more informative stochastic process than the population process commonly used in the literature. By exploiting the structure of the master equation governing this process, we develop a novel technique for calculating the exact solution of the master equation – up to a desired precision – in certain models of stochastic epidemiology. We demonstrate the potential of our method by solving the master equation associated with the stochastic SIR epidemic model. MATLAB software that implements the methods discussed in this paper is freely available as Supporting Information S1.

Highlights

Stochasticity can play an important role when studying a disease that spreads through a population of individuals [1,2,3]
A common approach to modeling this problem is by means of a Markov process, whose probability distribution satisfies a deterministic differential equation known as the master equation
We propose a promising numerical method to address this problem and illustrate its potential using a simple example based on the stochastic SIR epidemic model

Summary

Introduction

Stochasticity can play an important role when studying a disease that spreads through a population of individuals [1,2,3]. Accurate evaluation of the probability distribution of a Markov process requires a prohibitively large number of Monte Carlo samples for most systems of interest. Monte Carlo sampling is mostly used to estimate statistical summaries of the underlying stochastic population dynamics, such as means and variances. To evaluate the solution of the master equation, a number of approximation techniques have been proposed in the literature, such as the system-size expansion method of van Kampen [1,5,6]. The system-size expansion method for example can only produce a normal approximation to the solution of the master equation. If the probability distribution of the population process is bimodal, this method will produce erroneous results

Methods

Results

Conclusion