Modeling Influenza Virus Infection: A Roadmap for Influenza Research

Influenza A virus (IAV) infection represents a global threat causing seasonal outbreaks and pandemics. Additionally, secondary bacterial infections, caused mainly by Streptococcus pneumoniae, are one of the main complications and responsible for the enhanced morbidity and mortality associated with IAV infections. In spite of the significant advances in our knowledge of IAV infections, holistic comprehension of the interplay between IAV and the host immune response (IR) remains largely fragmented. During the last decade, mathematical modeling has been instrumental to explain and quantify IAV dynamics. In this paper, we review not only the state of the art of mathematical models of IAV infection but also the methodologies exploited for parameter estimation. We focus on the adaptive IR control of IAV infection and the possible mechanisms that could promote a secondary bacterial coinfection. To exemplify IAV dynamics and identifiability issues, a mathematical model to explain the interactions between adaptive IR and IAV infection is considered. Furthermore, in this paper we propose a roadmap for future influenza research. The development of a mathematical modeling framework with a secondary bacterial coinfection, immunosenescence, host genetic factors and responsiveness to vaccination will be pivotal to advance IAV infection understanding and treatment optimization.

Advanced knowledge of the immune system during IAV infection is a key step towards the establishment of improved prevention, control and treatment strategies for IAV infection. To this end, in recent years mathematical modeling has risen as a promising tool to quantify IAV infection as well as to identify potential therapeutic targets to control IAV. Mathematical models have been developed with different grades of detail providing quantitative insights into viral dynamics within a host [9][10][11]. Although the picture of IAV kinetics resulting from mathematical modeling and experimental studies has improved, knowledge of the mechanisms of viral infection control, immune modulation, and bacterial coinfection is still largely fragmented. Three excellent reviews in influenza mathematical modeling can be found in [9][10][11]. In this paper, we update the state of the art of the mathematical models for IAV infection, focusing mainly on parameter estimation issues. To this end, we go beyond the review presenting a mathematical model to exemplify different problems in mathematical biology. In the last section, we highlight relevant open questions, creating a roadmap for mathematical modeling towards a global understanding of influenza infection research.

IAV Pathogenesis
IAV infects mainly epithelial cells of the upper respiratory tract. However, IAV can infect cells of the lower respiratory tract in more severe cases [12][13][14]. Typical symptoms of IAV include fever, muscle aches, fatigue, runny nose and sore throat [15]. IAV infections are self-limited due to the development of a protective immune response (IR) [16]. An IAV infection starts when virions [17] enter the upper portion of the respiratory tract. Figure 1a shows major IAV infection events within a host. In the first phase, the binding and fusion capacity of the IAV hemagglutinin protein [18] facilitates virion adhesion to the receptors of the epithelial cell surface (binding). Approximately 20 min after infection, the virion is endocytosed inside the epithelial cell (endocytosis) [19]. Once inside the cell, the viral RNA is released and used as template to produce new virions, which are subsequently released into the extracellular environment (budding). The period between successful infection of the cell and the virus release is denominated the "eclipse phase" [20,21], which ranges between 5 and 12 hours post infection (hpi) [9,22]. After that, the new virions can start to infect other cells. The loss of virions can be attributed either to loss of virus infectivity or clearance by the immune system.  IAV replication peaks around 2-3 days post infection (dpi). Viral clearance is observed in most of the cases between 5 and 10 dpi [23]. In IAV infections, pathogen specific antibodies (Abs) and cytotoxic CD8 + T cells (CTLs) are detectable after 5 dpi and peak at 7 and 10 dpi respectively [24]. However, Immunoglobulin (Ig)G Abs peak later around 25 dpi [25]. Natural killer (NK) cells and interferon (IFN)-I are the main components of the innate immune system attributed to play a role in controlling IAV infections [26]. These components provide early immune containment of the IAV infection. Figure 1b depicts the dynamics of IAV infection and the principal elements of the immune system. Mathematical models capable of quantifying the main immune factors involved in the responses to IAV infections have the potential to provide critical information to advance our understanding of the immunity to this virus. Importantly, this information can aid in novel vaccine design and/or lead to strategies that limit complications such as bacterial coinfections. Additionally, the information derived from these models can suggest strategies to boost the IRs to IAV vaccines in populations that are prone to complications (e.g., the elderly).

In Vivo Systems
The first mathematical model to describe IAV dynamics was developed in 1976 by Larson et al. [27]. The model was fitted to viral titer data of mice infected with IAV (H3N2). After thirty years without modeling efforts, a work that described the IAV infection dynamics was presented by Baccam et al. [26], which adopted the well-known target cell model [28][29][30].
Several mathematical works have tried to model the eclipse phase in vivo [26] and in vitro [31,32]. These models have aimed at representing the time frame of the infection more adequately. This has resulted in an additional state, in which newly infected cells rest in a latent phase before becoming productively infected cells (I). Thus, the model in Equations (1)-(3) with the eclipse phase can be represented as follows:U where E represents the cells in the eclipse phase, which can become productively infected at rate k. In Holder and Beauchemin [32], the authors considered different time distributions for modeling the eclipse phase and viral release by infected cells. The different mathematical model formulations were fitted with in vitro data. The results showed that the time distribution forms of the eclipse phase and viral release directly affect the parameter estimation. Baccam et al. [26] fitted the model Equations (4)-(7) with data from human volunteers infected with IAV A/HongKong/123/77 (H1N1). The estimated biological parameters e.g., viral clearance and cell half-life provided quantitative means of the viral infection dynamics. However, as argued by the authors, the parameter values should be used with caution due to identifiability problems (parameters can not be estimated in a unique way from the respective experimental data). In a similar direction Handel et al. [33] applied the target cell model to human influenza data to assess the emergence of resistance to neuraminidase inhibitors. However, this study also raised identifiability issues and confidence intervals were not provided. Recently, Dobrovolny et al. [34] adopted a double target cell model [35] with the eclipse phase to predict the efficacy of neuraminidase inhibitors on uncomplicated and more severe infections inside a host. For this purpose, the authors considered two different types of target cells (default and secondary). The first, is the fraction of cells available to IAV infection; meanwhile the second type of cell is the one accessible for severe IAV infections. The model was capable of mimicking the dynamics of uncomplicated IAV infection, including the IR in the secondary cell population dynamics. In another study, Petrie et al. [36] investigated parameter uncertainty using the target cell model and explored the possible reduction of parameter uncertainty fitting by measuring both infectious (via tissue infection culture dose 50 (TCID 50 )) and total viral load (via reverse transcription polimerase chain reaction (RT-PCR)). In addition, Petrie et al. [36] revealed that the variation in TCID 50 assay sensitivity and calibration may affect the parameter estimation.

Mathematical Models Including the Immune Response
Although the target cell model can predict the dynamics of IAV infection without considering the host IR, control of the IR on the viral infection is fundamental to viral clearance [37]. To better understand the factors shaping the IAV infection course, a comprehensive model should incorporate the IR dynamics and its interactions with IAV. Numerous mathematical models have been introduced to investigate and predict the dynamics of the IR during IAV infections [11]. The first model was proposed by Bocharov and Romanyukha [38], which considered different mediators of the immune system. This study uncovered that CTLs and Abs are the main players controlling IAV infections. Hancioglu et al. [39] adopted the same modeling approach as Bocharov and Romanyukha [38], concluding that the initial viral load influenced disease severity. The target cell model was extended by Baccam et al. [26] to include the role of IFN-I. The IFN-I dynamics were modeled with the equation dF dt = sE (t − τ) − αF, where F is the level of the IFN-I, s is the secretion rate of IFN-I by infected cells, α is the IFN-I clearance rate and τ the lag period necessary to secret IFN-I. In this model, the influence of IFN-I is assumed to inhibit the viral replication rate through the fractional form p 1+ F , where p is the viral replication rate and is the effectiveness of IFN-I. Although, problems of parameter identifiability arose, the mathematical model including IFN-I dynamics was able to describe the double peaks present in some viral titer data. The role of IFN-I in preventing the infection of new cells was also investigated in [39][40][41][42][43]. In these works, the target cell model was extended taking into account another compartment where the cells remain refractory to IAV infection. In Hancioglu et al. [39], the dynamics of IFN-I were considered to promote the resistance of epithelial cells to IAV infection. An increased production rate of IFN-I and induction to resistance rate were fundamental to control disease duration and damage. Saenz et al. [41], modeled viral and IFN-I dynamics using data from equine IAV infection, revealing the extensive role of innate immunity in controlling the rapid peak of IAV infection. Using the same equine IAV infection data but with less parameters than [41], Pawelek et al. [40] included the IFN-I response in the target cell model. It was shown that the rapid and viral decline after the peak can be explained by the killing of infected cells mediated by IFN-I activated cells, such as NK cells. Hernandez-Vargas et al. [42] revealed that proinflammatory cytokines levels (IFN-α, IFN-γ, tumor necrosis factor (TNF)-α) could be responsible for slowing down viral growth in aged mice, limiting the activation of CTLs and causing the impaired IR to IAV infection reported in elderly.
Recently, the mechanisms of innate IR in IAV reinfection data were investigated in [43,44]. Cao et al. [43] compared different in-host reinfection models including different IFN-I control mechanisms of viral infection on in vivo data. The authors revealed that a model that included a cell in an IFN-induced state-of-resistance cannot explain the observed viral hierarchy. Moreover, the authors postulated that the dynamics of secondary infection strongly depend upon the inter-exposure time.
Another relevant component of the innate IR is driven by NK cells [45,46]. Canini and Carrat [45] introduced a mathematical model with NK activation induced by IFN-I. Activated NK cells (N) follow the equation dN dt = F − ρN, where ρ is the NK cell death rate and F the level of IFN-I. Fitting of the viral kinetic and symptom data showed a correlation between within-host parameters and illness course. In addition, the adaptive IR dynamics directed by CTLs and Abs have been investigated in different works, all of which have concluded that CTLs play a relevant role in the clearance of the viral infection [11,25,34,39,[47][48][49][50][51][52]. Lee et al. [52] developed a mathematical model quantifying the relations between viral replication and adaptive immunity. This mathematical model predicted that (i) CD4 + T cells have a prominent role in antibody persistence; (ii) antiviral drugs need to be administered within 2 dpi; and (iii) CTLs in the lung are as effective as neutralizing Abs for virus clearance, when present at the time of challenge. A relevant study dissecting the early and late IAV infection phases was conducted by Miao et al. [25]. In that study, the data of the IRs were fitted with a splines method [53], which was then used as inputs of the viral dynamics. The authors estimated important kinetic parameters (e.g., the half-life of infected cells, the infection rate of target cells) revealing the crucial role of CTLs in clearing the infection. In the same year, Handel et al. [47] identified two models with and without IR. Using a simplified Abs response, the authors argued that the data considered were not sufficient to discriminate the effect of the IRs. Tridane and Kuang [54] investigated viral clearance, modeling the interactions between IAV-infected epithelial cells and CTLs. The authors assumed two models for CTLs response (logistic and threshold growth). The systematic analysis showed that both model, fail to represent the mechanisms of viral clearance. An elegant work by Dobrovolny et al. [11] assessed the ability of mathematical models including different components of IR to predict the viral titer course. The authors compared these models with "knockout" experimental data, finding that no single model was able to explain the different viral kinetics. Recently a within-host mathematical model was developed by Price et al. [51], merging the innate and adaptive IR dynamics with inflammation. Authors showed a positive effect of controlled inflammation by IR.

In Vitro Systems
The complexity of the host immune system poses a significant challenge for the identification of critical factors involved in the interaction with IAV. This complexity has limited the exploration of some of their basic biological functions. An in vitro framework can simplify the system providing a small number of components to study. For example, key biological aspects can be inferred by analyzing viral infections of cell cultures (in vitro). In fact, a large array of mathematical models, focused on understanding the interactions between the virus and target cells, have been developed using in vitro generated data. For instance, Moehler et al. [55] applied the target cell model to include cell death as well as replenishment of target cells. This model described the viral kinetics of the equine IAV (A/Equi 2 (H3N8), Newmarket 1/93) in Madin-Darby canine kidney (MDCK) cells (microcarrier culture) and concluded that the number of available target cells is fundamental for viral growth.
Beauchemin et al. [31] explored different variants of the target cell model. In the experiments, the dynamics of MDCK cell infection by IAV (A/Albany/1/98 (H3N2)) were studied in the presence of amantadine and showed that the efficacy of this drug to block the infection of target cells was 56%-74%. In Schulze-Horsel et al. [56] infection of MDCK cells with different IAV strains was studied. The model described the dynamics of virus particle release (infectious virions and hemagglutinin content) and the parameters were estimated for different IAV strains. The results suggested that strains with slow initial infection dynamics but late induction of apoptosis resulted in higher viral yields. In Holder et al. [57], the authors compared fitness of the wild-type (WT) IAV (A/Brisbane/59/2007 (H1N1)) and a strain that acquired the oseltamivir-resistance mutation H275Y in its neuraminidase (NA) gene. Two different experimental settings were used: plaque assay and viral yields. The plaque assay indicated a higher fitness for the WT strain. In contrast, the viral yield assay suggested a higher fitness for the H275Y strain. Interestingly, in silico versions of these assays showed that the plaque assay was more sensitive for studying the duration of the eclipse phase, while the viral yield assay was sensitive for virus production. The different sensitivities of two assays could explain the conflicting results in terms of fitness. Later, Pinilla et al. [21] revealed that the principal effects of the H275Y substitution on the pandemic H1N1 (H1N1pdm09) strain were to lengthen in the mean eclipse phase of infected cells (from 6.6 to 9.1 h) and decrease (by 7-fold) the viral burst size. More recently, Paradis et al. [58] using the same model, compared the viral fitness of H275Y and I223V NA mutants with the pandemic H1N1 (H1N1pdm09) WT strain. The eclipse phase duration was longer for both mutants than the eclipse phase length of WT strain (by 2.5 and 3.6 h respectively). These works show the enormous benefits of mathematical modeling and how in vitro experimental data can support the development of a mathematical model and vice versa. Table 1 presents the most recent list of works concerning mathematical models in IAV infections. Quantitative results can be largely affected by the host and the experimental setting. Animal models can allow more controlled experiments [18]. IAV infection kinetics have been studied in different hosts, including avian species (wild bird and domestic poultry), mammals (swine, pony, ferret, mouse, guinea pig, non-human primates) [18]. The mouse model is widely used mainly because of its low cost, availability of different strains and "knockouts" of genes that lead to deficiencies of specific immune populations or products (i.e., T cells, cytokines, etc.) [59].  [54] However, it is worth mentioning that IAV transmission in mice is less efficient than in humans [18], and thus might not reflect well the infection in humans. Correspondingly, many mathematical models have derived the model parameters from murine experiments [10,25,27,42,47,52]. In vitro experimental data in cell lines (e.g., MDCK cells) also contributed valuable informations [55,57]. A major problem when comparing mathematical models is the difference between experimental settings in which the experiments are performed. For example, viral titers, as defined in many works, are provided in different units (e.g., RNA copies, plaque-forming units (PFU/mL), focus forming units (FFU/mL), TCID 50 or egg infection dose 50 (EID 50 /mL)). Moreover, for the use of different units, knowing what they actually measure is critical for constructing reliable models and obtaining accurate parameter estimates (e.g., RNA copies include both infectious and non-infectious virions whereas TCID 50 only indicates those infectious ones) [36]. However even in the case where experimental data present the same viral titer units, experiments conducted in different labs cannot be compared because the viral units are not standardized [58]. Moreover, Paradis et al. [58] suggested that even similar experiments performed in the same laboratory but at different time points cannot be compared. In addition to the noisy nature of biological measurement, the data used in mathematical models were not usually collected with a modeling purpose in mind. As a consequence, the data collected may not be optimal for identifying the model parameters. In particular, there are usually many measurements in the viral growth phase while measurements in the clearance phase are missing. These issues could be solved, if a close and iterative collaboration between biologists and modelers is established before the experimental design.

Parameter Estimation: A Continuous Challenge
Identification of infection kinetic parameters is one of the fundamental tasks to advance the knowledge of IAV infection mechanisms. These may help and allow experimentalists to test different scenarios. Nevertheless, estimating properly the parameters is a big challenge when the data are limited. In [9], the authors reviewed the parameters of IAV kinetics from several studies and showed a large variation of IAV infection kinetics parameters, which is the result of different experimental settings (e.g., in vivo, in vitro), different hosts and IAV strains.
Thus, parameter estimation procedures induce many questions on the robustness of the results. This highlights the problem in developing detailed mathematical models with a large number of parameters to be identified on a limited data set. Concurrently, another bottleneck is represented by the methodological tools used to perform parameter estimation. We overview briefly the mathematical modeling approach for model identification in Figure 3.
To generate a quantitative understanding and to explore particular hypotheses, the formulation of different mathematical models needs to be driven by the experimental data (step 1). Note that multiple models can provide the same fit with observed experimental data. Thus, it becomes necessary to choose between different models. The standard approach to model selection is first estimate all model parameters from the data, then select the model with the best-fit error and some penalty on model complexity (Akaike Information Criterion, Bayesian Information Criterion [65]). After the model selection process, a key obstacle in order to obtain quantitative knowledge is the identifiability (step 2) of model parameters. A mathematical model is said to be identifiable when the parameter set can be uniquely determined. This can be achieved from the mathematical model structure (structural) and from the experimental data (practical). The structural and practical identifiability are necessary in mathematical models to reach significant predictions [66][67][68][69][70][71]. A very reliable method to test both the structural and practical identifiability is the profile likelihood method proposed by Raue et al. [72]. The idea behind this approach is to explore the parameter values, requiring for each parameter the optimization procedure of the cost function with respect to all other parameters. In particular, for each parameter, a range of values centered at the optimized value is generated in an adaptive manner. Re-optimization of the cost function with respect to the other parameters is performed for each value of the parameters. The aim of this approach is to detect directions where the likelihood flattens out [72].
Once the identifiability analysis is completed, parameter uncertainty analysis (step 3) is necessary to assess the large variability usually encountered in the biological data. The most frequent approach for parameter estimation is the bootstrap method. Bootstrapping is a statistic method for assigning measures of accuracy to estimates [73]. The nonparametric bootstrap considers data to be independent and identically distributed, whereas the parametric bootstrap requires imposing on the data a distribution assumption which is usually unknown.
Bootstrap methods are frequently used as a conventional tool to take into account the uncertainty of the estimated parameter by calculating the confidence interval from bootstrap samples [10,25,42,57,74]. The bootstrap methods can be affected by the large variation of a few measurements or by the imposed distribution assumption that is usually unknown. As a result, the parameters confidence intervals can span a broad range [9]. Improvements of the bootstrap method in mathematical modeling have been proposed in [75][76][77]. Alternatively, the Bayesian approach could deal more efficiently with the parameters uncertainty, as well as the model prediction [78]. At the moment, there were only a few applications of the Bayesian methodology in mathematical modeling literature [79].
The outcome values of the bootstrap can serve to construct scatter plots. These are graphical sensitivity analysis methods, simple but useful tools to test the robustness of the results and highlight parameter dependence. Finally, once parameter distributions are inferred, one can test the predictive power of the mathematical model (step 4), generate new knowledge or hypotheses of the biological process of interest and guide the design of new experiments.
In summary, as presented in [80], future works on mathematical modeling of IAV infection should present the methodologies and parameter analysis to support the veracity of the estimations. In the next section, we exemplify how to conduct such mathematical modeling approach based on a mathematical model with the basic interactions between IR in IAV infections.

Case Study: Identification of a Mathematical Model of IAV Infection Including the Immune Response
Several studies investigated quantitatively the relevance of CTLs to clear IAV infections [25,[39][40][41][42]52]. Due to the identifiability limitations to estimate the parameters in the target cell model only with viral load data, we propose a minimalistic model able to fit IAV and IR dynamics ( Figure 4). In this example, we follow the mathematical modeling approach described in the previous section and schematized in Figure 3. Step 1: Mathematical Modeling

E V
We consider a mathematical model assuming that IAV (V) replication induces the proliferation of CTLs (E). We assume that the clearance of the infected cells by CTLs can be represented by the clearance of the virus. The dynamics of V and E can be described by the following equations: We define the CTLs replenishment rate s E = c e E 0 , where E 0 is the initial number of CTLs and c e is the half life of CTLs. The steady state should be satisfied to guarantee the CTLs homeostatic value E = s E /c e in the absence of viral infection (V = 0). The CTLs proliferate at a rate r. We assume the activation of CTLs proliferation by V follows a Michaelis Menten growth with half saturation constant k e . The V growth is modeled with a logistic function with maximum carrying capacity K v and growth rate p. V is cleared at a rate c v E.
Step 2: Identifiability Analysis Here, we consider the murine data from [81] to estimate the model parameters.
In this cross-sectional study, young and aged BALB/c mice were infected with 50 µL (50-100 PFU) of the IAV (H1N1) (PR8) . Lung viral titers were measured at different time points by plaque assays. Various components of the IR were also monitored. In this section, the viral and CTLs time course from the young mice are considered. To reduce the number of parameters to be identified, we fix the value of c e = 2 × 10 −2 (half life is approximately T = 34 days) as reported by [82] and K v equal to the maximum value of the viral titer data in [81]. Therefore, the profile likelihood analysis can be computed on the unknown model parameters, this is shown in Figure 5. Note, that all the parameters present a profile likelihood with a minimum, implying that parameters are identifiable [72]. Therefore, the results suggest the possibility to infer model parameters from the available experimental data set.
Step 3: Parameter Uncertainty After checking the identifiability properties of parameters, we proceed with the parameter uncertainty analysis. To this end, we consider the nonparametric bootstrap to obtain parameter estimates and 95% confidence intervals [76].   A thousand different sampling repetitions from the data in [81] are performed. In each repetition, we draw a random sample from each time point and estimate the parameter. Model fitting and parameter distributions from nonparametric bootstrap are shown in Figures 6 and 7, respectively. Parameter estimates and 95 % confidence intervals are shown in Table 2, respectively. The constraints for model parameter optimization, lower and upper bound, were broadly chosen to consider a large biological parameter space (see Table 2). The quality of fitting, as shown in Figure 6, indicates that the model is able to describe IAV infection and IR in murine data. The computed confidence intervals in Table 2 indicate a large interval only for the parameter k e , while they result in a low extent for the remaining parameters. The large confidence interval of the parameter k e may be attributed to the flat profile likelihood observed in Figure 5c. Quantitative comparisons of parameters in Table 2 with the target cell model parameters are not possible due to the differences in model structures and units. Nevertheless, we can compare the viral production rate p in Table 2 with the viral growth rate presented in Smith et al. [83]. Authors proposed an approximation for the viral growth which is equivalent to the Figure 8 in the growth phase. Our p estimates [3.43; 6.08] d −1 are in the same range of the estimate (6.59 d −1 ) in [83].  [26,40,42]. In order to explore possible parameter dependencies, we build scatter plots from the computed parameters values using the nonparametric bootstrap as it is shown in Figure 8.   The scatter plots in Figure 8 show inter-dependencies between p-c v , p-r and c v -r parameters. This implies that the estimation of one parameter influences the estimate of the others. For example, in the case of the parameters p-c v , increasing the value of the p increases the value of c v , whereas increasing the value of p decreases the value of r. These parameter correlations could lead to possible parameter estimation problems. Moreover, we evaluate the dependence of model parameter estimates p, c v , r, k e to viral carrying capacity K v (data not shown). The estimates of p and r showed no correlation with the parameter K v , while there is a positive correlation of c v and k e estimates with K v .
Ad hoc experiments to estimate parameter values could be a possible solution to alleviate parameter inter-dependencies and, therefore, estimate the remaining parameters more adequately. For example, in the previous work proposed by Pinilla et al. [21], virus titers were incubated without target cells and followed up to determine the remaining infectious titers. In this way, the approximate values of the viral clearance rate could be determined and provide a more accurate estimates for the whole set of kinetics parameters. The same procedure to reduce parameters inter-dependencies was followed in [31,58].
We want to close this section remarking how the mathematical modeling approach presented in Figure 3 can be useful to to advance our understanding of IAV and IR dynamics. We would like to point out that due to the nature of experimental data and the simple structure of our model further analysis are needed to provide more quantitative insights in the role of CTLs and IR to control IAV infection, as previous works pointed out [39][40][41][42]52].

Discussion and Future Perspectives
During last years, mathematical models of IAV have contributed to gain quantitative insights into viral infection mechanisms, within host dynamics and population scales [85]. These models considered different levels of complexity trying to explain the fundamental aspects of IAV infection and integrating most of the available biological knowledge. Next, we discuss what we consider are the main challenges for influenza research. These include: (i) bacterial coinfections; (ii) aging of the immune system and the role in IAV infection; (iii) challenges for influenza vaccination; (iv) host and IAV genetic factors (Figure 9). In particular, as illustrated in Table 1, for coinfection and aging only two mathematical modeling works have been carried out to date.

Bacterial Coinfection
Next to morbidity and mortality caused by IAV infection itself, enhanced susceptibility to bacterial coinfections represents a major risk. This was first recognized for the 1918/19 influenza outbreak (Spanish flu), as retrospective studies revealed that most of the several million of deaths were caused by severe complications due to secondary infections with bacterial pathogens [86]. A defined predisposition to life-threatening infections by bacteria (secondary infections) in IAV infected individuals was documented during epidemics in the 1950's and 1960's [87], as well as during the 2009 (H1N1) swine flu outbreak [88]. Altogether, epidemiological data clearly attest a high fatal synergism between IAV and many bacterial pathogens, such as Streptococcus pneumoniae (S. pneumoniae) [89].
For many decades, the straight-forward correlation between respiratory epithelial destruction caused by the viral infection and the facilitated adherence, growth and spread of bacterial pathogens remained the most accepted explanation for the increased susceptibility to secondary bacterial infections [90,91]. However, owing to a number of studies published over the last years, it is now largely understood that IAV infection also modulates the host's immune system leaving it unable to mount an adequate anti-bacterial defense [92,93]. The proposed underlying mechanisms include the desensitization of alveolar macrophages [94], as well as the IFN-γ mediated suppression of their phagocytic function [95,96]. The recruitment and function of neutrophils, NK cells, and gamma-delta T cells have also been shown to be modulated during IAV infections [97][98][99]. Additionally, IFN-I has been identified as a major upstream mediator of some of these effects [99,100]. With regard to the respiratory epithelium, which hosts viral replication, IAV-mediated suppression of repair and regeneration processes [101], as well as the CD200-CD200R axis [102] have also been identified as players in the synergism between IAV and secondary bacterial pathogens.
Even though there have been rapid advances in our understanding of immune modulation through IAV, many questions still remain open. Enhanced susceptibility to bacterial infections can persist long after the virus has been cleared from the respiratory tract [94], but it remains largely elusive how and through which cell types these long term-effects are mediated. Furthermore, little attention has been paid to the question of whether certain viral strains are more efficient in predisposing the host to bacterial superinfection and, in turn, whether certain bacterial strains are favored, e.g., S. pneumoniae [92,103]. In this context, it is important to notice that primary IAV infections influence all aspects of bacterial pathogenesis, such as acquisition, transmission, colonization, replication, and adherence [103][104][105]. Even the live attenuated influenza vaccine (LAIV) increases the duration and density of S. pneumoniae nasopharyngeal carriage, which is a crucial prerequisite for invasive pneumococcal disease [106].
Despite the extensive research on the IAV role to facilitate bacterial colonization, the effects of S. pneumoniae on IAV infection dynamics are poorly understood. In McCullers [92], it is argued that the pneumococcal infection may strengthen IAV infection. The most plausible mechanisms proposed involve bacterial virulence factors interfering with the IR, such as, bacterial proteases facilitating endocytosis and synergy between bacterial and viral neuraminidase [92]. Ultimately, it will be essential to answer these open questions in order to complete our understanding of the synergism between IAV and secondary bacterial pathogens. This will in turn aid in the development of more effective prophylactic strategies and/or novel treatments.
Mathematical models can play an important role to unravel the key processes operating during coinfection. To the best of our knowledge, the only mathematical work to represent IAV-pneumococcus coinfection and the synergy between the virus and bacteria was developed by Smith et al. [64]. This work also addresses the role of the innate IR (alveolar macrophages) to clear bacteria. The down-regulation of bacterial clearance was modeled in [64] with a combination of two Michaelis-Menten functions including four parameters, assuming macrophages dynamics constant. We would like to highlight that even in the simpler form of bacterial clearance (c B B), the parameter c b is not identifiable. Then, taking into account a more complex form modeling the mechanisms of IR in the bacterial clearance as proposed in [64] will require the design of new experimental data. Thus, parameter values in [64] should be interpreted with caution. Further mathematical modeling research should include the interactions between macrophages, IFN-γ and CTLs, which can be important to uncover the underlying mechanisms between the innate and adaptive IR in bacterial clearance.

Aging of the Immune System and the Role in IAV Infections
The United Nations population division has predicted that by 2050 the number of elderly persons will increase from the current 600 million to nearly 2 billion worldwide. In developed countries, due to the rise in average life expectancy and reduced birth rates, 25% of the population is expected to be older than 65 years by 2050, making the elderly the fastest growing segment of the population [107]. Immunosenescence refers to the biological aging process associated with a progressive decline in systemic immunity and increased susceptibility to diseases such as cancer, autoimmune pathologies, chronic diseases (e.g., type 2 diabetes), poor responses to vaccination, and increased vulnerability to common infectious diseases [108]. The immune system becomes compromised, both quantitatively and qualitatively. A reduction in some cell populations and/or impaired responses has been described in murine models, and humans [81,109,110]. Both innate and adaptive IRs are affected. For instance, alterations in NK cells, NK T cells, dendritic cells (DC), monocytes and macrophages, and neutrophils have been described within the innate immune system [111][112][113]. On the adaptive side, immunosenescence has been associated with decreased clonal diversity of naive T cells, increased memory T cells, decreased CTLs function, decline of CD28 expression, increased number of regulatory T cells, and altered T cell signaling [81,[114][115][116][117][118]. New insights into the changes induced by immunosenescence are coming from systems biology approaches [119,120]. Development of IRs and immunosenescence are complex biological processes and therefore not the result of isolated events, single proteins, chemicals, enzymes, or individual cell types. These changes are the coordinated effects of age-induced alterations in gene regulation and expression, signaling pathways within the cells, the interaction of various cell types, tissues, and biological networks. In sum, it is the synergistic effect of the individual parts that leads to impaired immunity and consequently immunosenescence caused by aging.
Mathematical models hold the potential to aid in elucidating correlates of protection for infectious diseases in the elderly. In Hernandez-Vargas et al. [42] age-related changes and their effects on influenza virus infection dynamics were examined. This work considered different IR components in young and aged mice: CTLs, NK cells, IFN-I, IFN-γ, and TNF-α. The fits of mathematical models to the experimental data from young and aged mice suggested that the increased levels of IFN-I, IFN-γ, and TNF-α (the "inflammaging" state) could promote the slower viral growth observed in aged mice, which consequently could limit the stimulation of immune cells and contributes to the reported impaired responses in the elderly. This study illustrated the advantages gained by the mathematical modeling approach to dissect quantitatively the role of the different IR mediators in young versus aged. However, more efforts are needed to uncover the age-related changes. It is important to consider that preexisting immunity to IAV, due to past infections/vaccination can directly alter the outcome described in [42].

Challenges for Influenza Vaccination
Although vaccination remains the cornerstone of influenza prophylaxis, the protection provided by seasonal vaccines is only partial, especially in high-risk groups such as the elderly [121,122]. Furthermore, new influenza virus strains arising from this zoonotic pathogen increase concerns about the threat of a new pandemia [123]. In this regard, the global spread of the (2009) (H1N1) pandemic virus within 2 months clearly demonstrated that current vaccine platforms cannot meet an emerging pandemia [124,125]. The global vaccine manufacturing capacity was inadequate and unable to supply enough vaccine doses promptly [126,127]. Thus, in order to increase vaccination efficacy huge efforts need to be taken to optimize vaccines (e.g., by co-admixing with adjuvants) which are already in the market and/or to develop new vaccination strategies [128] (Figure 10).
A possible strategy could be the development of a mucosal vaccine. Today, most of the licensed seasonal vaccines are administered via parenteral routes failing to promote mucosal immunity [129]. This is an important issue considering that the respiratory tract constitutes the port of entry for IAVs. In this context, mucosal immunization offers the potential to stimulate protective IRs at both the systemic level and the site of viral entry, e.g., by the induction of secretory IgA and CTLs activity [130][131][132]. However, current split and subunit vaccines do not elicit adequate IRs when administered via mucosal routes as they show reduced immunogenicity. On the other hand, the cold-adapted attenuated intranasal vaccine is only licensed for a very narrow subpopulation group. To overcome this constrain, adjuvants can be included in such mucosal vaccination strategies. Indeed, no mucosal adjuvant has been licensed worldwide so far. Mostly, this is related to either a lack of activity or safety concerns. For example, the widely used adjuvant aluminum hydroxide (alum) does not provide appropriate mucosal adjuvant activity, whereas the use of derivatives of the heat-labile enterotoxin of Escherichia coli administered by an intranasal route has been linked to Bell's palsy [133,134].  Other novel adjuvants such as Polyethyleneimine (PEI) [135] or the novel nanoemulsion W(80)5EC were already tested in volunteers and elicited both systemic and mucosal immunity following a single intranasal vaccination [136]. It was shown that these compounds may significantly contribute towards the development of innovative mucosal vaccines against influenza, especially for high-risk groups. Beside the need of adjuvants, novel vaccination strategies would also be needed to stimulate efficient cellular IRs. For example, a vaccine based on the ectodomain of influenza matrix protein 2 (M2e) has been applied with adjuvants and carriers, or fused with flagellin [137]. The results of preclinical studies in animal models have shown that M2e-specific Abs limit IAV replication and reduce mortality rate [138]. Other multimeric vaccines based on proteins were tested in clinical trials benefited from Montanide adjuvant formulation [139,140]. The drawback remains in the limited amount of antigen available with protein vaccines, and their lower capacity for inducing cytotoxic defences compared with replicating viruses.
Efficacious immune defence against influenza virus requires both humoral and cytotoxic arms, as witnessed by the limited protection offered by current seasonal vaccines which induces mainly humoral responses [141]. Nevertheless, replicating vaccines were shown to be able to enhance both humoral and cytotoxic IRs. Thus, DNA vaccines as well as adenovirus and vaccinia virus vector-based vaccines have been developed. These vaccines are able to stimulate cross-clade neutralising antibodies and cell-mediated immunity (CMI) [142][143][144]. Another promising approach consists of self-replicating replicon-RNA(RepRNA) which is readily generated by synthetic means, and incorporated into biodegradable delivery vehicles. RepRNA replication increases the number of RNA templates and expressed antigens, inducing long-lasting humoral and CMI defences [145,146]. However, several new approaches are under development aiming in the generation of universal vaccines [147].
To the best of authors knowledge, modeling IAV vaccination as well as in-host vaccine effect has not been explored. Until now, mathematical models have focused mainly on the effect of the vaccines on population-level spread of IAV [148,149]. However, quantitative efforts are necessary to gain comprehension into the immunological signatures of vaccination, evaluating the efficacy of different vaccines, adjuvants, delivery systems and immunization routes [150][151][152]. This will only be possible with the help of quantitative models of affinity maturation of B cells in germinal centers [153] giving rise to Abs with high affinity to the respective pathogen [154]. However, affinity is not the only important Ab property in this context. The quantity and the affinity of the generated Abs are normally contradicting properties, and it might be interesting to test strategies to overcome this limitation of the overall strength of an Ab response [155] in the context of vaccinations. Furthermore, Ab cross-reactivity might be essential to determine the efficiency of a vaccine. Specificity and cross-reactivity of Abs emerging in response to monovalent and polyvalent vaccinations were investigated for the example of malaria [156]. A similar comparison of different influenza vaccines would enable us to address the age-dependent response to influenza vaccination in terms of both safety and efficacy [147], which represents a main unsolved challenge in vaccinology.

Host and IAV Genetic Factors
In the past, molecular mechanisms focusing on the IAV pathogen have been broadly investigated in influenza research. The processes of adhesion, entry, replication and assembly with respect to the individual mutations of the IAV were studied intensively. However, not only the viral but also host factors are crucial for the outcome of an influenza infection. Recently, crucial host factors could be identified (see Table 3).
Together with the help of animal models the complexity can be reduced, and the factor of interest can be studied under controlled circumstances. While transmission of IAV is mainly studied in ferrets, mice are advantageous since their immune system and genome are well characterized and research tools are available. Importantly, several aspects of the murine immune system are similar to the human one. Mouse models demonstrated a clear genetic effect on susceptibility to influenza A infection [157]. In addition to the availability of inbred mouse strains, targeted modifications of the genome resulting in transgenic mice are easily possible. A novel mouse population, the Collaborative Cross, with a genotypic variation comparable to the human population generated from eight founder strains might provide an excellent model system in future [158]. Preliminary studies using this population for IAV verified the significance of this approach [159,160]. Using such a systems approach this study contributes to a better understanding of underlying biological networks under host genetic control during an IAV infection [159]. Table 3. Host Genetic factors. Host factors identified having a crucial role to determine severity of IAV infection in different hosts.

CPT2
Related complication as influenza-associated encephalopathy [164] TMPRSS2 Resistance to IAV infection [165][166][167] Due to the advent of modern "omics" techniques, a huge amount of data is becoming available in individuals. The gene network complexity can be only achieved with the help of systems-oriented approaches [168]. The deterministic framework with differential equations is the most widely used mathematical approach to describe key biological mechanisms. This is feasible only with high numbers of particles or species where the random fluctuations around the species concentrations are negligible [169]. Nevertheless, stochastic fluctuations are prominent in the presence of low mRNA copy numbers, becoming the "fingerprint" of gene regulatory network processes such as transcription and translation [170,171]. Thus, a stochastic modeling framework can play an important role to elucidate and dissect the host genetic mechanisms responsible for the IAV infections outcome (e.g., identify new drug targets, test the efficacy of treatments in a genetically heterogeneous model, identify susceptibility and resistance [169,170,172]).