A Likelihood Approach for Real-Time Calibration of Stochastic Compartmental Epidemic Models

Stochastic transmission dynamic models are especially useful for studying the early emergence of novel pathogens given the importance of chance events when the number of infectious individuals is small. However, methods for parameter estimation and prediction for these types of stochastic models remain limited. In this manuscript, we describe a calibration and prediction framework for stochastic compartmental transmission models of epidemics. The proposed method, Multiple Shooting for Stochastic systems (MSS), applies a linear noise approximation to describe the size of the fluctuations, and uses each new surveillance observation to update the belief about the true epidemic state. Using simulated outbreaks of a novel viral pathogen, we evaluate the accuracy of MSS for real-time parameter estimation and prediction during epidemics. We assume that weekly counts for the number of new diagnosed cases are available and serve as an imperfect proxy of incidence. We show that MSS produces accurate estimates of key epidemic parameters (i.e. mean duration of infectiousness, R0, and Reff) and can provide an accurate estimate of the unobserved number of infectious individuals during the course of an epidemic. MSS also allows for accurate prediction of the number and timing of future hospitalizations and the overall attack rate. We compare the performance of MSS to three state-of-the-art benchmark methods: 1) a likelihood approximation with an assumption of independent Poisson observations; 2) a particle filtering method; and 3) an ensemble Kalman filter method. We find that MSS significantly outperforms each of these three benchmark methods in the majority of epidemic scenarios tested. In summary, MSS is a promising method that may improve on current approaches for calibration and prediction using stochastic models of epidemics.

3. Determine the event that just realized: Select the eventq such that where u 2 is another random number drawn from a uniform distribution U ((0, 1]). 4. Update epidemic state given the realized eventr. 5. Until a desired time condition is satisfied, increment time t ← t + τ and move to Step 1.  These SITR trajectories are used to calculate the number of new diagnoses in week i 721 which is equal to T i + R i − T i−1 − R i−1 . If for a simulated trajectory the conditions on 722 attack rate (30% -50% for the mild scenario, 50%-70% for the severe scenario and 723 70%-100% for the extreme scenario) and timing of peak (between 10 and 20) are 724 fulfilled, we keep the trajectory for the respective scenarios or we otherwise we reject it. 725

726
Equations of the SEITR model 727 The Susceptible -Exposed -Infected -Treatment -Recovered model contains one more 728 compartment than the SITR model, namely Exposed. Disease transmission can be 729 modeled using the following ODE model: where k 1 = k 2 = 2θ 1 and θ 1 is the disease transmission rate, θ 2 is the rate of seeking 731 treatment while infectious, and θ 3 is the rate of recovering.

732
In our analysis, we assume that at time t = 0 one population member becomes 733 infected. The model initial condition can therefore be defined as R 0 and duration of infectiousness: For R 0 or the mean duration of infectiousness, 741 we have f Mi|θ (m|θ) = 1, as R 0 and the mean duration of infectiousness only depend on 742 the parameter vector sampled from the posterior π i .

743
R eff : For R eff as the performance target, M i is defined as M i = R 0 Si N . As S is the only state that the target depends on, the function f Mi|θ (m|θ) needs to yield the probability for this state. This is accomplished by integrating over the remaining states: Infection prevalence: If the target is the infection prevalence: i , the function f Mi|θ (m|θ) needs to yield the probability for this state and, hence, is calculated by integrating over the remaining states.
As the I.Poi benchmark does not calculate a belief state for time t i , we use the 744 simulation model to carry out simulations with parameters according to π i (θ|Y i ) and 745 consider these simulations as a sample from the belief state Π i .

746
All Predictions: In case of prediction for one week, three weeks (cumulative and specific) or the attack rate, it holds f Mi|θ (m|θ) = P (m|Y i ) with P as in Eq. (9).

747
Median integrated relative error (mIRE): Next, we will calculate a median relative error (mIRE) over the 50 trajectories as  Poi. We analyzed the computational time for the severe scenario with a population size 752 of 10000. Carrying out the 1000 stochastic simulations for benchmark I.Poi for the 1000 753 parameters takes about 100 minutes. The additional time for evaluating the likelihood 754 approximation is rather short with about only 1 second, see table 1 for numbers.

755
The computational times are shorter before the peak and longer after the peak,

759
We note that the implementation was not carried out in a way to maximize speed, so 760 there is opportunity for additional speed up with more efficient implementation.

761
However, these numbers demonstrate that all the algorithms can be used on a personal 762 computer in real time. Computational time in seconds for the calibration of one time course in average for the severe scenario in the population with 10000 individuals. * The 100 min are the computational time to carried out 1000 stochastic simulations for each parameter with the Gillespie algorithm.
(b) Choose the set Θ that includes the values of parameters θ for which the approximate likelihood function (10) should be calculated.
(b) Update the likelihood function: (c) Update the parameter posterior distribution: (b) Update the likelihood function: (c) Update the parameter posterior distribution: (c) Choose an initial set of particles Ψ. Its elements ψ ∈ Ψ contain values for the parameters and states: ψ = (θ, ν).

Calibration: For each observation
(b) Calculate the prior for the observations: obs,i as in Eq. (13). (d) Calculate the prior variance of the observed quantity as σ 2 prior = Variance ψ (y) .
(e) Calculate the prior co-variance of the observed quantity and each unobserved component of the particle ψ as σ m = co-variance ψ (y) , ψ (m) , where m denotes the unobserved component. ❼ ODE solution x(t, x 0 ; θ) for an ODE x ′ integrated for time t with initial value x 0 and parameter θ such 773 as in equation (5).
% and minimize by solving derivative equals 0.