Fitting stochastic epidemic models to gene genealogies using linear noise approximation

TL;DR

This paper introduces a Bayesian approach using linear noise approximation to efficiently fit stochastic epidemic models to gene genealogies, improving inference of disease dynamics from genetic data.

Contribution

It presents a novel Bayesian model that combines phylodynamic inference with stochastic epidemic models using LNA for computational efficiency.

Findings

Successfully applied to Ebola genealogies from 2014 epidemic

Provides joint posterior distributions of transmission parameters

Enhances computational tractability of stochastic epidemic inference

Abstract

Phylodynamics is a set of population genetics tools that aim at reconstructing demographic history of a population based on molecular sequences of individuals sampled from the population of interest. One important task in phylodynamics is to estimate changes in (effective) population size. When applied to infectious disease sequences such estimation of population size trajectories can provide information about changes in the number of infections. To model changes in the number of infected individuals, current phylodynamic methods use non-parametric approaches, parametric approaches, and stochastic modeling in conjunction with likelihood-free Bayesian methods. The first class of methods yields results that are hard-to-interpret epidemiologically. The second class of methods provides estimates of important epidemiological parameters, such as infection and removal/recovery rates, but ignores variation in the dynamics of infectious disease spread. The third class of methods is the most advantageous statistically, but relies on computationally intensive particle filtering techniques that limits its applications. We propose a Bayesian model that combines phylodynamic inference and stochastic epidemic models, and achieves computational tractability by using a linear noise approximation (LNA) --- a technique that allows us to approximate probability densities of stochastic epidemic model trajectories. LNA opens the door for using modern Markov chain Monte Carlo tools to approximate the joint posterior distribution of the disease transmission parameters and of high dimensional vectors describing unobserved changes in the stochastic epidemic model compartment sizes (e.g., numbers of infectious and susceptible individuals). We apply our estimation technique to Ebola genealogies estimated using viral genetic data from the 2014 epidemic in Sierra Leone and Liberia.

Tables

Table 1. Table A-1: Table for posterior median, 95% BCIs and ESSs for MCMC algorithms applied to Ebola data in Sierra Leone and Liberia.

	Sierra Leone				Liberia
		post med	95%BCI	ESS		post med	95%BCI	ESS
LNA	$I_0$	4.63	[1.28,14.41]	151	$I_0$	3.49	[1.03, 9.95]	1630
	$R_0$	1.69	[1.33.2.23]	167	$R_0$	1.67	[1.29,2.24]	942
	$\gamma$	32.47	[23.08,47.65]	345	$\gamma$	37.21	[25.98,53.13]	1704
	$\sigma$	14.71	[10.33,21.84]	141	$\sigma$	14.83	[10.41,20.70]	870
ODE	$I_0$	2.71	[1.10,6.37]	249	$I_0$	4.31	[1.89,9.27]	1236
	$R_0$	1.612	[1.30,2.09]	141	$R_0$	1.83	[1.41.2.44]	796
	$\gamma$	39.32	[26.63.55.82]	368	$\gamma$	38.31	[27.27.53.43]	1608
	$\sigma$	12.13	[8.61,16.98]	113	$\sigma$	13.67	[9.78,19.20]	879

Equations

$$\mathrm{Pr}(\mathbf{g})\propto \prod _{k=2}^n\lambda (c_{k-1})\exp \left(-{\int }_{c_{k-1}}^{c_k}\lambda (\tau )\mathrm{d}\tau \right),$$

$$\lambda (t)=\lambda (l(t),\beta (t),\mathbf{X}(t))=\left(\genfrac{}{}{0ex}{}{l(t)}{2}\right)\frac{2\beta (t)S(t)}{I(t)},$$

$$\mathrm{Pr}(\mathbf{g}\mid {\mathbf{X}}_{0:T},{\mathbf{\theta }}_{0:T})\propto \prod _{k=2}^n\left(\genfrac{}{}{0ex}{}{l(c_{k-1})}{2}\right)\frac{2\beta (c_{k-1})S(c_{k-1})}{I(c_{k-1})}\exp \left(-\sum _{i=0}^{i_k-1}{\int }_{{\mathcal{I}}_{i,k}}\left(\genfrac{}{}{0ex}{}{l_{i,k}}{2}\right)\frac{2\beta (\tau )S(\tau )}{I(\tau )}\mathrm{d}\tau \right).$$

$$\mathrm{Pr}\left({\mathbf{X}}_{0:T},{\mathbf{\theta }}_{0:T}\mid \mathbf{g}\right)\propto \mathrm{Pr}\left(\mathbf{g}\mid {\mathbf{X}}_{0:T},{\mathbf{\theta }}_{0:T}\right)\mathrm{Pr}\left({\mathbf{X}}_{1:T}\mid {\mathbf{X}}_0,{\mathbf{\theta }}_{0:T}\right)\mathrm{Pr}\left({\mathbf{\theta }}_{0:T}\right)\mathrm{Pr}\left({\mathbf{X}}_0\right),$$

$$\mathrm{Pr}({\mathbf{X}}_{1:T}\mid {\mathbf{X}}_0,{\mathbf{\theta }}_{0:T})=\prod _{i=1}^T\mathrm{Pr}({\mathbf{X}}_i\mid {\mathbf{X}}_{i-1},{\mathbf{\theta }}_{i-1}),$$

$$\mathbf{A}=\blockarray{ccl}\mbox{susceptible}&\mbox{infected}\\ \block{(cc)c}-11\mbox{infection}\\ 0-1\mbox{removal}.\\$$

$$\displaystyle\mathbf{h}(\mathbf{X}(t),\boldsymbol{\theta}(t))=\left(\begin{array}[]{c}\beta(t)S(t)I(t)\\ \gamma(t)I(t)\end{array}\right),\mathbf{H}=\mbox{diag}\left(\mathbf{h}(\mathbf{X}(t),\boldsymbol{\theta}(t))\right)=\left(\begin{array}[]{cc}\beta(t)S(t)I(t)&0\\ 0&\gamma(t)I(t)\end{array}\right).$$

$$\mathrm{d}\mathbf{\eta }(t)={\mathbf{A}}^T\mathbf{h}(\mathbf{\eta }(t),{\mathbf{\theta }}_{i-1})\mathrm{d}t,t\in [t_{i-1},t_i].$$

$$\mathrm{d}\mathbf{M}(t)=\mathbf{F}(\mathbf{\eta }(t),{\mathbf{\theta }}_{i-1})\mathbf{M}(t)\mathrm{d}t+\sqrt{{\mathbf{A}}^T\mathbf{H}(\mathbf{\eta }(t),{\mathbf{\theta }}_{i-1})\mathbf{A}}\mathrm{d}{\mathbf{W}}_t,t\in [t_{i-1},t_i],$$

$$\mathrm{d}\mathbf{m}(t)$$

$$\mathrm{d}\mathbf{\Phi }(t)$$

$${\mathbf{X}}_i\mid {\mathbf{X}}_{i-1},{\mathbf{\theta }}_{i-1}\sim \mathcal{N}\left(\mathbf{\eta }(t_i)+\mathbf{\mu }\left({\mathbf{X}}_{i-1}-\mathbf{\eta }(t_{i-1}),\Delta t_i,{\mathbf{\theta }}_{i-1}\right),\mathbf{\Phi }(t_i)\right).$$

$$R_0(t)=\sum _{i=1}^T{R}_{0_{i-1}}{1}_{[t_{i-1},t_i)}(t),$$

$$R_{0_i}=R_0(t,{\mathbf{\delta }}_{1:T},\sigma )=R_0\exp \left(\sum _{k=1}^i\sigma {\delta }_k\right),\text{mbox}fori=1,\dots ,T,$$

$${\mathbf{X}}_i=\mathbf{\eta }(t_i)+\mathbf{\mu }({\mathbf{X}}_{i-1}-\mathbf{\eta }(t_{i-1}),\Delta t_i,{\mathbf{\theta }}_{i-1})+{\mathbf{\Phi }}_i^{1/2}{\mathbf{\xi }}_i,$$

$$\mathrm{Pr}(I_0,R_0,\gamma ,{\mathbf{\delta }}_{1:T},{\mathbf{\xi }}_{1:T},\sigma |\mathbf{g})$$

$$\propto \mathrm{Pr}(\mathbf{g}|{\mathbf{X}}_{0:T},{\mathbf{\theta }}_{0:T})\mathrm{Pr}(I_0)\mathrm{Pr}(R_0)\mathrm{Pr}(\gamma )\mathrm{Pr}({\mathbf{\delta }}_{1:T})\mathrm{Pr}({\mathbf{\xi }}_{1:T})\mathrm{Pr}(\sigma ).$$

$$\text{mbox}MAE=\frac{1}{T+1}\sum _{i=0}^T|{\hat{R}}_{0_i}-R_0(t_i)|,$$

$$\text{mbox}MRAE=\frac{1}{T+1}\sum _{i=0}^T\frac{|{\hat{I}}_i-I(t_i)|}{I(t_i)+1},$$

$$\text{mbox}MCIW=\frac{1}{T+1}\sum _{i=0}^T\left[{\hat{R}}_{0_i}^{0.975}-{\hat{R}}_{0_i}^{0.025}\right],$$

$$\text{mbox}MRCIW=\frac{1}{T+1}\sum _{i=0}^T\frac{{\hat{I}}_i^{0.975}-{\hat{I}}_i^{0.025}}{I(t_i)+1},$$

$$\text{mbox}{ENV-R}_0=\frac{1}{T+1}\sum _{i=0}^T\text{mathbbm}1\left({\hat{R}}_{0_i}^{0.025}\le R_0(t_i)\le {\hat{R}}_{0_i}^{0.975}\right),\text{mbox}ENV-I=\frac{1}{T+1}\sum _{i=0}^T\text{mathbbm}1\left({\hat{I}}_i^{0.025}\le I(t_i)\le {\hat{I}}_i^{0.975}\right).$$

$$\hat{N}(t)=\hat{\beta }(t)\hat{S}(t)\hat{I}(t)\cdot \Delta t,$$

$${\tilde{a}}_{i1}{\mathcal{X}}_1+\cdots +{\tilde{a}}_{id}{\mathcal{X}}_d\stackrel{{\theta }_i}{\to }{\bar{a}}_{i1}{\mathcal{X}}_1+\dots +{\bar{a}}_{i1}{\mathcal{X}}_d,$$

$$\mathbf{A}:=\{{\tilde{a}}_{ij}-{\bar{a}}_{ij}\}\in {\mathbb{Z}}^{q\times d},$$

$$h_i={\theta }_i\prod _{j=1}^d\left(\genfrac{}{}{0ex}{}{X_j}{{\tilde{a}}_{ij}}\right).$$

$$\mathbf{h}(\mathbf{X},\mathbf{\theta })={\left(h_1,\dots ,h_q\right)}^T,\mathbf{H}(\mathbf{X},\mathbf{\theta })=\text{mbox}diag\left(\mathbf{h}(\mathbf{X},\mathbf{\theta })\right).$$

$$\mathrm{d}\mathbf{X}={\mathbf{A}}^T\mathbf{h}\left(\mathbf{X},\mathbf{\theta }\right)\mathrm{d}t,\mathbf{X}(0)={\mathbf{x}}_0,$$

$$\text{mbox}Susceptible+\text{mbox}Infected\stackrel{\beta }{\to }\text{mbox}Exposed+\text{mbox}Infected,$$

$$\text{mbox}Exposed\stackrel{\mu }{\to }\text{mbox}Infected,$$

$$\text{mbox}Infected\stackrel{\gamma }{\to }\text{mbox}Recovered.$$

$$\mathbf{A}=\blockarray{cccc}\mbox{Susceptible}&\mbox{Exposed}\mbox{Infected}\\ \block{(ccc)c}-110\mbox{reaction }\eqref{eq:reaction3.1}\\ 0-11\mbox{reaction }\eqref{eq:reaction3.2}\\ 00-1\mbox{reaction }\eqref{eq:reaction3.3}\\ ,$$

Full text

Key words: Coalescent, Susceptible-Infectious-Recovered model, state-space model, phylodynamics, Ebola virus

1 Introduction

Phylodynamics is an area at the intersection of phylogenetics and population genetics that studies how epidemiological, immunological, and evolutionary processes affect viral phylogenies constructed based on molecular sequences sampled from the population of interest [Grenfell et al., 2004, Volz et al., 2013]. Phylodynamics is especially useful in infectious disease modeling because genetic data provide a source of information that is complimentary to the traditional disease case count data. Here we are interested in inferring parameters governing infectious disease dynamics from the genealogy/phylogeny estimated from infectious disease agent molecular sequences collected during the disease outbreak. Working in a Bayesian framework, we develop an efficient Markov chain Monte Carlo (MCMC) algorithm that allows us to work with stochastic models of infectious disease dynamics, properly accounting for stochastic nature of the dynamics.

Currently, learning about population-level infectious disease dynamics from molecular sequences can be accomplished using three general strategies. The first strategy relies on the coalescent theory — a set of population genetics tools that specify probability models for genealogies relating individuals randomly sampled from the population of interest [Kingman, 1982, Griffiths and Tavare, 1994, Donnelly and Tavare, 1995]. Using a subset of these models [Griffiths and Tavaré, 1994], it is possible to estimate changes in effective population size — the number of breeding individuals in an idealized population that evolves according to a Wright-Fisher model [Wright, 1931]. Such reconstruction can be done assuming parametric [Kuhner et al., 1998, Drummond et al., 2002] or nonparametric [Drummond et al., 2002, 2005, Minin et al., 2008, Palacios and Minin, 2013, Gill et al., 2013] functional forms of the effective population size trajectory. In the context of infectious disease phylodynamics, nonparametric inference is the norm and the estimated effective population size is often interpreted as the effective number of infections or the effective number of infectious individuals. However, reconstructed effective population size trajectories are not easy to interpret and estimation of parameters of disease dynamics is difficult to accomplish if one wishes to maintain statistical rigor [Pybus et al., 2001, Frost and Volz, 2010].

Another way to learn about infectious disease dynamics from molecular sequences is to model explicitly events that occur during the infectious disease spread and to link these events to the genealogy/phylogeny of sampled individuals using birth-death processes. For example, a Susceptible-Infectious-Removed (SIR) model includes two possible events: infections and removals (e.g., recoveries and deaths), represented by births and deaths in the corresponding birth-death model [Stadler et al., 2013, Kühnert et al., 2014]. Other SIR-like models (e.g., SI and SIS models) differ by the number and types of the events that are needed to accurately describe natural history of the infectious disease [Leventhal et al., 2013]. Although these methods are more principled than post-hoc processing of nonparametrically estimated disease dynamics, they are not easy to scale to large datasets and/or high dimensional models. For example, in order to fit phylodynamic birth-death models to genomic and epidemiological data Vaughan et al. [2018] use particle filter MCMC. However, computational burden of particle filter MCMC methods is usually very high. Moreover, these methods often struggle with convergence when the dimensionality of statistical model parameters is even moderately high [Andrieu et al., 2010].

Structured coalescent models provide the third strategy of inferring parameters governing spread of an infectious disease [Volz et al., 2009, Volz, 2012, Dearlove and Wilson, 2013]. These models assume infectious disease agent genetic data have been obtained from a random sample of infected individuals, allowing for serial sampling over time. Although similar to the birth-death modeling framework, the structured coalescent models have two advantages. First, one does not have to keep track, analytically or computationally, of extinct and not sampled genetic lineages. Second, the density of the genealogy can be obtained given the population level information about status of individuals: for example, in the SIR model it is sufficient to know the numbers of susceptible, ($S(t)$), infectious, ($I(t)$), and recovered, ($R(t)$), individuals at each time point $t$. The second advantage comes with two caveats: 1) such densities can be obtained only approximately and 2) evaluating densities of genealogies is not straightforward and involves numerical solutions of differential equations. Even in cases when these caveats are manageable, the density of the assumed stochastic epidemic model population trajectory remains computationally intractable. One way around this intractability assumes a deterministic model of infectious disease dynamics [Volz et al., 2009, Volz, 2012, Volz and Pond, 2014], which potentially leads to overconfidence in estimation of model parameters. Particle filter MCMC offers another solution [Rasmussen et al., 2011, 2014], but, as we discussed already, these methods are difficult to use in practice, especially in high dimensional parameter spaces.

In this paper, we develop methods that allow us to bypass computationally unwieldy particle filter MCMC with the help of a linear noise approximation (LNA). LNA is a low order correction of the deterministic ordinary differential equation describing the asymptotic mean trajectories of compartmental models of population dynamics defined as Markov jump processes (e.g., chemical reaction models and SIR-like models of infectious disease dynamics) [Kurtz, 1970, 1971, Van Kampen and Reinhardt, 1983]. LNA can also be viewed as a first order Taylor approximation of Markov population dynamics models represented by stochastic differential equations [Giagos, 2010, Wallace, 2010]. A key feature of the LNA method is that it approximates the transition density of a stochastic population model with a Gaussian density [Komorowski et al., 2009].

Inspired by recent applications of LNA to analysis of Google Flu Trends data [Fearnhead et al., 2014] and disease case counts [Buckingham-Jeffery et al., 2018], we develop a Bayesian framework that combines LNA for stochastic models of infectious disease dynamics with structured coalescent models for genealogies of infectious disease agent genetic samples. Our approach yields a latent Gaussian Markov model that closely resembles a Gaussian state-space model. We use this resemblance to develop an efficient MCMC algorithm that combines high dimensional elliptical slice sampler updates [Murray et al., 2010] with low dimensional Metropolis-Hastings (MH) moves. Using simulations, we demonstrate that this algorithm can handle reasonably complex models, including an SIR model with a time-varying infection rate. We apply this SIR model to a recent Ebola outbreak in West Africa. Our analysis of data from Liberia and Sierra Leone illuminates significant changes in the Ebola infection rate over time, likely caused by the public health response measures and increased awareness of the outbreak in the population.

2 Methodology

2.1 Genealogy as data

We start with $n$ infectious disease agent molecular sequences obtained from infected individuals sampled uniformly at random from the total infected population. Further, we assume that a phylogenetic tree, or genealogy, $\mathbf{g}$ relating these sequences has been estimated in such a way that the tree branch lengths respect the known sequence sampling times. Such estimation can be performed with, for example, BEAST — a leading software package for Bayesian phylogenetic studies, particularly popular among molecular epidemiologists who collect and analyze viral genetic sequences [Suchard et al., 2018]. The genealogy is represented by a tree structure with its nodes containing two sources of temporal information: coalescent and sampling times. The coalescent times correspond to the internal nodes of the tree, which are defined as the times at which two lineages in the tree are merged into a common ancestor. The sampling times, corresponding to the tips of the tree, are the times at which molecular sequences were sampled. Note that sampling times are observed directly, while coalescent times are estimated from molecular sequences during phylogenetic reconstruction.

To perform inference about infectious disease dynamics using the above genealogy we need a probability model that relates the genealogy and infectious disease dynamics model parameters. Without too much loss of generality, we assume that the infectious disease is spreading through the population according to the SIR model — a compartmental model that at each time point $t$ tracks the number of susceptible individuals $S(t)$, number of infected/infectious individuals $I(t)$, and number of removed individuals $R(t)$ [Bailey, 1975, Anderson and May, 1992]. We assume that the population is closed so $S(t)+I(t)+R(t)=N$ for all times $t$, where $N$ is the population size that we assume to be known. This constraint implies that vector $\mathbf{X}(t)=(S(t),I(t))$ is sufficient to keep track of the population state at time $t$. We follow common practice and model $\mathbf{X}(t)$ as a Markov jump process (MJP) with allowable instantaneous jumps shown in Figure 1 [O’Neill and Roberts, 1999]. The assumed MJP process $\mathbf{X}(t)$ is inhomogeneous, because we allow the infection rate $\beta(t)$ and removal rate $\gamma(t)$ to be time-varying.

The structured coalescent models assume that only coalescent times $c_{1}<c_{2}<\cdots<c_{n-1}$ provide information about the population dynamics. These times are modeled as jumps of an inhomogeneous pure death process with rate $\lambda(t)$, where each “death” event corresponds to coalescence of two lineages and $\lambda(t)$ is called a coalescent rate. Then the density of the genealogy, which serves as a likelihood in our work, is written as

[TABLE]

where $c_{n}$ denotes the most recent sequence sampling time. The dependence of coalescent rate on the assumed population dynamics can be complicated and mathematically intractable, but luckily approximations exist for some specific cases. For the SIR model the approximate coalescent rate can be obtained via the following formula:

[TABLE]

where $l(t)$ is the number of lineages present at time $t$. Note that when the number of susceptibles is not changing significantly relative to the total population size (i.e., $S(t)\approx N$) and infection rate is constant (i.e., $\beta(t)=\beta$), the structured coalescent reduces to the classical Kingman’s coalescent, where we interpret $I(t)/(2\beta N)$ as the effective population size trajectory [Kingman, 1982]. It is possible to find approximate coalescence rate for general compartmental models, but closed form expressions exist only for a few models with a low number of compartments (e.g., SI, SIR) [Volz et al., 2009, Volz, 2012, Dearlove and Wilson, 2013].

Since we allow sequences to be sampled at different times $s_{1}<s_{2}<\cdots<s_{m}=c_{n}$, some inter-coalescent times are censored. To deal with this censoring algebraically, each inter-coalesecent interval $[c_{k-1},c_{k})$ is partitioned by the sampling events into $i_{k}$ sub-intervals: $\mathcal{I}_{0,k},\ldots,\mathcal{I}_{i_{k}-1,k}$. The intervals that end with a coalescent event are defined as $\mathcal{I}_{0,k}=[c_{k-1},\min\{c_{k},s_{j}\})$, for $s_{j}>c_{k-1}$ and $k=2,\ldots,n$. Let the number of lineages in each interval $\mathcal{I}_{i,k}$ be $l_{i,k}$. Then the number of lineages at each time point $t$ can be written as $l(t)=\sum_{k=2}^{n}\sum_{i=0}^{i_{k-1}}1_{\{t\in I_{i,k}\}}l_{i,k}$. If the interval $\mathcal{I}_{i,k}$ ends with a coalescent time, the number of lineages in the next interval will be decreased by $1$. If the interval ends with a sampling event $s_{i}$, then the number of lineages in the next interval is increased by $n_{i}$ — the number of sequences sampled at time $s_{i}$. Figure 2 shows an example of a genealogy with labeled coalescent times, sampling times, number of lineages, and the corresponding intervals.

We are now ready to connect the SIR model and a genealogy with serially sampled tips with the help of a structured coalescent density/likelihood. First we discretize the time interval between the time to most recent common ancestor $c_{1}$ (time corresponding to the root of the tree) and the most recent sampling time $s_{m}$ using a regular grid $t_{0}<t_{1}<\cdots<t_{T}$ ($t_{0}<c_{1}$ and $t_{T}>s_{m}$). Using this grid, we discretize the latent epidemic trajectory by assuming that $\mathbf{X}(t)=\sum_{j=1}^{T}\mathbf{X}_{j-1}\mathbf{1}_{[t_{j-1},t_{j})}(t)$, where $\mathbf{X}_{j}=(S_{j},I_{j})$ is a column vector. Similarly, we discretize the infectious disease dynamics parameter vector trajectory $\boldsymbol{\theta}(t)=(\beta(t),\gamma(t))$ so that $\boldsymbol{\theta}(t)=\sum_{j=1}^{T}\boldsymbol{\theta}_{j-1}\mathbf{1}_{[t_{j-1},t_{j})}(t)$, where $\boldsymbol{\theta}_{j}=(\beta_{j},\gamma_{j})$ is also a column vector. We collect latent variables $\mathbf{X}_{j}$s and parameters $\boldsymbol{\theta}_{j}$s into matrices $\mathbf{X}_{0:T}$ and $\boldsymbol{\theta}_{0:T}$ respectively. The SIR structured coalescent density/likelihood then becomes

[TABLE]

Since $S(t)$, $I(t)$, and $\beta(t)$ are piecewise constant functions, the integrals in the above formula are readily available in closed form and are fast to compute.

2.2 Bayesian data augmentation

2.2.1 Posterior distribution

Given genealogy $\mathbf{g}$, our goal is to infer the latent SIR population dynamic $\mathbf{X}_{0:T}$ and rate parameters $\boldsymbol{\theta}_{0:T}$ over time grid $t_{0}<t_{1}<\cdots<t_{T}$. Let $\Pr(\mathbf{X}_{0})$ and $\Pr(\boldsymbol{\theta}_{0:T})$ denote the prior densities for the initial compartment states and the SIR parameters respectively. The posterior distribution for the population trajectory $\mathbf{X}_{0:T}$ and parameters $\boldsymbol{\theta}_{0:T}$ given observed genealogy $\mathbf{g}$ is

[TABLE]

where $\Pr\left(\mathbf{g}\mid\mathbf{X}_{0:T},\boldsymbol{\theta}_{0:T}\right)$ is the structured coalescent likelihood introduced in Section 2.1 and $\Pr(\mathbf{X}_{1:T}\mid\mathbf{X}_{0},\boldsymbol{\theta}_{0:T})$ is the likelihood function for discrete observations of trajectory $\mathbf{X}_{1:T}$ given the initial value $\mathbf{X}_{0}$:

[TABLE]

where the factorization comes from the assumed Markov property of the disease dynamics. However, the SIR transition density $\Pr(\mathbf{X}_{i}\mid\mathbf{X}_{i-1},\boldsymbol{\theta}_{i-1})$ becomes intractable as population size $N$ grows large, making it difficult to perform likelihood-based inference for outbreaks in large populations.

2.2.2 Linear noise approximation

To furnish a feasible computation strategy for large populations, we use a linear noise approximation (LNA) method, in which the computationally intractable transition probability $\Pr\left(\mathbf{X}_{i}\mid\mathbf{X}_{i-1},\boldsymbol{\theta}_{i-1}\right)$ is approximated using a closed form Gaussian transition density. The LNA method replaces the MJP discrete state space with a continuous state space of $\mathbf{X}(t)$ to approximate the counts of at time $t$, under the following constraints: $S(t)>0,I(t)>0$ and $S(t)+I(t)\leq N$. To briefly explain how this approximation is obtained, we will need additional notation.

The SIR MJP instantaneous transitions, depicted in Figure 1, are encoded in an effect matrix

[TABLE]

Each row in matrix (5) represents a type of transition event and each column corresponds to a change in the susceptible and infected populations. Next, we define a rate vector $\mathbf{h}$ and a rate matrix $\mathbf{H}$:

[TABLE]

The above notation, as well as subsequent developments based on it, can be generalized to other epidemic models and, more generally, to a large class of density dependent stochastic processes, such as chemical reaction and gene regulation models [Wilkinson, 2011]. See Section A-1 in the Appendix for more details on this generalization.

Consider a transition from $\mathbf{X}_{i-1}$ at time $t_{i-1}$ to $\mathbf{X}_{i}$ at $t_{i}$. Recall that we assume that the SIR rates $\boldsymbol{\theta}(t)$ take constant values $\boldsymbol{\theta}_{i-1}$ in $[t_{i-1},t_{i})$. The LNA represents the value of the next state $\mathbf{X}_{i}$ as $\mathbf{X}_{i}=\boldsymbol{\eta}(t_{i})+\mathbf{M}(t_{i})$, where $\boldsymbol{\eta}(t_{i})$ is a deterministic component and $\mathbf{M}(t_{i})$ is a stochastic component. The deterministic component $\boldsymbol{\eta}(t_{i})$ can be obtained by solving the standard SIR ODE that in our notation can be written as

[TABLE]

The stochastic part $\mathbf{M}(t_{i})$ corresponds to the solution of the following SDE at time $t_{i}$:

[TABLE]

where \mathbf{F}(\boldsymbol{\eta}(t),\boldsymbol{\theta}_{i-1}):=\dfrac{\partial\mathbf{A}^{T}\mathbf{h}(\mathbf{X}(t),\boldsymbol{\theta}_{i-1})}{\partial\mathbf{X}}\Big{|}_{\mathbf{X}=\boldsymbol{\eta}(t)} is the Jacobian matrix of the deterministic part $\mathbf{A}^{T}\mathbf{h}(\mathbf{X}(t),\boldsymbol{\theta}_{i-1})$ in (11) evaluated at $\boldsymbol{\eta}(t)$. The solution of SDE (12), $\mathbf{M}(t)$, is a Gaussian process and can be recovered by solving two ordinary differential equations governing the mean function $\mathbf{m}(t):=\mathbf{E}[\mathbf{M}(t)]$ and covariance function $\boldsymbol{\Phi}(t):=\mathbf{Var}(\mathbf{M}(t))$:

[TABLE]

for $t\in[t_{i-1},t_{i}]$. A heuristic derivation of LNA, based on Wallace [2010], is given in Section A-2 of the Appendix. Let $\boldsymbol{\eta}_{t_{i-1}},\mathbf{m}_{t_{i-1}},\boldsymbol{\Phi}_{t_{i-1}}$ denote the initial values of $\boldsymbol{\eta}(t),\mathbf{m}(t),\boldsymbol{\Phi}(t)$ at time $t_{i-1}$ in differential equations (11), (13), and (14) respectively. There are two options for choosing these initial conditions: the non-restarting LNA of Komorowski et al. [2009] and the restarting LNA of Fearnhead et al. [2014]. In this paper, we will use the non-restarting LNA by Komorowski et al. [2009] with the following choice of initial conditions:

$\boldsymbol{\eta}_{t_{i-1}}=\boldsymbol{\eta}(t_{i-1})$, where $\boldsymbol{\eta}(t_{i-1})$ was obtained by solving the ODE (11) using parameter vector $\boldsymbol{\theta}_{i-2}$ over the interval $[t_{i-2},t_{i-1}]$, 2. 2.

$\mathbf{m}_{t_{i-1}}=\mathbf{X}_{i-1}-\boldsymbol{\eta}(t_{i-1})$, 3. 3.

$\boldsymbol{\Phi}_{t_{i-1}}=\mathbf{0}$.

Solving the system of ODEs (11), (13), (14), we obtain $\boldsymbol{\eta}(t_{i})$, $\mathbf{m}(t_{i})$, and $\boldsymbol{\Phi}(t_{i})$. The solution $\mathbf{m}(t_{i})$ will be a function of the initial value $\mathbf{X}_{i-1}-\boldsymbol{\eta}(t_{i-1})$, the interval length $\Delta t_{i}:=t_{i}-t_{i-1}$ and the SIR rates $\boldsymbol{\theta}_{i-1}$. To make this dependence explicit, we write $\mathbf{m}(t_{i}):=\boldsymbol{\mu}\left(\mathbf{X}_{i-1}-\boldsymbol{\eta}(t_{i-1}),\Delta t_{i},\boldsymbol{\theta}_{i-1}\right)$. Since (13) is a first order homogeneous linear ODE, the solution $\boldsymbol{\mu}\left(\mathbf{X}_{i-1}-\boldsymbol{\eta}(t_{i-1}),\Delta t_{i},\boldsymbol{\theta}_{i-1}\right)$ is a linear function of $\mathbf{X}_{i-1}-\boldsymbol{\eta}(t_{i-1})$. Hence, the transition from $\mathbf{X}_{i-1}$ to $\mathbf{X}_{i}$ follows the following Gaussian distribution:

[TABLE]

To summarize, the derived conditional Gaussian densities $\Pr(\mathbf{X}_{i}\mid\mathbf{X}_{i-1},\boldsymbol{\theta}_{i-1})$ allow us to compute the density of the latent SIR trajectory (4). As a result, our augmented posterior distribution of $\mathbf{X}_{0:T}$ and $\boldsymbol{\theta}_{0:T}$, shown in equation (3), can be computed up to proportionality constant and approximated via “standard” (not particle filter) MCMC approaches.

2.3 Reparameterization, priors, and MCMC algorithm

2.3.1 Reparameterizing SIR rates

We have experimented with multiple parameterizations of our inhomogeneous SIR model and found that the following parameterization works best with our proposed MCMC algorithm for approximating the posterior distribution (3). First, recall that we allow SIR rates to vary with time. Since it is much more likely for the infection rate to be time variable, we are going to assume a constant removal/recovery rate $\gamma$. This leaves us with the following parameters: infection rates on a grid $\boldsymbol{\beta}$, removal rate $\gamma$, and initial SIR state $\mathbf{X}_{0}=(S_{0},I_{0})$. Since we are interested in modeling an emerging infectious disease outbreak, we set the initial counts of susceptibles to $S_{0}=N-I_{0}$. Initial counts of infected individuals, $I_{0}$, is assumed to be low and treated as an unknown parameter with a lognormal prior distribution. Instead of the time-varying infection rate $\beta(t)$, we parameterize our SIR model with a time-varying basic reproduction number $R_{0}(t)=[\beta(t)N]/\gamma$. The reproduction number is interpreted as the average number of cases that one case generates over its infectious period in a completely susceptible population. Since our infection rate changes in a piecewise manner, the basic reproduction number varies over time in a piecewise manner too:

[TABLE]

where $R_{0_{i}}=[\beta_{i}N]/\gamma$ is the reproduction number corresponding to the time interval $[t_{i-1},t_{i})$. Let $R_{0}=R_{0_{0}}$ be the initial basic reproductive number and $\delta_{i}=\log\left(R_{0_{i}}/R_{0_{i-1}}\right)/\sigma$ be a normalized log ratio of $R_{0}(t)$ over two successive time intervals. Then, interval-specific basic reproduction numbers can be written as

[TABLE]

where we assume a priori that $\delta_{i}$s are independent standard normal random variables.

This construction implies that log-transformed piecewise constant reproduction numbers, $\log(R_{0_{i}})$s, a priori follow a first order Gaussian Markov random field (GMRF) with standard deviation $\sigma$ that controls the a priori smoothness of $R_{0}(t)$ trajectory [Rue, 2001, Rue and Held, 2005]. In addition to speeding MCMC convergence, working with $R_{0}(t)$ is convenient, because this trajectory is dimensionless and retains its interpretation when one changes the population size $N$. The initial $R_{0}$ is assigned a $\mbox{lognormal}(a_{1},b_{1})$ prior. We use a $\mbox{lognormal}(a_{2},b_{2})$ prior for the inverse of standard deviation $1/\sigma$.

2.3.2 Reparameterizing SIR latent trajectories

We reparameterize the latent SIR trajectory $\mathbf{X}_{1:T}$ with a sequence of independent Gaussian random variables $\boldsymbol{\xi}_{1:T}$, following a non-centered parameterization framework of Papaspiliopoulos et al. [2007]. According to formula (15), conditional on $\mathbf{X}_{i-1}$, $\mathbf{X}_{i}$ can be written as

[TABLE]

where $\boldsymbol{\xi}_{i}\overset{\text{iid}}{\sim}\mathcal{N}(0,\mathbf{I})$ for $i=1,\ldots,T$ and $\mathbf{I}$ is a $2\times 2$ identity matrix. In our parameterization, we will treat $\boldsymbol{\xi}_{1:T}$ as random latent variables and the SIR latent trajectory $\mathbf{X}_{1:T}$ as a deterministic transformation of $\boldsymbol{\xi}_{1:T}$. More details about our non-centered parameterization of $\mathbf{X}_{1:T}$ can be found in Section A-3 of the Appendix.

2.3.3 MCMC algorithm

Using our new parameterization, we are now interested in the posterior distribution of the initial number of infected individuals, $I_{0}$, removal rate, $\gamma$, the initial basic reproduction number, $R_{0}$, standardized vectors, $\boldsymbol{\delta}_{1:T}$ and $\boldsymbol{\xi}_{1:T}$, and GMRF standard deviation, $\sigma$:

[TABLE]

The latent variables $\mathbf{X}_{0:T}$ and parameter vector $\boldsymbol{\theta}_{0:T}$ are deterministic functions of new parameters $I_{0}$, $\gamma$, $R_{0}$, $\boldsymbol{\delta}_{1:T}$, $\boldsymbol{\xi}_{1:T}$, and $\sigma$. We use the following MCMC with block updates to approximate this posterior distribution. We update high dimensional vector $\mathbf{U}=(\log(R_{0}),\boldsymbol{\delta}_{1:T},\log(\sigma))$ using the efficient elliptical slice sampler [Murray et al., 2010]. Vector $\boldsymbol{\xi}_{1:T}$ is updated the same way in a separate step. Initial number of invected individuals $I_{0}$ and removal rate $\gamma$ are updated using univariate Metropolis steps. The full procedure is described in Algorithm 2, which together with details of the elliptical slice sampler can be found in Section A-4.1 of the Appendix. After MCMC is done, we report posterior summaries using natural parameterization. For example, we report posterior medians and 95% Bayesian credible intervals (BCIs) of the piecewise latent reproduction number trajectory, $R_{0_{i}}$ for $i=0,\ldots,T$, and latent trajectory $\mathbf{X}_{0:T}$.

2.3.4 Implementation

Our R package called LNAPhylodyn provides an implementation of our MCMC algorithm. The package code is publicly available at https://github.com/MingweiWilliamTang/LNAphyloDyn. This repository also contains scripts that should allow one to reproduce key numerical results in this manuscript.

3 Simulation experiments

3.1 Simulations based on single genealogy realizations

In this section, we use simulated genealogies to assess performance of our LNA-based method and to compare it with an ODE-based method, where we replace equation (18) with its simplified version: $\mathbf{X}_{i}=\boldsymbol{\eta}(t_{i})$. Under our assumption of a fixed and known genealogy and constant $R_{0}$, our ODE-based method closely resembles previously developed methods by Volz et al. [2009] and Volz and Siveroni [2018]. To compare ODE-based and LNA-based models in a Bayesian nonparametric setting, we equip the ODE model with the GMRF prior for time-varying $R_{0}(t)$, described in Section 2.3.1. We use the same MCMC algorithm for both LNA-based and ODE-based models, except we do not have a separate step to update latent vector $\boldsymbol{\xi}_{1:T}$ (equivalently, $\mathbf{X}_{0:T}$) in the ODE-based inference. See Algorithm 3 in the Appendix for a more detailed description of the ODE-based MCMC.

The simulation protocol consists of two steps. First, given the population size $N$ and pre-specified parameters $\gamma$, $I_{0}$, and $R_{0}(t)$, we simulate one realization of the SIR population trajectory based on the MJP using the Gillespie algorithm [Gillespie, 1977]. Next, we generate realistic lineage sampling times and simulate coalescent times from the distribution specified by density (2) using a thinning algorithm by Palacios and Minin [2013].

We test LNA-based and ODE-based methods under three “true” $R_{0}(t)$ trajectories over the time interval $[0,90]$:

Constant (CONST) $R_{0}(t)$. $R_{0}(t)=2.2$ for $t\in[0,90]$. Recovery rate $\gamma=0.2$. Initial counts of infected individuals $I_{0}=1$. Total population size is $N=$ 100,000. 2. 2.

Stepwise decreasing (SD) $R_{0}(t)$. $R_{0}(t)=2,t\in[0,30)$, $R_{0}(t)=1,t\in[30,60)$ and $R_{0}(t)=0.6,t\in[60,90]$. Recovery rate $\gamma=0.2$. Initial counts of infected individuals $I_{0}=1$. Total population size $N=$ 1,000,000. 3. 3.

Non-monotonic (NM) $R_{0}(t)$. $R_{0}(t)=1.4\times 1.015^{0.5t},t\in[0,30]$, $R_{0}(t)=1.750\times 0.975^{t-30},t\in[30,80]$ and $R_{0}(t)=0.4583,t\in[80,90]$. Recovery rate $\gamma=0.3$. Initial counts of infected individuals $I_{0}=3$. Total population size $N=$ 1,000,000.

For all simulations, we use $\mbox{lognormal}(1,1)$ prior for $I_{0}$. The parameters of the lognormal priors for the initial $R_{0}$ and inverse standard deviation $1/\sigma$ are set to $a_{1}=0.7,b_{1}=0.5$ and $a_{2}=3,b_{2}=0.2$ respectively, in such a way that a priori $R_{0}(t)$ trajectory stayed within a reasonable range of $[0,5]$ with 0.9 probability. We assign an informative prior for $\gamma$ in each simulation scenario, because prior information about this parameter is usually available: (1) CONST: $\gamma\sim\mbox{lognormal}(-1.7,0.1)$, (2) SD: $\gamma\sim\mbox{lognormal}(-1.7,0.1)$, (3) NM: $\gamma\sim\mbox{lognormal}(-1.2,0.1)$. We set the grid size to $T=36$, with $t_{i}-t_{i-1}=2.5$ for $i=1,\ldots,36$. For both LNA-based and ODE-based methods, we use 300,000 MCMC iterations. All MCMC chains appeared to converge (trace plots are shown in Section A-5.3.1 of the Appendix). The effective sample sizes of all unknown quantities were above 100.

The first row of Figure 3 shows point-wise posterior medians and 95% BCIs for the basic reproduction number trajectory, $R_{0}(t)$. Our LNA-based method performs well in capturing the continuous dynamics of $R_{0}(t)$. Though our approach may not perfectly catch the discontinuous changes in $R_{0}$ in the SD scenario, the method provides BCIs that are able to capture most of the $R_{0}(t)$ trajectory. The ODE-based method yields similar results in the CONST case and the SD case, but fails to capture the decreasing trends in the NM scenario.

The second row in Figure 3 shows posterior summaries of removal rate $\gamma$. Both LNA-based and ODE-based methods provide good estimates in the CONST scenario, with posterior modes centered at the true value and higher posterior densities at truth when compared with the prior. In the SD and NM scenarios with the time varying $R_{0}(t)$, the posterior estimates from the LNA-based method and ODE-based method, though still centered at the truth, do not differ much from the prior distribution.

Posterior summaries of $S(t)$ and $I(t)$ are depicted in the third and fourth rows of Figure 3. The two methods produce similar results in the CONST and SD scenario, as both of them have narrow BCIs covering the true trajectories. However, in the NM case, while the LNA-based method manages to recover the latent SIR trajectory trend, the BCIs from the ODE-based method fail to cover the true prevalence trajectory in the middle and at the end of the epidemic. Somewhat counterintuitively, LNA-based method produces BCIs for the latent trajectories, $S(t)$ and $I(t)$, that are narrower than its ODE counterparts. We suspect this is a result of the ODE-based method poor estimation of the basic reproduction number trajectory at the end of the epidemic.

3.2 Frequentist properties of posterior summaries

In this Section, we design a simulation study based on repeatedly simulating SIR trajectories using MJP with pre-specified parameters. The simulations are based on the non-monotonic $R_{0}(t)$ trajectory scenario in Section 3.1 with the same parameter setup, except the parameters of the lognormal prior for the initial $R_{0}$ are set to $a_{1}=0.7,b_{1}=0.3$. Simulating SIR dynamics under low initial number of infected individuals $I_{0}$ can end up with low prevalence trajectories that end at the beginning of the epidemic, or trajectories having unrealistically high prevalence, which are less likely to be observed during real infectious disease outbreaks. Therefore, while simulating SIR trajectories we reject such “unreasonable” realizations to arrive at 100 simulated trajectories. The details of the rejection criteria are given in Section A-5.2 of the Appendix. For each simulated SIR trajectory, a realization of a genealogy is generated using the structured coalescent process. We use both LNA-based and ODE-based model to approximate the posterior distribution of model parameters and latent variables for each genealogy.

We use three metrics to evaluate models based on their estimates of $R_{0}(t)$ and $I(t)$: average error of point estimates (posterior medians), width of credible intervals, and frequentist coverage of credible intervals. Since the value of $R_{0}(t)$ is greater than 0 and usually upper-bounded by 20 (i.e, it stays within the same order of magnitude), we will measure accuracy using an unnormalized mean absolute error (MAE):

[TABLE]

where $\hat{R}_{0_{i}}$ is the posterior median of $R_{0}(t_{i})$. In contrast, $I(t)$ varies from one at the beginning of the epidemic to thousands at the peak, so to evaluate accuracy of prevalence estimation we use the mean relative absolute error (MRAE):

[TABLE]

where $\hat{I}_{i}$ is the posterior median of I$(t_{i})$. We assess precision of $R_{0}(t)$ estimation based on the mean credible interval width (MCIW):

[TABLE]

where $\hat{R}^{0.025}_{0_{i}}$ and $\hat{R}^{0.975}_{0_{i}}$ denote the lower and upper bounds of the 95% BCI for $R_{0_{i}}$. Similar as our measure of accuracy, precision of $I(t)$ estimation is quantified via mean relative credible interval width (MRCIW):

[TABLE]

where $\hat{I}^{0.025}_{i}$ and $\hat{I}^{0.975}_{i}$ specify the lower and upper bounds of the $95\%$ BCI of $I(t_{i})$. In addition, we compute the “envelope” (ENV) — a measure of coverage of BCIs the true trajectory — for $R_{0}(t)$ and $I(t)$ as follows:

[TABLE]

Sampling distribution boxplots of $R_{0}(t)$ posterior summaries are depicted in the left three plots of Figure 4. The LNA-based method yields significantly lower MAE compared with the ODE-based method. As a trade-off, the MCIWs produced by the LNA-method are generally higher, as expected since the LNA-based method incorporates the stochasticity in the population dynamics. With less bias and wider BCIs, the LNA-based method BCIs result in better coverage than ODE-based BCIs, as shown by the envelope boxplots.

Sampling distribution boxplots of $I(t)$ posterior summaries, shown in Figure 4, are similar to the $R_{0}(t)$ results, with the LNA-based method generally having lower MRAEs, higher MRCIWs and a better coverage/envelope than the ODE-based method. Again, somewhat counterintuitively, the MRCIWs for the LNA-based method are smaller than the ODE counterparts. This is likely caused by significant bias in $R_{0}(t)$ estimation by the ODE-based method.

We also report the absolute error (AE) and 95% BCI widths for removal rate $\gamma$ in Figure 4. We note that an informative prior has been chosen for $\gamma$, because this parameter is weekly identifiably from genetic data alone. The LNA-based method yields a slightly higher AE than the ODE method. Both methods produce similar BCI widths.

4 Analysis of Ebola outbreak in West Africa

We apply our LNA-based method to the Ebola genealogies reconstructed from molecular data collected in Sierra Leone and Liberia during the 2014–2015 epidemic in West Africa [Dudas et al., 2017]. We use a Sierra Leone genealogy, depicted in the top left plot of Figure 5, which was estimated from 1010 Ebola virus full genomes sampled from 2014-05-25 to 2015-09-12 in 15 cities. The Liberia genealogy, shown in the top left plot of Figure 6, was estimated from a smaller number of samples: 205 Ebola virus full genomes sampled from 2014-06-20 to 2015-02-14. The original sequence data and the reconstructed genealogies are publicly available at https://github.com/ebov/space-time.

When Ebola virus infections were detected in West Africa in mid-Spring of 2014, various intervention measures were proposed and implemented to change behavior of individuals in the populations through which Ebola was spreading. Border closures, encouragement to reduce individual day-to-day mobility, and recommendations on changing burial practices were among the broad spectrum of interventions attempted by multiple countries. It is reasonable to expect that these intervention measures resulted in lowering the contact rates among members of the populations, which in turn reduced the infection rate, or equivalently the basic reproduction number.

When analyzing the Sierra Leone and Liberia genealogies, we rely on conclusions of Dudas et al. [2017] and assume the population in each country to be well mixed. Furthermore, we assume Ebola spread to follow SIR dynamics. For each country, the population size is specified based on its census population size in 2014, with $N=$ 7,000,000 for Sierra Leone and $N=$ 4,400,000 for Liberia. As in our simulation study, we use the lognormal prior for $R_{0}$ with $a_{1}=0.7$ and $b_{1}=0.5$ and the lognormal prior for the inverse standard deviation $1/\sigma$ with $a_{2}=3,b_{2}=0.2$. Recall that this prior setting ensures that a priori $R_{0}(t)$ stays within a reasonable range of $[0,5]$ with probability 0.9. For removal rate $\gamma$, we use an informative lognormal prior with mean 3.4 and variance 0.2 based on previous studies [Towers et al., 2014]. The parameter $1/\gamma$, interpreted as the length of the infectious period, is expected to be 8-18 days for each country a priori. The total time span for each genealogy is divided evenly into 40 intervals, which results in grid interval lengths, $\Delta t_{i}$s, to be 12.41 days for Sierra Leone and 6.9 days for Liberia. We run the MCMC algorithm in Section 2.3 for 3,000,000 iterations for Sierra Leone data and 750,000 iterations for Liberia data. The posterior samples are obtained by discarding the first 100,000 iterations and saving every 30th iteration afterward. The trace plots in Section A-5.3.2 of the Appendix indicate the MCMC algorithm has converged and achieved good mixing in each case.

Figures 5 and 6 show results for Sierra Leone and Liberia respectively, with intervention events mapped onto the calendar time on the x-axis. Our LNA-based method estimates the initial $R_{0}$ in Sierra Leone during 2014–2015 to be 1.68, with 95% BCI of $(1.33,2.23)$. Similarly, $R_{0}$ in Liberia during 2014–2015 has a point estimate 1.67 and a 95% BCI $(1.29,2.24)$. Our estimate of initial $R_{0}$ in Sierra Leone is consistent with the estimates of Stadler et al. [2014], who fitted multiple birth-death models to 72 sequences at the early stages of the outbreak. Volz and Pond [2014] used a susceptible-exposed-infectious-recovered (SEIR) model with a constant $R_{0}$ and estimated it to be 2.40 (CI: $(1.54,3.87)$). Althaus [2014] assumed an exponentially decaying $R_{0}(t)$ with an estimated initial $R_{0}$ of 2.52 (CI: $(2.41,2.67)$). The discrepancies between our and SEIR-based estimates are not unexpected, because SEIR models generally yield higher $R_{0}$ estimates than SIR models when applied to the same dataset [Wearing et al., 2005, Keeling and Rohani, 2011]. Our estimated $R_{0}$ for Liberia is in agreement with results of Althaus [2014], who fitted a SEIR model to incidence data and arrived at an estimated $R_{0}$ of 1.59 (CI: $(1.57,1.60)$).

The $R_{0}(t)$ dynamics in the two countries share a similar pattern: with (1) a decreasing trend that starts in Spring/Summer of 2014, (2) a stable/constant period until the end of September 2014 and (3) a final decrease below 1.0 (epidemic is contained) around November 2014. Since the number of susceptible individuals did not change significantly over the course of the epidemic, relative to the total population size, the basic and effective reproduction numbers, $R_{0}(t)=\beta(t)N/\gamma$ and $R_{\text{eff}}(t)=\beta(t)S(t)/\gamma$, are approximately equal. This allows us to compare our $R_{0}(t)$ estimation results with previously estimated changes in $R_{\text{eff}}(t)$. Our estimation of early $R_{0}(t)$ dynamics in Sierra Leone agrees with results of Stadler et al. [2013], who concluded that the effective reproduction number did not significantly decrease until mid June. Our estimated $R_{0}(t)$ trajectory suggests that later interventions, such as border closures and release of burial guides, may have been helpful in controlling the spread of the disease. The infectious period for Sierra Leone epidemic is estimated to be 11.2 days with a 95% BCI (7.6,16). For Liberia, the infection period has a point estimate of 9.8, with a 95% BCI $(6.87,14.05)$. The posterior median of the total number of infected individuals (final epidemic size) is 7,284 and its 95% BCI is $(3397,14870)$ for Sierra Leone, which is close to 8,706 total confirmed number of cases reported by Centers for Disease Control (CDC). Liberia had a smaller epidemic than Sierra Leone, with estimated total infected individuals being 2,842 and a 95% BCI of $(1296,6173)$. These results are also in agreement with 3,163 total confirmed cases from CDC.

We perform an out-of-sample validation by comparing our results with weekly reported confirmed incidence in Sierra Leone and Liberia from the World Health Organization (WHO). The posterior predictive weekly incidence at time $t$, denoted by $\hat{N}(t)$, is approximated by

[TABLE]

where $\hat{\beta}(t),\hat{S}(t)$ and $\hat{I}(t)$ are the posterior estimates of the infection rate, number of susceptible and number of infected individuals at time $t$ respectively, and $\Delta t:=7/365$ corresponds the time interval of one week. We plot the posterior predictive estimates of weekly incidence together with the corresponding weekly reported confirmed incidence. For both countries, our model-based incidence 95% BCIs cover the reported incidence counts from WHO, suggesting that our time varying SIR model can estimate incidence well from genetic data alone. We note that our estimated latent incidence should be greater than the reported incidence, because not all Ebola cases were reported and recorded. However, the discrepancy between latent and reported incidence should not be large, because Ebola reporting rate was high. For example, Scarpino et al. [2014] estimated that 83% of Ebola cases were reported.

We also report results from the ODE-based method and superimpose these results over LNA-based results on Figures 5 and 6. For the relatively small Liberia genealogy, the ODE-based and LNA-based methods yield similar parameter estimates. However, the larger Sierra Leone genealogy produces substantial differences between ODE-based and LNA-based estimates of the $R_{0}(t)$. The ODE-based method captures the decreasing trend of $R_{0}(t)$ in Spring and Summer of 2014, but provides narrow BCIs with unrealistic short term fluctuations in the basic reproduction number trajectory.

5 Discussion

In this paper, we propose a Bayesian phylodynamic inference method that can fit a stochastic epidemic model to an observed genealogy estimated from infectious disease genetic sequences sampled during an outbreak. Our statistical model can be viewed as semi-parametric: with (1) a parametric SIR model describing the infectious disease dynamics and (2) a non-parametric GMRF-based estimation of the time varying basic reproduction number. To the best of our knowledge, this is the first method combining a Bayesian nonparametric approach with a deterministic or stochastic SIR model for phylodynamic inference (although see [Xu et al., 2016] for a similar approach applied to more traditional epidemiological data). Our use of LNA allows us to devise an efficient MCMC algorithm to approximate high dimensional posterior distribution of model parameters and latent variables. Our LNA-based method produces posterior summaries with better frequentist properties than the state-of-the-art ODE-based method, underscoring the importance of working with stochastic models even in large populations. We showcase our method by applying it to the Ebola genealogies estimated from viral sequences collected in Sierra Leone and Liberia during the 2014–2015 outbreak. Our nonparametric estimates of $R_{0}(t)$ in Sierra Lione and Liberia suggest that the basic reproduction number decreased in two-stages, where the second stage brought it below 1.0 — a sign of epidemic containment. The second stage of $R_{0}(t)$ decrease closely follows in time implementation of interventions, pointing to their effectiveness.

Our method relies on the assumption that population dynamics follow a SIR model. However, it may be desirable to be able to relax this assumption. For example, in Ebola spread modeling some authors used a SEIR model that assumes a latent period during which an infected individual is not infectious [Althaus, 2014, Volz and Siveroni, 2018]. One future direction of this work is to generalize the LNA-based method to fit complicated compartmental epidemic models, including models with multi-stage infections like SEIR model and models with the population stratified by sex, age, geographic location, or other demographic variables. The structured coalescent likelihoods under these models may not have closed-form expressions. However, Volz [2012] and Müller et al. [2017] propose several strategies to approximate structured coalescent likelihoods. It would be interesting to combine these approximation strategies with LNA to estimate parameters of complex stochastic epidemic models from genealogies.

The experiments in Section 3.1 and Appendix Section A-6 indicate that one has to pay close attention to parameter identifiability when fitting SIR models to genealogies or to sequence data directly. Identifiability may not be a problem under an assumption of a constant $R_{0}(t)$. However, the removal rate tends to be only weakly identifiable in the scenarios with a time-varying basic reproduction number, in which the estimation can be sensitive to the choice of priors. In Section A-6 of the Appendix, we demonstrate that putting a weakly informative prior on the removal rate can cause bias not only in the estimation for removal rate, but also can lead to a failure in recovering the reproduction number and latent population dynamics. Therefore, successful inference of SIR model parameters using genealogical data should rely on a sound informative prior for the removal rate. This constraint is not a big shortcoming in practice, since prior information about the removal rate, or mean length of the infectious period, is usually readily available from patient hospitalization data [WHO Ebola Response Team, 2014].

Since parameter identifiability is a recurring problem in infectious disease modeling, integration of multiple sources of information is of great interest. Using particle filter MCMC, Rasmussen et al. [2011] demonstrated that jointly analyzing genealogy and incidence case counts considerably reduces the uncertainty in both estimation of latent population trajectory and SIR model parameters, compared with estimation based on a single source of information. We plan to use our LNA-based framework to perform similar integration of genealogical data and incidence time series. Another possible source of information is the distribution of genetic sequence sampling times. Karcher et al. [2016] proposed a preferential sampling approach that explicitly models dependence of the sampling times distribution on the effective population size. The authors demonstrated that accounting for preferential sampling helps decrease bias and results in more precise effective population size estimation. It would be interesting to incorporate preferential sampling into our LNA-based framework by assuming a probabilistic dependency between sampling times and latent prevalence $I(t)$.

Our method assumes a genealogy/phylogenetic tree is given to us. In reality, genealogies are not directly observed and need to be inferred from molecular sequences. Ideally, uncertainty in the genealogy should be handled by building a Bayesian hierarchical model and integrating over the space of genealogies using MCMC. In fact, implementations of such Bayesian hierarchical modeling already exist for nonparametric, birth-death, and ODE-based phylodynamic approaches [Drummond et al., 2005, Minin et al., 2008, Gill et al., 2013, Stadler et al., 2013, Volz and Siveroni, 2018]. Therefore, an important future direction will be to extend our LNA framework to fitting stochastic epidemic models to molecular sequences instead of genealogies. Similarly to the structured coalescent model implementation of Volz and Siveroni [2018], the easiest way to achieve this will be integration of our LNA MCMC algorithm into popular open source phylogenetic/phylodynamic software packages, such as BEAST, BEAST2, and RevBayes [Suchard et al., 2018, Bouckaert et al., 2014, Höhna et al., 2016].

Acknowledgements

We thank Jon Fintz for discussing the details of implementing Linear Noise Approximation. We are grateful to Michael Karcher and Julia Palacios for patiently answering questions about nonparametric phylodynamic inference and phylodyn package. M.T and V.N.M. were supported by the NIH grant R01 AI107034. M.T., T.B., and V.N.M. were supported by the NIH grant U54 GM111274. G.D. was supported by the Mahan postdoctoral fellowship from the Fred Hutchinson Cancer Research Center. This work was supported by NIH grant R35 GM119774-01 from the NIGMS to TB. TB is a Pew Biomedical Scholar.

A-1 A general framework for stochastic kenetic models

A-1.1 Stochastic model generalization

In Section 2, we provide an example of the linear noise approximation (LNA) for the SIR model. The LNA framework can be also generalized to other types of the stochastic kinetic models in Infectious Disease Epidemiology and in Systems Biology. Here, we give a general representation of the stochastic kinetic model by viewing it as a reaction network system. The notation is based on the work of Fearnhead et al. [2014].

Let’s start with a reaction system with $d$ reactants $\mathcal{X}_{1},\ldots,\mathcal{X}_{d}$ and $q$ reactions. Without loss of generality, each reaction is assumed to have a constant rate parameter $\theta_{i}$ for $i=1,\ldots,q$ and $\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{q})$ denotes the rate vector of the system (this framework can be extended to handle stochastic kinetic models with time-varying rates as in Section 2 of the main text). The transition event in the $i$th reaction ($i=1,\ldots,q$) has the following form:

[TABLE]

where $\tilde{a}_{ij}$ and $\bar{a}_{ij}$ are non-negative integers representing the number of $\mathcal{X}_{j}$ in the $i$th reaction equation. In a compartmental stochastic epidemic model, the coefficient $\tilde{a}_{ij}$ will be either 0 or 1. The transitions in the reaction system can be encoded in an effect matrix,

[TABLE]

with each row corresponding to a certain type of reaction event and each column representing the change in the counts of reactants. Let $X_{j}(t)$ denote counts/population of the $\mathcal{X}_{j}$ at $t$, and the population state at time t can be tracked by vector $\mathbf{X}(t):=(X_{1}(t),\ldots,X_{d}(t))$. Let $h_{i}$ denote the reaction rate of the $i$th reaction, where $h_{i}$ can be written as

[TABLE]

Hence, following the same notation as in Section 2.2.1 of the main text, the rate vector $\mathbf{h}$ and the rate matrix $\mathbf{H}$ can be defined as

[TABLE]

Given the above notation, the deterministic ordinary differential equation model of the reaction system can be written as

[TABLE]

where $\mathbf{x}_{0}$ is a vector of initial counts of reactants $\mathcal{X}_{1},\ldots,\mathcal{X}_{d}$.

A-1.1.1 Example: SEIR model

The above general representation of stochastic kinetic models can be directly applied to stochastic epidemic models. Here, we illustrate this on a Susceptible-Exposed-Infected-Recovery (SEIR) model. SEIR model is an extension of the SIR model that assumes a latent period called “Exposed”, in which an infected individual does not have the ability to infect others. The exposed individual will eventually become infectious with rate $\mu$. As in the SIR model, an infectious individual has removal/recovery rate $\gamma$. The transition events between different states of the SEIR model are depicted in Figure A-1.

Following the stochastic kinetic model representation, the SEIR model can be viewed as a reaction system of four reactants — susceptible, exposed, infected, and recovered individuals — and the following three “reactions”:

[TABLE]

Since the recovered population never interacts with individuals in other compartments, we will only keep track of the counts of susceptible, exposed, and infectious individuals at time $t$, denoted by $S(t)$, $E(t)$, $I(t)$ respectively. The effect matrix $\mathbf{A}$ for the SEIR model can be written as:

[TABLE]

with columns representing compartments and rows representing reactant changes during reaction events.

If we let $\mathbf{X}(t)=(S(t),E(t),I(t))$ denote the state vector at time $t$, then the rate vector $\mathbf{h}$ for the SEIR model is

[TABLE]

A-2 Derivation of the linear noise approximation

A-2.1 SDE approximation for MJP

A stochastic way to approximate the MJP model is to use the Stochastic Differential Equation (SDE) approximation, also known as the chemical Langevin equation (CLE) [Gillespie, 2000]. The SDE method can be viewed an approximation of the MJP at time $t$, obtained by applying a normal approximation to the Poisson distributed number of state transitions in a small interval of time $(t,t+\Delta t)$ [Gillespie, 2000, Wallace, 2010]. The deterministic part in SDE corresponds to the right hand side of ODE (11) and stochastic part is related to the variance of the system. The SDE for general stochastic kinetic models can be written as

[TABLE]

where $\mathbf{W}_{t}$ denote a $d$ dimensional Wiener process and the square root $\sqrt{\cdot}$ means the Cholesky triangle of the $d\times d$ covariance matrix.

A-2.2 LNA approximation of the SDE

Since in the main text we assume the rate $\boldsymbol{\theta}(t)$ varies in a piecewise constant way, without loss of generality, we use the notation $\boldsymbol{\theta}$ for the rate in a given time interval where the LNA is applied.

Theorem A-2.1** (Linear Noise Approximation for SDE).**

Let $\boldsymbol{\eta}(t)$ be the solution of ordinary differential equation (11) with initial value $\boldsymbol{\eta}_{0}$. Let $N$ be the system size, which is the total number of individuals in the system (In SIR model, $N$ will be the total population, i.e $N=S+I+R$), $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{q})$ denote the vector of rate parameters in $q$ reactions. Then the solution $\mathbf{X}(t)$ of the SDE (A-11) satisfies the following equation

[TABLE]

as $N\rightarrow+\infty$.

Proof.

The following derivation is based on [Wallace, 2010].

We rescale both the compartment size and reaction rates as follows:

[TABLE]

where $m_{i}=\sum_{j=1}^{d}\tilde{a}_{ij}$ is the sum of coeffcients in the left hand side of $i$-th reaction as in Section A-1. The transformed $\tilde{\mathbf{X}}(t)$ represents the proportion of individuals/reactants each compartment with respect to the total population size. Then we have $\mathbf{h}(\mathbf{X}(t),\boldsymbol{\theta})=N\mathbf{h}(\tilde{\mathbf{X}}(t),\tilde{\boldsymbol{\theta}})$ and $\mathbf{F}(\boldsymbol{\eta}(t),\boldsymbol{\theta})=\mathbf{F}(\tilde{\boldsymbol{\eta}}(t),\tilde{\boldsymbol{\theta}})$. Hence, the SDE (A-11) becomes

[TABLE]

Let $\tilde{\boldsymbol{\eta}}(t)$ be the solution of the ODE

[TABLE]

and we have $\boldsymbol{\eta}(t)=N\tilde{\boldsymbol{\eta}}(t)$, where $\boldsymbol{\eta}(t)$ is the solution of the ODE (11). $\tilde{\boldsymbol{\eta}}(t)$ can be viewed as a scaled version solution of ODE (11). Let $\boldsymbol{\xi}(t)=\sqrt{N}\left(\tilde{\mathbf{X}}(t)-\tilde{\boldsymbol{\eta}}(t)\right)=\dfrac{1}{\sqrt{N}}\left(\mathbf{X}(t)-\boldsymbol{\eta}(t)\right)$ denote the scaled residual, then the rescaled compartment size vector $\tilde{\mathbf{X}}(t)$ can be written as

[TABLE]

After using first order Taylor expansion of $\mathbf{h}(\tilde{\mathbf{X}}(t),\tilde{\boldsymbol{\theta}})$ and $\mathbf{H}(\tilde{\mathbf{X}}(t),\tilde{\boldsymbol{\theta}})$ around $\tilde{\mathbf{X}}=\tilde{\boldsymbol{\eta}}(t)$, the SDE (A-15) becomes

[TABLE]

where \mathbf{F}(\tilde{\boldsymbol{\eta}}(t),\boldsymbol{\theta}):=\dfrac{\partial\mathbf{A}^{T}\mathbf{h}(\tilde{\mathbf{X}}(t),\boldsymbol{\theta})}{\partial\tilde{\mathbf{X}}}\Big{|}_{\tilde{\mathbf{X}}=\tilde{\boldsymbol{\eta}}(t)} is the Jacobian matrix of the deterministic part $\mathbf{A}^{T}\mathbf{h}\left(\tilde{\mathbf{X}}(t),\boldsymbol{\theta}\right)$ in (11) at $\tilde{\boldsymbol{\eta}}(t)$. By subtracting (A-16) and multiplying by $\sqrt{N}$ on the two ends, the above equation becomes a differential equation with respect to $\boldsymbol{\xi}$:

[TABLE]

After multiplying by $\sqrt{N}$, the above equation gives us (A-12). ∎

Recall that $\mathbf{M}(t)$ is the solution of (12) with initial condition $\mathbf{M}(0)=\mathbf{X}_{0}-\boldsymbol{\eta}_{0}$. We can use $\boldsymbol{\eta}(t)+\mathbf{M}(t)$ as an approximation of $\mathbf{X}(t)$. Based on the local Lipschitz property of $\mathbf{F}(\boldsymbol{\eta}(t),\boldsymbol{\theta})$ with respect to $t$ and $\mathbf{A}^{T}\mathbf{H}(\boldsymbol{\eta}(t),\boldsymbol{\theta})$, $\mathbf{X}(t)$ can be approximated by $\boldsymbol{\eta}(t)+\mathbf{M}(t)$ with

[TABLE]

for fixed $t$ as system size $N\rightarrow+\infty$.

A-2.3 Derivation of equations (13) and (14) in the main text

Lemma A-2.2** (Solution of linear ODE system).**

Let $\mathbf{F}(t)\in\mathbb{R}^{d\times d}$ and $\mathbf{X}(t)\in\mathbb{R}^{d}$ be function of defined on $\{t:t\geq 0\}$ that satisfies the following linear ODE

[TABLE]

For $t\geq 0$, the solution of (A-20) can be represented as

[TABLE]

where $\boldsymbol{\Sigma}(t,0)$ is the solution of ordinary differential equation in $\mathbb{R}^{d\times d}$

[TABLE]

Lemma (A-2.2) gives the solution of linear ODE. Hence, the solution for the main text linear ODE 13 is on $[t_{i-1},t]$ will be

[TABLE]

where $\mathbf{m}_{i-1}$ is the initial state at $t_{i-1}$ and $\boldsymbol{\Sigma}(t,t_{i-1})$ is the transition matrix by

[TABLE]

and $\mathbf{m}_{i-1}$ is the initial value for $\mathbf{m}$ at time $t_{i-1}$.

Theorem A-2.3**.**

Let $\{\mathbf{M}(t)\}_{t\geq 0}\in\mathbb{R}^{d}$ be stochastic process that satisfies the following stochastic differential equation,

[TABLE]

Then the solution of (A-25) is the Gaussian process

[TABLE]

with mean process $\mathbf{m}(t):=\mathbf{E}\left[\mathbf{M}(t)|\mathbf{M}(t_{0})\right]$ satisfies (13) and variance process $\boldsymbol{\Phi}(t):=\mbox{Var}\left[\mathbf{M}(t)|\mathbf{M}(t_{0})\right]$ satisfies (14).

Proof.

Define matrix function $\boldsymbol{\Sigma}(t,t_{0})$ as (A-24). First we apply the linear transform $\tilde{\mathbf{M}}(t)=\boldsymbol{\Sigma}^{-1}(t;t_{0})\mathbf{M}(t)$. Based on Ito’s lemma, (A-25) can be simplified as a SDE of $\tilde{\mathbf{M}}(t)$:

[TABLE]

with solution.

[TABLE]

Then the solution of $\mathbf{M}(t)$ is

[TABLE]

$\boldsymbol{\Sigma}(t,t_{0})\mathbf{M}(t_{0})$ in (A-28) is a deterministic function of $t$. The integral $\int_{t_{0}}^{t}\boldsymbol{\Sigma}^{-1}(s,t_{0})\sqrt{\mathbf{A}^{T}\mathbf{H}\left(\eta(t)\right)\mathbf{A}}\mathrm{d}\mathbf{W}_{s}$ in (A-28) should be Gaussian random variable with mean $\mathbf{0}$ since it is a linear combination of the increments of Brownian motion with different variance. Hence, the $\mathbf{M}(t)$ should be a Gaussian process. By taking the expectation of (A-28), the mean of $\mathbf{m}(t)=\mathbf{E}[\mathbf{M}_{t}]$ satisfies

[TABLE]

which corresponds to the solution of ODE (13).

For the variance process, from (A-28),

[TABLE]

By differentiation with respect to $t$, (A-30) becomes

[TABLE]

which is the result in (14). ∎

A-2.4 Relationship between LNA and other methods

The SDE approach can be viewed as a normal approximation based on a $\tau$-leaping step for the MJP. The LNA can be derived either directly from Taylor expansion of the transition probability of the MJP or the Taylor expansion of the transition density of the SDE. The ODE solution can be considered as a limit of the mean trajectory of the MJP when system size $N$ goes to infinity. ODE solution can also be viewed as the deterministic part for SDE (A-11) and the mean process for LNA based on (A-52). Figure  A-2 depicts relationships between different dynamical system representations as a diagram.

A-3 Non-centered parameterization

In LNA, the latent trajectory $\mathbf{X}(t)$ is decomposed into the deterministic part $\boldsymbol{\eta}(t)$ plus a stochastic part $\mathbf{M}(t)$ that follows a multivariate Gaussian distribution with mean $\mathbf{0}$. However, the population size at the $i$-th time interval $\mathbf{X}_{i}$ depends on rate parameter $\boldsymbol{\theta}$ and is correlated with other population sizes $\mathbf{X}_{j}$s in the trajectory, leading to mixing issues for the MCMC chain, especially when we introduce multiple change points for reproduction number $R_{0}$.

Here we take the idea of non-centered parameterization from Papaspiliopoulos et al. [2007], Bernardo et al. [2003] and reparameterize the latent trajectrory in terms of residuals $\mathbf{X}_{i}-\boldsymbol{\eta}_{i}$ for $i=1,\dots,T$. Given rate parameters $\boldsymbol{\theta}_{i-1}$, ODE solution $\boldsymbol{\eta}_{0:T}$, fundamental matrix $\boldsymbol{\Sigma}(t_{i},t_{i-1})$ and variance matrix $\boldsymbol{\Phi}_{i}$ in (14), the trajectory $\mathbf{X}_{0:T}$ can be parameterized using standard Gaussian noise $\boldsymbol{\xi}_{1:T}$ based on the following iterative equations:

[TABLE]

Let $\mathbf{M}_{i}:=\mathbf{X}_{i}-\boldsymbol{\eta}_{i}$ denote the residual in grid cell $i$. Based on (A-33), the residual process satisfies

[TABLE]

Then $\mathbf{M}_{0:T}$ can be viewed as a Gaussian Markov random field with mean $\mathbf{0}$ that follows the Markov property on a chain graph. Let $\boldsymbol{\Sigma}_{i}$ be the abbreviated notation of $\boldsymbol{\Sigma}(t_{i},t_{i-1})$ and $\mathbf{P}_{i}=\boldsymbol{\Phi}_{i}^{1/2}$. From (A-35) , $\mathbf{M}_{i}$ can be written as

[TABLE]

Since $\boldsymbol{\Sigma}_{i}$ and $\mathbf{P}_{i}$ are governed by rate parameters $\boldsymbol{\theta}_{i-1}$ and initial value $\mathbf{X}_{0}$, then we define the transform matrix $\mathbf{L}(\mathbf{X}_{0},\boldsymbol{\theta}_{0:T})\in\mathbb{R}^{2T\times 2T}$,

[TABLE]

A linear relationship between $\mathbf{X}_{1:T}$ and the reparameterized noise $\boldsymbol{\xi}_{1:T}$ can be established with the help of the above transform matrix $\mathbf{L}$,

[TABLE]

Instead of directly updating $\mathbf{X}_{1:T}$, we will apply the above transform and update the Gaussian noise $\boldsymbol{\xi}_{1:T}$ instead. The MCMC approach will focus on sampling parameter $I_{0},R_{0},\gamma,\boldsymbol{\delta}_{1:T}$, $\boldsymbol{\xi}_{1:T},\sigma$ with the posterior likelihood

[TABLE]

In summary,the transformation that allows us to move from parameterization in terms of $\mathbf{X}_{0:T},\theta_{0:T}$ to the parameterization in terms of $I_{0},R_{0},\gamma,\boldsymbol{\delta}_{1:T}$, $\boldsymbol{\xi}_{1:T},\sigma$ are based on the following equations:

$R_{0_{i}}:=R_{0}(t_{i})=R_{0}\cdot\exp(\prod_{j=1}^{i}\delta_{j}\sigma)$ — a function of $R_{0}$, $\boldsymbol{\delta}_{1:i}$ and $\sigma$. 2. 2.

$\beta_{i}:=\beta(t_{i})=\dfrac{NR_{0}(t_{i})}{\gamma}$ — a function of $R_{0}$, $\boldsymbol{\delta}_{1:i}$, $\sigma$ and $\gamma$. 3. 3.

$\boldsymbol{\theta}_{i}=(\beta_{i},\gamma)$ — a function of $R_{0}$, $\boldsymbol{\delta}_{1:i}$, $\sigma$ and $\gamma$. 4. 4.

$\boldsymbol{\theta}_{0:T}$ — a function of $R_{0}$, $\boldsymbol{\delta}_{1:T}$, $\sigma$ and $\gamma$. 5. 5.

$\mathbf{X}_{0}=(N,I_{0})^{T}$. 6. 6.

$\mathbf{X}_{1:T}=\boldsymbol{\eta}_{1:T}+\mathbf{L}(\mathbf{X}_{0},\boldsymbol{\theta}_{0:T})\boldsymbol{\xi}_{1:T}$ — a function of $R_{0}$, $\boldsymbol{\delta}_{1:T}$, $\sigma$, $\gamma$, $I_{0}$ and $\boldsymbol{\xi}_{1:T}$.

A-4 MCMC details

A-4.1 Elliptical slice sampler

The elliptical slice sampler, proposed by Murray et al. [2010], aims at sampling from posterior distributions associated with probability models with a latent a priori zero-mean Gaussian random vector $\mathbf{X}\in\mathbb{R}^{d}$ with covariance $\Sigma(\boldsymbol{\theta})$, i.e., $\mathbf{X}\sim\mathcal{N}(\mathbf{0},\Sigma(\boldsymbol{\theta}))$. We use $L(\mathbf{Y}|\mathbf{X},\boldsymbol{\theta})$ to denote the likelihood function for observed data $\mathbf{Y}$ given latent variable $\mathbf{X}$ and parameter $\boldsymbol{\theta}$. Hence, the target posterior distribution for $\mathbf{X}$ given is

[TABLE]

where $\pi(\boldsymbol{\theta})$ is the prior distribution for parameter $\boldsymbol{\theta}$. The goal of elliptical slice sampler is to obtain posterior samples of latent variable $\mathbf{X}$ from $p(\mathbf{X}|\mathbf{Y},\boldsymbol{\theta})$. The proposal step in elliptical sampling consists of two parts: (1) proposing an auxiliary random vector $\mathbf{Z}\in\mathbb{R}^{d}$ from distribution $\mathcal{N}\left(\mathbf{0},\Sigma(\boldsymbol{\theta})\right)$, (2) proposing a variable $\alpha\in[0,2\pi]$ as an angle parameter. In elliptical slice sampler, a new state $(\mathbf{X}^{\prime},\mathbf{Z}^{\prime})$ is proposed by rotating the previous state $(\mathbf{X},\mathbf{Z})$ with angle $\alpha$,

[TABLE]

For any given $\alpha$, this transition leaves the joint prior probability invariant, i.e,

[TABLE]

Hence, $(\mathbf{X}^{\prime},\mathbf{Z}^{\prime})$ are considered at the proposed state and the ratio and the propose transition probability from $(\mathbf{X},\mathbf{Z})$ to $(\mathbf{X}^{\prime},\mathbf{Z}^{\prime})$ should equal that from $(\mathbf{X}^{\prime},\mathbf{Z}^{\prime})$ to $(\mathbf{X},\mathbf{Z})$, i.e

[TABLE]

The algorithm for elliptical slice sampler is given in Algorithm 1.

Notice that iterations will stop only when a new sample is accepted. Hence, the elliptical slice sampler has acceptance rate $1$, meaning that it will always update the latent random vector $\mathbf{X}$ at each MCMC iteration.

A-4.2 MCMC algorithm for the LNA-based SIR model

In this framework, the observed data are the genealogy $\mathbf{g}$ estimated from a sample of sequences from virus hosts. The sufficient statistics for SIR structured coalescent likelihood would be the coalescent times $\mathcal{T}$ and sampling times $\mathcal{S}$. The unknown parameters and the latent variables are

The initial number of infected individuals: $I_{0}$. The initial population is parameterized as $\mathbf{X}_{0}=(N,I_{0})$, suppressing that $S_{0}=N-I_{0}$ and that there are no recovered individuals at time 0. 2. 2.

The initial basic reproduction number $R_{0}$. 3. 3.

The removal rate $\gamma$. 4. 4.

The hyperparameter $\sigma$ that controls the smoothness of $R_{0}(t)$ trajectory. 5. 5.

The parameters modeling the first order differences in $\log(R_{0}(t)):\boldsymbol{\delta}_{1:T}$. Note that assuming $\delta_{0}=1$, the infection rate $\beta_{i}$ can be obtained as

[TABLE]

The parameter $\boldsymbol{\theta}_{0:T}$ can be obtained from $R_{0},\boldsymbol{\delta}_{1:T},\gamma$ and $\sigma$. 6. 6.

Random noise for the population trajectory at $t_{1},\ldots,t_{T}$, i.e. $\boldsymbol{\xi}_{1:T}$, with $\boldsymbol{\xi}_{i}\sim_{iid}\mathcal{N}(0,\mathbf{I})$ a priori. The latent SIR trajectories $\mathbf{X}_{0:T}$ can be recovered from $\boldsymbol{\theta}_{0:T}$, $\mathbf{X}_{0}$ and random noise $\boldsymbol{\xi}_{1:T}$.

The MCMC update for parameters and latent variables is given in Algorithm 2.

A-4.3 MCMC algorithm for the ODE-based model

The MCMC algorithm for ODE-based method is similar to the LNA-based MCMC except there is no need to update the Gaussian noise $\boldsymbol{\xi}_{1:T}$ in the population trajectory. The MCMC updates of parameters and latent variables is given in Algorithm 3.

A-5 Details of the simulation study

A-5.1 Simulation details for Section 3.1 of the main text

Here we provide details for the specified sequence/lineage sampling times and number of samples in each simulation scenario:

CONST $R_{0}(t)$: Sampling times: $\mathcal{S}=\{5,10,50,70,80,90\}$, number of samples at each time: $\left\{2,20,300,300,200,200\right\}$. 2. 2.

SD $R_{0}(t)$: Sampling times: $\mathcal{S}=\{5,10,50,70,80,90\}$, number of samples at each time: $\left\{2,20,200,80,20,20\right\}$. 3. 3.

NM $R_{0}(t)$: Sampling times: $\mathcal{S}=\{5,30,50,70,80,90\}$, number of samples at each time: $\left\{2,50,250,100,20,20\right\}$.

A-5.2 Simulation details for Section 3.2 of the main text

The $R_{0}(t)$ trajectory in the simulations is set to

[TABLE]

which is depicted in the left plot of Figure A-4. The initial number of infected individuals is $I_{0}=3$ and the removal rate is set to $\gamma=0.3$. The population size is fixed to $N=1,000,000$. Epidemic trajectories are simulated using the SIR Markov jump process (MJP) and are accepted/rejected based on the following criteria:

Reject the SIR trajectories that ends before time $90$. The number of infected individuals should never drop to [math] for $t\in[0,90]$, i.e. $\min_{t\in[0,90]}I(t)>0$. 2. 2.

Reject the SIR trajectories with extremely high maximum prevalence: the maximum prevalence should be less or equal than 12,000, i.e., $\max_{t\in[0,90]}I(t)\leq 12000$. 3. 3.

Reject SIR trajectories with extremely low maximum prevalence. The maximum prevalence should be greater or equal than 600, i.e., $\max_{t\in[0,90]}I(t)\geq 600$.

The 100 simulated SIR prevalence trajectories are shown in the right plot of Figure A-4. We also plot the ODE trajectory under the same parameter setting.

A-5.3 Trace plots

A-5.3.1 Trace plots for simulations from Section 3.1 of the main text

Figure A-5 shows the trace plots of the log-posterior for the LNA-based method and ODE-based method in the three simulation scenarios from Section 3.1. The effective sample sizes (ESSs) for all parameters range between 100 to 400.

A-5.3.2 Trace plots for Ebola data

Figures A-6 and A-7 show trace plots of parameters $R_{0},I_{0},\gamma,\sigma$ for the LNA-based model and ODE-based model respectively applied to the Sierra Leone genealogy. Figures A-8 and A-9 show the analogous trace plots for the analysis of the Liberia genealogy. We also list posterior medians, 95% BCIs, and ESSs for each parameter in the MCMC algorithm in Table A-1.

A-6 Prior sensitivity analysis

A-6.1 Simulations based on single genealogy realizations

In Section 3.1, we put informative priors on the removal rate $\gamma$ and explore three different simulation scenarios. Although our LNA-based model successfully recovers the $R(t)$ dynamics and SIR trajectories, the posterior density of the removal rate is not too different from its prior in the SD and NM scenarios. In this section, we investigate sensitivity of our inferences to changes in the prior of the removal rate $\gamma$. For the same genealogies and parameter settings as in Section 3.1, we assign weakly informative priors to the removal rate $\gamma$:

CONST $R_{0}(t)$ scenario: $\gamma\sim\mbox{lognormal}(-1.7,0.25)$, 2. 2.

SD $R_{0}(t)$ scenario: $\gamma\sim\mbox{lognormal}(-1.7,0.25)$, 3. 3.

NM $R_{0}(t)$ scenario: $\gamma\sim\mbox{lognormal}(-1.2,0.25)$.

For each scenario, we fit a LNA-based model using 300,000 MCMC iterations. The first row in Figure A-10 shows the point-wise posterior medians and 95% BCIs for the basic reproduction number trajectories, $R_{0}(t)$. Our LNA-based method performs well in the CONST and SD scenario. However, for NM scenario, the method fails to fully capture the increase and decrease trend at the beginning and the end of the epidemic. The second row in Figure A-10 depicts the prior and posterior densities of the removal rate $\gamma$. The LNA-based method estimates the removal rate with good precision in the CONST scenario. However, for SD and NM scenario, the removal rate posterior densities are similar to the prior densities, but shift to the right from the truth. Posterior summaries of $S(t)$ and $I(t)$ are given in the third and fourth row of Figure A-10. The LNA-based method performs well in recovering the truth in the CONST and SD scenarios. In the NM scenario, the true trajectories are still covered by the wide BCIs, but the model seems to underestimate the $S(t)$ and overestimate $I(t)$.

A-6.2 Frequentist properties of posterior summaries

In this section, we repeat the simulation study in Section 3.2 with a weakly informative prior distribution on recovery: $\gamma\sim\mbox{lognormal}(-1.2,0.25)$. We fit LNA-based models to approximate the posterior distribtuion of parameters and latent variables for each genealogy, and compare that with the estimation in Section 3.2 with informative prior on $\gamma$ ($\gamma\sim\mbox{lognormal}(-1.2,0.1)$). To evaluate the performance, we use same metric defined in Section 3.2 and generate posterior summary boxplots in Figure A-11.

Sampling distribution boxplots of $R_{0}(t)$ posterior summaries are depicted in the left three plots of Figure A-11. The LNA-based model with informative prior (IP) on $\gamma$ yields significantly lower MAE and MCIW than that with weakly informative prior (WIP). Compared with IP, the LNA-based model with WIP has really poor envelope for the $R_{0}(t)$ trajectory.

Sampling distribution boxplot of $I(t)$ posterior summaries, shown in Figure A-11, are similar to $R_{0}(t)$ results, with IP yields significantly lower MRAE and lower MRCIW. Somewhat counter intuitively, the WIP cases end up with higher coverage for $I(t)$ trajectory. This is likely caused by the wide BCI under WIP.

We also report the absolute error (AE) and 95% BCI widths for the removal rate $\gamma$ in Figure A-11. Though the WIP prior still centered at the truth of $\gamma$, we can see really large absolute error in the removal rate estimation. The 95% BCI coverage for $\gamma$ under IP is 1. Though WIP yields wider BCIs, the coverage for $\gamma$ is only 0.65.