Review of Mathematical Modeling of the Inflammatory Response in Lung Infections and Injuries

Lung inflammation may occur due to viral and bacterial infections, structural damage, or inhalation of dangerous particles. These injuries may be quickly resolved by the immune system, treated effectively through various interventions, become chronic problems, or lead to death. Mathematical modeling has been used to understand immune system dynamics during a number of pulmonary infections and injuries, identify key mechanisms, and provide important insights into new treatments. In this review, we present long-accepted modeling techniques and novel strategies for simulating various lung injuries to highlight the usefulness of mathematical modeling in addressing these life-threatening conditions. Advances in computational power have allowed for a diverse collection of models ranging from those using only Boolean operatorsto complex hybrid multi-scale models, each specifically designed to address relevant biological questions. To illustrate the findings from these mathematical approaches, we present detailed examples, summarize results, and consider future directions from modeling influenza, pneumonia, COVID-19, tuberculosis, anthrax, and other non-infectious injuries.


INTRODUCTION
Inflammation in the lungs can occur for many reasons, from bacterial infections to stretch by mechanical ventilation. The immune response to these insults has been widely studied, but there is still much unknown, including the roles of immune cells such as macrophages, neutrophils and T cells and the downstream effects of subcellular signaling. Mathematical modeling is a useful tool in trying to understand complex processes like the immune system (Figure 1). Mathematical analysis and simulations can reproduce known biology and predict response to interventions. This review highlights long-accepted and well-documented modeling techniques as well as novel strategies that provide insight into the immune response in the lungs and inspiration for future models.

Biological Background
The immune response is triggered to protect the host from invading pathogens or tissue damage. Type I and type II alveolar cells interact with a foreign agent or respond to stress, releasing pro-inflammatory signaling proteins, including cytokines and chemokines, to elicit an immune response [1]. The innate immune system responds to these signals and provides the first line of defense, subsequently triggering the adaptive immune response which consists of immune mechanisms specific to the insult. Most immune responses involve varying intensities of both [2][3][4][5]. FIGURE 1 | An overview of lung inflammation, which is the focus of this review, and the mathematical techniques used to model the immune response. ODE, ordinary differential equation; PDE, partial differential equation.
There are two stages of inflammation: pro-inflammatory and anti-inflammatory, classified by the cytokines that up-or downregulate the immune cells and those produced by the immune cells. Typically the first stage is a pro-inflammatory response during which damaged cells and pathogens are destroyed [1,6]. Innate immune cells such as neutrophils, eosinophils, and macrophages are recruited to the damaged or threatened site [7]. These cells work to remove the insult, often by phagocytizing the pathogen and sending pro-inflammatory signals to recruit other immune cells [8].
In response to pro-inflammatory roles of cytokines and chemokines, neutrophils remove damaged tissue and pathogens through a number of mechanisms, including phagocytosis and release of neutrophil extracellular traps. They also produce molecules such as serine proteases and reactive oxygen species that aid in the containment of pathogens. Though neutrophils can be very effective at eliminating foreign agents, an excessive neutrophil response or unprogrammed cell death, called necrosis, can lead to chronic inflammation [9][10][11]. Eosinophils release their own extracellular traps that capture microorganisms, and granules within eosinophils release proteins that participate in the removal of the inflammatory stimuli [7]. Mast cells are immune cells involved in the response to various diseases and play a role in wound healing and asthma. Mast cells are recruited from the blood and release various cytokines and chemokines such as TNF-α and IL-4 [12].
In addition to these immune cells, macrophages and dendritic cells (DCs) eliminate the threat, phagocytize debris and other immune cells, and send signals, but are also antigen-presenting cells (APCs) [3]. APCs, in general, present antigens to adaptive immune cells triggering a specific defense against the insult. Tissue-resident macrophages are already present and form a primary line of defense against pathogens. Circulating monocytes mature into macrophages, which are rapidly recruited to the site of inflammation upon damage. Macrophages have high plasticity, with the two main subtypes classified as classically activated, or M1, and alternatively activated, or M2, although these groups have overlapping features [13]. M1 macrophages release pro-inflammatory cytokines and phagocytize microorganisms, damaged cells, and apoptotic cells [14][15][16]. M2 macrophages generally promote regulation of inflammation and wound repair and play a critical role in resolution of allergic inflammation [17,18]. As with macrophages there are many types of DCs, including CD103+ respiratory DCs which deal well with foreign particles and CD11b hi DCs which produce significant amounts of chemokines [1,19].
APCs acquire antigens and migrate to the lymph nodes (LN) where adaptive immune cells are activated. APCs display the pathogen antigens on their surfaces, promoting naive T cells to develop an adaptive response specific to that pathogen [20]. T cells fall under two main categories, CD8+ helper and CD4+ cytotoxic. CD8+ T cells kill infected cells directly whereas CD4+ T cells use more indirect methods such as supporting the production of antibodies through activation of B cells [21]. Natural killer (NK) cells, from the same family as B and T cells, are involved in injuries such as influenza and tuberculosis. NK cells release high levels of cytokines, especially IFN-γ , regulate the function of other immune cells, and aid in the transition to an adaptive immune response [22].
The second stage of the immune response is the antiinflammatory response. Pro-inflammatory signaling is downregulated and anti-inflammatory mediators and the antiinflammatory functions of pleiotropic mediators, such as IL-10, promote repair and a return to normal health [20,23]. This second phase is orchestrated by both the innate and adaptive arms of the immune response [1]. An anti-inflammatory environment and phagocytosis of apoptotic neutrophils causes M1 macrophages to shift to an M2 phenotype. This shift is characterized by the decrease in production of pro-inflammatory mediators and an increase in that of anti-inflammatory and repair mediators [14,24]. DCs also utilize a process called efferocytosis to clear apoptotic cells, contributing to initiation of the antiinflammatory response. In addition to their role in the proinflammatory phase, T cells and DCs can also produce cytokines with anti-inflammatory roles [19,25]. Another important cell type is the fibroblast, appearing at the beginning of the antiinflammatory phase and recruited by chemoattractants such as IL-1β and TNF-α. Fibroblasts produce extracellular matrix components to support healing of damaged tissue [26].
As evidenced by the various and overlapping roles of immune cells, both innate and adaptive, and the cytokines and chemokines they produce, the immune response is complex and in many ways specific to the insult type. Furthermore, the proper shift from a pro-inflammatory response to an antiinflammatory response is extremely important, and an imbalance in either phase can cause severe complications, especially for immunocompromised individuals [27][28][29]. The preceding paragraphs are meant to give an overview of the immune response and details related to specific diseases that are addressed in this review will be highlighted in the following sections.
Inflammation specifically in the lungs is caused either by a foreign insult, by tissue damage due to stretch or strain, and in rare cases, an autoimmune disorder. Infections of the lungs are some of the deadliest infections in the world [30,31], some of which have been studied for decades such as tuberculosis, influenza, and pneumonia. On the other hand, the novel coronavirus, causing coronavirus disease 2019 (COVID- 19) has precipitated a global pandemic since its first case in December 2019 and many of its underlying mechanisms and behaviors have yet to be confirmed [32]. COVID-19 is caused by SARS-COV2 and primarily attacks the respiratory system [33]. Symptoms range from none apparent to severe and lifethreatening, resulting in the need for support such as mechanical ventilation [34,35]. Co-infection with pathogens or preexisting medical conditions can be especially life-threatening. Bacterial pneumonia complicated more than 95% of the 50 million or more deaths during the 1918 influenza pandemic [36].
This review explores the mathematical modeling of various aspects of the biology summarized here. The first section focuses on acute infections that can be completely eradicated from the body given the proper immune response. We give a brief overview of COVID-19 models, but since in silico studies of the immune response to this disease are still in their infancy at the time of this review, we focus on the two infections, pneumonia and influenza. In fact, much of the morbidity associated with influenza is due to co-infection with pneumonia [37]. The bacteria which causes pneumonia, Streptococcus pneumoniae (Sp), can mutate and become resistant to antibiotics and the virus causing influenza can similarly mutate and develop new strains. In fact, antibiotic resistance is becoming an increasing concern in the global community and a significant threat to public health [38]. The need to adapt to evolving pathogens is a significant motivation for the development of mathematical models [39,40] [for more information on the pathogenesis of these infections, see [41][42][43][44] for influenza and [45] for Sp.].
Unlike the pathogens causing influenza and pneumonia, those that cause tuberculosis and anthrax may not necessarily be fully eradicated in the body. Instead, levels can be so low that they are not acted upon by the immune system and thus rendered inactive. This is called a latent infection, and no symptoms are exhibited although an active infection and subsequent transmission may result at any point [46]. Section 3 focuses on these pathogens.
Anthrax is a serious infection caused by the bacterium Bacillus anthracis. Inhalation of anthrax occurs when anthrax spores enter the body via the respiratory tract rather than via the digestive tract or wound. Although inhaled anthrax is the rarest form of an anthrax infection [47], it is of the greatest presentday concern due to its potential as a bioterrorism threat [48]. The immune response to anthrax is not completely known; however, it is known that pathogenesis is significantly different from the immune responses to other bacterial infections and thus requires a unique model.
The bacterium that causes tuberculosis, Mycobacterium tuberculosis, is thought to infect around one-third of the world's population. Although only 5-10% of these cases actually develop symptoms [49], tuberculosis (TB) was responsible for about 1.5 million deaths in 2015 and even more in previous years, making it one of the top causes of death in the world [31]. Although many advances have been made in the prevention and treatment of both TB and anthrax, the reasons individuals experience different outcomes is unclear, precipitating modeling efforts to better understand and treat these pathogens.
In addition to viral and bacterial infections, various immune responses can be elicited from the inhalation of dangerous particles due to smoking, air pollution, and occupational exposure [50]. The third section focuses on these insults.
Smoking can lead to Chronic Obstructive Pulmonary Disease (COPD), which is a largely irreversible degenerative lung disease characterized by chronic inflammation [51,52]. Another example of chronic inflammation occurs through repeated asthma attacks over time, resulting in airway remodeling. The severity of these attacks varies from patient to patient, rendering the disease difficult to manage clinically [53]. The immune response to damage induced by mechanical ventilation (MV) has also recently received attention, especially since it can exacerbate complications from the original need for MV [54,55]. For these conditions, the dynamics of the immune response are still incompletely understood, particularly the roles of macrophages and other inflammatory cells in airway damage and the resolution of inflammation in asthmatic and allergic episodes [56,57].
In the following sections, we elucidate the techniques used to model these injuries and insults. Many open questions still exist to better treat infections, especially under the threat of antibiotic resistance. Furthermore, chronic inflammation is of particular interest regarding lung cancer, since there is evidence that chronic inflammation contributes to carcinogenesis [58,59]. In addition, connections between smoking, COPD, and lung cancer have been drawn [60,61]. In this work, we examine the different causes of lung inflammation as described above and how mathematical modeling helps understand the inflammatory response in each.

Mathematical Background
The immune system spans the organ level down to the molecular level with positive and negative feedback loops and highly non-linear behavior. Mathematical modeling addresses scales that include subcellular and molecular interactions, communication between individual cells, interaction between various cell types and other signals at the tissue level, up to host and population levels. Various types of interactions within the immune system can be modeled using mathematical techniques, providing predictions and simulations for real-life scenarios and suggesting interventions, including development and testing of new personalized treatments [62,63]. Figure 1 gives an overview of the immune response in the lungs and highlights the mathematical modeling methods used to research diseases that trigger lung inflammation. Differential equations are commonly used for modeling biological systems because of their ability to effectively capture non-linear behavior and demonstrate interactions between numerous variables, measuring continuous populations over time. Each state variable represents a component of a system, such as a cell, protein, or a quality like overall damage [64]. Terms in each equation may represent rates of interactions between cells or mechanisms and thus have biological significance, which is useful when analyzing the system.
Ordinary differential equation (ODE) models are used to create kinetic models describing average interactions between large numbers of cells, molecules and/or individuals. ODEs typically describe dynamics over a continuous time period, assume that the environment is well-mixed and, therefore, do not account for spatial dynamics. The model can be calibrated by fitting to time series data. A variety of mathematical and computational techniques have been developed for ODE models, including uncertainty and sensitivity analysis, optimal control methods, bifurcation analysis, and parameter estimation and identifiability [62,65].
Some spatial features can be accounted for with compartmental ODE model. These have variables associated with cell or molecule populations in each compartment and their movement in and/or out of the compartments. For example, Marino et al. incorporated a lung compartment and a lymph node (LN) compartment into their model. This paints a more biologically realistic picture, accounting for migration of CD8+ T cells to the site of inflammation from the lymph nodes [66]. An ODE model is limited by the amount of data available; large and detailed data sets are required for higher confidence in the model and its parameters. Furthermore, the ability of ODEs to account for spatial dynamics is limited and the assumption of an environment in which populations are well-mixed is not always the case.
Delay differential equations (DDEs) are similar to ODEs, but incorporate time delays into the interactions between components. For example, Krishnapriya et al. incorporates a time delay into their model of H1N1 influenza to replicate the latency period, the time between an individual's exposure to the virus and the time when the individual becomes infectious. DDEs may result in a better fit for systems without modeling all components of the process causing the delay. Sensitivity analysis for the model by Krishnapriya et al. showed that this length of time may be important for determining the course of action in treatment [67]. Compared to ODEs, DDEs are more difficult to simulate and analyze [65].
Partial differential equations (PDEs) account for continuous time and space. Because of their multi-dimensional nature, they often require more computation time. However, PDEs have been used to track populations throughout a spatial domain, such as macrophages and cytokines in a granuloma [68] or T cells in LPSinduced experimental asthma. Stability and sensitivity analysis can be performed for PDEs [69], though this is more commonly done with ODE models.
Currently most of the differential equation models for lung inflammation are deterministic and do not include stochastic events or noise-induced variations in response. Therefore, output of the model is entirely determined by the parameter values and initial conditions, and the same results are obtained every time the simulation is repeated. Stochastic models, on the other hand, incorporate biological randomness this causes simulations with the same parameters and initial conditions to have different model output, which can sometimes vary significantly. Stochastic differential equations (SDEs) are differential equations that include one or more terms incorporating inherent randomness into the model. Stochasticity is most often applied to ODEs, though used in PDEs and DDEs as well. SDEs are more computationally complex and more difficult to parameterize than their deterministic counterparts, and to date have not been commonly used for immunological models [70].
Other stochastic models used to model lung inflammation include cellular automata (CA) models and agent-based models (ABMs). These are used to account for spatial interactions between immune components, but in a discrete manner rather than a continuous one. In both CA models and ABMs, cells, molecules, and other components of the immune response are identified as "agents." These agents interact on a defined two-or three-dimensional space based on a set of rules informed by the behavior of the biological system. The stochastic component is a strength of these types of models, since deterministic models such as ODEs and PDEs may lose accuracy when the molecules or cells being modeled have low total counts and randomness effects their interaction rates [62]. Using a set of simple rules, highly complex behavior can emerge.
The main difference between CA models and ABMs is that in CA models, agents do not move but rather interact and are changed based on their neighborhood (often a Moore neighborhood [71]) and other processes around them. ABMs include some, if not all, agents that move around the space, interacting with other agents and making decisions. The purpose of CA models is generally to detect patterns or structures that emerge from the rules, whereas ABMs seek to explore how changes in individual agents or interactions impact overall system behavior [72]. One limitation of ABMs and CA models is that they are computationally intensive. ABMs can incorporate a variety of agent types and scales but are difficult to parameterize or analyze mathematically. Spatial domains or cell counts are reduced, in order to have a reasonable number of agents. This can lead to results that are often more qualitative, such as examining overall tissue damage [50], as opposed to the quantitative population results of a differential equations model.
Boolean networks are discrete models in which individual units are represented as nodes and the interactions between them are edges, defined by Boolean logic. By simplifying interactions in this way, large models can be constructed with high-throughput data to discover underlying system properties [73]. For example, a molecule may be represented by 0 for "unactivated/off, " then 1 for "activated/on, " when bound to its receptor. The combinations of ON, OFF, AND, OR, and NOT can be combined for more detailed interactions. These models can be developed to include pathways with hundreds of interactions and connections. Boolean models may be pruned to exclude nodes and interactions that are discovered to be unimportant by the data and can incorporate multiple scales. Boolean models do not include spatial dynamics and do not capture non-linearities well [62,74]. Methods have been developed to analyze these large signaling networks. In Abou Jaoude et al. the authors reviewed these methods focusing on model attractors, reachability properties, reducing large models and assessing the effects on dynamics when varying external signals [75].
With the rise in computational power, models have become increasingly complex, allowing for multiple scales and model types to be integrated. Hybrid models contain two or more model types whereas multiscale models contain two or more scales; often models include both [76]. For example, Cilfone et al. incorporate details of the cytokines TNF-α and IL-10 with macrophages and other cell types to explore granuloma formation using an agent-based model for cellular and subcellular interactions. Additionally, ODEs are used to describe secretion of TNF-α and IL-10 from individual cells, and PDEs model diffusion of cytokines and chemokines throughout the lung environment. Hybrid and multiscale models are highly computationally expensive and difficult to calibrate. Optimization and sensitivity analysis are uncommon for these model types [70].
The amount and quality of available experimental data is typically a limitation for mathematical modelers because it limits their ability to accurately estimate model parameters. To account for this sparsity, advanced parameter estimation and uncertainty analysis methods have been developed [77,78]. Alternatively, models may be simplified to capture the overall behavior of the immune system with a few parameters.
After developing a mathematical model, there are some common themes when discussing future work. Models can be made more realistic by adding more detail to the current work; for example, specific immune cell roles can be explicitly modeled rather than using parameters [79]. Other directions for future work may include refining parameter values using different estimation techniques or additional experimental data and combining with other existing models or linking across scales [80,81].
It is worth noting that aside from replicating the dynamics of the immune system, many models focus on replicating the mechanics of the lung injuries themselves. Mathematical approaches to other topics in immunology have also been extensively covered, including T cell turnover and diversity [82,83], intercellular communication via cytokines [84], cancer and tumor growth [70,85], viral dynamics [41], and general immune system function [86].
In this review, we focus on mathematical models that examine the immune response caused by lung injury or disease. Each of the following sections, grouped by insult type, begins with a summary of the mathematical models used. We then choose a model to describe in greater detail; these illustrative examples showcase the methods utilized by many other models in the section and are representative of the overall direction of modeling in that field. By choosing a different model type for each section, we also provide the reader with examples of how mathematical and computational methods are implemented to illustrate immune system dynamics. Each section also includes a comprehensive view of conclusions from the models in which we link mathematical approaches to biological knowledge and identify patterns in model results. -INDUCED INFLAMMATION:  PNEUMONIA, INFLUENZA, AND COVID-19 Pneumonia, influenza, and COVID-19 are acute infections that can benefit from mathematical and computational modeling, especially given the ability of the influenza virus to mutate, the fact that much is still unknown about the novel coronavirus, and the comorbidities associated with all three. The models examined in this section aim to contribute to a more detailed and robust knowledge of the effects of various strains of these pathogens and the development of new and personalized treatments for the most susceptible individuals. Aside from the illustrative example, we outline models in this section, often grouping several together based on similarities. Tables 1, 2 focus on models for influenza and pneumonia, respectively, stating the type of model, summary of the key variables and parameters, and main results. Most mathematical models associated with COVID-19 at the time of this review are population-level spread of disease models. To our knowledge, the only model published so far on the immune response is an ODE model by Du and Yuan [121]. Wang et al. presents a community-driven effort to model COVID-19 [122], and a few other preliminary reports related to modeling the immune response to SARS-COV2 are available thus far [123,124]. Given the limited models at this point we do not include a table for COVID-19 models.

Overview and Models From the Literature
For influenza and pneumonia models, a majority are ODE-based models, which can encompass the diverse needs of the modeler. Many ODE models, especially within the context of influenza, simplify dynamics by including the virion and epithelial cells as variables. Then the severity of the infection is quantified by the amount of epithelial cells destroyed by the pathogen, tracked by one of the ODEs. The three-dimensional set of equations modeling the virus, target cells, and infected cells is commonly known as a target cell-limited model, and many of these ODE models are variations on the target cell-limited model. A model may not represent immune components by a separate equation but instead include the effects of the innate and/or adaptive immune system as parameters or terms within the equations; this provides added simplicity and computational efficiency, especially when there is a lack of data [99,104,108,[115][116][117][118]. For example, Miao et al. [6] and Smith et al. [94] describe the death of infected cells by a constant rate instead of explicitly modeling apoptosis by phagocytes such as macrophages. On the other hand, for a more biologically realistic model, facets of the immune response can be included as state variables to track the levels or concentrations of immune components during response to the infection. In these models, emphasis is commonly put on infected and uninfected cell populations [90,93,95,115,117], T cells [96,98,102], and macrophages [100,101,114].
Delay differential equations, though more mathematically and numerically complicated, may be used to represent delays such as the latency period (the time between entry of the pathogen into a cell and the resulting release of new pathogens from the host cell) or the time lag between production of interferon (IFN) and its effect on pathogen production [67,79,87,88,92,104,113]. In place of explicit delay differential equations, some modelers choose to simulate a delay by including an additional equation in an ODE model to represent different states of host cells, such as two separate populations of infected epithelial cells [90].
A few other differential equation models have been proposed. Schelker et al. considers the effects of spatial interactions through a PDE to understand how influenza is delivered into the nucleus [111] and Handel et al. uses a stochastic version of an SIR model to consider resistant strains of the virus [103].
Aside from these most common model types, Mitchell et al. developed a CA model and a set of ODEs and compared their results [107]. Murillo et al. and Smith et al. review models including differential equations and ABMs, anticipating the shift toward more complex, multi-scale models [109] and better incorporation of data [120]. A Boolean model was developed by Anderson et al. to reveal important subcellular signaling pathways within a dendritic cell [73]. We focus on summarizing more recent findings, though noteworthy contributions have heretofore been made. We select an example to illustrate one of the ways that models are used to shed light on lung inflammation. Results from a variety of models will then be discussed, demonstrating the advantages and limitations of various mathematical modeling methods.
The COVID-19 model by Du and Yuan [121] is a target cell-limited model, where the innate and adaptive immune responses are accounted for by parameters in the equations. Results from the model show that the timing of the innate and adaptive responses in relation to the time of viral load peak could be important in distinguishing mild and severe cases. Although modeling of COVID-19 is still in its early stages, understanding of the differences between this disease and other respiratory infections as well as knowledge of mathematical and computational tools to study them will be important in developing treatments and antiviral drugs.

Illustrative Example: Manchanda et al.
As mentioned previously, mathematical models can range from simple with a few dimensions to complex with multiple model types and scales. Manchanda et al. employed a simplified approach to the host response to influenza through a threedimensional system of ODEs [97]. A condensed model is advantageous when there is limited data available; such is the case here.
The system of ODEs is shown in Equations (1-3). The first variable P, is the pathogenicity of the virus, i.e., its ability to cause an infection. D, the overall host immune defense, includes both the innate and adaptive responses. I is inflammation due to the pro-inflammatory response. Figure 2 gives the model schematics [93] Infected cell subpopulations Reproductive number, rate of transition from infectious and productive Two-phase approximate solution captures the two distinct phases of the immune response; small changes in the eclipse phase length could lead to more drastic changes in viral growth rate. [94] Cells resistant to infection, IFN Infected cell death rate, antiviral efficiency Innate immune system is responsible for initial decrease in viral population and adaptive immune response is necessary for complete clearance. [95]
governing these interactions.
Pathogenicity, as represented by Equation (1), is controlled by the virulence of the virus (first term) and suppressed by the immune response (second term). Equation (2) shows that the immune response is activated in the presence of the pathogen and decays once the virus is cleared. Inflammation from the pro-inflammatory response, Equation (3), is characterized by the effects of cytokines and chemokines, modeled by the first term and taken from Kumar et al. [125]. The second term represents the anti-inflammatory response which counteracts pro-inflammatory behavior. The authors also include a "clinical score, " adapted from Smith et al. [126]. The clinical score (also called symptom score), a measure of how sick the host is, sheds light on the characteristics of an infection by quantifying its severity. It is calculated by: These were sufficient to characterize the four different strains studied, two of which are shown in Figure 2B. This model is also able to link viral pathogenicity with clinical symptoms.
To compare to the model clinical score, mice were assigned a laboratory clinical score ranging from 0 to 7 based on observation to gauge severity of the infection. For mice with a modelgenerated clinical score peak that occurred between 3 and 6 days, such as the Jena/5555 strain in Figure 2B, little inflammation was observed. Oh the other hand, high inflammation was seen for mice with a biphasic course of symptoms such as those with the Jena/5258 strain, also shown in Figure 2B. The authors speculate that this second peak in the biphasic response could be in part due to inflammation and high values for α, δ, and ǫ.
A small-scale model can direct future mathematical and experimental efforts in the study of influenza; the main drawback is that the explanation of outcomes is reduced to differences in inflammation. As we have described, the inflammatory response is highly complex and regulated; additional work is required to parse out the specific mechanisms of inflammation that contribute to these differences.

General Results and Conclusions From Modeling Influenza and Pneumonia
The conclusions from models of influenza and pneumonia generally fall under two broad categories: (1) a further understanding of the dynamics of the immune response and (2) recommendations for antibiotic and antiviral treatments, including in the case of coinfection.
A common thread throughout many of the influenza models is an emphasis on the role of the adaptive immune response in clearing the virus [79,88,91,92,95,96,98,102]. Handel et al. found that both the innate and adaptive, instead of only the innate, immune responses were needed in their model to accurately explain the data [92]. Similarly, Pawelek et al. concluded that although the considerable drop in the viral concentration after its peak is due to the innate immune response in their model, the complete elimination of virus is achieved by the adaptive immune response [95]. Other models focused on CD8+ T cells, main players in the adaptive immune response, because of their regulatory role via production of antiinflammatory cytokines throughout the course of infection [127]. Numerous models identified the presence and effectiveness of adaptive immune cells such as antibodies and cytotoxic cells in their model as vital to recovery time and viral load [79,88,91,96,102].
Another use of mathematical modeling is estimating parameters or values that have biological relevance and are otherwise difficult or even impossible to determine [40,90]. These parameter values can then help identify disease phenotypes or strains that could become pandemic [97]. One commonly estimated parameter is the reproductive number of the bacteria/virus or a similar threshold condition separating eradication of the pathogen and infection [67,87,107]. Other calculated values are the incubation period (the period of time between entrance of the pathogen into the host and recognition by the immune system) [106] and parameters related to the life of the pathogen such as pathogen half-life and bacteria doubling time [94,104,116].
Identifying sensitive parameters can also suggest effective treatment possibilities. Many important factors in the immune response to influenza can be quantified with only a few parameters, supporting the development of simplified systems. Based on identifiable and highly sensitive parameters, the model by Manchanda et al. revealed that adaptive immunity and early pro-inflammatory responses are some of the most impactful factors in the dynamics of influenza [97]. Baccam et al. calculated the reproductive number R 0 which represents the number of new infections that an infected cell can produce [90]. Kamal et al. utilized a small-scale model, finding that treatments influencing the parameter representing viral infectivity may be effective [99]. A model developed by Smith et al. indicates that differences in outcome could be explained by small changes in the mechanism which controls bacterial-induced impairment of macrophages [115].
Often different pathogenic behaviors arise from models, and modelers categorize these behaviors to understand what may cause specific outcomes and how an individual might transition from a healthy state to an infected one. Some models delineate between the eradication of the pathogen and sustained growth [87,114,117,119] whereas others go into greater detail by identifying different characterizations of the infections such as low and high pathogenicity [80,101]. For example, Almocera et al. identifies three outcomes of pneumonia based on bifurcation parameters: dose-independent clearance, in which the infection is cleared by the immune system regardless of the bacterial dose size; dose-independent persistence, in which the immune system is unable to remove the bacteria regardless of dose; and dose-limited clearance, in which the bacteria can be cleared depending on the dose size [119]. Comparing these distinct outcomes can identify possible processes or parameter thresholds responsible for the different behaviors.
As previously mentioned, mathematical modeling provides the ideal platform to test optimal dosage and hypothesize treatments before expensive and time-consuming experiments are performed. Beauchemin et al. [104] and Canini et al. [110] consider viral resistance and its effect on treatment efficacy by including random drug-resistant mutations. Models also vary intervention time and duration [80,88,99,103,113] and determine which mechanisms the model predicts to be most useful for developing new treatments [88,108,112]. With regards to intervention time, some models recommend starting treatment no later than 2 days after symptoms appear for maximum efficacy [88,101] though exact windows for effective treatment depend on various factors, such as type of intervention.
For example, in Schirm et al. a sum of time-dependent pulse functions was used to model antibiotic injections [81]. The authors found that, beginning 24 h after infection with Sp, doses of the antibiotic Ampicillin every 12 h was effective, whereas every 48 h was insufficient. Furthermore, they found a threshold of the initial number of bacteria which determined whether or not the immune system by itself could eradicate the infection. This could be useful in determining the least amount of antibiotics necessary to recover from pneumonia, decreasing chances of antibiotic resistance. Focusing on the cells affected by the infection, the model by Domínguez-Hüttinger et al. used the Boolean concept of "switches, " described as either on or off, to represent a pathway as activated or not and determine by a threshold value whether the immune system can clear the pathogen without treatment [80]. Based on the different pathogenic behaviors, at least one switch, i.e., pathway, needs to be turned off to eliminate the pathogen. Results from the model suggested that treatment within 36 h is extremely important in preventing sepsis; this was consistent with experimental observation.
Coinfection is of particular interest because of the potential threat it poses when patients already infected with the influenza virus are also infected with a bacterial infection, especially during an epidemic or pandemic. Models explore the effects of antiviral therapy in preventing the secondary infection [98,113,115,117,118]. Smith et al. emphasize the role of alveolar macrophages and how their impairment by virions affects the host's ability to fight the secondary bacterial infection [115,117]. Shrestha et al. also examined the effects of antiviral treatment on growth of Sp. in an ODE model and found a 4-day window in which antiviral treatment can reduce the severity of a secondary infection [113].
Many authors stated that future work could involve including mechanisms not explicitly modeled in their current work [79,98,106,110,115,118,128], refining parameter values [80,81,91,97,114,128], and combining their models with other existing models [80,81]. Additional considerations that would be helpful in a greater understanding of these infections include the effects of aging [98] and spatial information [109].

LATENT VS. ACTIVE INFLAMMATION: TUBERCULOSIS AND ANTHRAX
Anthrax and tuberculosis both exhibit latent infections, in which levels of the pathogen are so low they go undetected by the immune system but can spontaneously resurface as an active infection. The main purpose of the models in this section is to further understand the complex dynamics of the infections, providing insight into the effectiveness of antibiotics and development of other treatments, especially under the threat of growing antibiotic resistance. Tables 3, 4 provide details and frameworks of numerous tuberculosis and anthrax models, respectively, with summaries of model trends and results in the following sections. The pathogenesis of tuberculosis has been an ongoing subject of computational and mathematical modeling for decades; several thorough reviews on modeling TB have been published [164,176,180,181] including recent ones from Kirschner et al. [76] and Cardona et al. [182].

Overview and Models From the Literature
TB modeling has benefited greatly from an increase in computational power, leading to highly complex hybrid models. ODE models of TB generally cover cellular and subcellular interactions, especially with macrophage phenotypes, T cells,  [154]. PDE-only models have also been developed, modeling the spatial structure of granulomas, collections of immune cells and bacteria vital in fighting infection [68,162]. For a different spatial approach, some modelers chose an ABM or CA model [141][142][143][144][145]147]. All of these incorporate dynamics on a cellular and subcellular scale, and many are concerned with granuloma formation. Several, including those from Ray et al. [142], Cilfone et al. [143], and Bru and Cardona [147], are particularly interested in how the granuloma is affected by various cytokines. A few other methods have been used, including Boolean for cellular and subcellular interactions [146] and signaling pathways within a cell [163]. Additionally, Pienaar and Lerm used the computational Markov Chain Monte Carlo (MCMC) method to simulate molecular level interactions in which variables were discrete and stochastic but did not include a spatial component [148].
Multi-scale hybrid modeling has gained traction particularly in the area of tuberculosis modeling. Combinations of ABMs, ODEs, and PDEs are the most common and range from granuloma formation [165-168, 171, 178, 179] to properties and efficacy of antibiotic drugs [170,173]. Generally a tissue-scale ABM is integrated with ODEs or PDEs modeling for molecular dynamics. The development of these multi-scale models will help to answer open questions about how changes on a cellular or subcellular level can impact population-level dynamics; the time scales range from minutes to decades, and the complexity of these models will introduce challenges but also great opportunity for new treatments [164].
Anthrax models are primarily ODEs, exploring the interactions between anthrax spores and immune cells [48, 133-135, 137-139, 184]. Diverse mathematical tools are used to model migration of spores, including the use of compartments [48,134] or an additional equation [136]. A few other approaches have been used: a stochastic adaptation of an existing ODE model [140], use of a previously-develop immune system-simulating CA framework [129], and a probabilistic model called a competing risk model [130,131,139]. See Table 3 for more detailed descriptions of these models.

Illustrative Example: Marino and Kirschner
In the previous illustrative example, we explored the methods and benefits of a models with a few variables. These are particularly [141]

Macrophage and T cell subpopulations
Bacterial growth rates, T cell movement, TNF-α production In the context of granuloma formation, TNF-α plays multiple roles in control of pathogen. [142] Macrophage and T cell subpopulations, IL-10, TNF-α

IL-10 synthesis, TNF-induced apoptosis
Multi-scale approach explored individual roles of cytokines, particularly IL-10, in granuloma formation and overall immune response.
[  Integration of model and experimental data shows that a change in inoculum size results in higher immune cell levels but had little effect on T cell priming. [159] Exposed and infected granulomas, healthy tissue Adaptive immunity sensitivity to antigen stimulation, bacteria dissemination Population-level study of granulomas predict mechanistic differences between latent and active infection. [160] Iron, lipid content, NO Iron overload inhibitory effect, bactericidal effect of nitric oxide Early host response depends on ability of pathogen to interact with iron and lipid pathways as well as oscillatory progression of disease.
[ [164] Macropahges, T cells, TNF, TNF-α receptors STAT1 activation of macrophages, bacterial growth rate, T cell recruitment TNF-α receptor 1 plays a critical role in control of infection through granuloma function. [165] T cell and macrophage subpopulations

APC delay during migration, trafficking of T cells in lymph node
Multi-organ model couples diffusion of cells between compartments with spatial dynamics of granuloma; underscores importance of APC movement. [166] Macrophage and T cell subpopulations, NFkB signaling pathway components NFkB-macrophage ratio, IKKK activation, A20 and IkBa mRNA synthesis Infection dynamics across molecular, cellular, and tissue scales is heavily mediated by TNF-induced NFkB signaling in macrophages; NFkB processes are recommended as treatment strategies. [167] TNF, macrophages, T cells TNF-mediated macrophage activation and cytokine secretion Explored dynamics of granuloma formation and effects of inhibitor treatments. [168] Macrophage and T cell subpopulations

TNF-α receptor internalization, antibiotic regimen
Review: models provide the ability to perform virtual experiments and clinial trials, reducing cost and risk. [169] Macrophage and T cell subpopulations, TNF, IL-10 Host cell recruitment, TNF-induced apoptosis, immune cell activation Granuloma formation model with antibiotic pharmacodynamics revealed that pre-treatment host factors greatly impact treatment outcomes. [170] CD4+ and CD8+ T cells, antigen presenting cells

Basal influx & efflux, differentiation rates
Emphasized pairing of mathematical and experimental methods to best discover TB biomarkers. [171] Oxygen levels, macrophage subpopulations, metabolic pathway components

Dissemination, oxygen & chemokine diffusion
Correlations between immune response, physiological response to hypoxia, and metabolic networks revealed specific mechanisms that classified different outcomes. [172]

T cells, DCs, peptide-MHC complexes
Binding probabilities, DC-T cell binding time Recommendations for vaccine formulations can be determined through modeling of T cell events and interactions. [173] ABM, PDE Macrophages, chemokines, TNF, IL-10 Synthesis rates of TNF-α and IL-10, granuloma size Specific processes describing cytokine synthesis, signaling, and spatial distribution control TNF-α & IL-10 levels and determine outcomes. [174] ABM, ODE, PDE

Macrophages, T cells Recruitment of macrophages and T-cells by chemokines, spatial movement of immune cells
Comparison of four different model types reveals importance of spatial dynamics and macrophage activation. [175] APCs, T cells, cytokines (TNF, IL-10) Switching time, cytokine production, receptor binding Review of computation models of host-pathogen interactions; importance of multi-scale approaches. [176] TNF-neutralizing drugs, macrophages, T cells Drug-induced cell death, drug binding to membrane-bound TNF-α Analysis of TNF/TNF-α receptor roles predicts mechanisms to improve host response during anti-TNF-α therapy, especially since drug binding to TNF-α can impair granuloma function. [177] Macrophages, T cells, TNF, IL-10 Rates in NFkB pathway, ratio of TNF-α to IL-10, memory cell initial conditions Emphasized importance in considering multiple scales when designing a treatment. [178] (Continued) Frontiers in Applied Mathematics and Statistics | www.frontiersin.org Incorporation of DCs into hybrid multi-scale model reveals importance of early events in immune response such as DCs recruitment and T cell proliferation. [179] Bacteria nutrients, cytokines, chemokines, antibacterial drug Antibiotic properties, ratios of immune cells Review of models of within-host response to TB highlighted integration of different types of data and testing of antibiotic regimens. [76] ABM, ODE Microbe, metabolites TNF-α & macrophage activation Review of systems biology approaches, including high-throughput methods, to understanding TB. [180] Boolean, ODE, ABM TNF-α, macrophages Immune system mechanisms, vaccine effects Compared various types of models, enumerating advantages and disadvantages for active and latent infections. helpful when little data is available. Increases in computational power have allowed for added complexity in models. This can be seen in TB modeling, in which models have built on one another over the years, increasing in detail and scale. Furthermore, TB has been widely studied and extensive data has been collected. Thus, larger models can be developed and calibrated based on this information. In this section, we highlight one such model. A hallmark of the immune response to TB is the development of granulomas, spherical concentrations of immune cells designed to isolate and kill the bacteria. The authors had previously developed an ABM, GranSim (http://malthus.micro. med.umich.edu/GranSim/), to model the immune response to TB, including the formation of granulomas and traffic to and from the lungs, bloodstream, and lymph nodes (LN). In this work they further expand the model to explicitly include dendritic cells, then calibrate the model to data and investigate the mechanisms driving behavior in granulomas [179].
The version of GranSim used in this illustrative example is a hybrid multi-compartmental model, shown in Figure 3A. The ABM models the formation of a single granuloma in the lung, including interactions between immune cells such as macrophages, CD4+ and CD8+ T cells, signaling molecules, and Mtb, the bacteria that causes TB. In previous versions of GranSim, DCs were represented by a more general APC agent; in this one, DCs have a role unique from other APCs in that DCs are immediately stimulated upon coming into direct or indirect contact with Mtb. Additionally, macrophages do not leave the grid whereas DCs migrate to the LN for T cell priming.
The authors then linked the single granuloma model to the host scale; this connection is important to understanding how cellular-scale dynamics can affect the overall outcome of the host. The number of granulomas in a host range from 13 to 97 at 4 weeks. A scaling factor, scalingMDC, is introduced which multiplies the number of DCs migrating from a single granuloma by the total number of granulomas in the host. The total number of DCs is then passed to the bloodstream and LN compartments, modeled by a set of non-linear ODEs. This number directs the priming of T cells and their recruitment back to the lungs. The model can also control the time it takes for the DCs to reach the LN.
The model was calibrated using data primarily from humans and non-human primates (NHPs). Data was collected in the blood for T cell dynamics and in the lung for bacteria levels. The use of detailed spatial and temporal data in multiple locations is critical for a more complex model such as this, and allows for greater confidence in model results. Figure 3B shows model simulations compared to data for colony-forming units (CFU), i.e., number of bacteria, per granuloma. The stochasticity of the ABM allows the model to cover the ranges observed experimentally over time. Figure 3C shows individual snapshots of two such granulomas simulated by the model.
For uncertainty and sensitivity analysis, the authors used Latin Hypercube Sampling (LHS) to uniformly and randomly sample the parameter space and Partial Rank Correlation Coefficient (PRCC) to calculate the association between parameters and model output. PRCCs can be calculated over time to examine whether some mechanisms are more important than others at specific times.
Analysis shows that T cell priming and proliferation and the timing of DC movement to the LN are critical functions. The parameter k 13 , representing CD8+ T cell proliferation, generates high PRCC values for all inflammation outcomes tested, including intra-and extracellular Mtb burden and total activated macrophages. The authors also found that a delay in DC movement to the LN had a negative impact on memory T cell roles and that higher levels of DCs pre-infection resulted in better outcomes.
The authors concluded that timing is key when it comes to both the innate and adaptive immune responses to Mtb. Early protection against the bacteria results in better outcomes, and specific recommendations for intervention were strengthening DC migration and CD8+ and CD4+ T cell trafficking. The inclusion of DCs into the model allowed for new conclusions to be drawn about mounting a timely and effective immune response. Future work mentioned included an ABM formulation of LN to replicate vaccine trials and test immunotherapy strategies. By using a multi-scale model, the authors were able to examine specific details of the immune response mechanisms driving changes in model output, connecting cell behavior to overall host response. Sufficient data over multiple scales, such as that obtained by Marino and Kirschner, is needed to make effective conclusions; this should be considered when one is choosing a modeling technique.
Granuloma formation has been widely studied in the context of TB as one of the main methods used by the immune system to eliminate bacteria. Through the use of multi-scale modeling, not only granulomas as a whole but also individuals macrophages and their functions have been studied at length [66, 68, 129, 135, 137, 138, 146, 149-151, 155, 175, 188, 189]. Because granulomas are localized clusters of cells, spatial distribution was found to be a significant factor in infection outcome [156,174]. To examine molecular interactions in granuloma, Fallahi-Sichani et al. explored the dynamics of the NFκB signaling pathway across molecular, cellular, and tissue scales in their ABM [167]. Results highlight the need for a balance between activated and infected macrophages, responding to and regulated by NFκB signaling. Several roles of TNFα within the granuloma were also found to be critical to outcome, including interaction with its receptor [165] and polarization of macrophages via TNF-α [66,141,142,144,156].
Additionally, activation of and phagocytosis by macrophages are among the most important mechanisms of host defense [49,68,138,146,149,150]. The spectrum of macrophage phenotypes from pro-inflammatory to anti-inflammatory has been studied more recently, including in TB research. Specifically, the "switching time, " the time point at which M1 (classically activated) macrophages becomes more dominant than M2 (alternatively activated) macrophages was a useful way to measure the M1-M2 balance. Day et al. found that an increase in the parameter k 16 , which controls the ability of infected AAM to inhibit bacterial growth, decreased bacterial loads but increased the switching time [49,157,162,176]. Marino et al. adapted a similar approach, in which AAM immediately become infected when they engulf a bacterium and CAM kill the bacterium [66]. Models also predicted results of antibiotics and other treatments' interaction with granuloma structures [170,177].
Coupling multiple scales can elucidate how signals and changes on a molecular scale can influence the immune response on a cellular or tissue scale [76,143,145,164,165,167,168,175,176,178,179,190]. Another widely-accepted practice is the consideration of several locations in the body affected by TB through a multi-compartmental model [48,134,160,171,178,180,189]. This is helpful in understanding the interaction between tissues and recruitment of immune cells. Results from bringing together various scales, organs, and tissues have elucidated the processes that regulate TNF-α and IL-10 concentration [143,165,167,174,178,180], the importance of infection rate in outcome [145,175,189], and the role of cell-to-cell interactions [176,179].
Since latent infection is a significant issue and can develop into active infection at any time, many models investigated this disease scenario. Results suggest the important role of T cells [155,186], the phagocytic capability of immune cells [149,150], the presence of TNF-α [154], and the identification of specific proteins as drug targets [163,181]. Other models provide suggestions for the treatment of latent infections such as TNF-α inhibitors or antimicrobial drugs [153,154]. For example, Agliari et al. found an exact solution to their differential equations model which allowed them to predict the number of latent infections that would become active for patients undergoing anti-TNF therapy [185].
Signaling via cytokines, especially the pro-inflammatory cytokine TNF-α, has been studied extensively in the context of tuberculosis and to some extent in anthrax [133,135]. With TNF inhibitors as a current treatment, models recommend adjusted dosage and timing for better results for both active and latent infections [154,170,177,178]. The balance between TNF-α and IL-10, an anti-inflammatory cytokine, was also predicted to be a significant regulator of Mtb [66,143,146,149,159,174,178]. TNF-α has been the subject of a host of other TB models [49, 66, 142, 143, 165, 168-170, 174, 176-181] because of experimental evidence surrounding its important role. Virtual depletions, deletions, and clinical trials have also been studied in several models as a cost-effective path toward more personalized treatment and unique intervention strategies [66,149,169].
Sensitivity analysis revealed that common important values were bacterial replication and decay rates for both TB [68,141,145,172] and anthrax models [135][136][137][138]189]. These rates could often determine in the model whether bacteria persisted or were eliminated by the immune response [68]. For example, the model by Segovia-Juarez et al. showed that a slower bacterial growth rate within infected cells is worse for the host [141]. Bacterial replication and decay rates as well as other model components identified distinct pathological states including early clearance, granuloma formation, and persistent infection [68,140,150,152,160,165,166,172,179,188].
Future directions in TB modeling include further in silico testing of antibiotics and other treatment strategies, study of antibiotic resistance, simultaneous modeling of multiple granulomas, and continuation of multi-scale modeling for a more complete and integrated understanding of TB [76,181]. In the less-developed field of anthrax models, future work revolves around adding complexity such as more detail in host immune function [134] and heterogeneity and stochasticity [137,191]. Recommendations for treatment as informed by the models include a faster initial response by the immune system, potentially achieved through vaccination. Response from classically activated macrophages and T cells, which reduce bacterial loads earlier in the infection, are especially significant. These conclusions were achieved through using different types of models, highlighting the importance of diversity in modeling techniques.

MODELS OF NON-INFECTIOUS INJURY
From smoking to COPD and in asthma, the complex interplay from molecular to tissue scales in non-infectious injuries is unclear. Many models for these types of injuries have been developed but primarily focus on biomechanics and general inflammation [54,[192][193][194][195][196][197]. However, there have been some models for non-infectious injury that have explicitly studied the immune response.

Overview and Models From the Literature
The wide variety of injuries covered in this section necessitates a number of modeling methods. Reviews of mathematical and computational models of immune system mechanisms involved in fibrosis, autoimmunity, COPD, and asthma highlight the need for multi-scale approaches and greater computational power [198][199][200]. ABMs and CA models lend themselves well to this end as many non-infectious diseases and disorders are associated with inhalation and diffusion of particles. Furthermore, they can be integrated with other types of mathematical models such as ODEs and PDEs [54,55,201,202].
Simpler models still hold significant explanatory power. ABMs have been used in particular to track particle inhalation [50,53,56] and CA models have incorporated strain and its effect on the immune response [54,203]. ODEs [52,204,205], and PDEs [69,206,207] track populations of immune cells and concentrations of molecules such as CAM, AAM, IFN-γ , and TNF-α. Summaries of models for the immune response to noninfectious lung injury are found in Table 5.

Illustrative Example: Brown et al.
The level of detail described and approaches used by the model in this section is both representative of the types of models utilized for non-infectious insults and also provides an alternative approach to the two models described in the previous sections. Many models for non-infectious injuries are ABMs or CA models and examine different interactions between immune cells and signaling molecules. Whereas, Manchanda et al. considered inflammation as a single variable and Marino and Kirschner described multiple subpopulations of immune cells, Brown et al. finds a level of complexity between these two to identify three distinct pathologies for particle inhalation, highlighting the need for personalized treatments.
The ABM by Brown et al. [50] is a system with abstracted populations of cells to develop a higher-level view of events instead of an exhaustive description of every interaction. The Review of inflammatory mechanisms of airway diseases including COPD, asthma, and cystic fibrosis highlights findings from mathematical and computational models. [199] ABM, agent-based model; AAM, alternatively-activated macrophages; CA, cellular automata; COPD, chronic obstructive pulmonary disease; CAM, classically-activated macrophages; EF, elastic fragments; LPS, lipopolysaccharides; ODE, ordinary differential equation; PDE, partial differential equation; Th, T-helper.
authors noted that in reality, multiple types of cells could fall under the category of one variable. For example, TNF-α and TGF-β1 represent all pro-and anti-inflammatory cytokines, respectively. The schematic is shown in Figure 4A.
A simulation begins when particles are "released" into the system, initiating the immune response. Some examples of rules in the model include: • Macrophages on the same patch as a particle become activated for 50 time steps. • Activated macrophages release pro-inflammatory cytokines and for the last 5 time steps of activation, they release antiinflammatory cytokines. • Fibroblasts move toward and heal damaged tissue by depositing collagens and increasing the "Tissue Life" value. A trademark of mathematical modeling is quantifying events that are difficult or even impossible to measure experimentally. The model in this example, as well as several described in previous sections, assess tissue damage as overall health. Here, damage is based on a value, "tissue life, " which measures health on a scale from 0 to 100, where 100 is perfect health with no damage. Each patch has a health value, decreased by the presence of proinflammatory cells on that patch, and increased by the presence of anti-inflammatory cells. The average over all patches is the overall tissue life. Because ABMs focus more on qualitative representations of a system, tissue life is an important variable to quantify results. On a 50-by-50 grid, Brown et al. performed simulations with exposures of 5, 10, 15, or 20 particles occurring at intervals of 50, 100, or 200 time steps. To model different frequencies, particles were added into the model 50, 100, or 200 times. These conditions were meant to test theories that both frequency and quantity of particle inhalation, such as how often an individual smokes cigarettes vs. the number of cigarettes smoked, play a role in the development of chronic lung disease. Risk factors may include smoking, air pollutants, and occupations such as coal mining; studies have examined the differences between these insults and how the lungs and immune system respond [208][209][210]. Based on the degree and duration of particle exposure, the model elicited three distinct states, shown in Figures 4B-D. The model suggested the most damage is due to increased particulate exposure at higher frequencies, even with short duration.
Categorization of outcomes is common in mathematical modeling and allows for quantification of states through specific parameter ranges or combinations. Brown et al. showed that this can be achieved using lumped variables, another frequentlyused technique. In addition, they remarked that a well-mixed model such as a system of ODEs may not have identified these states, pointing to the need for spatial effects to be considered in modeling. A limitation of this model is that it is not calibrated to data. Also, the inclusion of more detailed mechanisms and interactions may further shed light on inflammatory activity.

General Results and Conclusions From Models of Non-infectious Injury
Although many of the key players in the immune response, such as macrophages, T-cells, and neutrophils, appear regardless of injury type, the progression of damage and the mechanisms by which the immune system brings about healing is very different depending on the type of damage.
Chronic obstructive pulmonary disease (COPD) is a disabling, irreversible disease often brought on by smoking and characterized by excessive inflammation and limited airflow [211]. A focus on macrophages, including the M1-M2 spectrum of activation, is a common denominator among COPD models [52,57,201]. Models by Cox [52] and Ceresa et al. [201] both identified positive macrophage feedback loops that contributed to inflammation.
In asthma, memory T cells are developed in response to allergens and are the focus of many mathematical models of asthma [204,206,207,212]. Models also identified important regulators of the immune response to asthmatic episodes including the cytokines IL-13 [204] and IL-4 [207].
Pulmonary fibrosis is a lung disease caused by an overremodeling of lung tissue after damage. Ben Amar et al. and Warsinske et al. are in agreement that a multi-scale approach, including cellular and molecular scales, is an important step toward a more complete understanding of fibrosis, both in the lung and other tissues [200,202].
Mechanical ventilation provides support to patients with difficulty breathing but can cause stress and strain resulting in an immune response. To this end, models were developed to better understand how strain levels and distributions interact with the immune system. Results from the models show that strain distribution plays an important role in immune cell populations [54,203] but that strain rate was not as important as parameters related to the immune response [55].
Future research in modeling these types of injury includes multi-scale modeling for better treatments involving multiple targets [55,199,213] and combining data-driven and mechanistic models [214].

CONCLUSION AND THE FUTURE OF MODELING
Mathematical modeling is a method for obtaining a deeper understanding of the immune response to lung injuries, complementing and informing experimental and clinical strategies. By summarizing methods and providing specific examples in the literature, we show how unique additions to well-documented modeling practices and novel, computationally intensive methods revealed new information about the pathogenesis of infections, recovery from damage, and the multi-scale complexity of lung dynamics.
Research on bacterial and viral dynamics has progressed greatly due to mathematical and computational modeling efforts, elucidating immune system mechanisms for pathogen clearance and quantifying important parameters in the viral life cycle. Mathematical research of the effects of aging on the immune system, coinfections, and use of vaccinations could improve the treatment of influenza. Complex mechanisms such as activation of different phenotypes of immune cells like macrophages should be studied further, as well as the positive and negative feedback loops that, when dysregulated, can exacerbate inflammation or contribute to the reappearance of infection [62]. ABMs, "omics" data, systems-based approaches, and high-performance computing have been suggested for future modeling efficacy [98,109,112].
With the introduction of COVID-19 onto the global stage, the need for mechanistic modeling has never been greater. Smallscale models such as that by Du and Yuan [121] borrowed from influenza models, making use of limited data. Though computational tools such as machine learning algorithms have become useful in detecting patterns and important predictors, they require large amounts of data that are not yet available for COVID-19. Therefore, mathematical modeling of the immune response to SARS-COV2 to better understand the unique mechanisms involved will be vital in the accelerated path toward developing vaccines and treatments.
The complex, non-linear interplay between cells, cytokines, proteins, and other host factors in the immune response to an insult has led to significant development of multi-scale models over multiple biological and time scales [164,180,200]. Many successful complex models were built upon previous models and then calibrated with experimental data. However, many modeling efforts begin with models that have few variables and/or group similar cell types and/or mediators. This approach is often used when there is limited data available. The resulting models are typically small-scale models that have the advantage of being mathematically analyzable in order to gain insight in the biological problem using stability analysis and by finding bifurcations [117,186]. These models have been used as the foundations for building more complex models. Mathematical modelers are more and more pairing classic mathematical methods with computational approaches to develop models which are calibrated and validated with data. These efforts from various fields have led to open-source platforms [129,[215][216][217].
Mathematical modeling will provide frameworks to test new hypotheses, analyze phenomena, perform virtual clinical trials, inform experimental strategies, and recommend general and personalized therapies. Mathematical and computational modeling is playing an increasingly important role in science as a whole, and, as shown in this review, can contribute greatly to cutting-edge treatment of a wide variety of lung injuries.

AUTHOR CONTRIBUTIONS
SM, AR, and RH designed and developed the project. SM compiled the information and wrote the manuscript under the supervision of AR. RH consulted on the biological aspects of this work. AR and RH reviewed and edited manuscript.

FUNDING
This work was supported by the National Science Foundation via award CMMI-1351162 (RH).