A large-scale immuno-epidemiological simulation of influenza A epidemics

Background

Agent based models (ABM) are useful to explore population-level scenarios of disease spread and containment, but typically characterize infected individuals using simplified models of infection and symptoms dynamics. Adding more realistic models of individual infections and symptoms may help to create more realistic population level epidemic dynamics.

Methods

Using an equation-based, host-level mathematical model of influenza A virus infection, we develop a function that expresses the dependence of infectivity and symptoms of an infected individual on initial viral load, age, and viral strain phenotype. We incorporate this response function in a population-scale agent-based model of influenza A epidemic to create a hybrid multiscale modeling framework that reflects both population dynamics and individualized host response to infection.

Results

At the host level, we estimate parameter ranges using experimental data of H1N1 viral titers and symptoms measured in humans. By linearization of symptoms responses of the host-level model we obtain a map of the parameters of the model that characterizes clinical phenotypes of influenza infection and immune response variability over the population. At the population-level model, we analyze the effect of individualizing viral response in agent-based model by simulating epidemics across Allegheny County, Pennsylvania under both age-specific and age-independent severity assumptions.

Conclusions

We present a framework for multi-scale simulations of influenza epidemics that enables the study of population-level effects of individual differences in infections and symptoms, with minimal additional computational cost compared to the existing population-level simulations.

Models of infectious disease epidemics, such as due to influenza A virus, have proven useful for understanding dynamics of disease transmission, vaccination strategies, and are frequently used as a decision support tool by public health officials. Single scale models of influenza A virus have been developed both on the host-level and on the population level in order to study epidemiological, immune response and viral characteristics of influenza A infection [1].

On the population scale, continuous models based on SIR (susceptible/infected/removed) approaches have been used for large-scale epidemiological predictions [2]. Continuous SIR models assume homogeneous populations and any model assumption about heterogeneity (such as human characteristics, network structure, and environmental factors) increases the number of equations and parameters, necessitating changing the model structure for any additional assumption. An alternative discrete approach to the study of disease epidemics is through agent-based models (ABM), where individuals are represented as autonomous agents whose infectious status is followed in time. Population density, age-structure, travel patterns and inter-individual contact patterns are derived or inferred from available census and demographic data. The infection propagates in the population according to disease transmission and duration rules, where inter-individual variability is represented in the form of statistical distributions. Because ABMs evolve individual agents, simulations require considerable computational cost. Many large-scale collaborative networks have developed and use large platform ABMs to make epidemiological predictions [3–6] and the impact of a variety of mitigation strategies. ABMs have been used in studies of school closure strategies [7, 8], in determining the role of subway travel in an epidemic [9], in the role of presenteeism on disease transmission in the workplace [10], and in vaccine allocation [11, 12].

Immuno-epidemiological models explore how host-level immune response affects population-level epidemiological patterns [13] and allow much finer tuning of host-behavior, compared to population-level models. The nested modeling approach embeds a mechanistic model of host-pathogen disease dynamics into an epidemiological model of infection by linking epidemiological parameters such as transmission rate or duration of infectiousness [14–16]. Immuno-epidemiological models consider assumptions about acquired immunity to study spread of infection [17–19] or host-parasite co-evolution [14, 20]. To account for heterogeneity in immune response, SIR models may be refined to incorporate individual host response through different linking mechanisms. Understanding the biological mechanisms underlying these variations in host response is important for estimating disease parameters and for designing optimal prevention strategies [21, 22].

Symptoms of influenza infection include fever, sore throat, chills, and cough, and have been attributed to the immune response to viral infection [23–26]. Clinical markers related to disease and disease-severity, largely due to host-immunity factors and prior exposures, are referred to as clinical phenotypes. Mathematical models of within-host influenza A infection have been used in estimation of kinematic parameters such as infection and clearance rates of the virus, as well as in prediction of dynamics and drug therapy strategies for individual infection [27, 28], yet have not been used to provide a mechanistic basis for varying disease phenotypes. Immuno-epidemiological modeling could be used to account for variability in clinical phenotype and evaluate the impact of symptom-based, biologically-based, and potentially behavior-based mitigation strategies at the population level.

A key challenge in deriving a realistic description of within-host responses is lack of data at the host-level. Human data is available from volunteer challenge studies where volunteers are inoculated intranasally, tend to be young (ages 18-35), healthy, and pre-screened for existing immunity to the experimental IAV strain. For H3N2 it has been shown that illness following intranasal inoculation is milder than illness following aerosol inoculation, manifested by shorter duration of cough and fever [29]. Mechanisms of transmission of naturally acquired infection include contact transmission, droplet spray transmission, and aerosol transmission. Therefore, using a host-level model to predict individual response to influenza infection requires altering model parameters based on healthy volunteer data in a biologically meaningful way.

Morbidity and mortality resulting from influenza infection is highest among the elderly (>65) and children. In a study of pandemic H1N1 in Ontario Canada, recovery was faster among patients under 18 years of age [30]. To gain qualitative insight of how host-level dynamics are varied as a function of age, animal studies prove to be useful. A study between aged and adult mice infected with sublethal doses of influenza virus (A/Puerto Rico/8/1934) showed slower recovery and delay in immune system activation and virus clearance in aged mice [31]. Differences in mouse response as a function of age provide insight in how to vary baseline responses in a within-host model across different age groups, while synthetic population data sets provide reasonable age distributions across a population [32, 33].

Naive integration of a within-host model in a populationlevel ABM requires the evaluation of a nonlinear system of differential equations for each individual infected, which is computationally prohibitive for realistically large simulations. In this manuscript, we present a prototypical hybrid immuno-epidemiological ABM model linking an equation-based within-host model and providing mechanistically-based host variability to an agent-based population level model. For outputs of the within-host model we choose the infectivity and symptoms scores during the time course of infection, which reflect likelihood of transmission and stay-at-home behavior of an infected individual. In this fashion, we can replicate existing results from simulated epidemics, but we introduce additional flexibility in studying variation of host-response at no additional computational cost. To eliminate computational cost, we replace the exact response with an approximation obtained by analyzing the most significant response modes of the within-host model. Further, the within-host model is simple to program and implement, and can be used to represent within-host dynamics in a variety of population-level platforms. We apply the hybrid ABM model to study questions arising from individual heterogeneity that cannot be addressed using a single scale model. In particular, we examine the role of varying severity of disease as a function of age on the dynamics of simulated epidemics. We find that epidemiological estimates such as attack rate and incidence are reduced when differences in host response across age groups are considered. Use of the immuno-epidemiological ABM model opens up many possibilities regarding the realistic impact of influenza on populations with different contact patterns or on mitigation strategies which depend on host phenotype or behavior.

Human volunteer data

The data set consists of viral load and symptom data from five previously published human volunteer studies conducted between 1995 and 1999 for 84 total individuals infected with A/Texas/36/91 (H1N1) [23, 24, 34–36]. Of these studies, only one [35] (n =17) does not contain symptom data. The viral titers across the five studies are all reported as averages across individuals in the respective studies in units of l o g ₁₀ TCID ₅₀/mL. To adjust the data across studies, we used mixed-effects modeling to compute a vertical shift for each study (random-effect), assuming that the distance between each curve is minimized. We computed weighted means and standard deviations (fixed-effect) on the shifted data.

In all studies, symptom scores were measured twice a day, and averaged. Across studies, 65 of the 67 individuals with symptoms data experienced at least one symptom. While similar symptoms were scored, their measurements were on different scales across studies. As done in [37], we scaled each study curve to its maximum clinical score. The summary curves for both titer and symptoms were calculated as a weighted average with weights as the number of individuals considered in each study. Data is summarized in the Additional file 1.

Overview of the immuno-epidemiological ABM model

We describe the process in which we link an equation-based model of within-host influenza A dynamics to a population scale ABM. A diagram of the procedure is shown in Figure 1. A collection of agents representing the population interacts within the simulation environment, transmitting infection in the process. Knowledge of the within-host immune response is used to determine the response of each agent to infection. Within each host, the immune response is reduced using principal component analysis (PCA) to identify significant directions in parameter space that describe most of the variability in immune response. Response surface analysis (RSA) captures the relationship between input (now simplified to two principal components) and model response (infectivity and symptoms) as an algebraic equation. Randomization over the principal components is used to generate population response variability. With this framework in place, we simulate disease outbreaks to see how changes at the within-host level affect a simulated epidemic in Allegheny County, PA.

Schematic diagram of the hybrid model. (a) Naive combination of equation based model of virus and symptoms dynamics of influenza A infection with a population level model. Clinical phenotypes of infection are input into the host-level model, which is integrated numerically and returns symptoms and infectivity of an infected individual to use for simulation of an epidemic. (b) In the hybrid ABM model response surface replaces host-level ODE model.

Full size image

Estimation of parameters for virus/symptoms model

We enhanced a target-cell limited ordinary differential equation (ODE) model first proposed by Baccam et al.[38] to model viral and symptoms dynamics within a host with an additional term to incorporate symptoms dynamics:

A schematic of the model is shown in the left panel of Figure 1(a). An individual is infected by an initial viral load V(0)=V ₀, which infects target cells (H) of the respiratory tract with infection rate β. Virus replicates within latent phase cells (J), which after average rate of time k become infectious (I) and shed virus at a rate p. Infected cells die at rate δ, primarily through cell necrosis, which triggers productions of cytokines such as type I interferons (IFN), interleukin(IL)-6, tumor necrosis factor (TNF)- α and IL-1 [39]. The novel contribution of this model is the inclusion of a variable Ψ quantifying the intensity of systemic symptoms of influenza, such as malaise and fever. Of all the variables of the model, it is most reasonable to have systemic symptoms dependent on infected cells, since both flu-like symptoms and severity of systemic symptoms have been shown to be correlated to cytokine levels [40, 41], particularly IL-6 and IFN- α/β[23, 24, 42], and inflammatory cytokines are, in turn, produced by infected cells. It has also been noticed experimentally that the number of days of viral shedding is correlated with cytokine levels, severity of symptoms, and mucus weight [43, 44]. While resolution of symptoms is a complex process, we assume a constant decay rate a. Model trajectories include viral load and symptoms scores, which are used as determinants of infectivity and symptomaticity in the population level simulation and represent state variables corresponding to data. We use Bayesian inference to construct an ensemble of models consistent with the data in order to explore the parameter space associated with our model and data. We sample the Bayesian posterior density P(α|D) using a Markov Chain Monte Carlo (MCMC) method with parallel tempering [45–48]. Ensemble modeling results in a collection of trajectories consistent with the data, in which each trajectory corresponds to a parameter set from the posterior distribution.

Abstraction of the within-host model for population-level evaluation

The ensemble of model parameters α sampled from the posterior distribution P(α|D) is used to determine how the immune response depends on model parameters, and how the response varies over the population [49, 50]. As large ranges of parameter values may all fit the data, we seek to identify directions in parameter space in which the model behavior is tightly constrained in order to make predictions. The ensemble of parameters can be characterized using the covariance matrix Θ=c o v(α), or the Hessian matrix of the cost manifold at a global minimum in parameter space; these two matrices are related as H≈Θ ^-1. A singular value decomposition, or PCA, of this approximated Hessian provides principal directions in parameter space corresponding to maximum changes in the model with respect to the cost function, here represented by total deviance from empirical virus and symptom values at each day during infection. We use the principal directions to vary parameter values for model evaluation, and define the hyper-plane P=s p a n{d ₁,d ₂} as the plane in parameter space spanned by the two leading principal directions. The computational domain is a discretization of the box [-0.5,0.5]× [-0.5,0.5]⊂P into a 20 by 20 grid. We select a random sub-sample of 1,000 parameter sets from the ensemble, which characterizes the “origin” of the hyper-plane. For each point on the grid, parameters are deviated from baseline along two principal directions by a distance specified by the grid coordinates. The deviated values are simulated in the model (Equations 1-5) yielding viral loads and symptoms responses for 11 days. The simulated median surface is then parameterized by a quadratic polynomial, which has the form

Population level model

The within-host representation is linked to the Framework for Replication of Epidemiological Dynamics (FRED) ABM software, an open source population-level large scale modeling system developed by the University of Pittsburgh Models of infectious Disease Agents Study (MIDAS) Center of Excellence. We use data from the synthetic population database, which is freely available and based on 2005-2009 U.S. census data [32] from Allegheny County, Pennsylvania. The data from the 1,164,879 agents assigns characteristics and behaviors such as occupation and workplace or school [32, 33] to each individual. Each household is assigned a latitude/longitude coordinate, income, size, and the sex, race, and age of household occupants, representing the distribution of U.S. households. Behavior and daily routine of agents are based on a set of assumptions that describe contact patterns and movements of individuals. Individuals interact daily with schools or workplaces, household, and neighborhood networks. In each scenario, individuals contact infection based on a fixed mean number of people per day. For detailed descriptions on how agents interact and perform daily routines, the reader is referred to [4, 7, 9–12, 51].

At any time in FRED, each agent is in one of four states: susceptible (S), exposed (E), infectious (I) and recovered (R). We assume that individuals transition from the exposed state to the infectious state when the computed infectivity crosses a defined threshold. All agents are initially susceptible to infection, and at t=0 days, 100 random agents are initialized in the exposed state. Contact with an infectious person has an assigned probability of transmission, which is further modulated by the individual’s infectivity (See Additional file 1).

Model ensemble

The marginals and correlation plots for the posterior distribution P(α|D) are shown in Figure 2 for a sample of 1 million parameter values. Values for the best-fit parameter set are provided in Table 1. Calculation of correlation coefficients across the ensemble shows three parameter pairs to be significantly correlated. Viral production (p) and degradation rates (c) are highly correlated (with correlation coefficient R=0.90), as are symptoms onset (θ) and clearance (a) terms (R=-0.85). A parameter pair of interest are inoculum size (V ₀) and infection rate (β), with correlation R=-0.80. To fit the data, either a small initial viral load pairs with a fast infection rate for the viral trajectory to reach the first data point, or a larger initial viral load pairs with a slower infection rate.

Parameter distributions and correlations of model parameters. (a) One dimensional projections (marginal distributions) of posterior distribution along each (log) parameter axis. Biological bounds constrain values of V ₀, c, δ and k. Other parameters are normalized against their baseline value and shown in log scales. (b) Density plots of correlations of select parameter pairs of parameters.

Full size image

Using a random sub-sample of 15,000 parameter sets from the ensemble, statistical representation of trajectories showing median, 25-75, and 5-95 percentiles of variable values are plotted at each time point in Figure 3. The model fits data well, including the variability in the symptom response. There is variability across trajectories in viral load during the first 2 days, due to the different inoculum sizes V ₀ spanning several log scales and related to the correlation between V ₀ and the infection rate.

Ensemble model trajectories with median (black), 25-75 quartile ranges (dark gray) and 5-95 quartile ranges (light gray) across a random sub-sample of 15,000 parameter sets from the posterior distribution.

Full size image

Representing biological phenotypes

Singular value decomposition on the approximated Hessian of the posterior distribution was performed to identify candidate directions. Due to the strong correlation of the inoculum size V ₀ with the infection rate β revealed in analysis on the ensemble, it is only necessary to vary one of these two parameters (V ₀ and β) to achieve observable differences and the model is more sensitive to changes in infection rate. Therefore, we fix V ₀ associated with a given parameter set and reduced the parameter space to seven dimensions. Variation in inoculum size is still present due to the range of values sampled in the posterior distribution. Combined, directions one and two account for 89.2% of the variability, while the other five components account for the remaining. The leading two components give the following directions in parameter space:

Median trajectories and 69 percentile envelopes from 1,000 random parameter sets sampled from the ensemble varied in each of these directions is shown in Figure 4. The first direction d ₁=d ₁(α) (Figure 4b) corresponds to dynamics related to recovery from infection, and is dominated by δ, the infected cell decay rate. Slower clearance of infected cells results in prolonged viral load, and more severe symptoms that last longer. Changes along this direction resemble differences between adult and aged mice [31], supporting the biological significance of this direction.

Variability of host-level model responses along principal directions of the ensemble. (a) d ₁ and (b) d ₂ (Equations 7 and 8). Black represents the best fit baseline set, red and magenta colored curves represent positive scalar values, and blue and cyan curves represent negative scalar values. Curves are evaluated from a random sub-sample of 500 parameter sets in our posterior distribution. The dark center line of each trajectory is the median, with the 16-85’th percentiles shaded (representing the 69% center values).

Full size image

The second direction d ₂=d ₂(α) (Figure 4a) corresponds to onset dynamics which occur in the first two days of infection; virus that quickly reaches a high viral load corresponds to earlier symptoms than with those with a viral load that is slower to peak. Correlation between time of viral and symptoms peak was observed between the different experimental studies [23, 24, 34, 36], as well as between individual volunteers in an H3N2 study [52].

The response surfaces can be represented as arrays (see the Additional file 1) that provide a fast and efficient way to evaluate infection response. Figure 5 shows examples of the computed response surfaces, which are generated by fixing lines along the 4 edges and 2 centerlines of the response surface domain and evaluating Equation 6 (with coefficients taken from tables in Additional file 1: Methods) for points along each fixed line in Figure 5. Position on the response surface is described by coordinates (x ₁,x ₂) ∈ [ 0,1]× [ 0,1], where the point (0.5,0.5) corresponds to the nominal case (origin of the hyper-plane) associated with the data fitting. Along each line, we fix one scalar value x _j ,j=1,2 and select the other scalar value randomly from a uniform distribution. We repeat this for 35 random trajectories on each line, at x _j =0 (blue), x _j =0.5 (red) and x _j =1.0 (green) for j=1 (Figure 5a) and j=2 (Figure 5b).

Variability of host-level response across the response surface. The colored trajectories in (a) and (b) correspond to points on the response surface located on equally colored lines shown in (c). The surface is parametrized by the coordinates (x ₁,x ₂)∈[0,1]×[0,1]. (The point (0.5,0.5) corresponds to the nominal case). The solid black lines represent thresholds above which individuals are assumed to be infectious or symptomatic.

Full size image

Within-host epidemiological parameters derived in a population-level model

As a starting point for the full population-scale model, we assign two uniform random numbers (x _j ∈U(0,1),j=1,2) to every individual in the population that becomes infected. These numbers correspond to coordinates in the response plane and are used to compute trajectories describing infectivity and symptomaticity of that agent during infection. We assume infection is symptomatic only if the respective trajectory is above a scalar threshold value, which also provide a metric for estimating the duration of latent and incubation periods in the model. The threshold values can be varied according to the particular purpose of the model, as we illustrate in the age-severity study below. Sample thresholds are plotted in Figure 5 as solid black lines.

Symptoms are assigned values 0 (asymptomatic) if below the threshold value or 1 (symptomatic) if above. The symptoms threshold S ^∗=0.717 is selected so that 33% of individuals fall below it in order to maintain a percentage of asymptomatic cases consistent with literature estimates [37]. The viral load threshold is selected at V ^∗=1.35, which gives a mean duration of infectivity of 5.2 days. The mean onset to infectivity calculated from evaluating the response surfaces is 0.83 days. This threshold was selected so the model would be comparable with estimates in the literature, and a comparison is found in the Discussion [6, 30, 43, 53–59].

The length of onset and duration periods are estimated by interpolating where the evaluated trajectories crosses the threshold line and the means are reported in Table 2. Histograms representing the distributions of times of onset and duration of infectiousness and symptoms, as well as time of peak symptoms and infectivity and corresponding values are shown in Figure 6 for response surface evaluations corresponding to one million random numbers (infected individuals). We separate the asymptomatic cases (maximum symptoms score < 0.717) and show them in red, while symptomatic cases are in blue. The maximum infectivity over all infections is 3.71 (compared to the maximum infectivity data point of 2.88 log10 TCID ₅₀), which we use as a linear coefficient to normalize computed viral load into transmission probability in the agent-based simulation. We refer to the above described within-host model as the “baseline” model, and calibration of the baseline case is discussed in the Additional file 1.

Histograms showing distribution of onset, duration, and peak of infectivity and symptoms. Histograms represent distributions across 1 million random numbers (representing infected individuals), evaluated on the response surface (Equation 6). Asymptomatic individuals (maximum symptoms score < 0.717) are shown in blue, while symptomatic individuals are in red. (a) Onset and duration of infectivity and symptoms are computed from each response surface evaluation. (b) Time and value of infectivity and symptoms at their peak.

Full size image

Comparison of the baseline model with a reference ABM within-host model

We first show that the baseline within-host model produces similar results to a standard within-host model used in ABM studies, which we call the “reference” ABM model [4, 10]. In the reference within-host model in FRED, the mean duration of infectivity (4.1 days) and mean latent period (1.2 days) can be evaluated directly using the discrete probability distributions explicitly used in the model, and give an average total duration of infection as 5.3 days from inoculation to resolution [4]. Each infected agent is either symptomatic or not, based on a Bernoulli random variable.

In both cases in this comparison, all individuals who become infectious have an equal probability of contracting the same response, meaning that when an individual becomes infected they are assigned random numbers to represent illness independent of information such as a demographic or pre-existing conditions. Therefore the observed epidemiological differences are due only to the differences in the models of within-host response. The main difference between the two within-host models is that symptomatic individuals have a higher infectivity lasting over a shorter period and asymptomatic individuals have a lower infectivity lasting over a longer period in the reference ABM within-host model, compared to the new baseline within-host model.

To measure the differences resulting from the two within-host models in a population-level simulation, we run each epidemic scenario with 50 distinct random seeds, and compare epidemiological measurements in Table 3. We observe that the peak incidence corresponding to the baseline within-host model occurs earlier than the reference model (31.1 days compared to 25.1 days), which can be attributed to the shorter incubation period and reduced transmission rates in symptomatic individuals. Also observed is a larger peak incidence in the baseline case (36,811 compared to 41,883 infected), attributed to longer duration of infection. The reproductive rate R ₀ is estimated by the average number of secondary infections from the original 100 exposed individuals on day 0 computed on day 10. We observe that R ₀ is slightly higher in the baseline at 2.14+0.38 than default FRED at 1.98+0.37. We measure the average time to run one simulation over 100 days, and find that the Baseline case has an average run time of 79 seconds and in the standard ABM model has an average run time of 80 seconds on a laptop.

Age-severity study

Aging alters both the innate and adaptive branches of the immune system, influencing response to influenza virus infection [60]. We assume that age of an infected individual maps to a scalar input value along d ₁, reflecting longer and more severe infection and symptoms in older individuals. To illustrate the flexibility of our immune-epidemiological approach, we conduct a study where we compare the simulated epidemics between three different age-severity maps: a linear increasing function, and two functions based on mortality data from past epidemics in 1918 and 1957 [61]. Mortality patterns associated with epidemics traditionally form a “U” curve, reflecting higher mortality and more severe disease in the very young and the elderly. A distinctive feature of the 1918-19 influenza pandemic was a “W” shaped pattern in the mortality curve, showing a secondary peak in mortality rate for the 25-35 years age group [62, 63]. The age-severity model thus assumes mortality rates are indicators of severity of disease experienced by different age groups. Piece-wise linear functions map age to [0,1] for all three maps, and mortality data is scaled so the area under all three curves (for a g e∈[0,100]) is preserved. In this way the model considers two factors linking age with disease severity: deterioration of the immune system with age and immunological memory. Immunological memory refers to the adaptive immune response to pathogen the body has previously encountered, and contributes to explaining the “W” shape in the 1918 pandemic mortality curve [64].

The median age in the United States is 37.2 years, while Allegheny County has a median age of 41.3 [65], among the oldest in the nation. Franklin County, Ohio is comparable to Allegheny in population size, density and household size (see Table 4), but has a younger age distribution (median age 33.4). The age-severity model is also studied on Franklin County in order to compare how different age distributions affect the overall epidemic. Plots of the three “age-to-phenotype” maps superimposed over age distribution histograms for both counties are presented Figure 7(a)-(b).

Age-severity study results. Panels (a) and (b) show age-distribution of population in two counties (left axis). Superimposed are three age-severity models showing the assigned value of x ₁ coordinate as a function of age. (c) Attack rates for each scenario. (d) Fold change from the baseline model for each age-severity model.

Full size image

Epidemiological estimates averaged across 50 random simulations are reported in Table 5. In the baseline scenario, the higher attack rate in Franklin County illustrates the role age distribution has on the outcome of an epidemic (Figure 7(c)). Younger individuals, particularly school aged, have higher contact probabilities than older individuals. This is supported by decreased attack rates in both counties when the age-severity models are imposed. Across the three age-severity models, the W-curve severity model has the highest attack rate, illustrating how more severe and prolonged illness in the 25-35 age group can increase epidemic severity compared to either the linear or U-shaped model. The differences between the linear and the U-shaped models are less dramatic, particularly in Allegheny County.

We plot the fold difference in the attack rates in the age-severity models to their baselines in order to make a direct comparison between the models (Figure 7(d)). We observe a larger fold change for the W-model in Allegheny County, showing that for these counties, the decreased illness severity in the 45-65 ranges has a larger impact than the increased illness severity in the 25-35 age range. The larger fold change in the linear and U-shaped models of Franklin County implies a more dramatic change in epidemic outcome when age-severity assumptions are made in the younger age distribution scenario. The model assumes symptomatic cases by determining threshold maximum symptoms score. Less severe illness is more likely to correspond to asymptomatic cases, so the age-severity model also has tendencies to have an age dependence on symptomatic cases. We report the percentage of symptomatic patients in each scenario in Table 5.

The strategy for incorporating the host-level model described here could be employed with other deterministic population level models as well as other well-established ABM simulation platforms such as EpiSimS [3]. Our aim here is to develop a methodology to bridge the gap between within-host models of virus/symptoms dynamics and population-level models that may also capture human movement and demographics. FRED provides a ready to use framework for population-level modeling using the synthetic population data set, allowing easy extension of the within-host model to simulations of populations with different demographic descriptions (here by distribution of age). We note that integration of our within-host model into a deterministic population level model may have similar results with less complexity, and we provide the coefficients for the response surfaces in the supplemental for straightforward use.

Our model implicitly uses age-structured contact patterns in the synthetic data set and ABM, and the results from the age-severity study reflect this relationship between distribution of age and epidemic severity. Previous studies have derived location specific contact matrices, either through use of synthetic populations with ABM network structures [67, 68], or from surveys of European communities [69]. In these cases, two mixing patterns emerge: a strong contact structure between people of similar age and a weaker structure between children and middle-aged adults, most likely due to parent-child contact. Use of contact matrices based on social networks better capture age-specific infection patterns of infectious diseases which helps in estimation of the basic reproductive number R ₀. In our age-severity study, we provide a method to further fine-tune these estimates depending on host-response to a circulating strain. The computed values of R ₀ in this study are derived from assumptions about one particular experimental strain, and results may vary depending on characteristics unique to different circulating strains.

When integrating our derived within-host model into the population-level framework of FRED, further assumptions were required on scaling and thresholds. The symptoms threshold was set so approximately 1/3 of individuals in the baseline scenario are asymptomatic, consistent with accepted public health estimates [37]. However, for the infectivity threshold I ^∗, there was flexibility. Different values of I ^∗, correspond to different latent periods, impacting the duration of time an individual is in the exposed category. The threshold (1.35) was selected such that the average baseline latent periods are 0.83 days, and mean duration of infectivity is 5.17 days. Estimates for latent period are experimentally unknown due to difficulty in measuring, but a mathematical model predicts 0.4-1.5 days [53]. Duration of infectivity has been estimated at 4.8 days [37], and various studies have reported estimates at 3.38 days [30], 6.6 days [54], 5.0 days [55], and 3.1 days [43]. Our mean value of 4.98 days thus is consistent with existing estimates. This assumption reflects that in the baseline model, the incubation period exceeds the latent period, so there is a period in which an individual is not symptomatic and shedding virus. In one modeling study, a lag of 1.9 days between latent and incubation periods was assumed [6].

In our study, the incubation period is estimated at 2.14 days, duration of symptoms about 4 days for the symptomatic cases, in which both estimates rely on our choice of symptoms threshold. Incubation period has been estimated anywhere from 1.4 days (1.3 - 1.5 95% CI) [56], 1.5-2 days [57], 2-3 days [58], to 4 days [30]. Estimates for duration of symptoms vary from 4.5-5 days [37], 6 days [59], 5.6 days [55]. We remark that lowering the symptoms threshold would coincide with a shorter incubation period. One estimate that is well agreed upon is that viral load peaks at around 2 days post inoculation, and symptoms peak at around 3 days p.i. [37], which is consistent with our simulation measurements as well.

We selected the mathematical model of Baccam et al.[38] to describe viral and symptoms dynamics for three reasons. First, structural identifiability analysis showed parameter values could be estimated given the data. Second, it is the simplest model that generated the required inputs for a population level simulation (infectivity and symptoms score), and last, the model fit the data. Our method also provides a method for varying host-level responses from experimental baseline measurements in order to predict various illness scenarios under different host assumptions. As a starting point, we used a nominal response calibrated to one non-pandemic seasonal human influenza strain (A/Texas/36/91). As more data allows more subtle description of disease phenotypes characterized by more complex within-host models, our method could be easily scaled to such models as, although computation of the response-surfaces is more complex, it is done a priori and thus requires no additional resources for the population-level simulation.

Community based studies also give insight corresponding to our assumption that infection in the elderly may be characterized by longer viral load, and longer and higher symptoms scores. One study suggests a median duration of symptoms of 7 days for the entire community, compared to 8 days for older people [30]. Another study estimates that duration of infectiousness/viral shedding is 6.3 days for children, and 6.7 days for adults [54]. Using a linear mapping from age to infection severity is a simplistic assumption without biological motivation, which is why three different mappings were studied.

Our model estimates of duration and intensity of infectivity and symptoms show asymptomatic infections correspond to reduced viral load, which is consistent with previous studies of asymptomatic disease and transmission, but the degree is unknown [70]. In a community-based study of 2009 H1N1 in Hong Kong, peak viral load averaged over the asymptomatic subjects was 3.2e3 copies/mL compared to 3.6e7 copies/mL averaged over the symptomatic subjects [57]. In a recent human volunteer study of H3N2/Wisconsin, 50% of asymptomatic subjects had evident viral shedding and a reported a significantly reduced viral load in the asymptomatic cases [52]. However, the exact mechanism contributing to asymptomatic infection is not incorporated in this model, and a study of these mechanisms and how they influence an epidemic are important future considerations.

This study makes the assumption that infectivity scales linearly with contact rate and is proportional to the log10 of the viral load as measured by nasopharyngeal swabs. Although this is plausible, the exact relationship between viral load and infectivity on disease transmission is unknown. One model, which studies the relationship between virus/host dynamics and a population-level model, multiplied viral load by the amount of nasal discharge as the estimate of infectivity [71]. Another model of measles used the area under the infected cell curve to model transmission strength by 'amount’ of infection which was also used as a determinant for symptomatic or asymptomatic infection [72]. In a household study of influenza transmission in Hong Kong, models of molecular viral shedding and log10 molecular viral shedding were studied to compare viral load with transmission [57]. The case in which viral load is used as an indicator of transmission leads to nearly all transmission occurring in the first 1-2 days after inoculation, implying that intervention strategies must be rapid. Using the log10 model of viral shedding to transmission corresponds to significant transmission lasting strongly up to 3-4 days after viral inoculation. Therefore, knowing the exact relationship between infectivity and disease transmission would improve modeling and have significant implication for mitigation strategies.

In the mathematical model describing the host immune response, symptoms scores are generated from the biological mechanisms described by the ODE model, which greatly simplifies complex immune system dynamics. The immune system is composed of innate and adaptive components, which respond to viral infection through inflammation and creation of virus-specific antibodies and effector cells to help resolve the infection [73]. Previous models have been studied within the context of improved understanding of population-level dynamics have incorporated interactions of the innate immune system [53, 74]. Incorporating features such as inflammation and virus-specific antibody production would yield more realistic within-host dynamics, and could explicitly account for factors such as preexisting immunity. Importantly, more complex models could be attempted using response surface methodology at little additional computational cost for large-scale simulations. Further, the ensemble methodology allows the computation of the thickness of the response surface, which would add an additional stochastic element better representing our incomplete knowledge of individual responses.

The added flexibility of a mechanistic within-host to population-level models might be particularly relevant when evaluating the impact of host-level containment measures and additional factors known to demonstrate significant within-host effect depending on immune status. The ability to include this flexibility at minimal computational cost also offers exiting possibilities when simulating new variants, co-infection with bacterial pathogens, for providing a biological basis for morbidity, undoubtedly a major driver of agent behavior and cost of disease, and for mortality related to disease, which is currently lacking in existing population-level simulators.

We have developed a multi-scale immuno-epidemiological ABM model linking an equation-based within-host model to an agent-based population level model. The key innovation is the use of response surface analysis to simplify the computation of output trajectories of the equation-based within-host model, which makes the immuno-epidemiological model just as efficient and scalable as the underlying population-level ABM. The hybrid ABM model replicates results of simulated epidemics with simpler within-host assumptions such as the reference FRED within-host model, and provides additional flexibility to explore questions relevant to host heterogeneity at no extra computational cost.

Research reported in this publication was supported by the National Institute Of General Medical Sciences of the National Institutes of Health under Award Number U54 GM088491, R01-GM83602, and NSF RTG DMS0739261. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

The authors wish to acknowledge David Galloway and Molly Eggleston at the Public Health Dynamics Laboratory in the University of Pittsburgh Graduate School of Public Health for their help during the preparation of this manuscript.

Affiliations

1. Department of Mathematics, University of Pittsburgh, Pittsburgh, PA, USA

   • Sarah Lukens
   •  & David Swigon

2. Pittsburgh Supercomputing Center, Carnegie Mellon University, Pittsburgh, PA, USA

   • Jay DePasse
   •  & Shawn T Brown

3. School of Computer Science, Carnegie Mellon University, Pittsburgh, PA, USA

   • Roni Rosenfeld

4. Center for Vaccine Research, Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, PA, USA

   • Elodie Ghedin

5. Department of Critical Care Medicine, University of Pittsburgh, Pittsburgh, PA, USA

   • Gilles Clermont

6. Department of Biostatistics, University of Pittsburgh Graduate School of Public Health, Pittsburgh, PA, USA

   • John Grefenstette

7. Department of Computational and Systems Biology, University of Pittsburgh, Pittsburgh, PA, USA

   • Ericka Mochan

8. Department of Epidemiology, University of Pittsburgh Graduate School of Public Health, USA

   • Donald S Burke

9. Department of Biological Sciences, University of Notre Dame, South Bend, IN, USA

   • Sarah Lukens

Corresponding author

Correspondence to Sarah Lukens.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

SL conceived of and designed the study, analyzed and interpreted data, ran models, and drafted the manuscript. JDP made essential computational contributions, particularly with design of host class and implementation of the model at both scales. RR helped conceive study at all levels. EG helped conceive study and edited the manuscript. EM participated in study design and coordination, edited the manuscript. SB contributed to study design, computational programming, and edited the manuscript. JG contributed to study design, computational programming and edited the manuscript. DB helped conceive study and contributed to study design, and provided funding for the study. DS helped conceive study, participated in its design and coordination, edited the manuscript, provided funding for the study. GC conceived of the study, and participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.

Below are the links to the authors’ original submitted files for images.

Reprints and Permissions