Skip to main content
Advertisement
  • Loading metrics

Sequential infection experiments for quantifying innate and adaptive immunity during influenza infection

  • Ada W. C. Yan ,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

    a.yan@imperial.ac.uk

    Affiliations School of Mathematics and Statistics, The University of Melbourne, Parkville, Victoria, Australia, MRC Centre for Global Infectious Disease Analysis, Department of Infectious Disease Epidemiology, School of Public Health, Imperial College London, London, United Kingdom

  • Sophie G. Zaloumis,

    Roles Conceptualization, Methodology, Supervision, Writing – review & editing

    Affiliation Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Parkville, Victoria, Australia

  • Julie A. Simpson,

    Roles Conceptualization, Funding acquisition, Writing – review & editing

    Affiliation Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Parkville, Victoria, Australia

  • James M. McCaw

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

    Affiliations School of Mathematics and Statistics, The University of Melbourne, Parkville, Victoria, Australia, Centre for Epidemiology and Biostatistics, Melbourne School of Population and Global Health, The University of Melbourne, Parkville, Victoria, Australia, Modelling and Simulation, Infection and Immunity Theme, Murdoch Childrens Research Institute, The Royal Children’s Hospital, Parkville, Victoria, Australia

Abstract

Laboratory models are often used to understand the interaction of related pathogens via host immunity. For example, recent experiments where ferrets were exposed to two influenza strains within a short period of time have shown how the effects of cross-immunity vary with the time between exposures and the specific strains used. On the other hand, studies of the workings of different arms of the immune response, and their relative importance, typically use experiments involving a single infection. However, inferring the relative importance of different immune components from this type of data is challenging. Using simulations and mathematical modelling, here we investigate whether the sequential infection experiment design can be used not only to determine immune components contributing to cross-protection, but also to gain insight into the immune response during a single infection. We show that virological data from sequential infection experiments can be used to accurately extract the timing and extent of cross-protection. Moreover, the broad immune components responsible for such cross-protection can be determined. Such data can also be used to infer the timing and strength of some immune components in controlling a primary infection, even in the absence of serological data. By contrast, single infection data cannot be used to reliably recover this information. Hence, sequential infection data enhances our understanding of the mechanisms underlying the control and resolution of infection, and generates new insight into how previous exposure influences the time course of a subsequent infection.

Author summary

The resolution of an influenza infection requires different components of the immune response to work together. Despite recent advances in mathematical modelling, we do not well understand how much each broad immune component contributes to this process at a given time. In this study, we show that laboratory data on the amount of virus over the course of a single infection is insufficient for inferring the contribution of each broad immune component. However, if the animals are exposed to two different virus strains with only days separating exposures, then the timing and strength of protection provided by the first infection against the second provides crucial additional information. We show how mathematical models can be used to recover the timing and strength of each immune component, thus enhancing our understanding of how an infection is controlled, and how a previous exposure changes the time course of a subsequent infection.

Introduction

The influenza virus infects epithelial cells in the respiratory tract, causing respiratory symptoms such as coughing and sneezing, and systemic symptoms such as fever. Three main components of the immune response—innate, humoral adaptive and cellular adaptive immunity—work together to control an infection. Experiments have revealed the contribution of each major immune component to resolution of an infection, by suppressing each immune component in turn [15]. However, current mathematical models do not agree on how each major immune component contributes quantitatively.

A study by Dobrovolny et al. [1] highlights these discrepancies. The study showed that eight existing viral dynamics models [613] made different qualitative predictions when different components of the immune response were removed. Each model failed to reproduce the effect of removing at least one of the three components discussed above. The discrepancies arose because many models were only fitted to viral load data from a single infection.

It has been shown that many models can fit the viral load for a single infection well, including models without a time-dependent immune response which are thought to be less biologically realistic [6]; however, if data for multiple initial conditions are available, the viral load may have more features to distinguish between competing models [1416]. One way of altering the initial conditions is through a previous or ongoing infection. We previously conducted a series of experiments where ferrets were sequentially infected with two influenza strains [17, 18]. When a short time interval (1–14 days) separated exposures, a primary infection protected against a subsequent infection. This protection likely arose through cross-immunity, whereby the immune response stimulated by one strain also protects against infection with another.

While immune markers indicated the approximate timing of each arm of the immune response [19], the strength of cross-protection due to each component was difficult to measure experimentally. We hypothesised that mathematical models can be used to gain further insight from these types of experiments. Few existing models include interactions between influenza strains on short timescales; hence, we constructed viral dynamics models to reproduce the qualitative observations of these experiments [20, 21]. The models also reproduce observations from a range of experiments where immune components were suppressed [22].

Here, we use simulations to show that these mathematical models allow us to extract the timing and strength of cross-protection from sequential infection data. By attributing cross-protection to specific immune components, the models lead to new insight into how previous exposure influences the time course of a subsequent infection. Moreover, we find that compared to single infection experiments, sequential infection experiments provide richer information on host immunity, and thus are potentially a powerful tool to study immune-mediated control of a primary infection.

Results

Synthetic data

As a first step to compare the information made available by sequential infection versus single infection experiments, we generated synthetic datasets for each scenario. Mimicking the experimental procedure of Laurie et al. [17], we generated a sequential infection dataset where ferrets were exposed to two influenza strains, with intervals of 1, 3, 5, 7, 10 and 14 days between exposures; and a single infection dataset where ferrets were exposed once only. Details are given in the Materials and Methods section. Using synthetic data means that we know the ‘true’ contribution of each immune component in resolving a single infection, and the ‘true’ extent of cross-protection between infections.

Fig 1 shows a subset of the synthetic data. For a single infection, the viral load trajectory can be split into exponential growth, plateau and decay phases. For short inter-exposure intervals (1–5 days), infection with the challenge virus was delayed; for long inter-exposure intervals (7–14 days), infection with the challenge virus was unaffected. These features of the synthetic data match the qualitative results of Laurie et al. [17] for infection with influenza A followed by influenza B, or vice versa. The parameter values were chosen such that the delay was due to innate immunity. This choice was made because experimentally, innate immune markers such as type I interferon were observed to be elevated 1–5 days after a primary infection [19], and our previous mathematical model incorporating the innate immune response made predictions consistent with the observed temporary immunity [20]. The full set of synthetic data is provided in S1 Fig.

thumbnail
Fig 1. A subset of the synthetic data.

(a) The line shows the simulated ‘true’ viral load for a single infection, with the arrow showing the time of exposure. The simulated viral load with noise is shown as crosses. The horizontal line indicates the observation threshold (10 RNA copy no./100μL); observations below this threshold are plotted below this line. Values below the observation threshold were treated as censored. (b—c) For sequential infections with the labelled inter-exposure interval, the dashed and dotted lines show the simulated ‘true’ viral load for a primary and challenge infection respectively; the arrows show the times of the primary and challenge exposures. The simulated viral load with noise is shown as crosses.

https://doi.org/10.1371/journal.pcbi.1006568.g001

Verification of the fitting procedure

In this section, we first verify that we could recover the simulated ‘true’ viral load by fitting our model to the data. In the next section, we will sample from the joint posterior distributions thus obtained to extract the contribution of each immune component.

Fig 2a shows that the simulated ‘true’ viral load was recovered accurately when fitting the model to either sequential infection or single infection data. The blue and green (overlapping) areas are 95% credible intervals predicted by the models fitted to the sequential infection and single infection data respectively. Both shaded areas included the simulated ‘true’ viral load, shown as the dotted line. This consistency indicates that the fitting procedure accurately recovered the viral load.

thumbnail
Fig 2. Verification that the fitting procedure recovered the viral load.

(a) For a single infection, the blue and green areas are the 95% credible intervals for the viral load (in the absence of noise), as predicted by the models fitted to the sequential infection and single infection data respectively. (b—c) For sequential infections with the labelled inter-exposure interval, the grey and blue areas show the 95% credible intervals for the primary and challenge viral load respectively, predicted by the model fitted to sequential infection data. The other elements of the figure are identical to Fig 1: the dashed and dotted lines show the simulated ‘true’ viral load for a primary and challenge infection respectively; the arrows show the times of the primary and challenge exposures; and the horizontal line indicates the observation threshold.

https://doi.org/10.1371/journal.pcbi.1006568.g002

Fig 2b and 2c confirm that fitting to sequential infection data accurately recovered the viral load for different inter-exposure intervals. The grey and blue areas show the 95% credible intervals for the primary and challenge viral load respectively.

Comparing the immunological information in each dataset

Next, we compared the behaviour of the fitted models to the behaviour of the ‘true’ parameters, to determine the information in each dataset on

  • the effect of each immune component in controlling a single infection;
  • cross-protection between strains; and
  • each immune component’s contribution to cross-protection.

The effect of each immune component in controlling a single infection.

In Fig 3, we removed various immune components from the model. We then compared predictions of the viral load for a single infection by the models fitted to the two datasets.

thumbnail
Fig 3. Predicting the viral load for a single infection when various immune components were absent.

The vertical lines indicate, for the ‘true’ parameter values, the times at which the immune components labelled under each panel took effect. These times were determined by when the viral load for the baseline model (red dotted line) deviated from the viral load when the immune components were absent (black dashed line). These times were recovered using sequential infection data in all of the panels (95% prediction intervals for the viral load in blue), while the timing of adaptive immunity in (a) was recovered using single infection data (intervals in green). In addition, the viral load when adaptive immunity was suppressed was accurately predicted using sequential infection data (a). However, the viral load was not accurately predicted using either dataset in the remaining scenarios (b—e). Prediction intervals were constructed without measurement noise.

https://doi.org/10.1371/journal.pcbi.1006568.g003

First, we showed how the viral load trajectory for the ‘true’ parameters changed when adaptive immunity was suppressed. We defined the ‘baseline’ as the viral load when all immune components were present (red dotted lines in Fig 3, which are the same as the black lines in Figs 1a and 2a). Suppressing adaptive immunity prevented resolution of the infection (black dashed line in Fig 3a), which was consistent with findings of a previous experiment [5]. The viral load deviated from the baseline trajectory at 4 days post-exposure (vertical line), indicating that this was the time at which adaptive immunity took effect.

We then asked whether the models fitted to the two datasets predicted this change. Chronic infection in the absence of adaptive immunity was only predicted using sequential infection data (Fig 3a). Single infection data did not enable consistent prediction of this outcome, as indicated by the broadening prediction interval. However, both datasets enabled recovery of the time at which the viral loads in the presence and absence of adaptive immunity deviated (the vertical line in Fig 3a). Hence, the timing of adaptive immunity was accurately estimated using either dataset.

In Fig 3b–3e, we repeated this procedure, suppressing (b) innate immunity, (c) all immunity, (d) humoral adaptive immunity, or (e) cellular adaptive immunity. When innate immunity was suppressed, the peak viral load increased and the recovery time decreased. The increase in peak viral load was consistent with previous studies where innate immunity was suppressed [2, 23]. The decrease in recovery time was not observed in these studies, and may be an artefact of modelling adaptive immunity as independent of innate immunity. When both innate and adaptive immunity were suppressed, the peak viral load increased, and resolution of the infection was delayed. These changes were consistent with a previous experiment where innate immunity was suppressed [2].

Fig 3b shows that sequential infection data enabled accurate inference of when the viral loads in the presence and absence of innate immunity deviated, hence recovering the timing of innate immunity. By contrast, the model fitted to single infection data predicted that the viral loads could deviate much earlier. Neither model accurately predicted how the infection resolved in the absence of innate immunity; however, the prediction intervals for the model fitted to sequential infection data were tighter, and the peak viral load was consistently predicted to be higher than for the baseline model. Similarly, when both innate and adaptive immunity were absent, the model fitted to sequential infection data recovered the timing of overall immunity, but could not predict the viral load in the absence of immunity (Fig 3c).

Without innate immunity, the viral load peaks due to target cell depletion, and without any immune response, the infection resolves due to target cell depletion. The lack of predictive ability indicates that both datasets lack information on how target cells would hypothetically become depleted, and how this depletion would affect the viral load, in the absence of the immune response. One is thus cautioned against using parameter values from a model fitted to data in immunocompetent hosts to make predictions in situations where target cells may become severely depleted, such as if individuals are immunocompromised.

Fig 3d and 3e show that neither dataset enabled prediction of how the viral load changed when (d) the humoral adaptive immune response or (e) the cellular adaptive immune response was removed. This implies that sequential infection data (of the type reported in Laurie et al. [17]) cannot be used to distinguish the contributions of antibodies and cellular adaptive immunity to resolution of infection. In detail, the ‘true’ parameters predicted that when humoral adaptive immunity was disabled, the viral load rebounded instead of continuing to decrease (black dashed line in Fig 3d). When cellular adaptive immunity was disabled, resolution of the infection was delayed (black dashed line in Fig 3e). The fitted model’s predictions ranged from no delay to a chronic infection.

Cross-protection between strains.

Given the above mixed results, we then tested whether sequential infection data accurately captured the timing and extent of cross-protection, by simulating the viral load for inter-exposure intervals other than those where data was provided. We selected new inter-exposure intervals of 2 and 20 days; the former lay between inter-exposure intervals included in the original data (1, 3, 5, 7, 10 and 14 days), while the latter lay outside this range. Then, using the models fitted to the original data (that is, not re-fitting to the new data), we predicted the challenge viral load for these new inter-exposure intervals. Because a primary infection could greatly affect a challenge infection, but not vice versa, we focused on the behaviour of the challenge infection.

Fig 4 shows that predictions by the model fitted to sequential infection data (blue areas) were accurate. By contrast, predictions using single infection data (green) did not agree well with the ‘true’ viral load. Note that to predict cross-protection using single infection data, we used the model assumptions that innate immunity was completely non-specific and antibodies were completely strain-specific, and considered the optimistic scenario where we had independent, perfect information about the proportion of cellular adaptive immunity that was cross-reactive (details in the Materials and methods section). Even then, single infection data did not accurately capture cross-protection.

thumbnail
Fig 4. Predicting the outcomes of further sequential infection experiments.

Sequential infection data, but not single infection data, enabled prediction of further sequential infection experiment outcomes. The lines show the simulated ‘true’ viral loads for inter-exposure intervals of (a) 2 and (b) 20 days. The shaded areas show the 95% prediction intervals for the challenge viral load.

https://doi.org/10.1371/journal.pcbi.1006568.g004

Each immune component’s contribution to cross-protection.

Having shown that the sequential infection data captures the timing and extent of cross-protection between strains, we then asked whether such cross-protection could be attributed to the ‘correct’ mechanisms (the same mechanisms as given by the ‘true’ parameters). These mechanisms are

  • target cell depletion due to the infection and subsequent death of cells;
  • innate immunity; and
  • cellular adaptive immunity.

In our model, antibodies are strain-specific and thus do not contribute to cross-protection.

Before analysing the behaviour of the fitted models, we quantified how each immune component contributed to cross-protection for the ‘true’ parameters. In Fig 5, for a one-day inter-exposure interval, we plotted in red the challenge viral load for the baseline model (the original model fitted to the data, where all three of the above immune components can mediate cross-protection). We observed that the challenge infection was delayed relative to a primary infection. We then modified the baseline model such that only a subset of immune components mediates cross-protection, as detailed in the Materials and Methods section. We used the modified model to predict the viral load (in black), and compared it with the baseline viral load. (The blue areas will be discussed shortly).

thumbnail
Fig 5. Predictions of the challenge viral load for a one-day inter-exposure interval when the mechanisms mediating cross-protection were restricted.

The black solid lines show the challenge viral load for the ‘true’ parameter values when the mechanisms mediating cross-protection were restricted using (a) model XC, (b) model XIT, (c) model XI, or (d) model XT. The red dotted lines show the viral load for the baseline model. Comparing the two sets of lines revealed that innate immunity mediated cross-protection, whereas cellular adaptive immunity and target cell depletion did little to mediate cross-protection. The model fitted to sequential infection data accurately predicted the challenge outcomes for models XC and XIT, but not model XI or model XT (95% prediction intervals shown). It thus correctly attributed cross-protection to target cell depletion and/or innate immunity, but could not definitively distinguish between the two. For clarity, the viral load for the primary infection is not presented in this figure.

https://doi.org/10.1371/journal.pcbi.1006568.g005

For example, in Fig 5a, we modified the baseline model such that only cellular adaptive immunity, and not target cell depletion or innate immunity, can mediate cross-protection. We denoted this modified model ‘model XC’. Unlike the baseline model (red dotted line), the challenge viral load for model XC was not delayed (black solid line); in fact, it closely resembled that for a single infection. Comparing the two simulations led to the conclusion that cellular adaptive immunity did not play a major part in cross-protection.

We then modified the baseline model such that both target cell depletion and innate immunity can mediate cross-protection, but cellular adaptive immunity cannot do so. We denoted this model ‘model XIT’. The challenge viral loads according to model XIT and the baseline model were very similar (overlapping lines in Fig 5b). Hence, for the ‘true’ parameters, cross-protection was mediated by innate immunity and/or target cell depletion.

To distinguish between these two mechanisms, we constructed model XI, where only innate immunity, and not target cell depletion or cellular adaptive immunity, can mediate cross-protection. Once again, the challenge viral load was very similar to the baseline model (overlapping lines in Fig 5c). We also constructed model XT, where only target cell depletion, and not innate immunity or cellular adaptive immunity, can mediate cross-protection. The challenge viral load for model XT was not delayed, and resembled that for a single infection (Fig 5d). We concluded that the cross-protection was largely mediated by innate immunity.

Having demonstrated the utility of the modified models, we returned to the original question of whether sequential infection data could be used to distinguish between mechanisms for cross-protection. We sampled parameter sets from the joint posterior distributions obtained by fitting the baseline model to sequential infection data, and used them as inputs for models XC, XIT, XI and XT respectively, to generate the blue areas in Fig 5. If the fitted parameters and the ‘true’ parameters predict the same infection outcomes under the modified models, then the fitted model attributes cross-protection to the ‘correct’ mechanisms.

Models XC and XIT made the same predictions using the fitted parameters (shaded area) and the ‘true’ parameters (black line), so sequential infection data enabled us to accurately attribute cross-protection to target cell depletion and/or innate immunity, rather than cellular adaptive immunity (Fig 5a and 5b). On the other hand, the fitted parameters did not consistently predict the challenge outcome for models XI and XT (Fig 5c and 5d). Hence, we could not use sequential infection data to consistently rule out the possibility that cross-protection occurred due to both target cell depletion and innate immunity. However, only a very small proportion of trajectories sampled from the joint posterior distribution incorrectly attributed the delay to both target cell depletion and innate immunity (S2 Fig). Moreover, the fitted model was able to rule out the possibility that target cell depletion alone was responsible for cross-protection, as the 95% credible interval for model XT does not include the 95% credible interval for the baseline viral load.

Similarly, we were unable to disentangle different mechanisms of innate immunity from the sequential infection data alone (S3 Fig and S1 File).

For infection with heterologous influenza A strains, rather than the influenza A and B strains discussed thus far, we hypothesise that innate and cellular adaptive immunity contribute to cross-protection at different inter-exposure intervals [24]. S2 File presents the same analysis for this scenario, where we were able to unravel these different contributions using sequential infection data.

Our findings are robust to the ‘true’ parameters chosen, given that the parameters capture the qualitative observations by Laurie et al. [17]. S5 File presents the same analysis for a different set of ‘true’ parameters where the degree of cross-reactivity in the cellular adaptive immune response is low. For the most part, the same qualitative results were obtained: sequential infection data enabled inference of the timing of innate and adaptive immunity, prediction of the viral load for further experiments with different inter-exposure intervals, and prediction of the viral load for models XC and XIT. The main difference was that the model fitted to sequential infection data was only able to capture the timing of adaptive immunity, and not how the infection would resolve if adaptive immunity were removed. A possible explanation for the different result is that for the parameters in the main text, the viral load showed a clear plateau while innate but not adaptive immunity was active, enabling the fitted model to predict that the viral load would stay at that plateau in the absence of adaptive immunity. By contrast, the viral load for the different set of parameters in the supplementary material did not show a clear plateau during the innate immunity phase, possibly reducing the fitted model’s ability to infer the viral load in the absence of adaptive immunity. Another minor difference was that the immune components whose timing could be inferred using single infection data were different from in the main text, although the overall finding—that sequential infection data enabled inference of the timing of more immune components—was the same.

For a given set of ‘true’ parameters, repeating the entire study with different simulated noisy datasets did not change our findings (S6 File).

Because our model does not capture every biological detail of the experimental system, we also tested whether our findings were robust to model misspecification. We generated data using a model from a study by Zarnitsyna et al. [25], modified to include a variable degree of cross-reactivity between strains. We then fitted the model in the present study to the generated data. Even though the data was generated using a different model, we were still able to infer the timing of innate and adaptive immunity, predict the outcome of further experiments with different inter-exposure intervals, and distinguish the contributions of different components of the immune response to cross-protection. This analysis is presented in S7 File.

In summary, the synthetic sequential infection data enabled accurate inference of the contribution of cellular adaptive immunity to cross-protection, as well as the combined contributions of target cell depletion and innate immunity. However, using this data alone, we could not conclusively distinguish the contributions of innate immunity and target cell depletion to cross-protection, or distinguish the contributions of different innate immune mechanisms.

Discussion

Advantages of sequential infection experiments

In this study, we have simulated experiments which investigate the interaction of influenza strains through sequential infections, then explored how mathematical models could be applied to the data to gain insight into immune mechanisms. Our analysis has shown that the sequential infection study design, compared to the single infection study design, provides richer information for inferring the timing and strength of each immune component.

We have identified three main advantages of sequential infection data. The first advantage is in inferring how each immune component helps to resolve a single infection. We found that the synthetic sequential infection data captures the timing of innate and adaptive immunity during a single infection, and thus enables accurate prediction of the outcomes of some in silico experiments where immune components were removed. In contrast, we could not consistently infer the timing and strength of innate immunity from single infection data. Moreover, single infection data contains information only on the timing of adaptive immunity, but not the effects of suppressing adaptive immunity.

The second advantage is that sequential infection data contains more information on the effects of cross-protection. We were able to use the model fitted to the sequential infection data to precisely predict outcomes of further such experiments using the same strains but different inter-exposure intervals. Using the model fitted to the single infection data greatly reduced predictive power.

The third advantage is in inferring the contribution of each immune component to this cross-protection. For the dataset in the main text, we were able to infer that cellular adaptive immunity played little role in cross-protection, and that innate immunity and/or target cell depletion led to the observed cross-protection. We also showed that target cell depletion alone could not explain this cross-protection.

Collectively, the above findings strongly suggest that analysing real sequential infection data using mathematical models will help infer the timing and strength of host immunity, which are difficult to measure directly in laboratory experiments. Such mathematical models will not only have the ability to explain observed experimental outcomes, but the ability to predict outcomes of new experiments, which can then be tested in the laboratory. These findings are particularly important as sequential infection experiments are increasingly being used to study the role of the immune response during infection with influenza and other respiratory pathogens [18, 26].

Limitations of sequential infection experiments

This study has highlighted some limitations of quantifying the immune response using virological data from sequential infection experiments alone.

Firstly, using the synthetic sequential infection data, we could not discriminate between the effects of cellular and humoral adaptive immunity in controlling a primary infection. If the effects of cellular and humoral adaptive immunity need to be distinguished, such as to predict the effects of vaccines that boost these components separately, quantities other than the viral load may need to be measured.

Secondly, we could not definitively rule out the possibility that target cell depletion contributed significantly to cross-protection. We were also unable to distinguish the roles of different innate immune mechanisms in cross-protection. Some modelling applications may require the strengths of different innate immune mechanisms to be known separately. An example of such an application is modelling the effect of treatments that modulate the innate immune response, such as the toll-like receptor-2 agonist Pam2Cys which has been shown to stimulate innate immune signals and reduce influenza-associated mortality and morbidity in animal studies [27].

In this simulation-based study, we were able to compare inferred quantities to a ‘ground truth’, to understand which quantities were inferred accurately, and which inferences may need to be treated with caution. For example, in S4 File, we show that the marginal posterior distributions for some parameters exhibit bias, such as that for the basic reproduction number R0. These apparent biases reflect correlations in the joint posterior distribution, and would be difficult to identify without a simulation-based study. This approach is thus crucial for sound interpretation of future studies fitting models to experimental data.

In addition to total viral load data, the study by Laurie et al. [17] also included infectious viral load measurements for single infection ferrets, and serological responses and cytokine levels at limited time points. Inclusion of this data could help to address the above limitations; the utility of this additional data can be assessed by further simulation-based studies.

New experiments could also be conducted to improve parameter estimates, leading to more accurate inference of the timing and strength of immune components. Previous studies have measured viral decay rates in vitro and incorporated these estimates into model fitting [28, 29]. In vitro studies can also directly measure the time course of those immune mechanisms which are active in vitro [30].

Future work and concluding summary

Now that we have shown how mathematical models can increase the utility of sequential infection experiments, fitting the model to the experimental data by Laurie et al. [17] is a priority. A simulation-estimation study alone cannot validate the mathematical model used, or infer the effects of host immunity against the pathogens in the experiments. However, this simulation-based study ensures that results will be interpreted appropriately when the models are fitted to data.

We have demonstrated that compared to single infection experiments, the sequential infection study design helps us to better understand cross-protection on short timescales. Further, data from sequential infection experiments helps to discriminate between existing models for a primary infection, leading to an improved understanding of the control and resolution of infection.

Materials and methods

The model

Viral dynamics.

The viral dynamics model is based on a model we previously published [24]. It incorporates three major components of the immune response—innate, humoral adaptive and cellular adaptive.

Fig 6 shows a compartmental diagram of the model for a single strain. The system is described by a coupled set of ordinary differential equations (Eqs 14). (1)

thumbnail
Fig 6. The within-host influenza model for a single strain.

(Top) Viral dynamics and innate immune response; (middle) humoral adaptive immune response; (bottom) cellular adaptive immune response. Solid arrows indicate transitions between compartments or death (shown only for immune-enhanced death processes); dashed arrows indicate production; plus signs indicate an increased transition rate due to the indicated compartment.

https://doi.org/10.1371/journal.pcbi.1006568.g006

Eq 1 describes the dynamics of target cells (T), infected cells (I) and virions (Vinf and Vtot for infectious and total virions respectively). Virions (Vinf) bind to target cells (T) to infect them; infected cells (I) produce virions; and infected cells and virions both decay at a constant rate. Target cells also regrow, with an imposed carrying capacity. Immunity is mediated by the compartments R (resistant cells), F (type I interferon), A (antibodies) and E (effector CD8+ T cells), the dynamics of which will be described shortly. Descriptions of model parameters are given in S1S4 Tables, and in our previous publication [24].

The compartment Vinf refers to the number of infectious virions in the host; however, an infected cell produces both infectious and non-infectious virions, the latter of which arise due to defects introduced during the viral replication process [31, 32]. Moreover, in the experiments conducted by Laurie et al. [17], the total viral load, rather than the number of infectious virions, was measured. An additional complication is that the concentration of total nasal wash, rather than the absolute number of virions, was measured. Hence, we include an equation for the total virion concentration Vtot.

Innate immunity is mediated through type I interferon (F), the production of which is stimulated by infected cells. Three effects of type I interferon are modelled through Eqs 1 and 2:

  • rendering target cells temporarily resistant to infection (TR);
  • decreasing the production rate of virions from infected cells; and
  • increasing the decay rate of infected cells.
(2)

The humoral adaptive immune response is mediated by antibodies (A), which bind to virions and neutralise them, rendering them non-infectious. Naive B cells (B0) are stimulated by virus to proliferate and differentiate into plasma cells (P), which produce antibodies. Eq 3 describes these processes. (3)

The cellular adaptive immune response is mediated by effector CD8+ T cells (E). Infected cells stimulate the differentiation of effector CD8+ T cells from their naive counterparts (C); effector CD8+ T cells then increase the death rate of infected cells. Some effector CD8+ T cells remain after a primary infection as memory CD8+ T cells. After a refractory period (represented by the M stage), they are modelled as having the same dynamics as naive cells, and can be re-stimulated to become effector CD8+ T cells upon challenge. Eq 4 describes these processes. (4) When more than one strain co-infects the host, the strains interact in three ways:

  • competition for target cells, which become depleted due to the infection and subsequent death of cells;
  • innate immunity, which acts across all strains; and
  • cellular adaptive immunity, which can be partly cross-reactive.

Activation of each of these mechanisms by the primary virus lowers the effective reproduction number of the challenge strain, but to different extents depending on parameter values. Note that because the model includes target cell regrowth, infection with the challenge virus can become established despite target cell depletion due to the death of infected cells. Each naive CD8+ T cell pool can be stimulated by one or more virus strains, depending on model parameters; cross-reactivity arises when a T cell pool can be stimulated by more than one virus strain. The clearance of infected cells by effector CD8+ T cells is similarly strain-specific. The antibody response is modelled as completely strain-specific, with no cross-reactivity between strains. It is thus unnecessary to include long-term humoral adaptive immunity. Extension of the model to include the potential effects of antibody-mediated cross-protection (as reviewed by [33]) is the subject of future work.

S4 Fig illustrates the model for two strains and three T cell pools; the equations are given in S1 File. Three T cell pools is a parsimonious choice, to allow for one pool to be cross-reactive between strains and two pools to be strain-specific, one for each strain.

Observation model.

Observations were simulated from the ‘true’ viral load by adding lognormal noise and imposing a detection threshold. Mathematically, the measured viral load yqfk for each virus q = 1, 2, …,Q, ferret f = 1, 2, …, F and measuring time point tqfk where k = 1, 2, …, Kqf is given by (5) β is a vector of parameter values, uf is the inter-exposure interval for ferret f, eqfk is the measurement error, and Θ is the detection threshold. Vtotq(tqfk, uf, β) is the solution to the two-strain version of Eqs 14 for the Vtotq compartment at time tqfk for the given parameter values and inter-exposure intervals. Θ takes the value 10 RNA copies/100μL in the experiments by Laurie et al. [17]. A measured viral load of 0 denotes that the viral load is below the detection threshold.

Therefore the likelihood of the model given the data is (6) In the second line of Eq 6, the likelihood when the data is below the detection threshold is obtained by integrating the probability density function from 0 to the detection threshold, i.e. treating the data below the threshold as censored [34]. The vector θ contains the parameters β, the inter-exposure intervals uf, the time points tqfk, and the standard deviation σ of the measurement error.

Simulated experiments

The model and the chosen ‘true’ parameters were used to generate synthetic data akin to that in Laurie et al. [17]. For six ferrets, intervals of 1, 3, 5, 7, 10 and 14 days separated exposures to two influenza strains. In addition, thirteen ferrets were exposed to a single influenza strain only. The sequential infection dataset consists of the viral load for the six sequential infection ferrets and one single infection ferret; the single infection dataset consists of the viral load for the thirteen single infection ferrets. The number of single infection ferrets was chosen such that the number of exposures to influenza virus is the same in each dataset, and so the number of data points is roughly the same.

Selection of model parameters to generate synthetic data

The ‘true’ parameter values chosen to generate the synthetic data are given in S1S4 Tables. The parameters were assumed to be identical between the two strains, except for the parameters governing cross-reactivity in the cellular adaptive immune response. In addition to the criteria discussed in the Results section, the parameters were chosen to reproduce qualitative behaviour for a single infection when immune components are suppressed:

  • when the innate adaptive immune response is absent (F → 0), the peak viral load increases [2];
  • when the humoral adaptive immune response is absent (A → 0), the viral load rebounds [3];
  • when the cellular adaptive immune response is absent (E → 0), resolution of the infection is delayed [4]; and
  • when both arms of the adaptive immune response are absent (A, E → 0), chronic infection ensues [5].

For an extensive evaluation of a very similar model’s behaviour under these types of conditions (for single infection events), see [22].

In addition to parameter values, initial values were required when simulating infections. For a single infection, the initial values for all compartments in Eqs 14 except T, Vinf, Vtot, C and B0 were zero. The initial values of C and B0 (naive T and B cells respectively) were normalised to 1. The initial values of T and Vinf (the number of target cells and the concentration of infectious virions respectively) were estimated parameters. The initial concentration of total virus was then Vtot(0) = γαVinf(0), where γ and α were conversion parameters described in S1 Table. For sequential infections, the conditions at the time of the primary exposure were as above; the system was integrated until the time of the challenge exposure, at which Vinf,2(0) infectious virions for the challenge strain was added to the system, and the total concentration of the challenge strain was set to Vtot,2(0).

Model fitting

Parameters to be estimated.

All model parameters were estimated, except for the following parameters which were either fixed or a function of other estimated parameters. We fixed two parameters—the number of plasmablast division cycles (nB) and the number of effector CD8+ T cell division cycles (nE)—to be 5 [35, 36] and 20 [37] respectively. In addition, when fitting the model to single infection data, we considered the optimistic scenario where we had independent, perfect information about the proportion of cellular adaptive immunity during a primary infection that was cross-reactive with the challenge strain. As one T cell pool was cross-reactive between strains and two pools were strain-specific, this amounted to fixing the proportion of cellular adaptive immunity attributed to each T cell pool. We did so by fixing the numbers of infected cells for half-maximal stimulation of naive CD8+ T cells kCj1 to their ‘true’ values for each T cell pool j. Then when we extended the model to two strains, we set kCjq to these same ‘true’ values. We then calculated the clearance rates of infected cells by effector CD8+ T cells κEjq by taking the fitted value of κE11, and applying the formula κEjq = κE11 kC11/kCjq (see S1 File).

Instead of fitting the infectivity (β) and the production rate of infectious virions from an infected cell (pVinf), we fitted the basic reproduction number R0 (Eq 7) and the initial viral load growth rate r (Eq 9), as we hypothesised that these were more intimately linked to features of the viral load curve. Practically speaking, we proposed a new value for R0 (or r), calculated the corresponding values of β and pVinf, solved the model equations, calculated the likelihood of the data given the parameters, and accepted or rejected the new value for R0 (or r).

The basic reproduction number R0 is the mean number of secondary infected cells due to (the virions produced by) a single infected cell. The expression for R0 is (7) and is the same as that for a model without a time-dependent immune response [38].

The viral load during early infection can be approximated by (8) Arenas et al. [39] showed using a simulation-estimation study that this parameter was well estimated even when only viral load data was available.

The expression for r, derived by linearising Eq 1 around the disease-free equilibrium [40], is (9)

Prior distributions.

We began with a uniform distribution in parameter space whose bounds along each dimension are given in S1S4 Tables. Note that parameter estimation was performed in a parameter space where all parameters except for the standard deviation of the measurement error σ were log transformed. Then, we excluded regions of parameter space where the parameters log10 β and log10 pVinf, which were not directly estimated but were instead recovered from Eqs 7 and 9, were outside the bounds given in S1 Table.

The priors were deliberately chosen to be wide because previous parameter estimates came from a range of experimental systems, and parameters with similar physical definitions could vary in value depending on the model used. The bounds for viral replication parameters were based on those by Petrie et al. [41] where the equivalent parameters exist. Otherwise, where multiple estimates existed in the literature (as cited in the tables), the bounds were chosen to encompass all of them. Where we could only find a single estimate, bounds spanning at least an order of magnitude were chosen (unless the parameter is a pure rate parameter, as discussed shortly). Where no estimate was found, we assigned very wide bounds spanning much more than one order of magnitude. In general, the bounds for pure rate parameters (those with units day−1 only) were chosen to be narrower as their order of magnitude was known, whereas bounds for parameters such as R0 were much wider.

Furthermore, for computational efficiency, some minimal constraints on the behaviour of the viral load and timing of various immune components were incorporated into the prior distribution. These constraints were imposed because parameter sets that generate ‘unreasonable’ viral load trajectories for a single infection caused large delays in numerical integration of the two-strain differential equations. The inclusion criteria were that for a single infection,

  • the total viral load rises by at least one order of magnitude during infection;
  • the total viral load peaks 0–7 days post-exposure;
  • the humoral adaptive immune response is not active too early (five days post-exposure, the neutralisation rate of virus by antibodies, κAA, does not exceed 103 day−1); and
  • the cellular adaptive immune response is not active too early (five days post-exposure, the clearance rate of infected cells by effector CD8+ T cells, , does not exceed 103 day−1).

If the viral load trajectory (in the absence of measurement noise) predicted by a parameter set does not fulfil all of these conditions, the value of the prior distribution is zero at that point in parameter space.

Model fitting algorithm.

We fitted the model using the Metropolis algorithm [42, 43] embedded within a Gibbs sampler structure [44], implemented in Octave 3.8.2 [45]. To evaluate the likelihood, Eqs 14 were solved using the CVODE solvers developed by [46], implemented in MATLAB [47]. Of the available solvers, a backward differentiation formula method in variable order, variable step, fixed leading coefficient form was chosen. Extinction was enforced by defining an infection to have resolved if both the number of infected cells and virions was below 0.1.

To assess convergence, three chains were run in parallel using different starting parameter values θ0 drawn from the prior distribution. The procedure for determining the number of iterations for which to run the chains is detailed in S1 Text. For efficient mixing, the proposal distributions were tuned such that the proportion of accepted proposals was not too low or too high, as detailed in S1 Text. For each of the three chains, parallel tempering (as developed by [48] and reviewed by [49]) was implemented to improve exploration of parameter space. The number of iterations before testing whether to swap chains in the parallel tempering process was set to 10. During the calibration period for the proposal distributions, the temperatures were also calibrated [50], as detailed in S1 Text. Once convergence was reached, the effective sample size was calculated for each chain (using the iterations that were kept following the burn-in process) using the effectiveSize function in the coda [51] package in R [52]. Convergence diagnostics for the chains are shown in S3 File.

The marginal posterior distributions in this study are plotted in S4 File using all samples from the chains (after burn-in), without thinning. When using the joint posterior distribution to make predictions, we used 104 parameter sets corresponding to uniformly spaced iterations in each of the chains.

Results were visualised using MATLAB R2015b [53]. Code to reproduce all figures is provided at https://bitbucket.org/ada_yan/sim_est_immunity.

Model predictions

First, to determine whether the fitted model captured the timing and strength of each immune component during a primary infection, we used parameter sets from the joint posterior distribution to simulate the viral load during a single infection, using a modified model where either

  • adaptive immunity is suppressed;
  • innate immunity is suppressed;
  • innate and adaptive immunity are suppressed;
  • cellular adaptive immunity is suppressed; or
  • humoral adaptive immunity is suppressed.
95% prediction intervals were constructed using these simulations.

Second, to determine whether the fitted model captured cross-protection between strains, we used parameter sets from the joint posterior distribution to simulate different inter-exposure intervals.

Third, to determine whether the fitted model captured the contribution of each immune component to cross-protection between strains, we used parameter sets from the joint posterior distribution to simulate the viral load during sequential infection, using a modified model where either

  • cross-protection is only mediated by cellular adaptive immunity, and not target cell depletion or innate immunity (model XC);
  • cross-protection is mediated by innate immunity, but not target cell depletion or cellular adaptive immunity (model XI);
  • cross-protection is mediated by target cell depletion, but not innate immunity or cellular adaptive immunity (model XT); or
  • cross-protection is mediated by target cell depletion and/or innate immunity, but not cellular adaptive immunity (model XIT).

Details of models XC, XI, XT and XIT are provided in S1 File.

Table 1 summarises the model modifications in this section.

Supporting information

S1 Fig. The full set of synthetic data.

(a) The line shows the simulated ‘true’ viral load for a single infection, with the arrow showing the time of exposure. The simulated viral loads with noise for the thirteen single infection ferrets are shown as crosses. The horizontal line indicates the observation threshold (10 RNA copy no./100μL); observations below this threshold are plotted below this line. Values below the observation threshold were treated as censored. (b—g) For sequential infections with the labelled inter-exposure interval, the dashed and dotted lines show the simulated ‘true’ viral load for a primary and challenge infection respectively; the arrows show the times of the primary and challenge exposures. The simulated viral load with noise is shown as crosses. The sequential infection dataset consists of the viral load for the six sequential infection ferrets and one single infection ferret; the single infection dataset consists of the viral load for the thirteen single infection ferrets.

https://doi.org/10.1371/journal.pcbi.1006568.s001

(PDF)

S2 Fig. Trajectories for models (a) XI and (b) XT generated using 100 uniformly sampled parameter sets from the MCMC chains after burn-in, for the model fitted to sequential infection data.

The green trajectory incorrectly attributed the delay observed in the baseline model to both target cell depletion and innate immunity.

https://doi.org/10.1371/journal.pcbi.1006568.s002

(PDF)

S3 Fig. Sequential infection data did not enable accurate prediction of the challenge viral load for a modified model where only one of the three innate immune mechanisms mediates cross-protection.

The challenge viral load for the ‘true’ parameter values and a modified model where cross-protection is mediated by only one innate immune mechanism (models XI1–XI3, red line) was compared to the viral load for the baseline model (black line). (a—c) show results for models XI1–XI3 respectively. At a one-day inter-exposure interval, the delay in the baseline model occurred due to a combination of innate immune mechanisms 2 and 3. Prediction intervals for the viral load for models XI1–XI3 according to the model fitted to sequential infection data (blue areas) did not accurately recover the viral load according to the ‘true’ parameters. Hence, the fitted model did not attribute cross-immunity to the correct mechanisms of the innate immune response.

https://doi.org/10.1371/journal.pcbi.1006568.s003

(PDF)

S4 Fig. Compartmental diagram for two strains and three T cell pools.

Cells infected with influenza strain 1 stimulate naive CD8+ T cells in pools 1 and 3, and are cleared by effector CD8+ T cells in these pools. Cells infected with influenza strain 2 stimulate naive CD8+ T cells in pools 2 and 3, and are cleared by effector CD8+ T cells in these pools.

https://doi.org/10.1371/journal.pcbi.1006568.s004

(PDF)

S1 Text. More details on the model fitting procedure.

https://doi.org/10.1371/journal.pcbi.1006568.s005

(PDF)

S2 Text. Notes on biologically plausible ranges for the parameters pVratio, α and γ, as provided in S1 Table.

https://doi.org/10.1371/journal.pcbi.1006568.s006

(PDF)

S1 File. Two-strain model equations for the baseline and modified models, model modifications from our previous study, and compartmental diagrams for the modified models.

https://doi.org/10.1371/journal.pcbi.1006568.s007

(PDF)

S2 File. Results for an additional set of parameters where the degree of cross-reactivity in the cellular adaptive immune response is high.

https://doi.org/10.1371/journal.pcbi.1006568.s008

(PDF)

S3 File. Convergence diagnostics for the MCMC chains.

https://doi.org/10.1371/journal.pcbi.1006568.s009

(PDF)

S4 File. Marginal posterior distributions for the model parameters.

https://doi.org/10.1371/journal.pcbi.1006568.s010

(PDF)

S5 File. Results for a different set of ‘true’ parameters where the degree of cross-reactivity in the cellular adaptive immune response is low.

https://doi.org/10.1371/journal.pcbi.1006568.s011

(PDF)

S6 File. Results for different noisy datasets with the same ‘true’ parameters.

https://doi.org/10.1371/journal.pcbi.1006568.s012

(PDF)

S7 File. Results for a data set generated using a different model.

https://doi.org/10.1371/journal.pcbi.1006568.s013

(PDF)

S1 Table. Viral replication parameter values and prior bounds.

Note that the values and prior bounds are given in logarithmic space. For example, the value of log10 R0 was log104.9 and the prior bounds were [0, 3]. Hence, the value of R0 was 4.9 and the prior bounds of R0 were [1, 1000]. β and pVinf were not directly fitted, but their values as recovered from Eqs 7 and 9 could not exceed the bounds given. Because total virions include infectious virions, the total virion decay rate should be slower than the infectious virion decay rate. Hence, the difference between the infectious and total virion decay rates δVinfδVtot, rather than the infectious virion decay rate δVinf, was fitted to ensure that the former quantity was positive. Notes on biologically plausible ranges for the parameters ptot, α and γ are given in S2 Text.

https://doi.org/10.1371/journal.pcbi.1006568.s014

(PDF)

S2 Table. Innate immune response parameter values and prior bounds.

https://doi.org/10.1371/journal.pcbi.1006568.s015

(PDF)

S3 Table. Values and prior bounds for the cross-reactivity parameters in the cellular adaptive immune response.

The number of infected cells for half-maximal stimulation of naive/memory CD8+ T cells kCjq and the clearance rate of infected cells by effector CD8+ T cells κE11.

https://doi.org/10.1371/journal.pcbi.1006568.s016

(PDF)

S4 Table. Adaptive immune response and observation model parameter values and prior bounds.

https://doi.org/10.1371/journal.pcbi.1006568.s017

(PDF)

Acknowledgments

We thank Karen L. Laurie for designing and performing the experiments modelled, and providing virological insight. We also acknowledge helpful discussions with Patricia T. Campbell, Pengxing Cao, Steven Riley and Alexander E. Zarebski.

References

  1. 1. Dobrovolny HM, Reddy MB, Kamal MA, Rayner CR, Beauchemin CAA. Assessing mathematical models of influenza infections using features of the immune response. PLoS ONE. 2013;8(2):e57088. pmid:23468916
  2. 2. Seo SH, Hoffmann E, Webster RG. Lethal H5N1 influenza viruses escape host anti-viral cytokine responses. Nat Med. 2002;8(9):950–4. pmid:12195436
  3. 3. Iwasaki T, Nozima T. Defense mechanisms against primary influenza virus infection in mice I. The roles of interferon and neutralizing antibodies and thymus dependence of interferon and antibody production. J Immunol. 1977;118(1):256–263. pmid:401512
  4. 4. Yap KL, Ada GL. Cytotoxic T cells in the lungs of mice infected with an influenza A virus. Scand J Immunol. 1978;7(1):73–80. pmid:305613
  5. 5. Kris RM, Yetter RA, Cogliano R, Ramphal R, Small PA. Passive serum antibody causes temporary recovery from influenza virus infection of the nose, trachea and lung of nude mice. Immunology. 1988;63(3):349–53. pmid:2832312
  6. 6. Baccam P, Beauchemin C, Macken CA, Hayden FG, Perelson AS. Kinetics of influenza A virus infection in humans. J Virol. 2006;80(15):7590–7599. pmid:16840338
  7. 7. Handel A, Longini IM Jr, Antia R. Towards a quantitative understanding of the within-host dynamics of influenza A infections. J R Soc Interface. 2010;7(42):35–47. pmid:19474085
  8. 8. Pawelek KA, Huynh GT, Quinlivan M, Cullinane A, Rong L, Perelson AS. Modeling within-host dynamics of influenza virus infection including immune responses. PLoS Comput Biol. 2012;8(6):e1002588. pmid:22761567
  9. 9. Saenz RA, Quinlivan M, Elton D, MacRae S, Blunden AS, Mumford JA, et al. Dynamics of influenza virus infection and pathology. J Virol. 2010;84(8):3974–3983. pmid:20130053
  10. 10. Hancioglu B, Swigon D, Clermont G. A dynamical model of human immune response to influenza A virus infection. J Theor Biol. 2007;246(1):70–86. pmid:17266989
  11. 11. Bocharov GA, Romanyukha AA. Mathematical model of antiviral immune response III. Influenza A virus infection. J Theor Biol. 1994;167(4):323–360. pmid:7516024
  12. 12. Miao H, Hollenbaugh JA, Zand MS, Holden-Wiltse J, Mosmann TR, Perelson AS, et al. Quantifying the early immune response and adaptive immune response kinetics in mice infected with influenza A virus. J Virol. 2010;84(13):6687–6698. pmid:20410284
  13. 13. Lee HY, Topham DJ, Park SY, Hollenbaugh J, Treanor J, Mosmann TR, et al. Simulation and prediction of the adaptive immune response to influenza A virus infection. J Virol. 2009;83(14):7151–7165. pmid:19439465
  14. 14. Li Y, Handel A. Modeling inoculum dose dependent patterns of acute virus infections. J Theor Biol. 2014;347:63–73. pmid:24440713
  15. 15. Ahmed H, Moore J, Manicassamy B, Garcia-Sastre A, Handel A, Antia R. Mathematical analysis of a mouse experiment suggests little role for resource depletion in controlling influenza infection within host. ArXiv e-prints. 2017;.
  16. 16. Handel A, Liao LE, Beauchemin CAA. Progress and trends in mathematical modelling of influenza A virus infections. Curr Opin Syst Biol. 2018; https://doi.org/10.1016/j.coisb.2018.08.009.
  17. 17. Laurie KL, Guarnaccia TA, Carolan LA, Yan AWC, Aban M, Petrie S, et al. Interval between infections and viral hierarchy are determinants of viral interference following influenza virus infection in a ferret model. J Infect Dis. 2015;212(11):1701–1710. pmid:25943206
  18. 18. Laurie KL, Horman W, Carolan LA, Chan KF, Layton D, Bean A, et al. Evidence for viral interference and cross-reactive protective immunity between influenza B virus lineages. J Infect Dis. 2018;217(4):548–559. pmid:29325138
  19. 19. Carolan LA, Rockman S, Borg K, Guarnaccia T, Reading P, Mosse J, et al. Characterization of the localized immune response in the respiratory tract of ferrets following infection with influenza A and B viruses. J Virol. 2016;90(6):2838–2848.
  20. 20. Cao P, Yan AWC, Heffernan JM, Petrie S, Moss RG, Carolan LA, et al. Innate immunity and the inter-exposure interval determine the dynamics of secondary influenza virus infection and explain observed viral hierarchies. PLoS Comput Biol. 2015;11(8):e1004334. pmid:26284917
  21. 21. Yan AWC, Cao P, McCaw JM. On the extinction probability in models of within-host infection: the role of latency and immunity. J Math Biol. 2016;73(4):787–813. pmid:26748917
  22. 22. Cao P, Wang Z, Yan AWC, McVernon J, Xu J, Heffernan JM, et al. On the role of CD8+ T cells in determining recovery time from influenza virus infection. Front Immunol. 2016;7:611. pmid:28066421
  23. 23. Hoshino A, Takenaka H, Mizukoshi O, Imanishi J, Kishida T, Tovey MG. Effect of anti-interferon serum of influenza virus infection in mice. Antiviral Research. 1983;3(1):59–65. pmid:6191656
  24. 24. Yan AWC, Cao P, Heffernan JM, McVernon J, Quinn KM, La Gruta NL, et al. Modelling cross-reactivity and memory in the cellular adaptive immune response to influenza infection in the host. J Theor Biol. 2017;413:34–49. pmid:27856216
  25. 25. Zarnitsyna VI, Handel A, McMaster SR, Hayward SL, Kohlmeier JE, Antia R. Mathematical model reveals the role of memory CD8 T cell populations in recall responses to influenza. Front Immunol. 2016;7(May):165. pmid:27242779
  26. 26. Chan KF, Carolan LA, Korenkov D, Druce J, McCaw J, Reading PC, et al. Investigating viral interference between influenza A Virus and human respiratory syncytial virus in a ferret model of infection. J Infect Dis. 2018; p. jiy184.
  27. 27. Tan ACL, Mifsud EJ, Zeng W, Edenborough K, McVernon J, Brown LE, et al. Intranasal administration of the TLR2 agonist Pam2Cys provides rapid protection against influenza in mice. Mol Pharm. 2012;9(9):2710–2718. pmid:22823162
  28. 28. Pinilla LT, Holder BP, Abed Y, Boivin G, Beauchemin CAA. The H275Y neuraminidase mutation of the pandemic A/H1N1 influenza virus lengthens the eclipse phase and reduces viral output of infected cells, potentially compromising fitness in ferrets. J Virol. 2012;86(19):10651–10660. pmid:22837199
  29. 29. Paradis EG, Pinilla LT, Holder BP, Abed Y, Boivin G, Beauchemin CAA. Impact of the H275Y and I223V mutations in the neuraminidase of the 2009 pandemic influenza virus in vitro and evaluating experimental reproducibility. PLoS ONE. 2015;10(5):e0126115. pmid:25992792
  30. 30. Mitchell H, Levin D, Forrest S, Beauchemin CAA, Tipper J, Knight J, et al. Higher level of replication efficiency of 2009 (H1N1) pandemic influenza virus than those of seasonal and avian strains: kinetics from epithelial cell culture and computational modeling. J Virol. 2011;85(2):1125–35. pmid:21068247
  31. 31. Nayak DP, Chambers TM, Akkina RK. In: Cooper M, Eisen H, Goebel W, Hofschneider PH, Koprowski H, Melchers F, et al., editors. Defective-interfering (DI) RNAs of influenza viruses: origin, atructure, expression, and interference. Berlin, Heidelberg: Springer Berlin Heidelberg; 1985. p. 103–151.
  32. 32. Marriott AC, Dimmock NJ. Defective interfering viruses and their potential as antiviral agents. Rev Med Virol. 2010;20(1):51–62. pmid:20041441
  33. 33. Ekiert DC, Wilson IA. Broadly neutralizing antibodies against influenza virus and prospects for universal therapies. Curr Opin Virol. 2012;2(2):134–141. pmid:22482710
  34. 34. Ahn JE, Karlsson MO, Dunne A, Ludden TM. Likelihood based approaches to handling data below the quantification limit using NONMEM VI. J Pharmacokinet Pharmacodyn. 2008;35(4):401–421. pmid:18686017
  35. 35. Marchuk GI, Petrov RV, Romanyukha AA, Bocharov GA. Mathematical model of antiviral immune response. I. Data analysis, generalized picture construction and parameters evaluation for hepatitis B. J Theor Biol. 1991;151(1):1–40. pmid:1943135
  36. 36. Sze DMY, Toellner KM, García de Vinuesa C, Taylor DR, MacLennan ICM. Intrinsic constraint on plasmablast growth and extrinsic limits of plasma cell survival. J Exp Med. 2000;192(6):813–821. pmid:10993912
  37. 37. van Stipdonk MJB, Lemmens EE, Schoenberger SP. Naïve CTLs require a single brief period of antigenic stimulation for clonal expansion and differentiation. Nat Immunol. 2001;2(5):423–429. pmid:11323696
  38. 38. Beauchemin CAA, McSharry JJ, Drusano GL, Nguyen JT, Went GT, Ribeiro RM, et al. Modeling amantadine treatment of influenza A virus in vitro. J Theor Biol. 2008;254(2):439–451. pmid:18653201
  39. 39. Arenas AR, Thackar NB, Haskell EC. The logistic growth model as an approximating model for viral load measurements of influenza A virus. Math Comput Simul. 2017;133:206–222.
  40. 40. Nowak MA, Lloyd AL, Vasquez GM, Wiltrout TA, Wahl LM, Bischofberger N, et al. Viral dynamics of primary viremia and antiretroviral therapy in simian immunodeficiency virus infection. J Virol. 1997;71(10):7518–25. pmid:9311831
  41. 41. Petrie SM, Butler J, Barr IG, McVernon J, Hurt AC, McCaw JM. Quantifying relative within-host replication fitness in influenza virus competition experiments. J Theor Biol. 2015;382:259–271. pmid:26188087
  42. 42. Metropolis N, Ulam S. The Monte Carlo method. JASA. 1949;44(247):335–341. pmid:18139350
  43. 43. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. J Chem Phys. 1953;21(6):1087–1092.
  44. 44. Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell. 1984;PAMI-6(6):721–741.
  45. 45. Eaton John W H S Bateman David, Wehbring R. GNU Octave version 3.8.1 manual: a high-level interactive language for numerical computations. CreateSpace Independent Publishing Platform; 2014. Available from: http://www.gnu.org/software/octave/doc/interpreter.
  46. 46. Cohen SD, Hindmarsh AC, Dubois PF. CVODE, a stiff/nonstiff ODE solver in C. Computers in Physics. 1996;10:138–148.
  47. 47. Vanlier J, Tiemann CA, Hilbers PAJ, van Riel NAW. A Bayesian approach to targeted experiment design. Bioinformatics. 2012;28(8):1136–1142. pmid:22368245
  48. 48. Geyer CJ. Markov chain Monte Carlo maximum likelihood. In: Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface. Interface Foundation of North America; 1991. p. 156.
  49. 49. Earl DJ, Deem MW. Parallel tempering: theory, applications, and new perspectives. Phys Chem Chem Phys. 2005;7:3910–3916. pmid:19810318
  50. 50. Kone A, Kofke DA. Selection of temperature intervals for parallel-tempering simulations. J Chem Phys. 2005;122(20):206101. pmid:15945778
  51. 51. Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6(1):7–11.
  52. 52. R Core Team. R: A Language and Environment for Statistical Computing; 2016. Available from: https://www.R-project.org/.
  53. 53. MATLAB. R2015b. Natick, Massachusetts: The MathWorks Inc.; 2015.