Computational and Theoretical Analysis of Human Diseases Associated with Infectious Pathogens
View this Special IssueResearch Article  Open Access
Zhenzhen Shi, ChihHang J. Wu, David BenArieh, Steven Q. Simpson, "Mathematical Model of Innate and Adaptive Immunity of Sepsis: A Modeling and Simulation Study of Infectious Disease", BioMed Research International, vol. 2015, Article ID 504259, 31 pages, 2015. https://doi.org/10.1155/2015/504259
Mathematical Model of Innate and Adaptive Immunity of Sepsis: A Modeling and Simulation Study of Infectious Disease
Abstract
Sepsis is a systemic inflammatory response (SIR) to infection. In this work, a system dynamics mathematical model (SDMM) is examined to describe the basic components of SIR and sepsis progression. Both innate and adaptive immunities are included, and simulated results in silico have shown that adaptive immunity has significant impacts on the outcomes of sepsis progression. Further investigation has found that the intervention timing, intensity of antiinflammatory cytokines, and initial pathogen load are highly predictive of outcomes of a sepsis episode. Sensitivity and stability analysis were carried out using bifurcation analysis to explore system stability with various initial and boundary conditions. The stability analysis suggested that the system could diverge at an unstable equilibrium after perturbations if (maximum release rate of Tumor Necrosis Factor (TNF) α by neutrophil) falls below a certain level. This finding conforms to clinical findings and existing literature regarding the lack of efficacy of antiTNF antibody therapy.
1. Introduction
Sepsis, currently defined as a systemic inflammatory response (SIR) to an infectious agent or trauma, is increasingly considered an exaggerated, poorly regulated innate immune response to microbial products [1, 2]. Under health conditions, intruding pathogens are eliminated by immune cells in the immune system. If overwhelming immune response occurs, unbalanced responses between immune cells lead to unexpected harmful patient outcomes such as high fevers, flushed skin, and elevated heart rate, resulting in sepsis. Possible progression to severe sepsis is marked by generalized hypotension, tissue hypoxia, and Organ Dysfunction [1]. Severe sepsis can further develop into septic shock under longlasting severe hypotension [3], ultimately leading to death.
Severe sepsis and septic shock during an infection are the primary causes of death in intensive care settings [4]. On average, sepsis causes 250,000 deaths per year in the United States [5]. Among patients in intensive care units (ICUs), sepsis is the second highest cause of mortality [6] and the 10th leading cause of death overall in the United States [7]. An average of 750,000 sepsis cases occur annually, and this number continues to increase each year [6]. Care of patients with sepsis can cost as much as $60,000 per patient, resulting in a significant healthcare burden of nearly $17 billion annually in the United States [8, 9]. Sepsis development in a hospitalized patient can lead to extended hospital stays and consequently increase financial burdens. Cross and Opal [10] pointed out the lack of rapid, reliable assays that could be used to identify the stage or severity of sepsis and to monitor the use of immunomodulatory therapy. However, no such assays are available because complexity of the inflammatory response and the unpredictable nature of septic shock in individual patients render the effect of targeting isolated components of inflammation with supportive therapy difficult to predict [10, 11].
Development of a nonbiased, predictive model and modelderived policies that prevent patients from experiencing severe consequences of sepsis (e.g., Organ Dysfunction) is critical for improving ICU patient care. As studies of mechanisms leading to sepsis development significantly progress due to discoveries of new inflammatory proteins and increased knowledge of the interaction of host cells and pathogens, mathematical models have been developed as dynamic knowledge representation of complicated biological processes. Specifically, the models have been used to simulate dynamic patterns of selected essential indicators in disease progression by integrating cellular and molecular pathways in an immune system. These mathematical models offer potential for understanding complex dynamic systems and, therefore, are used by researchers from various fields to simulate immune response to specific disease [12–15]. Development of modeling techniques could allow novel strategies for disease treatment, oriented at compromising harmful effects of inflammatory responses, to be proposed or tested in model simulations.
In order to construct a mathematical model of sepsis, we searched literatures and found two representative system dynamics mathematical models (SDMMs) of Acute Inflammatory Response (AIR) in previous studies. In 2004, Kumar et al. [12] presented a simplified 3equation SDMM to describe mathematical relationships between pathogen, early proinflammatory mediators, and late proinflammatory mediators in sepsis progression. In 2006, Reynolds et al. [13] proposed a mathematical model for AIR that included a timedependent, antiinflammatory response in order to provide insights into a variety of clinically relevant scenarios associated with inflammatory response to infection.
2. System Dynamics Mathematical Model Development
Existing mathematical models focused on inflammation in the literature proved that mathematical modeling is a valid approach for simulating disease progression [12–16]. However, the number of variables used, the limited control of system parameters, and the inclusion of many variables involved in real immune response were not modeled in detail. Therefore, oversimplification in AIR models [13, 15] limited AIR behaviors and biological relevance of simulated results. For example, simulated results from AIR models [13, 15] failed to capture a dampened oscillated infection in AIR progression. In addition, existing mathematical models are incomplete representations of sepsis because simulated AIR in both mathematical models [13, 15] is considered to be an initial stage of sepsis progression. Therefore, to improve on current models, we developed an 18equation SDMM to incorporate the most influential variables for septic response development during innate immune response and adaptive immune response. We included equations to represent pathogen load, phagocyte (including neutrophil and monocyte) activation, early and late proinflammatory cytokine release, tissue damage, antiinflammatory cytokine release, CD4+ T cell activation, CD8+ T cell activation, B cell activation, and antibody release. Indicator selection was based on knowledge of cellular and molecular pathways of sepsis from experts in the field and extensive literature review [4, 17–28]. We chose Salmonella as a “targeted” pathogen strain in our mathematical model and simulated immune responses to Salmonella in the liver of mice. Immune responses to Salmonella infections have been investigated extensively in [20, 29–33]; therefore, an abundance of data exists for accurate incorporation of relationships among variables to support our SDMM. We used a series of known and hypothesized kinetics of biological system components, including conventional logistics function, law of mass action, and MichaelisMenten kinetics to build SDMM from subsystems and mimic interactions between indicators. We combined these formulated but generalized dynamic modeling techniques into a comprehensive SDMM framework to describe sepsis progression, by measuring the steady state of various components in inflammatory responses. In the following seven subsections, we present a detailed description of mathematical construction for each subsystem in a mouse hepatic inflammatory response during SDMM development.
2.1. Process Description
AIR typically occurs when immune cells, such as tissue macrophage, detect intruding pathogens or existing tissue damage and emit a signal to resting phagocytes, such as neutrophil and monocyte (two types of immune cells), in the blood vessels near the infected tissue. Resting phagocytes are activated and migrate towards the site of pathogens or damaged tissue that has recognizable proteins on surface similar to proteins of immune cells. Once activated phagocytes reach the infection site, they engulf and consume the pathogens. Meanwhile, these activated phagocyte cells release proinflammatory cytokines such as Tumor Necrosis Factorα (TNFα), Interleukin1 (IL1), Interleukin6 (IL6), Interleukin8 (IL8), and HighMobility Group Protein B1 (HMGB1) that activate and recruit additional resting phagocytes from circulation to the infection site. All activated phagocytes eliminate pathogens and secrete substances that accelerate the killing of healthy cells and induce inflammation in the initial process of inflammatory response. In the later stage of AIR progression, several types of antiinflammatory mediators, such as Interleukin10 (IL10), are released by activated phagocytes (primarily monocytederivedmacrophage). These antiinflammatory cytokines inhibit the production of proinflammatory cytokines, consequently inhibiting further recruitment of resting phagocytes. We translated biological processes of AIR to a logical chart, as shown in Figure 1. An explicit description for each biological process is presented in the following six subsections.
2.2. Step 1: Kupffer Cell Local Response Model
Macrophages, one of the innate host’s first lines of defense against bacterial pathogens, are antimicrobicidal cells that often determine outcomes of an infection [21]. Furthermore, hepatic macrophages (also known as Kupffer Cells or resident liver macrophages) constitute 80%–90% of tissue resident macrophages in the body and significantly influence propagation of liver inflammation [34, 35]. Kupffer Cells within the liver trap and eliminate a majority of bacteria that enter the blood stream [22]. During the initial stage of AIR, Kupffer Cells eliminate pathogens, specifically Salmonella, during local immune responses.
We developed a Kupffer Cell local response model, defined as interactions between the pathogen and Kupffer Cell [34], consisting of the following:
In (1), denotes pathogen load, represents a constant growth rate for pathogens, and represents maximum carrying capacity of the pathogen. Parameter represents phagocytosis (killing) rate of Kupffer Cells when Kupffer Cells begin to kill pathogens. Although phagocytosis rate is dependent on time in a slowSshape curve [23], the phagocytosis rate does not change if the phagocytosis rate versus time is assumed to be linear. Therefore, we relaxed the condition that phagocytosis rate is constant in our model and assumed was constant [23]. Equation (2) represents changes of Kupffer Cells over a unit of time, and denotes the amount of Kupffer Cells in the liver that is available for pathogen binding. Parameter term represents a constant proliferation (replenishment) rate for Kupffer Cell population, and represents maximum carrying capacity of Kupffer Cells in the liver of mice. Parameter term represents the unbinding rate of binding Kupffer Cells and represents the killing rate of free Kupffer Cells induced by binding to intruding pathogens.
A standard logistic function was used to model pathogen population growth with limited maximum carrying capacity, identified as the first term in (1) [36]. The second term of (1) models local Kupffer Cell responses or decrease in pathogen population phagocytized by initial tissue resident macrophages (Kupffer Cells). The process of phagocytosis includes two steps: pathogenligand binding to receptors of Kupffer Cells and phagocytosis by Kupffer Cells. We used a Hilltype function and receptorligand kinetics to model the two basic steps [21, 30, 34, 37–39]. First, we defined the rate of pathogen binding to Kupffer Cells as a Hilltype function () in which represents a strong affinity of pathogen binding to Kupffer Cells and is Kupffer Cell concentration that phagocytoses half the pathogens. Second, we modeled pathogen to Kupffer Cell receptors using receptorligand kinetics (), where represents pathogen concentration. We determined pathogen concentration using the number of pathogens divided by maximum carrying capacity of the pathogen (10^{8} cells in the liver of mouse [32]). The final variable to determine pathogen decrease was the phagocytosis rate of pathogens by Kupffer Cells (represented by ) times the portion of pathogens binding to Kupffer Cells ).
We assumed that Kupffer Cells population growth followed a standard logistic growth pattern with a constant proliferation (replenishment) rate, denoted as , and a maximum carrying limit, K_{∞}, represented by the first term in (2). Because pathogen binding did not preclude phagocytosis of additional pathogens after completion of phagocytosis, we used receptorligand kinetics to model the release of Kupffer Cells from the bindingcomplex, represented by the second term in (2); represents the rate by which Salmonella are phagocytosed by free Kupffer Cells and made available for additional interactions with Salmonella. The decreasing number of free Kupffer Cells was due to free Kupffer Cells binding to pathogen, described by the third term , and the natural decay of free Kupffer Cells represented by the fourth term in (2). Free Kupffer Cells become binding Kupffer Cells once they bind to pathogen, as described by the first term in (3). The second term in (3) measures decreasing (releasing) portion of binding Kupffer Cells. The definition of parameters and corresponding experimental data for each system parameter in Kupffer Cell local response model are summarized in Table 1 (refer to Appendix).

2.3. Step 2: Neutrophil Immune Response Model
Simulated results (data not shown) from our Kupffer Cell local response model indicated that Kupffer Cells may not sufficiently eliminate infection, especially when the local infection is overwhelming. Furthermore, pieces of evidence in biological studies have shown that recruitment of neutrophils (one type of immune cells) from circulation to the infection site significantly contributes to AIR progression because neutrophil is able to kill pathogens. Neutrophils accumulation is induced by a proinflammatory cytokine called “TNFα” that is released by Kupffer Cells or activated neutrophils in the tissue. Release of cytokines follows trafficking machinery, and cytokines are released via proteinprotein interactions initiated by ligand binding to receptors [61, 62]. The mechanism of cytokine release is depicted in Figure 2.
We modeled a proteinprotein interaction as MichaelisMenten kinetics [63] and derived our neutrophil immune response model as follows:
Equation (4) was further derived from (1) in the Kupffer local immune response by incorporating phagocytic effects of neutrophils, represented by term . Equations (5) and (6) are cited from (2) and (3).
Equation (7) represents changes of proinflammatory cytokines (denoted by T), such as TNFα, released by binding tissue resident Kupffer Cells () and binding activated neutrophils (). Because TNFα was released after pathogens bound to receptors of tissue resident Kupffer Cells or activated neutrophils, we modeled the process of TNFα release as a combination of MichaelisMenten kinetics and receptorligand kinetics [64]. In (7), the release of TNFα from Kupffer Cells was initiated by receptorligand kinetics, followed by enzymatic kinetics (MichaelisMenten), represented by the term , where represents the maximum production rate of TNFα by binding Kupffer Cells. The release of TNFα is a combined effect of receptorligand kinetics and enzymatic kinetics; therefore, we incorporated both terms in the model to represent combined effects of TNFα releasing processes. Similarly, we used receptorligand kinetics and MichaelisMenten kinetics to model the release of TNFα from binding activated neutrophils in the second term of (7). The third term in (7), , measures degradation of TNFα, with representing the degradation rate of TNFα per hour.
In (8), the first term is a standard logistic function to measure increase in the number of resting neutrophils per time unit (hour), represented by the influx of neutrophils into blood vessels per hour. The second term indicates that the decrease in number of resting neutrophils per time unit is due to the neutrophils activation process promoted by pathogen and proinflammatory cytokine TNFα, where denotes concentration of TNFα and denotes concentration of pathogens [24, 65, 66]. The third term in (8), , represents the natural decay of resting neutrophils; is defined as the apoptotic rate of resting neutrophils per time unit in hours. In (9), the first term exactly equals the second term in (8) because the increased population of activated neutrophils directly resulted from activation of the population of resting neutrophils. The second term of (9) used mass action kinetics () to model the release of activated neutrophils from the bindingcomplex and make activated neutrophils available for additional interaction with pathogens, where represents the bindingcomplex and represents the rate of activated neutrophils released from the bindingcomplex. Similar to the third term in (8), the third term of (9) models natural apoptosis of activated neutrophils. Equation (10) is similar to the derivation of (3) in the Kupffer local response model. We used a hyperbolic tangent function in (11) to represent a slowsaturation influx rate of neutrophils into hepatic parenchyma, thereby representing the rate of activated resting neutrophils. The definition and corresponding experimental data for newly added system parameters in the neutrophil immune response model are summarized in Table 2 (refer to Appendix).

2.4. Step 3: Damaged Tissue Model
Complexity in AIR progression is due to multiple effects induced by inflammatory cells. Recruitment of neutrophils helps clear local pathogen levels; however, those inflammatory cells are harmful because they release toxic molecules such as reactive oxygen species (ROS), which can damage host tissue [24, 66]. Recent experimental results have shown that neutrophils’ integrins adhere to ICAM1 receptors of hepatocytes and accelerate the killing process of distressed hepatocytes [67].
We assumed the binding process of neutrophils to hepatocytes (healthy liver cells) also followed receptorligand kinetics; therefore, we derived the following damaged tissue model:
In (12), D denotes the number of apoptotic hepatocytes or dead hepatocytes, represents the rate of apoptotic hepatocytes killed by activated neutrophils, and represents the recovery rate of apoptotic hepatocytes. Receptorligand kinetics represents the amount of apoptotic hepatocytes that bind to activated neutrophils, with binding rate modeled as a Hilltype function . Activated neutrophils have recently been found to kill apoptotic hepatocytes [67]. After neutrophils adhere to apoptotic hepatocytes, neutrophils release harmful chemical substances such as reactive oxygen species and proteases that accelerate death of apoptotic hepatocytes [67, 68]. When multiplying by , the entire first term in (12) represents the number of apoptotic hepatocytes killed by activated neutrophils per hour, which is the total number of dead hepatocytes per hour. The maximum number of apoptotic or dead hepatocytes does not exceed the total number of hepatocytes in the liver (represented by ). In addition, represents the recovery rate of apoptotic hepatocytes, and the second term in (12) is defined as the amount of recovering apoptotic hepatocytes. The definition of parameters and corresponding experimental data for newly added system parameters in damaged tissue model are summarized in Table 3 (refer to Appendix).

2.5. Step 4: Monocyte Immune Response Model
Recent biological experiments from the literature [69, 70] have shown that monocyte, recruited by the presence of HMGB1, significantly impacts liver inflammation and liver fibrosis. Upon liver injury, inflammatory Ly6cC (Gr1C) monocyte subset, as precursors of tissue macrophages in blood vessels near the infected site, is attracted and recruited to the injured liver via CCR2dependent bone marrow egress. The chemokine receptor CCR2 and its ligand MCP1/CCL2 promote monocyte subset infiltration upon liver injury and further promote the progression of liver fibrosis [26, 67]. Because evidence has shown that Tumor Necrosis Factorα (TNFα) induces a marked increase in CCL2/MCP1 production in dose and timedependent manners [71], we assumed the influx of monocytes from blood vessels to liver is induced by effects of HMGB1 and TNFα. Therefore, we modeled the influx of monocytes similarly to kinetics of neutrophils influx. According to existing literature, HMGB1 is released by necrotic cells and activated monocytes [22, 71, 72]. Therefore, we modeled the release of HMGB1 using receptorligand kinetics and enzymatic kinetics, similar to the release of TNFα, by incorporating effects of necrotic cells and activated monocytes:
In (13), we incorporate the effect of phagocytosis by monocytes into (4) because monocytes phagocytose pathogen by a CD14dependent mechanism [73]. We recalled the Hilltype function equation () to represent receptorligand binding kinetics between pathogens and activated monocytes. Because binding activated neutrophils are engulfed by infiltrating monocytes [27], we used to calibrate the killing process of binding activated neutrophils by activated monocytes, thereby modifying (10) to (14). Equations (15), (16), and (17) describe activation and migration of resting monocytes from blood vessels to infected tissue. In (15), (16), and (17), , , and represent resting monocytes, free activated monocytes, and binding activated monocytes, respectively. Principles used to build those three equations are similar to the principle used to build (8), (9), and (10) for the neutrophil immune response model. Equation (18) calibrates the release of HMGB1 per hour by activated monocytes and apoptotic hepatocytes. The process of releasing HMGB1 is similar to the process of releasing TNFα. The definition of parameters and corresponding experimental data for newly added system parameters in the monocyte immune response model are summarized in Table 4 (refer to Appendix).

2.6. Step 5: SDMM of Innate Immunity
As one type of antiinflammatory cytokines, IL10 was found to prevent subsequent tissue damage by inhibiting activation of phagocytes, including neutrophils and monocytes [74]. This antiinflammatory mediator, produced by macrophages, dendritic cells (DC), B cells, and various subsets of CD4+ and CD8+ T cells [75], follows the same mechanism as proinflammatory (TNFα and HMGB1) release. Because our main focus in this paper was to model innate immune responses, we ignored the release of IL10 by B cells and T cells during adaptive immune responses; therefore, we modeled the release of IL10 similarly to proinflammatory cytokine release:
In (20), represents the number of antiinflammatory cytokine (IL10) during AIR, and represents the release rate of antiinflammatory cytokine (IL10) by activated monocytes, derived from enzymatic kinetics. The first term in (20) calibrates the increase in the number of antiinflammatory cytokines every hour and the second term calibrates the decrease in the number of antiinflammatory cytokines every hour due to natural degradation. Corresponding parameters and their values are defined in Table 5 (refer to Appendix). After incorporating , the inhibition function of IL10, we derived a comprehensive mathematical model for innate immunity of AIR as follows. represents the dissociation rate of IL10 with initial estimated value equivalent to 0.02. Consider

In this 14equation SDMM, variables , , , , , , , , , , , , , and represent levels of pathogen, free Kupffer Cell, bound Kupffer Cell, TNFα, resting neutrophil, free activated neutrophil, bound activated neutrophil, rate of resting neutrophil activated under infection, damaged tissue, resting monocyte, free activated monocytes, bound activated monocytes, HMGB1, and IL10, respectively. These variables are identified and selected as essential indicators in AIR. All system parameters ( and so on), which reflect the strength of the host’s immune system, are adjustable during model simulation. Detailed description of system parameters is presented in the Appendix.
2.7. Step 6: SDMM Incorporated with Adaptive Immunity
Innate immunity plays a significant role in regulating pathogen clearance through multiple types of cell interactions, providing the first line of defense during early stages of inflammation. Compared to innate immunity, adaptive immunity is typically recognized as a late stage of immune response to infection activated by antigen presenting cells (APCs) [76]. The nature of adaptive immune response is more complicated than innate immune responses and involves numerous interactions among cells and cytokines. To simplify adaptive immunity, we selected four representative cells, including CD4+ T cells, CD8+ T cells, B cells, and antibodies, to simulate a series of immune responses during pathogenic inflammation. The 18equation SDMM incorporated with adaptive immunity is presented as follows:
Equation (48) describes the recruiting process of CD4+ T cells during adaptive immunity. The first term in (48) is a standard logistic function to describe the natural migration process of CD4+ T cells to the site of infection, and is a constant parameter to define the recruitment rate of CD4+ T cells from lymph node to the site of infection under undefined mechanisms in our SDMM. Activated monocytes that are phagocytizing pathogens were recognized as one type of APCs; APCs display major histocompatibility complex class II (MHCII) peptide on the surface available for binding to T cell antigenspecific receptor (TCR) [77]. APCs also activate the TCR on CD4+ T cells and enhance CD4+ T cell migration to the site of infection through a TCRMCHII receptorligand response [76], represented by the second term, . Similar to the receptorligand response we modeled in innate immunity, we used a Hilltype function to model the binding rate of activated monocytes to CD4+ T cells. Receptorligand kinetics represent the amount of CD4+ T cells activated by activated monocytes. Our model assumes that T cells become activated under TCRMCHII receptorligand response; however, we recognize that the activation process of T cells is much more complicated than we modeled because T cell activation requires at least two signals in order to become fully activated [77–80]. CD4+ T cells that undergo apoptosis are phagocytized by activated monocytes [81], represented by the third term in (48). We assume that free activated monocytes and binding activated monocytes phagocytize binding CD4+ T cells, represented by a receptorligand response , with the binding rate equal to and the phagocytosis rate equal to . The fourth term, , in (48) describes a natural apoptosis process of CD4+ T cell during migration and activation processes.
Similar to (48), (49) describes the recruitment process of CD8+ T cells during adaptive immunity. The activation process of CD8+ T cells through a major histocompatibility complex class I peptide (MHCI) TCR mechanism follows similar receptorligand kinetics of CD4+ T cells, represented by the second term, , in (49). The activation process of CD4+ T cells and CD8+ T cells is depicted in Figure 3.
CD4+ T and CD8+ T cells mediate the host response to sepsis in various ways. Experimental studies have shown that 1 effector cells proliferated by CD4+ T cells can improve the phagocytosis rate of Kupffer Cells, activated neutrophils, and activated monocytes through a receptorligand response [82]. To simplify our SDMM, we used CD4+ T cell population to substitute for effector cell population, and we measured a decrease in the amount of pathogens via CD4+ T celldependent interactions using receptorligand kinetics, represented by the sixth term in (34). CD8+ T cells are cytotoxic cells because their primary function is to kill infected target cells [82, 83]. Therefore, we incorporated receptorligand kinetics into the third term in (36), the fourth term in (40), and the third term in (45) to measure the decrease in binding Kupffer Cells, binding activated neutrophils, and binding activated monocytes. In SDMM, we used the population of binding Kupffer Cells, binding activated neutrophils, and binding activated monocytes to represent the population of infected cells under the assumption that binding cells bind to pathogens. Therefore, the population of binding cells was also used to represent the population of APCs in our SDMM.
Macrophage activation is related to IFNgamma released by T cells [84–86]. Because we did not calibrate IFNgamma in our SDMM, we calculated the monocyte activation process using CD4+ T cell and CD8+ T cell populations instead of interferongamma (IFNgamma) population for simplicity. Under this assumption, we revised the second term in (43) and the first term in (44) to . The newly revised term, , incorporates the CD4+ T cell and CD8+ T cell populations to reflect the role of CD4+ T cells and CD8+ T cells in the resting monocyte activation process.
or effector cells activate B cells to release antibodies [76]. Equation (50) describes the activation process of B cells by the CD4+ T cell population under the assumption that the CD4+ T cell population can represent and effector cell populations due to model simplification. The first term in (50) measures the migration process of B cells from lymph nodes to the site of infection, which is derived from a standard logistic function. Derivation of the second term, , in (50) is similar to derivation of the second terms in (48) and (49), following a receptorligand kinetics. Decrease in B cell population was induced by natural apoptosis, represented by the third term, , in (50). Plasma cells secrete antibodies [76], but we did not incorporate this specific mechanism into our SDMM. Instead, we modeled that antibodies were released by B cells. In (51), the release of antibodies from B cells is represented by the first term, , following receptorligand kinetics and enzymatic kinetics (MichaelisMenten), similar to TNFα, HMGB1, and IL10 release process described in innate immunity. The second term, , in (51) describes the natural catabolism of antibodies. When antibodies are released from plasmas cells, cells define the isotype of the antibody [76]; we did not model specific isotype of antibodies in our model. Antibodies can opsonize pathogen and contribute to further pathogen clearance at the late stage of inflammation [76, 82], as represented by the fifth term, , in (34). The definition and corresponding experimental data for newly added system parameters in SDMM incorporated with adaptive immunity are summarized in Table 6 (refer to the Appendix).

3. Simulated Results
Using SDMM, we identified three distinct dynamic patterns of indicators that represent three states of AIR progression: Healing Process, Persistent Infection, and Organ Dysfunction. Based on our computed results, we concluded that a Healing Process occurs when the level of pathogens, level of phagocytic cells (neutrophils and monocytes), and level of inflammatory cytokines (TNFα, HMGB1, and IL10) oscillate below threshold during infection. We recognized that a Persistent Infection occurs if inflammatory responses are active (damaged tissue oscillates above threshold during infection). We also recognized that Organ Dysfunction occurs if an overwhelming load of bacteria is observed. Computed results are shown in Figure 4.
In order to initially validate our SDMM, model behaviors were compared to results from experimental designs under specific parameter settings. If results did not match, model reconfiguration was implemented by adjusting the relationship between components (indicators) or finetuning parameter values. We compared our simulated results to experimental results [87] and simulated results from a latest version of an AIR progression mathematical model [13]. We observed that our simulated results had better agreement with experimental results compared to simulated results from the previous mathematical model because our simulated results captured a dampened oscillated infection. We recognized that this improvement of simulation accuracy is a result of additional cellular and molecular pathways of AIR progression incorporated into our SDMM compared to previous mathematical models [12, 13]. For example, we simulated the effect of monocytes in our SDMM by incorporating interactions of monocytes with other cells and cytokines. In contrast, previous mathematical models simulated the combined effect of neutrophils and monocytes with the limitation of oversimplification of AIR progression. Our simulated results indicated that time required for peak levels of TNFα, HMGB1, and IL10 is approximately 12 hrs, 18 hrs, and 24 hrs, respectively. These results are consistent with results from clinical trials [28], as shown in Figure 5.
We also explored the impact of pathogen initial load on phagocytic cells, inflammatory cytokines, and damaged tissue at low, medium, and high levels during AIR progression. We found that dynamic patterns of AIR progression were identified as “Healing Process” if the initial number of pathogens was set below 3.2 (result was transformed to a base10 logarithm) in simulation; dynamic patterns of AIR progression were identified as “Persistent Infection” if the initial number of pathogens was set between 3.2 and 5.9 (result was transformed to a base10 logarithm) in simulation; and dynamic patterns of AIR progression were identified as “Organ Dysfunction” if the initial number of pathogens was set above 5.9 (result was transformed to a base10 logarithm) in simulation. During some simulation replications, our findings are inconsistent with pieces of evidence found from experimental studies [88–90] that indicated outcomes of AIR progression are more likely to lead to a healthy state with a lowdose of pathogens, which will be further illustrated in the Discussion.
By incorporating adaptive immunity to SDMM, we generated dynamic patterns of pathogen count, dead hepatocyte count, activated neutrophil count, activated monocyte count, TNFα, HMGB1, IL10, CD4+ T cell, CD+ 8 T cell, B cell, and antibodies using Mathematica (Wolfram Mathematica 9.0). Computed results are shown in Figures 6 and 7.
Based on our computed results, we observed pathogen count converged toward 0 at approximately 14 days (336 hrs) after infection during a Persistent Infection when the effect of adaptive immunity was incorporated into the full model. Compared to Persistent Infection observed in innate immunity (shown in Figure 4(b)), the activated neutrophil count and HMGB1 count converged toward 0 at approximately 25 days (600 hrs) after infection. Convergence in TNFα count occurred at approximately 14 days after infection, earlier than convergence in HMGB1 count in innate immunity. The peak level of activated monocytes increased to 26000, which was 2 times higher than the peak level of activated monocytes observed in innate immunity. No additional dead hepatocytes were observed after 25 days (600 hrs) after infection because cells (activated neutrophils and activated monocytes) and cytokines (TNFα, HMGB1, and IL10) associated with further tissue damage converged toward 0, indicating adaptive immunity positively impacted outcomes of sepsis progression.
By incorporating CD4+ T cells, CD8+ T cells, B cells, and antibodies into innate immunity, we observed that elevated pathogen count during Organ Dysfunction began to drop at approximately 20 days after infection (500 hrs), and the process of pathogen clearance induced by adaptive immunity persisted approximately 5 days after infection. Pathogen count returned to 0 at 25 days after infection (720 hrs). Cells (activated neutrophils and activated monocytes) and cytokines (TNFα, HMGB1, and IL10) associated with innate immunity dropped significantly during simulation, but CD4+ T cells, CD8+ T cells, and B cells persistently elevated after 500 hrs after infection, indicating adaptive immunity’s contribution to pathogen clearance during the late stage of sepsis progression. A mice model infected with a high dose of Escherichia coli [52] showed that the number of CD4+ T cells, CD8+ T cells, and B cells persisted throughout 7 days, thereby conforming to dynamic patterns of CD4+ T cells, CD8+ T cells, and B cells observed in our SDMM.
4. Stability Analysis
In order to study model behaviors under various parameter settings and initial conditions, bifurcation diagrams were used to conduct stability analysis for each subsystem during model construction. The objective of stability analysis was to identify key parameters or key processes in sepsis episodes. Numerical analysis that we used is similar to the previous study [16].
We started with stability analysis by calculating equilibrium points in Kupffer Cell local response model. The equilibrium points were derived by setting equations in Kupffer Cell local response model free of the time (time is denoted by in equations), which imply that
To solve (52), (53), and (54), we firstly added (53) to (54), which eliminate the (53) and (54) to (55):
By solving (52) and (55) together, we could obtain the following feasible equilibrium points.
If ,
The above equilibrium points are valid if the following conditions are satisfied:
From the derived feasible equilibrium points, we obtained two diseasefree equilibrium points given as
We further calculated the associated Jacobian matrix to determine stability of the diseasefree equilibrium points; the Jacobian matrix was given as follows:where .
Replacing the first diseasefree equilibrium point into the Jacobian matrix above (59), we can further derive the following Jacobian matrix:
In order to find the associated eigenvalues with (60), we solved the following equation:
Using Mathematica (Wolfram Mathematica 9.0), we obtained the eigenvalues of (61) as follows:
Thus, we concluded that the first diseasefree equilibrium point is stable if and only if the following conditions are satisfied:
Following a similar procedure above, we replaced the second diseasefree equilibrium point into the Jacobian matrix in (59).
The Jacobian matrix associated with the second diseasefree equilibrium point was revised to
Again, by solving (65),
We obtained the eigenvalues associated with the second diseasefree equilibrium point, and the eigenvalues were expressed as follows:
Thus, the stability of the second diseasefree equilibrium can be achieved if and only if the following conditions are satisfied:
Because (the growth rate of pathogen) was assumed to be always larger than 0, we concluded that the diseasefree equilibrium points for Kupffer Cell local response model are always unstable.
In order to verify our conclusion, we did a numerical study on the second diseasefree equilibrium point . We found the diseasefree equilibrium point changed if pathogen load was changed from 0 to 2 at equilibria (a small perturbation was given); the simulated results of change in the diseasefree equilibrium point are shown in Figure 8.
We also analyzed stability of the pathogen saturation equilibrium point . By numerical analysis, we concluded that the pathogen saturation equilibrium point is stable if the following conditions are satisfied:
When > 0.5, the pathogen saturation equilibrium point became unstable. Simulated results of change in the pathogen saturation equilibrium point are shown in Figure 9.
Stability analysis of equilibrium points in Kupffer Cell local response model indicated that Kupffer Cell local response model is not a stable system. The diseasefree equilibrium point changed when the second infection occurred (P was changed from 0 to 2). However, recruiting more Kupffer Cells positively contributed to the pathogen clearance after a saturated infection (P = ), as shown in Figure 9.
Bifurcation diagrams are graphical tools to visualize dynamic system behavior changes with parameters. In this paper, we used Matcont to generate bifurcation diagrams. Matcont, a Matlab continuation package with a graphic user interface (GUI) for interactive numerical study of parameterized nonlinear ordinary differential equations (ODEs), computes curves of equilibria, limit points, Hopf point, limit cycles, fold, torus, and branch point bifurcation of limit cycles [91].
In bifurcation diagrams, axis represents equilibria of state variable and axis represents value of system parameter that generates equilibria. Therefore, bifurcation diagrams reflect change in equilibria of dynamic systems (change in number of equilibria or change in numerical value of equilibria) in relation to change in numerical value of system parameters. We analyzed stability of dynamic systems by identifying types of bifurcation points in bifurcation diagrams because bifurcation points are defined as points at which stability changes from stable to unstable. Two typical bifurcation points were evident in our bifurcation diagrams: limit point (marked as “LP” in Matcont) and Hopf point (marked as “H” in Matcont). Neutral saddle point was marked as “NS” in the bifurcation diagram, but it is not a bifurcation point for equilibrium because it is identified as a hyperbolic saddle. Figure 10 shows that change in equilibria of state variable pathogen is related to change in system parameters in the neutrophil immune response model.
(a)
(b)
(c)
(d)
LPs in bifurcation diagrams of neutrophil immune response model appeared when two equilibria merged into one equilibrium; the number of equilibria of dynamic systems changed when LPs were detected. LPs are also turning points at which dynamic systems change from stability to instability. In Figure 10(a), stable equilibria of pathogen are observed when system parameter increases from 0 to 4.93. When equals 4.93, LP is identified and unstable equilibria of pathogen are generated as decreases from 4.93 to 0. Therefore, equilibria of pathogen of our neutrophil immune response model are bistable when ranges from 0 to 4.93. Similarly, equilibria of pathogen in Figure 10(b) are bistable when system parameter ranges from 25 to 200. In Figure 10(c), equilibria of pathogen are bistable when ranges from 0 to 0.21.
A Hopf bifurcation, identified in Figure 10(d), is a periodic bifurcation in which a new limit cycle is born from a stationary solution. Hopf point, a turning point for periodic orbits, is detected when system parameter changes. The detected Hopf point in Figure 10(d) begins a limit cycle continuation in which two cycles collide and disappear. Because the first Lyapunov coefficient [92] is positive, an unstable limit cycle exists, bifurcating from this equilibrium. Figures 11(a) and 11(b) show the family of limit cycles bifurcating from detected Hopf point in Figure 10(d). The family of limit cycles is represented using limit cycle planes, such as TNFapathogen plane and pathogen plane. Figure 11(c) shows a limit cycle sphere represented by a TNFa, , and pathogen plane. Figure 11(d) indicates that two limit cycles occur when equals 5495.64 or 6265.00.
(a)
(b)
(c)
(d)
In Figure 11(c), the first family of limit cycle (small red cycle in the center of the sphere) spirals outward as system parameter decreases, and the second family of limit cycle appears when decreases to 5495.64 (a red cycle line appears). As increases from 5495.64, the second family of limit cycle spirals outward again. When increases to 6265.00, an unstable equilibrium is detected, as depicted in Figure 12(a). If value of is between 5495.64 and 6265.00, equilibria of the neutrophil immune response model are stable and converged, as shown in Figure 12(b). This finding infers either a high release rate of TNFα ( is above 6265.00) or a low release rate of TNFα ( is below 5495.64), thereby inducing generation of unstable equilibria in the neutrophil immune response model. From a biological response perspective, high release rate of TNFα indicates overproduction of proinflammatory cytokines related to overwhelming proinflammation; low release rate of TNFα leads to failure to recruit a sufficient amount of neutrophils related to infection clearance. Based on our stability analysis, we found that the release rate of TNFα can positively or negatively influence outcomes of AIR progression, thereby conforming to experimental perturbation findings regarding effectiveness of antiTNFα therapies [93–95].
Continued stability analysis on the monocyte immune response model indicated that change in system parameters , , and induces bistability of the monocyte immune response model. We observed that the monocyte immune response model was bistable if at least one of the following three conditions was met: was between 0 and 0.32, was between 0 and 0.28, or was between 0 and 0.21. Specifically, we observed that (maximum release rate of TNFα by activated neutrophil) and m_{t2} (number of activated neutrophils at which the reaction rate is half of the maximum production rate) are essential for oscillated monocyte immune response model. Similar to the neutrophil immune response model, limit cycles bifurcate from Hopf point. Therefore, we conclude that the oscillated infection is dependent on the amount of released TNFa and recruited neutrophils in AIR progression. However, released monocytes and associated cytokines such as HMGB1 do not contribute to oscillation in AIR progression.
Building upon the monocyte immune response model, we incorporated the effect of antiinflammatory cytokine (IL10) into the full model. We observed that Hopf point was detected when increased to 128000 because antiinflammatory cytokine inhibited activation of phagocytic cells (neutrophils and monocytes). This trend indicates that infection oscillation (harmful outcomes) requires additional proinflammation activated by neutrophils in the full model, compared to monocyte immune response model without including the effect of antiinflammatory cytokine. Therefore, our simulated results demonstrated that AIR progression is more likely to end with Healing Process if the effect of antiinflammatory cytokine is incorporated.
Strengthened (increased and m_{t2}) proinflammatory immune responses could also induce stable or unstable equilibria, leading to a dampened oscillated infection or diverged infection, similar to our observations in Figure 12. However, we observed that if high effect of antiinflammatory cytokine was incorporated (dissociation rate equal to a base10 logarithm 8) at the beginning of infection, AIR progression resulted in an unstable overwhelming pathogen load at equilibria (refer to Figure 12(a)). However, a stable dampened oscillated pathogen load at equilibria (refer to Figure 12(b)) was observed if medium effect of antiinflammatory cytokine (dissociation rate equal to a base10 logarithm 5) was incorporated. These observations confirmed that effects of antiinflammatory cytokine can be positive or negative to AIR progression depending on levels of antiinflammatory cytokine.
We conducted bifurcation analysis for the model incorporated with adaptive immunity, similar to bifurcation analysis we conducted in the neutrophil subsystem, monocytes subsystem, and full model. We selected four bifurcation diagrams as shown in Figure 13.
(a)
(b)
(c)
(d)
As shown in Figure 13(a), two Hopf bifurcations were detected at = 2.8 and = 4.1. Similarly, Hopf bifurcations were also detected in Figures 13(b) and 13(c) when = 17, = 38, or = 0.047. Compared to innate immunity, incorporation of adaptive immunity induced a further stabilized limit cycles, bifurcation from the equilibrium. Our stability analysis shown in Figure 13(d) illustrates that the Hopf bifurcation moves to lower value compared to Hopf bifurcation detected in innate immunity. The change in bifurcations indicated the contribution of adaptive immunity to sepsis progression.
In Figures 14(a), 14(b), and 14(c), the first family of limit cycle (small red cycle in the center of the sphere, marked as LPC) spiral outward as system parameter decreases, and the second family of limit cycle appears when decreases to 2.4 (a red cycle line appears). As increases from 2.4, the second family of limit cycle spirals outward again. A period doubling is detected when increases to 3, marked as PD in Figure 14. Because the first Lyapunov coefficient is negative, limit cycle bifurcations from the equilibrium are stable compared to unstable limit cycles detected in the neutrophil subsystem.
(a)
(b)
(c)
(d)
5. Discussion
Experimental results in literature have suggested that antiinflammatory mediator inhibits activation of phagocytes and reduces the ability of activatedphagocytes to attack pathogen [96], consequently related to mortality and severity of infection in sepsis [97, 98]. However, other experimental studies have shown that antiinflammatory cytokine downregulates production of secreted cytokines by inhibiting various behaviors of activated immune cells, thereby reducing the risk of tissue damage [28, 99, 100]. Our computed results from SDMM suggested that the effect of antiinflammatory cytokines could be a “doubleedged sword” for AIR because antiinflammatory cytokine would either decrease mortality associated with tissue damage or increase mortality associated with high load of bacteria. With a low effect of antiinflammatory cytokine (dissociation rate equal to a base10 logarithm 2), our computed results showed that antiinflammatory cytokine fails to inhibit the release of activated immune cells (activated neutrophils and activated monocytes) and subsequent cytokine production. Levels of damaged tissue significantly accumulated during the first 500 hours (approximately 20 days) of simulation. With the high effect of antiinflammatory cytokine (dissociation rate equal to a base10 logarithm 8), our simulated results and stability analysis demonstrated that sepsis progression leads to increased chance of death caused by overwhelming pathogen load at the end of simulation.
To further investigate effects of antiinflammatory cytokines, we simulated a medium effect of antiinflammatory cytokine (dissociation rate equal to a base10 logarithm 5) and compared simulated results to high effect of antiinflammatory cytokine and low effect of antiinflammatory cytokine. Our computed results showed that pathogen load decreases during the first 100 hours of infection in combination with the total amount of dead hepatocytes. Furthermore, we observed that production of activated neutrophils and activated monocytes declined to baseline near 0 at the end of simulation, indicating a positive trend of sepsis progression to a healthy pattern. Therefore, we conclude that the level of antiinflammatory cytokines significantly impacts direction of sepsis progression. We also conclude that levels of antiinflammatory cytokine and time of intervention of antiinflammatory cytokines determine outcomes of AIR under specific system configuration. Based on simulated results from our SDMM, we inferred that the survival rate of the host (chance of ending with a Healing Process) could be improved if a medium level of IL10 injection was set between 3 hrs and 6 hrs after infection.
We assert that care must be taken when applying simulated results to clinics before the implementer fully understands the underlying setting of the simulation. Because it was difficult to simultaneously incorporate every intermediate biological process of inflammatory response, reasonable assumptions must be made when building a mathematical model. In our SDMM, we did not model Salmonella replicating within neutrophils. However, experimental study [58] asserted that neutrophils and macrophages were the primary sites for Salmonella proliferation in a mouse. Therefore, Salmonella replication could be considered in the future model if additional literature supported this fact. Various T cell types were reported to be able to express IL10 under various conditions [101]. Therefore, IL10 production estimation is difficult because IL10 levels produced by T cells were various due to stimuli type or the strength of stimuli. In our model, we did not differentiate helper T cells from specific types that are identified in biological process. Plasma cells secrete antibodies [76], but we did not incorporate this specific mechanism in our SDMM. Instead, we modeled that B cells released antibodies. When antibodies are released from plasmas cells, cells define the isotype of the antibody [76]; however, we also did not model specific isotype of antibody in our model. Furthermore, we ignored the fact that antibody opsonization induces stimulation of cytokine release when they are phagocytized by inflammatory cells. We ignored the fact that antibody opsonization induces stimulation of the release of various cytokines from neutrophils and macrophages [76]. Also, we ignored effects of other proinflammatory cytokines such as IL1, IL12, and IL8 in our SDMM. Biological immune responses responding to infection are recognized as a series of complex processes including intracellular transductions (transfer of DNA) and intercellular pathways between cells. These biological processes will be developed with evolved understanding and continued investigation of cellular and molecular mechanisms [34], which could be further research interests in the field. In our SDMM, we used numerical count, or the number of indicators in the simulation, as the estimate of cell or cytokine number in AIR. In practice, physicians must translate data to measurable units of indicators similar to how we translated clinical data to simulation data. Furthermore, our conclusion regarding IL10 was drawn based on specific simulation settings including setting system parameters and initial loads of indicators. Initial system setting must be fully understood before considering application of IL10 level for preclinic experiments.
Based on our simulated results regarding antiinflammatory cytokine, we propose a hypothesis testing: If medium levels of antiinflammatory antibody were injected into the host with sepsis between 3 hrs and 6 hrs would survival rates of the host improve under hyperinflammation? The purpose of this hypothesis testing is to detect effective zones of the antiinflammatory antibody related to Healing Process of AIR in order to help develop therapeutic agents in preclinical trials.
According to our simulation study, we found that initial levels of pathogen significantly impact dynamic p