Data-driven prediction and origin identification of epidemics in population networks

Effective intervention strategies for epidemics rely on the identification of their origin and on the robustness of the predictions made by network disease models. We introduce a Bayesian uncertainty quantification framework to infer model parameters for a disease spreading on a network of communities from limited, noisy observations; the state-of-the-art computational framework compensates for the model complexity by exploiting massively parallel computing architectures. Using noisy, synthetic data, we show the potential of the approach to perform robust model fitting and additionally demonstrate that we can effectively identify the disease origin via Bayesian model selection. As disease-related data are increasingly available, the proposed framework has broad practical relevance for the prediction and management of epidemics.


Introduction
Robust prediction of the spread of an epidemic is critical to monitoring and halting its progress.The reliability of these predictions, which have high clinical and societal significance, hinges on the underlying mathematical models which quantify the spread and virulence of diseases.Several models have been proposed for predicting the spread of epidemics in real-world populations, allowing for the development of strategies for effectively managing disease spread via organized intervention.Perhaps the most well-known approach is Kermack and McKendrick's compartmental SIR model (and its extensions, such as SIRS and SEIR), a differential equation model which divides populations into groups corresponding to their relation with the disease (e.g., susceptible or recovered), which is widely studied due to its simplicity and predictiveness for several common diseases [1,2].More recent work has also incorporated the topological aspect of network structure by modeling explicit (or random) population networks through which diseases propagate [3,4,5,6], working toward a more holistic view of disease modeling which can include aspects such as demography, land use, and climate change [7].The predictions made by these mathematical models have been used to study a diverse set of historical and modern epidemics, including HIV [8,9], malaria [10,11,12], polio [13], and tuberculosis [14], by using a wide range of data assimilation techniques [15,16].
Many aspects of these models have also been analyzed from a more abstract mathematical perspective.Local bifurcation analysis has been performed on the man-environment-man and SIR models [17,18], while other work has used Lyapunov functions to determine endemic equilibria for SIRS and SEIR models [19,20] or has considered a mean-field approach [21,6,22,23].Recent work has addressed how the network structure influences the spread of disease via the initial conditions and network topologies [4,5] and the ways in which epidemics spread on random networks [3].Analytical results have also been obtained for the case of two competing (or promoting) diseases on a network [18,24].Moreover, many of the models in question have been used to design intervention policies or allocate vaccines via optimal control [20,25] or randomized interventions [26].
In this work, we introduce data-driven reverse engineering of models for the spread of an epidemic through a population network.The model structure and parameters are inferred from noisy observations using a Bayesian framework for uncertainty quantification (UQ); Bayesian inference enables robust predictions and the rational selection of the best among competing models using data-based evidence.At the same time, Bayesian inference involves sampling of the (potentially high-dimensional) parameter space, requiring repeated evaluations of the forward model.As such, in cases where the forward model is complex and computationally intensive (e.g., a large network with intricate connectivity), the Bayesian approach may be prohibitively expensive.Yet the Bayesian setting is of considerable practical interest given its potential applications to real-world data collection: efficient parameter estimation would enable calibration of the models with actual observations, and models could be compared based on their degree of fit.
In what follows, we apply our Bayesian uncertainty quantification framework Π4U to an extension of the SIR model to graphs.Π4U is an efficient parallel implementation of the transitional Markov chain Monte Carlo (TMCMC) algorithm, which offsets the complexity of UQ approaches by making use of modern parallel computation to run many copies of the model simultaneously [27,28].Using noisy, simulated data, we show that our method is able to efficiently estimate values of the model parameters and their underlying uncertainties.The Π4U framework is also convenient for Bayesian model selection, which we use to identify the origin of the epidemic (see, e.g., [29]) by considering each possible starting location as a separate model.Bayesian approaches thus show significant potential to aid in real-world epidemic modeling and mitigation.

SIR model
The SIR epidemic model decomposes a population into three eponymous groups: hosts who are susceptible to the disease, hosts who are infected and contagious (the infective group), and hosts who are neither susceptible nor infected, either via gained immunity from recovery or due to a vaccine, quarantine policies, or disease-related death (the removed group).

Single population model
Let S(t), I(t), and R(t) denote the size of the susceptible, infective, and removed groups, respectively, as functions of a continuous time t.The SIR model is based on three main assumptions: first, since the timescale on which the disease evolves is assumed to be much shorter than the timescale on which the population may evolve via, e.g., births or natural deaths, the population Y is assumed constant, and so for all t (note that individuals killed by the disease are considered part of the removed group).Second, members of the population are assumed to come into contact uniformly at random and at a constant rate β -this parameter governs the rate at which an infection can spread.Finally, the infective population recovers (or is otherwise removed from the infectives via, e.g., death) at a constant rate γ.These assumptions can thus be visualized as yielding the following set of ordinary differential equations: Namely, at a particular time t, S(t) susceptibles and I(t) infectives come into contact at a rate β, yielding βIS transitions from susceptible to infective (implicitly assuming that contact with an infective immediately infects a susceptible -if this assumption is not desired, the chance of disease transfer can be incorporated in β).Meanwhile, I(t) infectives are removed at a rate γ, yielding γI transitions from infective to removed.

Epidemic model on graphs
The SIR model is readily generalized to a directed graph with N vertices, a mathematical construct which can be thought of as modeling a collection of N distinct communities.Namely, let each node (community) be a distinct population whose dynamics evolve according to Eq. 3; the directed edges (connections between communities) are a convenient framework to dictate transfer between populations.Since each population itself has three groups (susceptible, infective, removed), three quantities are needed to describe movement.Here, we use λ i,j , η i,j , and g i,j to describe the rate of movement from node i to node j on the susceptible, infective, and removed groups, respectively; identifying each transition rate as the weight of the edge connecting i to j, these rates are naturally written as weighted adjacency matrices, here denoted Λ, H, and G.The SIR model on a network, now a system of N models corresponding to each population i, can then be written as or more succinctly in matrix form: (Here, F = (1, 1, . . ., 1) T is a vector of ones which simplifies the notation.)It should be emphasized that S, I, and R are N ×1 vectors whose i th element corresponds to the i th population.Note that if λ i,j = η i,j = g i,j = 0 for all i, j, i.e., there is no movement between populations, each model reduces to the single population model (Eq.3).

Bayesian methodology
The network model described by Eq. 5 is a predictive model for tracking the spread of an epidemic through a population network.Here, we introduce a Bayesian approach to the inverse problem, i.e., reverse-engineering aspects of the model itself using observed outputs.In real-world scenarios, these observed outputs (e.g., the number of infected patients at a particular set of community health centers) are noisy due both to observational noise (e.g., not all infections are reported) and to model error -epidemic models are mathematical equations introduced to represent the real system, and so will not exactly predict the noise-free measurements.
The inverse problem for epidemics is of considerable practical interest: accurate estimation of model parameters would allow for the identification of a disease's underlying characteristics (e.g., its infectivity β or the host recovery rate γ).Moreover, by defining a class of models corresponding to the disease having originated in different communities, Bayesian model selection could be used to probabilistically determine the initial outbreak location [29,30,31], thereby aiding in identification and mitigation of the vector of infection.We note that stochastic optimization techniques such as CMA-ES [32] may also be used to infer optimal model parameters from noisy data [33]; however, such techniques neither enable robust predictions nor provide a framework for model selection as does the following Bayesian framework.

Bayesian uncertainty quantification
Denote as θ ∈ R n the set of parameters corresponding to the model M of interest.Here, the model M is given by the SIR network model (Eq.5); its parameters include the infectivity β and recovery rate γ.By evolving the system forward in time, we can generate deterministic predictions for the system at a future time T -for example, the number of infected individuals present at a certain subset of nodes.
We first consider the problem of parameter estimation: suppose we observe a subset of noisy predictions from this model and wish to estimate the parameters θ ∈ R n which generated them.In particular, we will assume the observed data D ∈ R m obey the model-prediction equation where g(θ|M ) : R n → R m denotes the deterministic mapping of parameters to outputs and e is an additive error term.The posterior distribution of the parameters given the data is then given by Bayes' Theorem as p(θ|D, M ) = p(D|θ, M )π(θ|M ) ρ(D|M ) (7) in terms of the prior π(θ|M ), likelihood p(D|θ, M ), and evidence ρ(D|M ) of the model class, given by the multi-dimensional integral This scenario can be extended to the case where the model M is one of many models in a parameterized class M; the probability that the observed data were generated by any particular model M i is also given by Bayes' Theorem: In particular, under the assumption of a uniform prior on models, P r(M i |D) is directly proportional to the evidence ρ(D|M i ), and so model selection is "free" when the evidence is already calculated for parameter estimation [34,35].
In order to calculate the likelihood p(D|θ, M ) needed for Eq. 7, we need to postulate a probability model for the error term e.Here, we assume the model error e is normally distributed with zero mean and covariance matrix Σ.Assuming that errors at different nodes are uncorrelated, the covariance matrix becomes Σ = σI, where I is the m × m identity matrix.
If e is Gaussian, it follows that D is also Gaussian, and so the likelihood p(D|θ, M ) of the observed data is given as where is the weighted measure of fit between the model predictions and the measured data, | • | denotes determinant, and the parameter set θ is augmented to include parameters that are involved in the structure of the covariance matrix Σ (here, the noise level σ).
The main computational barrier in calculating the posterior distribution of parameters given by Eq. 7 is the complex forward problem g (the epidemic network model) which appears in the fitness J(θ, D|M ).The Π4U software [28] has two advantages in this respect: first, it approximately samples the posterior via transitional Markov chain Monte Carlo [36], described below, which is massively parallelizable; and second, it leverages an efficient parallel architecture for task sharing (see Appendix).

7:
Resample based on samples available in stage j in order to generate samples for stage j + 1 and compute likelihood p(D|θ j+1,k , M ) for each.

8:
if q j+1 > 1 then goto loop. 13: Recent literature suggests that q j+1 , which determines how smoothly the intermediate distributions transition to the posterior, should be taken to make the covariance of the plausibility weights at stage j smaller than a tolerance covariance value, often 1.0 [36,28].
Next, the algorithm calculates the average S j of the plausibility weights, the normalized plausibility weights w(θ j,k ), and the scaled covariance Σ j of the samples θ j,k , which is used to produce the next generation of samples θ j+1,k : Σ j is calculated using the sample mean µ j of the samples and a scaling factor b, usually 0.2 [36,28].
The algorithm then generates N j+1 samples θj+1,k by randomly selecting from the previous generations of samples {θ j,k } such that θj+1, = θ j,k with probability w(θ j,k ).These samples are selected independently at random, so any parameter can be selected multiple times -call n j+1,k the number of times θ j,k is selected.Each unique sample is used as the starting point of an independent Markov chain of length n j+1,k generated using the Metropolis algorithm with target distribution f j and a Gaussian proposal distribution with covariance Σ j centered at the current value.
Finally, the samples θ j+1,k are generated for the Markov chains, with n j+1,k samples are drawn from the chain starting at θ j,k , yielding N j+1 total samples.Then the algorithm either moves forward to generation j + 1 or terminates if q j+1 > 1.

Results
In the following results, we apply our high-performance implementation of Bayesian uncertainty quantification Π4U to a case study with two example network structures.In each case, we fix particular values of the system parameters (infectivity β, recovery rate γ, and time of observation T ) and use Eq. 5 to evolve the network forward in time via a 4th-order Runge-Kutta method.At the observation time T , the infective populations (and sometimes also the recovered populations) from a selected subset of communities, corrupted by additive Gaussian noise with noise level σ, are output as the noisy observed data.Namely, observed data D k are generated as where each deterministic population datum p k generated from the reference model is added to a zero-mean, unit-variance Gaussian k scaled by the noise level σ to yield the observed noisy datum D k .In order that the signal-to-noise ratio be high enough for meaningful estimation, we choose σ to be a fraction σ = 0.01α (or sometimes σ = 0.05α) of the standard deviation α of all model outputs p k .
We then use our method to approximately solve the inference problem by generating 10 4 samples from the posterior distribution of the model parameters given these noisy outputs, checking the validity of our approach by comparing the resulting distributions to the known reference values.While we use synthetic, model-generated noisy data, we note that the framework is readily extended to the incorporation of real-world data.
For ease of comparison, we present numerical results in terms of the rescaled parameters (θ β , θ γ , θ T , θ σ ), given by θ β = β/β 0 , i.e., the ratio between the estimated value and the true reference value.Accurate estimation thus results in scaled parameters close to 1.As the Π4U approach does not rely on the choice of a particular prior, we use a simple uniform distribution on [0.01, 2] × [0.5, 2] × [0.02, 5] × [0.5, 10] in this scaled parameter space.

Network 1: the 20-barbell graph
We first consider a network with two distinct populations, each with many highly interacting subcommunities.The two populations mix via a single route, modeled by a single connecting edge.This "barbell graph" -two complete 20-node graphs connected by a single edge -is illustrated in Fig 1 .In this case, we impose uniform transition rates between adjacent vertices of 0.02, 0.3, and 0.05 for the susceptible, infective, and recovered populations, respectively.The infection begins at node 1 with the configuration S 1 (0) = 5, I 1 (0) = 95, R 1 (0) = 0, and all other nodes are fully susceptible with configuration S i (0) = 100, We consider in particular the case of having information only from a limited subset of nodes; by placing the "sensors" at different locations (i.e., observing different subsets of nodes), we can test how the sensor configuration influences the parameter estimation procedure and the corresponding uncertainties.In particular, we assume observations of both the infective and recovered populations at the sensor locations.β and γ are positively correlated, i.e., similar outputs can be achieved by simultaneously raising both the infection rate and the recovery rate.Intuitively, a faster-spreading disease must be counteracted by quicker recovery in order for the dynamics to remain consistent.Similarly, both β and γ are negatively correlated with T ; a more infectious disease or quicker recovery would increase the speed of the system dynamics, meaning similar outputs would be observed earlier.
The recovered noise standard deviation σ has significantly larger uncertainty than do the system parameters β, γ, and T .Despite this comparatively large uncertainty, the reference noise value σ is recovered to within two standard deviations for all three sensor configurations, though the posterior means θ σ consistently overestimate the true magnitude.Uncertainties in other parameters depend strongly on sensor location; observing nodes on opposite sides of the graph yields much smaller uncertainties in β, γ, and T than observing nodes on the same side of the graph.Numerical values for T = 3 and T = 7 appear in the first and third block of Table 1.For both observation times, all sensor configurations recover β, γ, and T to within one standard deviation.The experiment with sensors at nodes 3 and 24 continues to be both the most accurate (in terms of posterior means) and precise (in terms of uncertainties), highlighting the importance of sensor placement in leveraging information from different timescales.
To test the robustness of the parameter estimation to increased noise, we next reconsider the time T = 5 case with noise level σ = 0.05α, five times the previously-used σ = 0.01α.Results appear in the right column of Fig 4 and in the fourth block of Table 1.Again, the experiment using simulated data from nodes 3 and 24 (Fig 4E ) recovers the parameters with comparatively lower uncertainty and greater accuracy.For all three sensor configurations, parameters are recovered to within one standard deviation, and so we conclude that the Bayesian uncertainty quantification approach to SIR models on the 20-barbell graph has significant robustness to observational noise.
A final experiment tests the approach in the context of model misspecification.We perturb the network by introducing an extra edge between nodes 2 and 40; while observations are generated from a network which includes this edge, we perform parameter estimation using the original model (without the extra edge).Results for observing nodes 3 and 24 for the T = 3, σ = 0.01α case with two different perturbation strengths appear in Figure 5 and Table 2.As compared to the original case (second row of Table 1 and Fig 2E), the most notable change is the significant increase in the estimated noise level -disagreements between the original model (now misspecified) and the observed data are reconciled by assuming a much higher magnitude of noise.Nonetheless, estimation of system parameters (β, γ, T ) for both perturbations is reasonably successful, with all reference values recovered to within one standard deviation, though uncertainties in estimation are 3 -4 times larger than without the perturbation.

Network 2: the three-group network
The second network considered is a 44-node graph comprising three large sub-networks with limited interaction (Fig 6).Each sub-network has a distinct topological structure and set of nonuniform transition rates (explicit values appear in the Appendix).We again consider three sensor config-     3. Sensors at the 7-node subset recover β, γ, and T to within one standard deviation, while larger subsets recover all parameters with comparable accuracy and much greater precision: the 7-node 95% confidence interval for the scaled infection rate θ β is [0.8627, 1.1937], which narrows significantly to [1.0020, 1.0256] in the 23-node case.Increasing the number of observed nodes has the additional effect of accurately recovering the scaled noise level θ σ , which is significantly overestimated when observing only the 7-node set.
Many real epidemic data sets include multiple observations over time; it is worth verifying that our approach can reasonably incorporate such observations.T (time from infection to first sample) remains the only unknown time in this setting, as the relative timing between samples is known.Blocks 2-4 of Table 3, the right column of .Comparatively, additional samples beyond the second provide only marginal benefit (when observing the 7-node set, moving from 2 to 4 samples reduces uncertainties for β, γ, T , and σ by factors of 1.35, 1.99, 1.24, and 1.69, respectively).We additionally note that observing Experiment Observed Set   Observed Set As with the 20-barbell graph, we also consider model misspecification via perturbation of the network.The final row of Table 3 presents parameter estimation results when observations at the 35-node subset are generated using an additional edge connecting nodes 16 and 44 (and thus Groups I and III); transition rates along the additional edge are chosen to match those of the edges connecting Groups II and III.As before, uncertainties for system parameter estimates (β, γ, and T ) are larger than without the perturbation (by factors of 1.33, 1.39, and 1.33, respectively), though recovered parameter values are nonetheless accurate (scaled values off by less than 0.03), and so parameter estimation in this setting remains successful.σ is again overestimated, as it must additionally account for disagreement between the model used to generate the data and the model used to perform parameter estimation.
We next augment the parameter set θ with the initial population vector S 0 , I 0 , and R 0 of the initially infected node, but take as known the observation time T = 5.In order that the reference values of all parameters be positive, the initial population vector at node 34 is altered to S 34 (0) = 5, I 34 (0) = 90, R 34 (0) = 5.The scaled parameter set (θ β , θ γ , θ S 0 , θ I 0 , θ R 0 , θ σ ) uses a uniform prior on  4. Compared to estimation of β, γ, and T , correlations between parameters are generally weaker in this context, though there do exist clear relationships (e.g., larger I 0 necessitates smaller β for the infection to spread at the same absolute rate).All sensor configurations recover the disease parameters β and γ accurately, with scaled values off by less than 10 −2 in all but one case.Conversely, uncertainty in the recovered distributions for S 0 and R 0 is significantly higher than for other parameters, owing to the small ground-truth values for these parameters relative to the total population at the origin (S 0 = R 0 = 5 out of 100 individuals at time t = 0) and their correspondingly minor effect on the behavior of the initial outbreak.

Origin of disease identification
Finally, we introduce a method for probabilistically identifying the origin of the epidemic (see, e.g., [29]) using the Bayesian model selection framework described in the Bayesian methodology section; recall that all observations are at the future time T , and so the origin may not be identified with certainty even when included in the set of observed nodes.We initialize the disease at node 1 with the standard initial configuration S 1 (0) = 5, I 1 (0) = 95 and R 1 (0) = 0, with all other nodes fully  Model Log Evidence P r(M j |D) susceptible with S i (0) = 100, I i (0) = R i (0) = 0, and use corrupted observations of infective and recovered populations from the 35-node subset.Other parameter values are identical to those used in the previous section; in this case, T = 5 is assumed to be known, with β and γ estimated from observations.Defining the model class M j as the model under which the disease originated from node j with the given initial vector, the log evidence for each model is generated from the model selection framework (see Eq. 9 in Methods).
A representative selection of results appear in Table 5. Model M 1 , corresponding to the correct origin of the disease at node 1, has a significantly larger log evidence than all other models considered.Models which place the origin at increasingly distant points generate increasingly less accurate and more uncertain results; models M 32 and M 43 , which originate the disease in Group III, find the estimated noise σ to be two orders of magnitude larger than the reference value.Out of the models shown, if the correct model M 1 is not considered, the probabilities shift to 0.003 for M 8 , 0.997 for M 14 , and all other probabilities ≈ 10 −13 or smaller, suggesting that topographic proximity to the true origin is the dominant factor in the evidence.
The effect of topographic proximity on the model evidence appears visually in Fig 9 .Models M i , i = 1, . . ., 44 correspond to the epidemic beginning at node i; each node in the graph is colored by the estimated log evidence log P r(M i |D), where D are the noisy data obtained from a reference simulation beginning at node 1.In order that the fine detail be more visible, the color mapping is nonuniform such that node 1 (log evidence −8.15) is the only node in its color bin; other bins from log evidence −360 to −250 capture the range of behavior in the remaining nodes.Evidence decays with topographic distance within the first group and becomes negligible in the second and third groups, highlighting the improbability of the epidemic starting at these distant nodes and producing the noisy observations D corresponding to an epidemic outbreak at node 1.
Lastly, we attempt origin identification via the same model selection approach for the perturbed model which generates data using an additional edge connecting nodes 16 and 44; results appear in Table 6.Despite the misspecification, the correct model M 1 is selected with near certainty, though its log evidence is many orders of magnitude smaller than in the no-perturbation experiment (−193.3884 vs. −8.1478)owing to the high level of noise required to explain disagreements with the observed data.Models M 32 and M 43 , both located in Group III, are also subject to significant drops in log evidence.A second notable difference is the reduced level of estimated noise for incorrect models M j -the additional edge allows the disease to quickly reach nodes which are distant under the original model, thereby requiring less noise to explain discrepancies with the observed data.Overall, selection among origin models appears robust to perturbation.As a specific example, the TMCMC method within Π4U is able to achieve an overall parallel efficiency of more than 90% on 1024 compute nodes of Swiss supercomputer Piz Daint running hybrid MPI+GPU molecular simulation codes with highly variable time-to-solution between simulations with different interaction parameters.

Transition matrices for three-group network
Transition rates were chosen to vary among groups in the three-group network in order to test the robustness of our approach to nonuniform rates.Populations moved between Group I and Group II (via the edge connecting node 14 to node 21) at a rate of 0.4, while the transition rate between Group II and Group III (via the edges connecting node 31 to nodes 25 and 26) was 0.2.Rates (λ i,j , η i,j , g i,j ) within Group I were selected randomly to be either (0,2,0.05,0.1)or (0.1,0.1,0.2),those for Group II were selected randomly to be either (0.15,0.2,0.1) or (0.3,0.1,0.1), and Group III had uniform rates of 0.05.

Figure 1 :
Figure 1: The 20-barbell graph.Two complete 20-node graphs are connected by a single edge.

Fig 3
shows the deterministic populations of the susceptible, infected, and recovered groups as a function of time for a selection of nodes involved in the experiments (since, e.g., node 12 is identical to node 3, only one is shown).Nodes 21 and 27, both contained in the initially susceptible subgraph, reach peak infective population around time t ≈ 5. Nodes on the side of the infection origin, conversely, achieve peak infective population at t ≈ 3. The results of parameter estimation at the reference observation time T = 5 thus suggest that observing nodes around the time when the infective population peaks improves the accuracy of the recovered parameters.This effect can also be seen in Fig 2 in the joint marginals of the system parameters β, γ, and T : nodes 3 and 12, which peak at t ≈ 3, give much more smeared marginals in Fig 2A, obtained using observations at (A) Sensors at nodes 3 and 12, T = 5 (B) Sensors at nodes 3 and 24, T = 5 (C) Sensors at nodes 24 and 27, T = 5 (D) Sensors at nodes 3 and 12, T = 3 (E) Sensors at nodes 3 and 24, T = 3 (F) Sensors at nodes 24 and 27, T = 3

Figure 2 :
Figure 2: Parameter estimation results for the 20-barbell graph with infection rate β = 0.02, recovery rate γ = 0.3, noise level σ = 0.01α, and time step ∆t = 0.005.In each experiment, noisy data from two nodes were used to track the epidemic.For each pair of nodes, the experiment was run for time of observation T = 5 (A, B, C), T = 3 (D, E, F), and T = 7 (shown in Fig 4).Histograms for each parameter are displayed along the main diagonal of the figure.Subfigures below the diagonal show the marginal joint density functions for each pair of parameters, while subfigures above the diagonal show the samples used in the final stage of TMCMC.Colors correspond to probabilities, with yellow likely and blue unlikely.

(A )
Susceptible population (B) Infective population (C) Removed population (D) Total population

Figure 3 :
Figure 3: Time evolution of the susceptible, infective, recovered, and total populations at nodes 3, 20, 21, and 27 of the 20-barbell graph.Nodes from different complete subgraphs have different trends and peak times.

time T = 5 ,
than in Fig 2D, corresponding to observations at T = 3; meanwhile nodes 24 and 27, which peak at t ≈ 5, have sharper distributions in Fig 2C, with observations taken at T = 5, than in Fig 2F, which has observations at T = 3.The parameter distributions obtained for data from nodes 3 and 24, shown in Fig 2B and Fig 2E, are very similar in both cases, suggesting that placing one sensor in each subgraph leverages information from both timescales.The parameter estimation results at T = 7, shown in the left column of Fig 4, corroborate this conclusion.Though all nodes in the graph are well past peak infective population at this time, using information on two different timescales (the two subgraphs) yields much sharper joint marginals (see, e.g., the joint distribution of β and T in Fig 4B as compared to Fig 4A and Fig 4C).

Figure 5 :
Figure 5: Parameter estimation results for the perturbed 20-barbell graph with infection rate β = 0.02, recovery rate γ = 0.3, time step ∆t = 0.005, and noise level σ = 0.01α at time T = 3 using observations from nodes 3 and 24.Data are generated with an additional edge between nodes 2 and 40, which uses (A) 10% or (B) 50% of the transition rate of other edges; parameter estimation uses the original model.See Fig 2 for description of subfigures.

Fig 7 ,
and the left column of Fig 8 present results when considering additional samples taken at evenly spaced intervals of length 1.Including a second observation yields a significant reduction in estimated uncertainties, especially when observing only 7 nodes (uncertainties for β, γ, T , and σ are reduced by factors of 23.51, 17.22, 19.89, and 9.70, respectively)

Figure 7 :
Figure 7: Parameter estimation results for the three group network.(A, B, C) Parameter estimation for β, γ, T based on information from 7-node, 23-node and 35-node subsets, respectively.Reference values are β = 0.02, γ = 0.3, and T = 5, with time step ∆t = 0.0005.(D, E, F) Same experiment using two observations from each sensor separated by an interval t = 1.See Fig 2 for description of subfigures.

Figure 9 :
Figure 9: Origin identification model evidence by node, three-group network.Color map of model selection log evidence for models initializing the epidemic at each node in the three-group network using noisy data from an epidemic simulation beginning at node 1. Nonuniform color mapping (right label) emphasizes differences among incorrect models (log evidence −360 to −250); log evidence for the correct origin (node 1) is notably larger (log evidence −8.15).

Figure 10 :
Figure 10: Task graph of the TMCMC algorithm (left) and parallel architecture of the TORC library (right).

Table 1 :
Numerical results for estimation of β, γ, T and σ for the 20-barbell graph.Blocks of the table are different times of observation or noise levels.Parameters are reported in terms of scaled values θ; accurate estimation thus results in values close to 1. Uncertainties, e.g., u β , are the ratio of a parameter's standard deviation to its mean.

Table 2 :
Numerical results for estimation of β, γ, T and σ for the perturbed 20-barbell graph.Data are generated with an additional edge between nodes 2 and 40, which uses the specified fraction of the transition rate of other edges; parameter estimation uses the original model.See Table1for description of values.

Table 4 :
Numerical results for estimation of β, γ, S 0 , I 0 , and R 0 for the three group network with σ = 0.01α.See Table1for description of values.multiple samples over time can reduce the correlation between parameter uncertainties (see, e.g., the highly-correlated posterior of β and γ in Fig 7A as compared to the nearly uncorrelated ellipse of Fig 8A).

Table 5 :
Subset of model selection results for the three-group network.Recovered scaled parameters θ and uncertainties u appear with the estimated log evidence for the model.

Table 6 :
). Subset of model selection results for the perturbed three-group network.Data are generated with an additional edge between nodes 16 and 44; parameter estimation uses the original model.Recovered scaled parameters θ and uncertainties u appear with the estimated log evidence for the model.