Abstract
Although the influenza A virus has been extensively studied, a quantitative understanding of the infection dynamics is still lacking. To make progress in this direction, we designed several mathematical models and compared them with data from influenza A infections of mice. We find that the immune response (IR) plays an important part in the infection dynamics. Both an innate and an adaptive IR are required to provide adequate explanation of the data. In contrast, regrowth of epithelial cells did not seem to be an important mechanism on the time scale of the infection. We also find that different model variants for both innate and adaptive responses fit the data well, indicating the need for additional data to allow further model discrimination.
1. Introduction
A quantitative understanding of infection dynamics is important for improved control of infectious diseases. To that end, mathematical models, combined with experimental data, can provide valuable insights. While some pathogens have been studied quantitatively (Nowak & May 2001; Perelson 2002; Dixit et al. 2004; Kirschner & Marino 2005; Asquith & Bangham 2007; Davenport et al. 2007), influenza has received little attention.
Seasonal influenza usually causes uncomplicated and transient infections in humans, with virus replication localized to the upper respiratory tract (URT). Two recent studies used viral load data from human volunteers infected with influenza and combined them with mathematical models to quantify the infection dynamics (Baccam et al. 2006; Handel et al. 2007). These studies showed that it was possible to describe the infection dynamics without the need to consider the immune response (IR). Instead, the decline of viral load after a few days could be attributed solely to the depletion of target cells, which are primarily epithelial cells lining the URT. It is quite possible that the virus dynamics in the URT is driven mainly by depletion of target cells. For instance, Francis & StuartHarris (1938) found that ferrets infected intranasally with a sublethal dose of influenza virus developed desquamation of the tracheal area by day 2 with complete destruction of the epithelium. The animals survived and fully regenerated the epithelial tissue within a few weeks.
On the other hand, reports from immunocompromised humans who shed influenza virus for prolonged periods suggest that the IR plays an important role in clearing the infection, or at least in preventing it from becoming chronic and potentially lethal (Rocha et al. 1991; Klimov et al. 1995; Boivin et al. 2002; Weinstock et al. 2003). The IR is likely to be especially important in more severe influenza infections of the lower respiratory tract (LRT). Such infections can lead to viral pneumonia and in the worst cases to death. Humans infected with the H5N1 avian strain often show such LRT infections (Tran et al. 2004; de Jong et al. 2006). Similarly, some of the hosts that died during the devastating 1918 pandemic seem to have succumbed to a viral pneumonia (the majority likely died due to secondary bacterial pneumonia; Morens & Fauci 2007; Morens et al. 2008). Autopsies show that severe influenza infections often involve significant damage of the LRT, which presumably contributes to the host's death (Giles & Shuttleworth 1957; Hers et al. 1958). Therefore, in LRT influenza infections, an IR that can quickly suppress the virus, without causing too much immunopathology, seems crucial.
Unfortunately, there is little kinetic data beyond virus load available from infected humans, and often the data are reported in a form that make them unsuitable for detailed quantitative studies. In fact, the mathematical studies for human URT influenza infections mentioned above (Baccam et al. 2006; Handel et al. 2007) showed that models that included either an innate (Baccam et al. 2006) or adaptive (Handel et al. 2007) IR component could also explain the observed virus dynamics; the available data were not sufficient to properly discriminate between models with and without an IR. In contrast, animal studies can usually provide more data. For instance, mice infected with influenza provide a good model for the more severe form of influenza pneumonia and have been extensively studied. Most of these studies point towards the importance of a functioning IR. However, the relative roles and contributions of the different components of the IR are less clear (Swain et al. 2004; Tamura & Kurata 2004; Doherty et al. 2006; Thomas et al. 2006). To make further progress towards a quantitative understanding of influenza A infection dynamics, we decided to fit mathematical models to some of the data from influenza A infections in mice. We used data from two experimental studies (Iwasaki & Nozima 1977; Kris et al. 1988) and fitted them to several mathematical models that describe the withinhost dynamics of the virus, target cells and different aspects of the IR. We find that the IR plays an important part in the infection dynamics. Both an innate and an adaptive IR are required to provide adequate explanation of the data.
2. Material and methods
2.1. Experimental data
We use data obtained from two experimental studies. The first experimental study by Kris et al. (1988) contains viral load data for primary infections with the H3N2 influenza strain A/Port Chalmers/1/73. In this study, both wildtype BALB/c mice with a functioning IR, as well as nude, athymic (nu/nu) mice were infected. The lack of a thymus in the nude mice leads to missing Tcells, which also leads to a largely nonfunctioning Bcell/antibody response. Kris et al. (1988) reported levels of virus load, given in 50 per cent egg infectious doses (in the following abbreviated as EID) for both the wildtype and nude mice.
The second study by Iwasaki & Nozima (1977; in the following abbreviated as IN) contains data from naive C3H/He mice infected with the H1N1 influenza strain A/PR8. The investigators reported viral titre, measured in 50 per cent plaque forming units (in the following abbreviated as PFU). In addtion, IN report cell damage as measured by lung lesions, levels of IgA, IgM and IgG antibody titres and interferon (IFN) levels (type I and II) for both serum and tracheobronchial washings. We make use of the tracheobronchial data from the experiments with immunecompetent mice and mice depleted of IgM antibodies (figures 1, 2, 7 and 8 in Iwasaki & Nozima (1977)). For additional details on the experimental procedures and data, we refer to the original studies (Iwasaki & Nozima 1977; Kris et al. 1988). Data were extracted from the published papers using Engauge Digitizer (digitizer.sourceforge.net).
2.2. Mathematical model
We employ mathematical models based on ordinary differential equations. The simplest such model is the one that describes the dynamics of uninfected and infected target (epithelial) cells and virus. Such models have been used extensively to study HIV (Perelson 2002) and more recently for influenza A (Baccam et al. 2006; Handel et al. 2007; Beauchemin et al. 2008). More detailed descriptions of these types of models can be found in these references. Briefly, uninfected epithelial cells, U, can become infected by virus at a rate b. This leads to latently infected cells, E, which after a mean time of 1/g hours start to produce virions, V, at rate p. Virus producing, infected cells, I, die after some time, 1/d. We ignore turnover of uninfected epithelial cells. However, we include the replacement of dead cells, D, by new susceptible epithelial cells at a rate λ. Free virus is cleared by nonspecific mechanisms (e.g. mucosal transport) at rate c.
In addition, we add equations for the IR. After infection, the innate IR is the first to respond. The innate response consists of both inflammatory cytokines and populations of cells such as neutrophils and macrophages (Tamura & Kurata 2004; Maines et al. 2008). We assume that such an innate response, F, is triggered upon infection and increases proportionally to free virus at rate w, with a firstorder clearance rate, δ. The impact of the innate IR on the virus is likely multifactorial. For instance, it is known that inflammatory cytokines render cells less susceptible to infection or reduce the production of virus in already infected cells. Alternatively, innate cellular components can increase the death of infected cells or clearance of infectious virions. In the main text, we assume that the innate IR reduces the production of infectious virions; some alternatives are discussed in appendix A.
While the innate IR is triggered rapidly, the adaptive response takes longer to reach high levels. Both the humoral component (Bcells, antibodies) and the cellular component (Tcells) of the adaptive IR have been found to play a role during influenza infections. Here, we focus on the humoral component; we briefly discuss Tcellbased responses in appendix A. The exact dynamics of Bcells/antibodies following infection is not well known. Activation of Bcells likely depends on viral load (antigen), while some of the expansion dynamics is probably independent of virus load, as has been found for Tcells (McHeyzerWilliams et al. 2006; Harty & Badovinac 2008; Batista & Harwood 2009). We model this by assuming that the adaptive response, X, is activated proportional to free virus load at rate f. Activation is followed by antigenindependent clonal expansion at an effective growth rate r. We only consider the expansion phase of the IR. This is justified for acute viral infections, since the peak and following contraction of the adaptive IR usually occurs days or weeks after the virus has been cleared, and we are only interested in the infection dynamics up to the point where the virus is cleared. When fitting to data, we consider X to represent antibodies, which we assume are proportional to the number of responding Bcells. The main mode of action for antiinfluenza antibodies is to neutralize the virus, thereby removing free infectious virions (Outlaw & Dimmock 1991). (Note that the experimental protocols only measure such infectious virions, consistent with our model.) We model neutralization as a massaction clearance term of free virions at rate k. We discuss several different model variants in appendix A.
With these assumptions, the model equations are given byTable 1 summarizes the model parameters. The factor γ in the equation for virions is needed for the conversion from infectious virions, the units of V in our model, to EID or PFU as measured by the experimental data. A value of γ = 1 assumes that one virion corresponds to one EID or PFU, γ < 1 allows for the possibility that more than one virion is required to produce one EID or PFU.
An important quantity for models such as this one is the basic reproductive number. This number, given by2.1describes the average number of infectious progeny virions produced by one virion at the start of the infection. A value of R_{0} > 1 leads to virus growth while a value lower than 1 will not.
2.3. Parameter bounds
Since initial fitting of the models to the available data indicated that the number of free parameters leads to degeneracies, we decided to fix several model parameters. We mention briefly at the appropriate places how results change if we fit these parameters instead of fixing them. Further, to ensure that our models are in agreement with known biology of the infection process, we imposed bounds on some of the fitted model parameters, based on estimates from the literature.
Experimental data suggest that production of virions starts as early as around 3 h after a cell becomes infected (Cairns & Destgroth 1957; White et al. 1965) and seems to be well under way at around 8 h postinfection (Henle et al. 1947; Cairns et al. 1952; Rimmelzwaan et al. 1998; Sidorenko & Reichl 2004). We therefore decided to fix the duration of the latent phase at 6 h (g = 4 per day), which also agrees with the best fit value from Baccam et al. (2006). Similarly, several studies indicate that cells die around 12–24 h after infection (Price et al. 1997; Brydon et al. 2003; Beauchemin et al. 2008). We therefore fixed the lifespan of productively infected cells at 12 h (d = 2 per day), which gives a total lifespan for infected cells (latent + productively infected phase) of 18 h.
Previous in vitro studies have reported the lifespan of virions to be several hours (Horsfall 1954; Choppin 1963), with a recent study suggesting a halflife of about 6 h (Beauchemin et al. 2008). It is not clear how that translates to the in vivo situation. On the one hand, virions could be optimized for the in vivo environment and survive longer. On the other hand, nonspecific clearance mechanisms, such as mucosal transport, likely reduce the effective halflife of virions significantly. We decided to use a value for the virion clearance rate of c = 10 per day, in line with previous estimates.
Several studies have shown that regeneration of the epithelium of the URT starts a few days after the onset of infection (Ramphal et al. 1979). It is not clear if the same is true for LRT epithelial cells, and if the newly generated epithelial cells are susceptible to virus infection right away (Francis & StuartHarris 1938; StuartHarris & Francis 1938). For those models where we include regeneration of epithelium, we bound the rate for λ between 0.5 and 0 per day.
We set no bounds on b and p since these are hard to estimate, especially given that the virus load is expressed in units of PFU (Iwasaki & Nozima 1977) or EID (Kris et al. 1988), which leads to b and p also having those units. We also do not enforce bounds on any of the parameters describing the IRs, since those are not well known.
Additionally, we need to specify or fit the initial number of target cells, U_{0}, and viral load, V_{0} (all other variables are initially zero). We estimate the number of target cells to be U_{0} = 7 × 10^{9}. This number is based on an estimated surface area of the LRT of mice of about 2175 cm^{2} (Ito et al. 2003) and an estimated epithelial cell surface area of about 3 × 10^{−7} cm^{2} (Farmer & Hay 1991). It is worth mentioning that a recent study based on an alternative approach (estimation of cell numbers based on measured protein levels) led to an estimate for the target cell population that is about 1000fold lower than our estimate (P. Thomas, P. Dash and C. Sanders 2009, personal communication). While a more precise measurement of the target cell population would be valuable, it fortunately does not matter much for our study. All results shown below are also valid for different number of target cells, the only difference would be a rescaling of the parameters p and γ.
For the viral load, only the study by IN (1977) reports the amount of virus in the inoculum. Further, it has been shown that the ‘take’ can be lower and depends on the details of the inoculation procedure (Yetter et al. 1980). We therefore treat the inoculum, V_{0}, as a free parameter and fit it to the data. However, we impose a lower bound of hundred infectious virions and upper bounds as described in the sections below.
2.4. Fitting procedure
We fitted the data using nonlinear least square fit routines, i.e. we minimized ∑_{i}(y_{i}^{d} − y_{i}^{m})^{2}, where d and m stand for data and model. Data for viral load, IFN and antibody titre are expressed in units of log_{10}, lung lesions/dead cells are expressed as a fraction between 0 and 1. Errors are assumed to be normally/lognormally distributed, with equal variance for a given dataset. Since for some of the fits, several different quantities are fitted simultaneously, it was necessary to give similar weights to the different datasets. We did so by rescaling the data for each different variable (log virus load, log antibody titre, fraction dead cells, log IFN) by its maximum value, to ensure that datasets for each variable have values between zero and one.
Data were fitted to the mathematical models using several different least square fit routines (lsqnonlin, nlinfit, fmincon) provided by Matlab R2007a (The Mathworks). More details on the solver routines can be found in the Matlab documentation, available at www.mathworks.com. The scripts used to produce the results are available from the authors upon request. We used more than one least square solver since we noticed that occasionally a certain fitting routine would return a local minimum; however, when this minimum was used as a starting value, a different solver routine was able to improve on it. We therefore switched between the different solvers until none of them was able to further improve the solution. We also varied the initial conditions for the parameters over a wide range to ensure that we obtained the globally best fit. In addition, we tried a few different global optimization algorithms, namely simulated annealing, genetic algorithms and direct search methods, all part of the Matlab Genetic Algorithm and Direct Search Toolbox. None of these algorithms led to different (better) results. In general, the optimization landscape was found to be rather rugged, which required occasional ‘manual intervention’ through switching solvers, tweaking initial guesses and restarting the fit several times. Because of this, we decided not to perform bootstrapping to obtain confidence intervals, since such a procedure does not allow manual interventions and we did not trust an ‘automated’ routine to properly converge to the right minimum every time, which would have led to false estimates for the confidence intervals.
While we assess the quality of model fits mainly based on biological relevance and not statistical criteria, we also computed the Akaike Information Criterion (AIC) as an additional way to assess the goodnessoffit for the different models. A lower AIC value means that the model describes the data better. We employ the corrected AIC, which accounts for low numbers of data points, given by2.2where m is the number of parameters plus one, n is the number of data points and SSR is the sum of square residuals obtained from the fitting routine, with the assumption that the variance is equal to SSR/n (Burnham and Anderson 2002).
3. Results
3.1. Fitting the first dataset
We begin by fitting the data from Kris et al. (1988). The data consist of viral load for influenza infection in wildtype BALBc mice with a functioning immune system, as well as nude mice, which do not have a thymus and therefore no Tcells, which also leads to a severely impaired Bcell response (symbols in figure 1). As can be seen, in the absence of well functioning B and Tcell responses, the virus persists at a high level; eventually the mice succumb to the infection. In contrast, if a functioning adaptive IR is present, the infection is cleared relatively rapidly. This shows that the adaptive IR is important and will have to account for the different virus dynamics between the two datasets.
Both Bcells/antibodies and Tcells have been suggested as the main drivers of virus clearance (Palladino et al. 1995; Mozdzanowska et al. 2000; Hasegawa et al. 2002; Lee et al. 2005). Here, we focus on an antibodymediated adaptive IR. Since no antibody data are available, the antibody response is essentially on an arbitrary scale and, without loss of generality, we can reduce the number of parameters by setting the antibodymediated virion clearance rate to k = 1. We fit the data for wildtype and nude mice simultaneously, with an antibody response present in the wildtype mice and absent for the nude mice. We further start with the assumption that regrowth of epithelium can be neglected (λ = 0). We also initially ignore a possible role of the innate IR (w = 0).
Figure 1a shows that a decent fit can be obtained. However, crucial features of the virus dynamics, such as the initial growth phase, are not well reproduced. This is captured by the unrealistically low basic reproductive number of R_{0} = 1.05. Along with this, the bestfit inoculum dose is V_{0} = 3.2 × 10^{6}, which is higher than the virus load data for day 1 in the BALBc mice. Unfortunately, Kris et al. do not report the virus inoculum used for their experiments. An earlier study from the same laboratory used a similar protocol (Kris et al. 1985), and in that study the inoculum is 10^{6} EID. However, their laboratory also showed in yet another study that the ‘take’ of virus (measured as the fraction of inoculum that could be recovered 10 min postinfection) can be very variable and much lower than the inoculum dose (Yetter et al. 1980). Given this uncertainty, we decided to set upper bounds on V_{0} and investigate how the results change. For V_{0} = 10^{5}, the fit is still reasonable; however, once V_{0} is bound by V_{0} = 10^{4}, which is just below the viral load for day one in the wildtype mice, the fit becomes very poor (figure 1b). (Allowing the fixed parameters c, d and g to vary leads to a somewhat lower SSR but cannot remedy the overall poor fit.) With the information we have available, we cannot categorically rule out that the initial effective inoculum is indeed at levels around 10^{6}. However, the increase in influenza virus titre is usually rapid, suggesting an R_{0} that is well above 1 (Baccam et al. 2006). A high R_{0} leads to significant cell destruction, which in turn leads to the virus running out of target cells and therefore a decline in viral load. This is the reason why the model cannot properly reproduce both a strong initial virus growth phase and a prolonged almost constant level of virus load, as seen in the nude mice.
One way a steady state could be maintained is if destroyed epithelial cells are replaced by new susceptible cells (similar to the steadystate dynamics observed in HIV, where the CD4^{+} Tcells, which are the virus's main target, are constantly replenished). To test this idea, we next fitted the model with λ ≠ 0. It is now possible to obtain a decent fit even with the bound V_{0} ≤ 10^{4} (figure 2). A higher value of R_{0} is obtained for this fit (table 2), which is more in line with the rapid increase in influenza viral titre usually observed. This model predicts a relatively high level of destruction of lung epithelial cells in both wildtype and nude mice. To maintain an almost constant level of virus in the immunocompromised mice, the effective reproductive number needs to be R ≈ 1. Since R_{0} is much larger, the only way to achieve R ≈ 1 is by depletion of most epithelial cells and then, through regrowth of epithelium at rate λ, a steady state of virus, infected cells, and a low level of uninfected target cells can be maintained. It is not clear if such a high level of lung epithelium destruction is reasonable. Presumably, this would lead to rapid death of the mice.
Instead of assuming that the reduction of R to a value of about 1 is mainly because of the depletion of target cells, it could also be possible that it is to a large extent owing to innate IR mechanisms. Despite the lack of a functioning adaptive IR, the nude mice still have a largely wellfunctioning innate response. Components of the innate IR, such as various cytokines, macrophages and neutrophils, have all been shown to play a role in influenza virus control (Tumpey et al. 2005; Koerner et al. 2007; Kim et al. 2008; Tate et al. 2008). We tested a model that includes an innate response, as described in §2. Since we again do not have data for the innate response, we can reduce the parameters by setting w = 1 without loss of generality. As figure 3 shows, we find a good fit, even in the absence of epithelial regrowth (λ = 0). Allowing for epithelial regrowth does not further improve the results (not shown). We assumed the innate IR to be the same for the wildtype and nude mice. However, because of the many feedback loops that exist between different components of the IR, it is possible that the nude mice also have an altered innate response. Allowing the parameters for the innate response to vary between the wildtype and nude mice better captures the differences in virus load between wildtype and nude mice during the first few days, but the added parameters are not justified on statistical grounds (not shown).
We can conclude from this section that both regrowth of epithelial cells and a combination of innate and adaptive IR can explain the data. The latter model has an additional parameter compared with the former model and is not justified on statistical grounds (table 2). However, we believe it to be biologically more plausible. First, a study in ferrets found that newly grown epithelial cells were initially resistant to infections (StuartHarris & Francis 1938), which could indicate that even though epithelial cell regrowth a few days postinfection is observed in mice (Ramphal et al. 1979), these might not be available for reinfection. A more recent modeling study for URT influenza infections in humans also concluded that the regrowth of epithelial cells could be neglected (Baccam et al. 2006). Second, in the presence of an innate response, the level of epithelial cell destruction is lower, which seems more plausible (see the next section). Lastly, many studies have shown that the innate response plays an important role during influenza infection (Tumpey et al. 2005; Koerner et al. 2007; Kim et al. 2008; Maines et al. 2008; Tate et al. 2008). Nevertheless, based solely on the data analysed in this section, we cannot draw a definite conclusion.
3.2. Fitting the second dataset
The previous dataset only provided virus load, albeit for wildtype and nude mice. This allowed us to draw some conclusions, but we could not completely resolve the importance of the innate IR. To make further progress, we turn to the study by IN (1977), which reports data not only on viral load, but also on percent lung lesions, IFN levels and antibody titres. In addition to infection of fully immunocompetent mice, IN performed studies in which they depleted IgM, IgG or IgA antibodies. They found that while IgG depletion did not change the infection dynamics and the mice survived (as did the wildtype mice), depletion of either IgM or IgA led to a change in virus dynamics and death of the mice. This again indicates the importance of the adaptive IR, and more specifically, the presence of certain types of antibodies. We focus here on the data for the IgMdepletion experiment, since these lead to almost complete abrogation of all three types of antibodies. Therefore, as in the previous dataset, we can assume that the difference between the immunocompetent and the antiIgMtreated mice is the presence or absence of the Bcell/antibody IR. (A possible Tcell response is discussed in appendix A.) The data for the wildtype and antiIgM experiments are reported in figures 1, 2, 7 and 8 in Iwasaki & Nozima (1977) (see symbols in figure 5).
First, we consider a model that allows for regrowth of epithelial cells and includes an adaptive IR, but we ignore the potential role of the innate response. The IN study reports the amount of virus inoculum the mice were exposed to (V_{0} = 4 × 10^{4} PFU), which we use as upper bound for the viral load in our fits. Figure 4 shows the best fit for the model, which is clearly rather poor.
Next, we tried to see if including an innate response can lead to a better fit. The IN study reports data on one important component of the innate IR, the IFN system (GarciaSastre 2001, 2002), on which we focus in the following. The IN study reports IFN levels for both the situation with an intact IR and the antiIgM scenario. While the increase and decrease in the IFN response in the immunocompetent mice can be well described with our innate IR model used in the previous section, this is not possible for the antiIgM scenario. Specifically, in this situation, after the initial virus increase and decline, the virus increases again around day 6 and reaches very high levels by day 12; this coincides with increases in lung pathology and is followed by death of the animals (figure 5a,b). Curiously, there is no concurrent secondary increase of IFN levels, at least not before day 8, which is the last day at which IN report IFN levels (figure 5d). We could accommodate the fact that IFN does not increase again, despite increase in virus load, by for instance ‘turning off’ growth (setting w = 0) at day 5 in our innate IR equation. However, this seems rather ad hoc and we are not aware of any biological mechanism that would justify such a model. We therefore decided that instead of modelling IFN, we fit it with a simple nonmechanistic equation (Regoes et al. 2004), namely exponential IFN increase until the peak at day 5, followed by exponential decay. This essentially reproduces an equation that would have an arbitrary turnoff for growth at day 5, but does not try to cast the IFN dynamics in the form of a (contrived) mechanistic model. This approach has the added benefit that we can include a time lag (see below) without having to deal with delay differential equations.
With the IFN levels determined (figure 5d), we next fit the available data for viral load, lung lesions and antibody titre to the model. We assume that the antibody response is present for the wildtype and absent for the antiIgM mice. Note that we fit the model to antibody data, in effect assuming that Bcell numbers and antibody levels follow an equivalent dynamics. Figure 5 shows the best fit for the model.
While the agreement for the antibody level and fraction of dead cells is decent, the fit to the viral load data is not. One problem could be our assumption that the reported data on lung lesions correspond directly to dead cells in our model. It is unclear what the exact correspondence between lung lesions and dead cells in our model should be. We therefore introduced a scaling factor between the lung lesion data and the dead cells in our model. However, allowing for such a scaling factor did not improve the results significantly (not shown) and the best fit was obtained for a value of this scaling factor of about 1. Another possibility is that the cytotoxic T lymphocytes (CTL)based IR is present and plays a role for both the immunocompetent and antibody depleted mice. However, extending the model to include such a CTLbased IR (with dynamics described by the same model as that for the antibody response) also did not improve the fit (not shown). Another possibly important feature could be that there is a delay between a given level of IFN and its action of reducing virus production. We therefore included a delay, such that F(t − τ) was acting on the reduction of virus production. This produced a better fit (figure 6), with a delay time of about τ = 1.3 days. To check the role of targetcell replenishment, we forced λ = 0 and found that it leads to an almost identical fit (not shown).
Interestingly, the model still predicts a decline in virus load at days 11 and 12 for the antiIgM mice. A somewhat better fit with a less steep virus decline on these days can be achieved for lower values of c, e.g. for c = 1 we find a fit with an SSR = 0.08. However, a halflife of virions of 24 h seems somewhat too long. Alternatively, increasing the lifespan of infected cells does also lead to a less steep virus decline on days 11 and 12 but the overall fit becomes worse.
We can conclude from this section that the data are best described by a model that includes an adaptive IR and an innate IR. Using IFN levels as the proxy for the innate response, we found that a time lag between levels of IFN and its action on the reduction in virus production improved the fit. Targetcell replenishment was not found to be important. Even the bestfitting model is not completely satisfactory; specifically the virus dynamics on days 11 and 12 for the antiIgM mice is not well captured.
4. Discussion
Two previous studies showed that viral load data of URT influenza infections from humans could be fitted with a simple resource limitation model, without the need to invoke an IR (Baccam et al. 2006; Handel et al. 2007). However, both studies also indicated that models that included either an innate response (Baccam et al. 2006) or an adaptive response (Handel et al. 2007) were consistent with the data. Unfortunately, not enough data are available to quantify the role of the IR in human infections. However, for influenza infections in mice, more data are available. Here, we studied influenza infections of the LRT of mice, which is a model for human influenza pneumonia—and might also partly apply to uncomplicated URT infections, though this needs to be investigated further. The availability of additional data for immunocompromised animals, as well as data for dead cells and IR components allowed us to make further progress in discriminating between different possible models for the infection dynamics. We find that for the experimental data investigated, both an innate and an adaptive IR are required to properly describe the infection dynamics. Our results additionally suggest that regrowth of epithelium is not an important driver of the infection dynamics.
When fitting the second dataset, we took IFN as a proxy for the whole innate response. Some previous studies in mice suggested that IFN did not significantly affect the infection dynamics (Graham et al. 1993; Price et al. 2000). One reason for this could be the fact that many laboratory mice strains do not have a functioning Mx1 gene. The Mx1 gene has been shown to be important for the induction of the IFNmediated antiviral state (Staeheli & Haller 1987; Staeheli et al. 1988; Grimm et al. 2007). Both mice strains used in the two experimental studies are Mx1 deficient. While IFN can also act through other pathways (GarciaSastre 2001, 2002), the relative importance of each of these pathways is not well understood. Additionally, it has been shown that components of the innate response other than IFN play an important role (Kim et al. 2008; Tate et al. 2008). Therefore, IFN might be a good proxy for the dynamics of the innate response, but is probably not the only important component that influences the infection dynamics. The role of IFN as a proxy for the innate response dynamics deserves further study. For the reasons mentioned above, we did not use a mechanistic model to fit the IFN response. It seems curious that the IFN levels continue to decline, despite a rebound of the virus. It would be interesting to see if this is a general feature. Unfortunately, we are only aware of data for IFN levels in immunocompetent hosts, where virus levels increase up to day 2 or 3, followed by virus decline and clearance. Not surprisingly, IFN follows that dynamics. It would be interesting to see additional data for the IFN dynamics for a situation such as the nude mice from the first dataset. Unfortunately, to our knowledge apart from the IN study, no other such studies/data exist in the literature.
Our results improved by including a lag between IFN levels and action of a given level of IFN. The importance of a lag was previously found for influenza infections in humans (Baccam et al. 2006). An earlier study of hepatitis C virus infections also reported a similar lag (Neumann et al. 1998). Since IFN acts in many different ways (Stark et al. 1998; GarciaSastre 2001, 2002), it is hard to suggest a specific mechanism for such a delay. But in general terms, IFN is upstream of a signalling cascade that eventually can lead to the activation of effectors that induce an antiviral state in cells (Stark et al. 1998; GarciaSastre 2001). It is therefore likely that a certain delay between a given IFN level and action on the virus/cells occurs.
We used data from two specific experimental studies to fit our models. It is well known that there can be large variability between different experiments, owing to experimental procedures that differ between laboratories, the use of different influenza strains, and the use of different mice strains. Indeed, the importance of certain types of antibodies, such as IgA, IgM and IgG, continues to be studied, with results that do not always agree with each other (Kris et al. 1985). We therefore do not believe that the ability of a model with only an antibodybased adaptive IR and no Tcellbased component means the latter plays no role. Indeed, as briefly described in appendix A, models that replace the antibody response with a CTL response can fit the data equally well. This implies that of course models with both an antibody and Tcell response will also be able to fit the data. However, the available data do not seem to justify building more such complicated models, since model discrimination would not be possible.
In summary, we have shown that both the innate and adaptive IRs are important to describe the dynamics of influenza A virus in LRT infections. Our findings also suggest that regrowth of epithelial cells is unimportant on the time scale of the infections considered here. However, the first dataset led to inconclusive results in this regard and this deserves further study.
Acknowledgements
AH was supported by NIH grant K25AI072193. IML and RA were supported by NIH grant U01GM070749. We thank Alan Perelson for discussions during the early stages of this work.
Appendix A
In the main text, we modelled an antibodybased adaptive IR and the IFN response in a way that is biologically reasonable but by no means the only possible model formulation. We considered a number of alternative models, which we briefly describe here.
A.1. Additional innate IR models
Since the innate response consists of both inflammatory cytokines as well as populations of cells such as neutrophils and macrophages (Maines et al. 2008), it can impact the virus dynamics in multiple ways. In the main text, we assumed that the action of the innate response is to reduce the production of infectious virions. Alternatively, inflammatory cytokines might render cells less susceptible to infection, while innate cellular components (e.g. neutrophils, macrophages) can increase death of infected cells or clearance of infectious virions. We therefore considered several alternative mechanisms for the action of the innate response.

F1: reduction of susceptibility to infection, b → b/(1 + κF).

F2: killing of infected cells, −κIF term added to the equation for infected cells.

F3: removal of infectious virions, −κVF term added to the equation for virions.
A.2. Additional adaptive IR models
We investigated several other models for the adaptive IR dynamics.

X1: virusindependent clonal expansion at a fixed rate, given by Ẋ = rX.

X2: increase in IR proportional to viral load (mass action), Ẋ = rXV.

X3: increase in IR according to mass action at low virus densities and virusindependent clonal expansion at high virus densities, Ẋ = rXV/(V + s).
As an alternative to a humoral, Bcell/antibody response, we investigated a CD8^{+}, CTL response, which has also been found to play a role during influenza infections (Mackenzie et al. 1989; Eichelberger et al. 1991; Bender et al. 1992). We modelled a CTL response by adding a massaction killing term, −kXI, to the equation for infected cells (and a term +kXI to the equation for dead cells), while removing the antibodybased massaction clearance term (−kXV) from the equation for the free virions. The CTL dynamics was assumed to follow the same equations as the humoral response.
In addition to altering the dynamics of the adaptive IR, we also investigated an alternative to the massaction killing term (−kXV), namely we considered killing that saturates at high virus load (Pilyugin & Antia 2000; Handel et al. 2009), modelled with a term (−kXV/(V + s)).
Finally, for both the adaptive and innate responses, we considered models in which the expansion/activation was proportional to the number of infected cells instead of the virus load (e.g. Ẋ = fI + rX instead of Ẋ = fV + rX for the IR dynamics discussed in the main text).
A.3. Additional results for the first dataset
We found results to be similar for almost all of the different innate and adaptive model variants described in the previous two sections. (In the following, when we say ‘similar’ we mean a few percent difference in the SSR.) For all combinations tested, models without innate IR and no epithelial regeneration could not fit the data well if we bounded V_{0} as described in the main text. Including epithelial regeneration led to fits very similar to those shown in the main text, i.e. almost complete epithelial destruction. Once an innate IR was included, results were similar to those shown in figure 3, independent of the presence or absence of epithelial regrowth. The different mechanisms for the innate response produced equally good fits, while for the different adaptive IR models the pure massaction model (X2) performed worse compared with the other model variants.
Replacing the antibodybased response with a CTLbased response (i.e. switching from killing free virions to killing infected cells) led to very similar fits. Again, the choice for the innate mechanism did not matter and the massaction model for the adaptive IR (X2) did not fare as well as the other three models.
Further, allowing for saturation in the adaptive IR killing did not led to significant improvements for any of the models. Equivalently, models in which we assumed that the innate or adaptive IR expansion depended on the number of infected cells led to results similar to the ones where we assumed dependence on viral load.
The IN (1977) study provided data for the antibody response, with a bestfit value for the expansion rate, r, which was lower compared with the result from the first dataset (table 2). We decided to check how the fits for the first dataset would change if we set an upper bound at r = 0.5, in line with the rvalue obtained for the IN data. This leads to a qualitatively equally good fit, but with an SSR that is about twice as large as for the model shown in the main text.
A.4. Additional results for the second dataset
For the second dataset, we again considered various combinations for the dynamics and mode of action of the innate and adaptive responses as described in the previous sections.
For all considered model combinations, models that did not include both the innate response and a time delay for the innate response produced poor fits similar to those shown in figures 4 and 5. Including the time delay leads to fits for the different model combinations that were of similar quality as the one shown in figure 6, no matter how we chose the adaptive response (all different expansion dynamics models, antibody or CTL response, with saturation in killing or without). We also found that innate IRs F1 and F2 led to worse fits, whereas F3 gave results similar to the one shown in the main text.
We also considered an alternative to the exponential growth and decline function we used to fit the IFN response for dataset 2. It turns out that a simple Gaussian function can provide a good fit to the logtransformed IFN data. The difference between a Gaussian function and the triangular function used in the main text is that IFN starts at a base level, then increases following the data, and declines back to a base level. While we do not know if such a shape is biologically more reasonable, we tried this alternative for the IFN response. We did not find much different results. Most of the time, the fits were slightly worse compared with the IFN function used in the main text.
A.5. Additional discussion
Every model needs to strike a balance between simplicity and realism. The right balance should be guided by the questions one wants to address. Equally important, mathematical models that are used to fit data need to be tailored to the available data. Obviously, the models we use to describe the IR dynamics are very simple. The many different combinations for specific aspects of the dynamics that we considered lead to a rather large number of models, most of which we only discussed very briefly and in a summary fashion in this appendix. Many more biologically realistic and complex models could be considered. However, given the fact that many of the arguably simplest models we used here were able to fit the data well suggests that more complicated models would be able to fit the data as well. Model discrimination (i.e. rejection of models that do not fit) would not be possible, which would prevent any further insights. While it will certainly be important to advance our models from the very simple ones we used here towards more realistic and detailed ones, this endeavour can only be successful if it is closely tied to data—data that are powerful enough to allow discrimination between alternative models. Currently, such data for influenza do not seem to be available, suggesting that future modelling efforts need to be combined with the production of new data. Additional experimental studies along the lines of those analysed here, namely the measurement of kinetic data for both the virus and IR for both wildtype infections and infections for which specific aspects of the IR have been knocked out, should be able to produce the data necessary for future model discrimination. A collaborative effort between experimentalists and modellers will hopefully move us further towards a quantitative understanding of the dynamics of this important pathogen.
 Received February 19, 2009.
 Accepted April 14, 2009.
 © 2009 The Royal Society