Computing the extinction path for epidemic models

In infectious disease modelling, the expected time from endemicity to extinction (of infection) may be analysed via WKB approximation, a method with origins in mathematical physics. The method is very general, but its uptake to date may have been limited by the practical difficulties of implementation. It is necessary to compute a trajectory of a (high dimensional) dynamical system, the ‘extinction path’, and this trajectory is maximally sensitive to small perturbations, making numerical computation challenging. The purpose of this paper is to make this methodology more accessible. Our method to achieve this is to present four computational algorithms, with associated Matlab code, together with discussion of various ways in which the algorithms may be tuned to achieve satisfactory convergence. One of the four algorithms is standard in this context, although we are able to somewhat enhance previously available code; the use of the three other algorithms in this context is novel. We illustrate our methods using three standard infectious disease models. Our results demonstrate that for each such model, our algorithms are able to improve upon previously available results.