The development of a fully-integrated immune response model (FIRM) simulator of the immune response through integration of multiple subset models
© Palsson et al.; licensee BioMed Central Ltd. 2013
Received: 15 November 2012
Accepted: 21 August 2013
Published: 28 September 2013
The complexity and multiscale nature of the mammalian immune response provides an excellent test bed for the potential of mathematical modeling and simulation to facilitate mechanistic understanding. Historically, mathematical models of the immune response focused on subsets of the immune system and/or specific aspects of the response. Mathematical models have been developed for the humoral side of the immune response, or for the cellular side, or for cytokine kinetics, but rarely have they been proposed to encompass the overall system complexity. We propose here a framework for integration of subset models, based on a system biology approach.
A dynamic simulator, the Fully-integrated Immune Response Model (FIRM), was built in a stepwise fashion by integrating published subset models and adding novel features. The approach used to build the model includes the formulation of the network of interacting species and the subsequent introduction of rate laws to describe each biological process. The resulting model represents a multi-organ structure, comprised of the target organ where the immune response takes place, circulating blood, lymphoid T, and lymphoid B tissue. The cell types accounted for include macrophages, a few T-cell lineages (cytotoxic, regulatory, helper 1, and helper 2), and B-cell activation to plasma cells. Four different cytokines were accounted for: IFN-γ, IL-4, IL-10 and IL-12. In addition, generic inflammatory signals are used to represent the kinetics of IL-1, IL-2, and TGF-β. Cell recruitment, differentiation, replication, apoptosis and migration are described as appropriate for the different cell types. The model is a hybrid structure containing information from several mammalian species. The structure of the network was built to be physiologically and biochemically consistent. Rate laws for all the cellular fate processes, growth factor production rates and half-lives, together with antibody production rates and half-lives, are provided. The results demonstrate how this framework can be used to integrate mathematical models of the immune response from several published sources and describe qualitative predictions of global immune system response arising from the integrated, hybrid model. In addition, we show how the model can be expanded to include novel biological findings. Case studies were carried out to simulate TB infection, tumor rejection, response to a blood borne pathogen and the consequences of accounting for regulatory T-cells.
The final result of this work is a postulated and increasingly comprehensive representation of the mammalian immune system, based on physiological knowledge and susceptible to further experimental testing and validation. We believe that the integrated nature of FIRM has the potential to simulate a range of responses under a variety of conditions, from modeling of immune responses after tuberculosis (TB) infection to tumor formation in tissues. FIRM also has the flexibility to be expanded to include both complex and novel immunological response features as our knowledge of the immune system advances.
KeywordsBiological networks Immune response Mathematical modeling Ordinary differential equation systems Systems biology
Mathematical models are a natural approach to improve our understanding of complex biological systems, and ultimately enabling us to predict their behavior and control them . In particular, the intricacies and nonlinear nature of the mammalian immune response have attracted considerable attention over the years , in no small part due to the role of the immune response in a variety of relevant human conditions. In addition, the existence of a mathematical model allows one to explore the known differences in immunity development between human and non-human species  by altering or excluding specific pathways, as dictated by experimental findings. The basic components of interest would in principle include: cellular or cytotoxic responses (i.e., the development of cells from the T-lineage that attack antigens directly), humoral responses (i.e. the endogenous production of antibodies from cells of the B-lineage that bind to the antigen receptor and hasten its removal) and signaling features (mainly, but not exclusively, through the cytokine network) .
In general, the development of quantitative models is often based on the selection of features of interest and their description in mathematical form, followed by their functional integration into a model that can be interrogated and/or used to predict features of interest. Such features can then be compared to experimental data. Similar procedures are followed for immune response models, but due to the system’s complexity, modeling and simulation efforts have focused on specific subsets of the system, such as the cellular responses [5–7], humoral responses [8, 9] and/or cytokine networks [10, 11], while sometimes excluding or simplifying other components from the model. As a general consideration, the development of comprehensive models is difficult and has to contend with the integrated network nature of the system, where the addition of one novel component necessarily requires defining the interactions of the new item with the remainder of the network.
Several modeling formalisms have been used in developing models of the immune system. Historically, these have been mostly categorized as differential equation models or agent-based models. Agent-based models or cellular automata models of the immune response have attracted great interest  since very early studies  and have been refined and proposed over the years in a manner that is responsive to new knowledge [14, 15]. Their greatest strength is their flexibility and relative ease of use , which makes them suitable to model very complex systems without having to mechanistically specify the known component interactions. Instead, the system is defined in terms of “computer agents”, which are sets of rules by which individual actors (i.e. populations of cells, or even individual cells) are created, interact and are destroyed. The modeling effort then focuses on monitoring the interactions among agents, which gives rise to complex, sometimes emergent behaviors that, depending on the rule base, can provide a striking similarity to the temporal evolution of the system being represented. Such models can then be used to develop answers to complex problems, including therapy optimization [17, 18]. As others have pointed out , despite their power, challenges remain with agent-based models, including the availability of widely accepted software and of model checking and goodness of fit strategies that resemble those commonly used for differential equation models. Differential equation models have provided tremendous insight in the dynamics of complex immunological networks [19, 20] and are still widely used, relatively easier to communicate and more readily shared than agent-based models. Some of these models can achieve remarkable degrees of complexity and realism . In addition, differential equations form the backbone of translational pharmacokinetic-pharmacodynamic (PK-PD) models [22, 23], the class of models that describe how drug dose influences response through quantitatively linking the drug dose to exposure (pharmacokinetic [24, 25]) and the exposure to response (pharmacodynamic ) in a living system. These historically are the models of choice in drug research and development.
The integration of multi-scale, realistic models of physiology with pharmacotherapeutic models is a desirable goal that would allow more mechanistic, predictive and overall useful models for drug research and development , in addition to enhancing collaborative efforts between biology and modeling. This effort is receiving renewed attention through the area of “systems pharmacology” , as explored in two successful interdisciplinary workshops hosted by the National Institutes of Health in 2008 (http://meetings.nigms.nih.gov/?ID=3447) and 2010 (http://meetings.nigms.nih.gov/?ID=8316). Issues related to efficient model sharing and model construction are also the purview of the Interagency Modeling and Analysis Group (IMAG) (http://www.imagwiki.nibib.nih.gov). Building these models accurately and efficiently represents a significant challenge. This has prompted the development of sophisticated software to facilitate integration of separate submodels  and parallel computation .
The recent availability of computational environments to functionally connect submodels without having to write ad hoc computer code is well complemented by the development of approaches to supervised “monolithic”  model integration. In our case, a differential equation framework was chosen for the development of an integrated immune response simulator, coupled with a useful framework, found in systems biology, for integrating multiple subset models into a coherent whole. In this framework, connectivity matrices are built to describe the global network structure, followed by introduction of rate laws to seamlessly integrate multiple biological processes [31, 32].
The Fully-integrated Immune Response Modeling (FIRM) simulator is a differential equation based integration of multiple existing models of the immune system [5, 8, 10]. It accounts for both the humoral and cellular immune response systems and attempts to parsimoniously represent the spatially distributed nature of the system. The goal of this integrated model is to specify antigen exposure over time and calculate predicted antibody levels and cell concentrations following biological perturbations such as immunization or infection. To develop FIRM, we used a pharmacokinetic / pharmacodynamic modeling approach to combine previously published individual models of humoral and cellular response with antigen exposure. FIRM includes both the antigen-specific antibodies and cell populations, and accounts for cytokines and adjuvant components as needed. It is a hybrid construct incorporating structures and parameter values from published models in multiple mammalian species.
This report outlines the stepwise integration of networks describing the cellular dynamics for both T and B-cell responses to bacterial infections and to tumor growth in a target organ. In addition, to illustrate the incorporation of novel mechanisms, we propose and integrate within the framework a new hypothesized model of regulatory T-cell kinetics accounting for immunoevasion.
Model formulation and content
The process of building an integrated simulator starts with the definition of the underlying physiological structure. This preliminarily defines the existing interrelations among all the variables of interest as a “superset” of cellular and molecular populations and reactions. Second, all the cellular and molecular state variables are identified and the interrelationships (transitions) between them determined. The structure of the networks is thereby specified. Third, the mathematical forms of the equations that describe the fluxes are then formulated and their numerical values determined (from literature or existing data). Usually, the first two steps involve the determination and selection of existing relationships that have a physiological basis. As such, they are somewhat easier than the third step, where such relationships need to be made specific and quantitative. The availability of plausible numerical values is a well-known rate limiting step in the definition and assembly of kinetic models, and in the rest of this section we will outline the approach we followed to inform FIRM’s parameterization.
where x is the vector of state variables (concentrations of various cell types and molecules) and v is the vector of fluxes from one state to the next (i.e. transport processes, reaction rates, cellular fate processes, etc., expressed in concentrations per unit time). S is a matrix that describes the structure of the network and its topology. Every column in S represents a flux and every row represents a state variable. The vector k contains the numerical values of the kinetic and physical constants (often, but not necessarily, expressed in units of inverse time). In general, the vector of fluxes v is a function of the state variables and the kinetic and physical constants characterizing transport and reaction processes.
The matrix S, the state variables and the fluxes corresponding to the final structure of the FIRM model are found in supplemental Additional file 1: Tables S1, S2 and S3. The mathematical form of all the flux variables are given in Additional file 1: Table S4, and the numerical values and their literature sources are found in Additional file 1: Table S5 and S6.
Spatial distribution features
To account for known features of the spatial distribution of the immune response components in our simulations in a parsimonious manner, the FIRM model has five separate tissue spaces where the cell and cytokine populations can travel: lung (assumed to have a volume of 1000 mL), blood (4500 mL), lymphoid tissues relevant to the cellular (10 mL) and humoral response (150 mL) and the sites of immune recognition (4500 mL). As populations of cells and molecules travel between biological spaces, their concentrations are multiplied by the respective volumes of distribution so as to maintain mass balance. This was particularly important for the population of infected macrophages, which changes dynamically as bacterial infection progresses and turns out to have a time-dependent variable volume of distribution whose features needed to be accounted for in the simulation. Additional file 1: Table S7 in the Supplementary Material contains the various tissue space volumes.
It is worth noting that this is a parsimonious representation since, to accurately represent spatially differentiated behaviors, one would have to define biological spaces for each spatially (and functionally) separate component of each organ or tissue that has a distinct pattern or behavior from other parts of the organ or tissue under consideration. These spaces would be defined so as to have different volumes, rate constants for accessibility, etc., to reflect their heterogeneous physical structure and, given the number of parameters required, would require detailed experimental information at the cellular and molecular level. In this sense, FIRM is a parsimonious model and does not reach this level of granularity, although the model structure is amenable to be extended and incorporate such considerations where warranted, required by the purpose of the modeling exercise and supported by the data.
Simulation platform and integration procedure
Results and discussion
(1) Procedure for integrating multiple dynamic models from different sources
where x is the vector of state variables (cell types and molecule concentrations) and v is the vector of fluxes from one state to the next (i.e., transport processes, reaction rates, cellular fate processes, etc.). S is a matrix that describes the structure of the network and its topology. Every column in S represents a flux and every row represents a state variable. The vector k contains the numerical values of the kinetic and physical constants. The contents of all these mathematical objects are found in Additional file 1: Tables S1 (Stoichiometric matrix), Additional file 1: Table S2 (variable list), Additional file 1: Table S3 (flux list), Additional file 1: Table S4 (rate laws), Additional file 1: Table S5 (kinetic constant values), and Additional file 1: Table S6 (miscellaneous parameter values).
The reconstructed network is shown in Figure 1A. It includes all the major T and B-cell types, the pathogens (antigens), as well as tumor and its debris. The regulatory effects of the major inflammatory growth factors and cytokines are shown in Figure 1B. These graphs highlight the boundaries of FIRM. The processes that we considered for inclusion have been described multiple times (see e.g. figure in  for a description) and relate to humoral and adaptive immunity. Briefly, the humoral side of the system describes the activity of B-cells that, when in contact with an antigen/pathogen, secrete antibodies which are essentially a secreted version of their receptor compatible with the antigen). The antibodies bind to the antigens and neutralize them. The link between the humoral and the cellular components of the system is provided by helper T-cells, which activate B-cells. This function is currently not included in the model since it was not a part of the constituent submodels. The cellular side starts with naïve T-cells recognizing antigen epitopes through antigen presentation via dendritic cells and macrophages and subsequently developing a T-cell driven response to the antigen. This process is particularly important when the pathogen is a cell, such as in cancer and bacterial immunity, which are both described in FIRM. T-cells can also differentiate to regulatory T-cells, which essentially mute features of the immune response.
The FIRM simulator was built in a step-by-step fashion, summarized in Figure 2, from both constituent models that have appeared in the literature [5, 8, 10] and novel mechanisms. A summary of those steps is below.
Overall summary of integration issues
Integration issues resolved with the mapping of the MK model [Ref.] onto the FIRM network structure
• Added basal state levels of resting macrophages and IDC using MK latency parameter values:
• Removed “HTL” (TH1) from activation of macrophages.
• BM conversion to B cells has the same rate as the death of BM.
- MR = 5 * 108 cells
• Using DB value for MA half-life.
• Antibody produced by BP, BA in blood, and BA in lymphoid B.
- IDC  = 5 * 107 cells
- η13 = 1 day-1
• Expansion of the BL model due to the relaxation of equilibrium assumptions required the creation of variables x41-x52 and fluxes v87-v100.
• Combined bursting (v3) and natural death (v14) of infected macrophages into one flux (v3). The reaction rate will be the summation of the two individual reaction rates.
• Created constant recruitment of TCP and THP in the lymphoid T, in fluxes v50 and v21, respectively, analogous to I1 and I2 from the DB model.
• Created “antigen” variable in blood and Site of recognition.
• Introduced MI (Infected Macrophages) as a separate space with a variable volume:
- I1 → ρ50 , I2 → ρ21
- “Antigen” can be either a tumor debris or a bacteria cell, they become a part of the “antigen” pool once they enter the blood
- VolumeMI = 8*10-9 * x3
• Accounted for TCP presence in the blood, created a separate variable x25.
- Permeates from lung to blood, and from blood to Site
• Fixed bacteria accounting issues:
• Temporarily changed HTL (TH1) in FACTOR to HTLP (THP).
- In the MK model, half of the amount of bacteria released during bursting was required to infect one macrophage. These two amounts have been made independent, but they are currently set to 25 bacteria for infection and 50 bacteria for bursting. These are the same values used in the MK model, but they can be changed easily.
• Using DB value for MR half-life.
• Created receptor sites on select B cell populations.
- Bursting (v3) is based on the ratio of intracellular bacteria to infected macrophages (x5/x3). The bursting will occur at a greater rate as the ratio approaches a set number (the macrophage’s capacity). This capacity is currently set to 50.
- η10 = 0.05 day-1
- x16, x17, and x18 have receptor sites
- Bursting (v3) releases x5/x3 extracellular bacteria into the system, instead of a fixed number.
• Redefined FACTOR with HTL (TH1).
- Receptor sites have 2 states: antigen-bound and free
- T-cell induced apoptosis (v4) releases x5/x3 extracellular bacteria into the system.
• Added TH1 (HTL) proliferation from the DB model as a negative term to the death flux v28.
• Expanded antigen-B cell interaction to include receptor sites and binding events.
- It is important to note that x5/x3 is a time-dependent ratio.
• Modeled differentiation of naïve T cells to TCP (v50) to mirror the action of v21 from the MK model.
- All antigen-receptor binding events occur at the same rate
• Combined naïve T cell death and recirculation from the MK model into one clearance flux (v20). The reaction rate will be the summation of the two individual reaction rates.
• Modified MAPC from the DB model. MAPC and its corresponding fluxes (v57, v58) will remain inactive and undefined. The functionality of MAPC described in the DB model, using the variable INFLAM, will be merged with the dendritic cells:
- The receptor-antigen binding event is a reversible reaction
• Added basal state levels of resting macrophages and IDC using MK latency parameter values but using the new clearance flux of naïve T cells (v20):
- Added term to recruitment of IDC cells (v15) using INFLAM as a trigger.
- x20 is assumed to have the same receptor state as x18
- T = 98,039 cells
■ (ρ21 + ρ50)/2 * INFLAM
• Expanded antigen-antibody interactions to include dynamic single- and double-bound states.
• Modified rate law v22. The MK formulation allowed for negative proliferation.
- Added term to migration/maturation of IDCs (v17) using INFLAM as a trigger.
- All antigen-antibody binding events occur at the same rate
• Accounted for THP presence in the blood, created a separate variable x13.
■ Used term from MK stimulation, but replaced x4 (PE) with INFLAM
- All bound antibody states are cleared at the same rate
• Used volume ratios to properly account for cell migrations across tissue space borders.
• Cut off an INFLAM feedback loop by globally redefining INFLAM without HTL (TH1) when substituting in FACTOR. Now, the only variable that determines INFLAM is tumor burden. The basic premise of the INFLAM loop is an increase in INFLAM causes dendritic cells to produce more helper T cells, and the creation of these helper T cells caused FACTOR to increase, which in turn caused INFLAM to increase.
- Binding events occur in both the blood (with “antigen”) and the lung (with extracellular bacteria)
• Eliminated flux v30. The migration flux of TH2 to the blood (v31) that was to be added with the B cell response will take its place. v31 will take the death rate of v30 (0.3333 day-1) as its reaction rate. Having two fluxes drain the TH2 population was leaving the TH2 levels in the lung much too low and ineffective.
• Added new fluxes to FIRM structure:
• Defined initial conditions with analytical solutions for B cells and B cell free receptors sites.
• Added basal state levels of IL-12 in the lung, produced by the basal levels of MR.
- v84 → death of THP in the blood
• Permeation of tumor debris to blood is turned off.
IL-12 = 5*108 (q78a/η79)
- v85 → death of TCP in the blood
- Tumor-antibody interaction lacks definition at this time
- v86 → proliferation of TH1 in the lung (removed negative term from v28)
• Defined initial conditions with analytical solutions for:
- MR, IDC, T, THP in the lymphoid T, THP in the blood, TCP in the lymphoid T, TCP in the blood, IL-12 in the lung
Immune response to tumor formation: Cellular response to tumor challenge. The second subset model identified for inclusion in the integrated model (by DeBoer et al., hereafter the DB model) described the inflammation response to the presence of a tumor. Its components were the growth of tumor mass, the activation of macrophages in response to the tumor cells, the proliferation and differentiation of cytotoxic and T1 helper cells, and the killing of tumor cells creating tumor debris. The original formulation of this subset model had all tracked populations in a single biological space; therefore, the cell populations described in the model had to be mapped to their appropriate organs. The scope of this model can be seen in Figure 3, together with its overlap and points of contact with the MK model. Once again, the state variables and fluxes associated with the content of this subset model were added to the model network.
The inclusion of the DB model marks the first integration in the FIRM system of two separately developed and reported kinetic models (MK and DB). The integration of two kinetic models resulted in some complexities, reflecting in turn: the state of biological knowledge revealed by the models, the assumptions made, the structure of the network, and the detail of the quantitative information. The process of integrating models from various sources and built for different mammalian species requires explicit and sometimes implicit biological and structural assumptions. We summarized those choices as “integration issues” and they are reported in detail in the Table 1. The resolution of integration issues is critical to the construction of an integrated model such as FIRM. For example, the DB model included processes and parameters which tended to be descriptive as opposed to mechanistic, reflecting the knowledge of immunology at the time. Therefore, the issues involved with the addition of the DB model included integrating phenomenological parameters into the FIRM state variable and rate law structure (and sometimes modifying them), selecting a value for kinetic constants that appear in both models, and including new network fluxes, mostly cellular, as necessary to integrate the models to ultimately ensure balance of mass.
B-cell response to antigen: Humoral response to antigen challenge. The third subset model that was included in FIRM (by Bell, hereafter the BL model) details the B-cell and antigen response to the presence of an antigen in the system. Included in the BL model is the exposure of naïve B-cells to antigen, the activation of B-cells that migrate to the lymphoid B via the blood, the differentiation of activated B-cells to plasma and memory B-cells, and the production of antibodies by said B-cells. Lastly, the antibodies work to eliminate the antigen from the blood and the target organ. This model’s components are again shown in Figure 3, which also highlights commonalities with the MK and DB models.
It should be mentioned at this point that the constituent models composing FIRM were not necessarily tailored to a particular animal species. The BL model description mentions data being obtained in rabbits, but it also states that the model recapitulates essential features of the mammalian immune response. The DB model is based on data in mice, while the MK model was built to be applicable to humans. This makes FIRM a hybrid model containing features of a few mammalian species. That being said, validation of models against data is a tool of paramount importance in model development. FIRM provides predictions of timing and extent of immune respone, which can be compared against experimental data. While we propose here the FIRM structure and have strived to maintain consistency of units across the models, we make no attempt at validation and later we propose approaches to test this important question.
Addition of Tregand TGF-β: With the integration of the MK, DB, and BL models, we have created a platform with which basic immunological simulations can be performed. The FIRM platform is flexible enough to easily introduce new information and physiological responses. This section demonstrates that principle through addition of regulatory T-cells (Treg). Treg play an important role in the immune system’s response to a tumor. Treg are a rather new discovery in the immunological field, and considered to be of great importance. Therefore, Treg were selected to be the first addition when expanding the FIRM platform beyond the original publications.
The basal state: As with all kinetic models, FIRM has a basal (unperturbed) state. In this particular case, the basal state represents the immune system components’ baselines when there is no antigen present. Such a basal state reflects steady state homeostasis and was calculated using the nominal parameters and running the model to stable steady state conditions. Calculated basal state cell populations include: 108 resting macrophages in the target organ, 5 × 107 dendritic cells in the target organ, ~8 × 104 naïve T-cells in the lymphoid T, ~103 T helper precursor cells in the lymphoid T, 3 × 103 T helper precursor t cells in the blood, ~2 × 104 naïve B-cells in the site of recognition, ~103 cytotoxic T precursors in the lymphoid T, 5 × 104 cytotoxic T precursors in the blood, ~3 × 102 molecules of IL-12 and ~3 × 102 of TGF-Beta in the target organ, 104 Treg in the blood, 2 × 102 in the lymphoid T and 106 in the target organ, and ~2 × 107 free receptor sites on the naïve B-cells in the site of recognition. The Treg in the blood account for about 9.1% of all T cells circulating in the blood . In the basal state, the system is free of antigens, and, therefore, antibodies. This basal state can then be perturbed by exposure to antigen reflecting various stimuli or pathological situations. The response to one or many antigens can be simulated.
The final version of the FIRM simulator has 55 nodes (cells and antibodies), 107 distinct processes and 171 parameters. As defined and with the postulated interactions and parameter values described above, the FIRM simulator is an initial step towards a simulator of the immune response capable of representing a variety of putative challenges. To demonstrate the use of FIRM we present four case studies that illustrate its different features and capabilities.
(2) Use of the simulator – case studies
Blood borne pathogen infection;
Spontaneous tumor rejection;
Influence of Treg on tumor rejection.
The first two case studies represent confirmatory simulations for comparison with the original MK and BL model results.
To simulate a pure TB infection, 100,000 extracellular bacteria were introduced into the lung (target organ) and flux v87 was deactivated. Flux v87 represents the permeation of extracellular bacteria between the target organ and the blood (Figure 1A). Flux v87 was thus deactivated to prevent the bacteria from migrating away from the lung. The TB infection is known to stay local in the lung and not permeate into the blood. FIRM was then simulated with its nominal parameters.
Blood borne pathogen
The humoral (B-cell) response was used by FIRM to eliminate the infection. In fact, in the case study, as expected, the cellular response is virtually non-existent. The antigens never come in contact with the dendritic cells in the tissue. Therefore, there is little antigen presentation to the T-cells in the lymphoid T to drive differentiation of the naïve T-cells. The graphs in Figure 8 show a fast humoral response flooding the blood vessels with antibodies. The antigens are quickly bound to the antibodies in the blood and are removed from the system as antibody-antigen complexes. This simulation resembles the behavior of the original BL model, confirming the performance of the integrated FIRM against the individual submodels.
Since no tumor-antibody interaction is defined in FIRM at this point (all the individual models, MK, DB or BL, lacked a mechanistic description of antibody-mediated cell kill), simulations of tumors will not include migration of tumor antigens (debris in the DB model) into the blood from the target organ. Therefore, the flux v88 was deactivated for this case study. Additionally, we changed the half-life of all T helper cells to 0.02 day-1 from 0.3333 day-1 (the remaining reaction rates were left at the nominal parameter values). The updated half-life was derived from the DB model, while the previous half-life was derived from the MK model. As a justification, through successive model runs we observed that an elevated half-life for the T-cells prevented them from encouraging proliferation of cytotoxic T-cells and therefore prevented tumor kill. This shows that tumor removal by the immune system may be influenced by the kinetics of these cells and also points to features of the model that greatly influence its predictions.
What truly happens in the mammalian immune system following challenge with a single proliferating tumor cell is unknown. However, it has been assumed that the immune system regularly removes tumor cells that may arise spontaneously (the concept of immunosurveillance). If anything, this simulation highlights the role played by the T helper population interaction with the tumor. In our hands, a change in this population’s reported half life to 0.02 day-1 from 0.3333 day-1 was sufficient to produce a more realistic output; however, the model structure may also have been incorrectly specified in the original models. Overall, these considerations point to experimental and research areas of focus. In addition, other parameters in the model may have similar strong influences on the model predictions.
Influence of Treg on tumor rejection
The Treg addition had a profound effect on how the comprehensive system reacts to the presence of a tumor. As such, it is a good example of how the FIRM simulator can be expanded to include additional, more realistic characteristics of immunological responses. This addition reduces to practice how FIRM could be expanded towards the ultimate goal of building a more comprehensive mathematical representation of immune system features.
Inflammation and immune response are thought to be a common denominator in human disease. A comprehensive simulator of the immune response in human tissues is thus a needed tool for a variety of applications. FIRM was undertaken as an initial step towards meeting this need. The development and use of FIRM showed that: (1) already developed and proven methods from systems biology could be used to facilitate the construction of a platform for the integration of subset models each focusing on a feature of the immune response; (2) these methods can be successfully implemented to generate an integrated simulator; and (3) that FIRM can account for a variety of challenges to the human immune system through its multi-scale characteristics. A structure such as FIRM can be used prospectively for iterative model building, as the scope of the simulator grows and as new discoveries are made and integrated in the framework; at the same time, the predictions from FIRM can be compared against experimental data, to improve the model and, consequently, mechanistic understanding of its underlying biology.
As molecular systems biology has developed over the past decade, large-scale and even genome-scale models have been formulated . These models are formed from reconstructed networks based on biochemical, genetic and genomic data, and thus have been able to compute a variety of phenotypic functions. This process has been particularly successful for metabolic models . Since such models are mostly based on the universal principles of flux or mass balance, this process can be applied to systems analysis of complex dynamic systems, such as the immune response. We were able to build FIRM based on a reconstructed network of the main cellular and molecular components involved in the immune response. The key challenge over the individual component models is that FIRM is built on multiple tissue spaces of different volumes, even varying volumes in some cases. Additionally, the model equations are formulated in terms of total counts of each variable in a given tissue space. This allows for complete and accurate accounting and balancing of all state variables. Concentrations of cellular and molecular species can then be computed by simply dividing the total amount at each time point with the biological volume of the space, thus allowing for direct comparison to experimental measurements.
The reconstructed network was then populated with information from published models that describe subsets of the immune response. These published models contained detailed information about the mathematical form of the flux equations. In addition, numerical values for all the parameters are provided in these subset models that are typically obtained from measurements or the available experimental literature.
The putative reconstructed network allows the mapping of multiple subset models and their ready integration under a unified format. In principle this is a simple process, but in practice it has been implemented as a stepwise procedure that reveals subtle integration challenges as additional subset cellular and regulatory processes are added. One of the major integration problems arises when discrepancies in the numerical values for the same parameter appear in different subset models. For example, in some cases, assumptions about the physiological processes accounted for in a subset model were not needed in the comprehensive network reconstruction and thus had to be relaxed or otherwise modified (just as an example, eliminating “HTL” (TH1) from macrophage activation in the DB model). This was relatively easy to accomplish since the network reconstruction has a manageable amount of biological and biochemical detail.
Integration of three T-cell responses, the full B-cell response and the regulatory action of key growth factors was performed. All these subsets of the immune response form a postulated, coherent whole as described by FIRM, which by its integrated nature is more comprehensive than any of the individual subsets and can be further expanded for additional realism. FIRM is thus capable of simulating both tumor and pathogen challenges to the immune system, either separately, or simultaneously. When FIRM is applied to simulate such challenges, it displays appealing multi-scale (time, cellular events, etc.) characteristics. This was demonstrated through case studies representing the formation and eradication of a tumor, and the response to TB infection. The challenges encountered in integrating the FIRM network and its component models are actually representative of the obstacles that can be encountered when attempting to synthesize published mathematical models in a cohesive whole. All biological system models are usually delimited by the boundaries of the system being studied – they are not comprehensive because they cannot be. Their scope is usually limited to the original question they were designed to answer. The biological modeling community is currently entering a new stage: it already moved from the development of individual models to the definition of databases and repositories for these models to be shared; the next evolution is to provide robust tools for on-demand component model integration. This task is addressed in this work. When performing model integration, care must be taken not to reach over the scope of the original component models, especially by trying to incorporate or match aspects that cannot be generalized. This is best decided on a case by case basis. FIRM is a starting point for the modeling community to consider these issues.
Given its current state and the iterative model building enabled by the systems biology approach used in its development, FIRM can and should be expanded to account for additional components of the immune response as new knowledge becomes available. An important area where FIRM could be expanded is the innate response, including natural killer (NK) cells. Models have been recently developed  to characterize this response component. It would also be useful to examine the behavior of FIRM in the context of reinfection, through refinement of the memory B-cell models and kinetics. In addition, a mechanism for antibody-dependent cellular cytotoxicity (ADCC) would provide an unambiguous link between the DB and the BL model and would enrich the FIRM structure. Perhaps the most important missing feature from FIRM is the activation of the B-cell cascade by the T-helper 2 cells, which is a known, essential mechanism for humoral response initiation. All these remain important areas for further development, whose investigation is made somewhat easier by FIRM’s modular structure.
Since it is possible that a range of parameter values could provide similar unperturbed state values, another area of future investigation would involve sensitivity analysis of the model. To explore this, the model could be expanded to include stochastic behavior so as to account for random variation of the state variables around a physiological set point. Additionally, the model could incorporate expected daily fluctuations, such as e.g. circadian changes, to account for the fact that basal states exhibit a range of behaviors around a baseline. Such analyses might be best performed when the model is used in conjunction with experimental studies, such that development of FIRM (or similar models) can be supported by joint simulations and experimentation.
This work exemplified an approach to constructing large-scale integrated simulators of complex dynamic structures such as the immune system. Ours is a hybrid approach bringing together conventional differential equation-based models with methods developed in systems biology over the past decade. The implementation of this approach leads to a postulated representation of the immune system that incorporates the underlying cellular processes and cytokine regulation based on elements of mammalian physiology. The simulator provided by FIRM is well suited to go through an integrated and interactive model building process with experimental validation to reach an increasing state of completion, similar to what as has been accomplished for genome-scale models of metabolism [37–39].
Initial discussions with Kenneth Luu on model development and the support of Scott Fountain are gratefully acknowledged. This work was partially presented at the American Conference of Pharmacometrics (ACoP) held in San Diego, CA, April 3–6, 2011 and was published in the conference proceedings in abstract form.
- Anderson ARA, Quaranta V: Integrative mathematical oncology. Nat Rev Cancer. 2008, 8 (3): 227-234. 10.1038/nrc2329.PubMedView ArticleGoogle Scholar
- Bauer AL, Beauchemin CAA, Perelson AS: Agent-based modeling of host-pathogen systems: The successes and challenges. Inform Sci. 2009, 179 (10): 1379-1389. 10.1016/j.ins.2008.11.012.View ArticleGoogle Scholar
- Mestas J, Hughes CCW: Of Mice and Not Men: Differences between Mouse and Human Immunology. J Immunol. 2004, 172 (5): 2731-2738.PubMedView ArticleGoogle Scholar
- Perelson AS, Weisbuch G: Immunology for physicists. Rev Mod Phys. 1997, 69 (4): 1219-1268. 10.1103/RevModPhys.69.1219.View ArticleGoogle Scholar
- De Boer RJ: Macrophage T lymphocyte interactions in the anti-tumor immune response: a mathematical model. J Immunol. 1985, 134 (4): 2748-2758.PubMedGoogle Scholar
- De Pillis LG, Radunskaya AE, Wiseman CL: A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 2005, 65 (17): 7950-7958.PubMedGoogle Scholar
- Kirschner D, Panetta JC: Modeling immunotherapy of the tumor–immune interaction. J Math Biol. 1998, 37 (3): 235-252. 10.1007/s002850050127.PubMedView ArticleGoogle Scholar
- Bell GI: Mathematical model of clonal selection and antibody production. J Theor Biol. 1970, 29 (2): 191-232. 10.1016/0022-5193(70)90019-6.PubMedView ArticleGoogle Scholar
- Kolev M: Mathematical modeling of the competition between acquired immunity and cancer. Int J Appl Math Comput Sci. 2003, 13 (3): 289-296.Google Scholar
- Marino S, Kirschner DE: The human immune response to Mycobacterium tuberculosis in lung and lymph node. J Theor Biol. 2004, 227 (4): 463-486. 10.1016/j.jtbi.2003.11.023.PubMedView ArticleGoogle Scholar
- Burke MA: Modeling the Proliferative Response of T Cells to IL-2 and IL-4. Cell Immunol. 1997, 178 (1): 42-52. 10.1006/cimm.1997.1125.PubMedView ArticleGoogle Scholar
- Cheng Y: A discrete computer model of the immune system reveals competitive interactions between the humoral and cellular branch and between cross-reacting memory and naïve responses. Vaccine. 2009, 27 (6): 833-845. 10.1016/j.vaccine.2008.11.109.PubMedPubMed CentralView ArticleGoogle Scholar
- Celada F, Seiden PE: A computer model of cellular interactions in the immune system. Immunol Today. 1992, 13 (2): 56-62. 10.1016/0167-5699(92)90135-T.PubMedView ArticleGoogle Scholar
- Castiglione F: A Network of Cellular Automata for the Simulation of the Immune System. INTERNATIONAL JOURNAL OF MODERN PHYSICS C. 1999, 10: 677-686. 10.1142/S0129183199000516.View ArticleGoogle Scholar
- Puzone R: IMMSIM, a flexible model for in machina experiments on immune system responses. Futur Gener Comput Syst. 2002, 18 (7): 961-972. 10.1016/S0167-739X(02)00075-4.View ArticleGoogle Scholar
- Kohler B: A systematic approach to vaccine complexity using an automaton model of the cellular and humoral immune system: I. Viral characteristics and polarized responses. Vaccine. 2000, 19 (7–8): 862-876.PubMedGoogle Scholar
- Cappuccio A, Castiglione F, Piccoli B: Determination of the optimal therapeutic protocols in cancer immunotherapy. Math Biosci. 2007, 209 (1): 1-13. 10.1016/j.mbs.2007.02.009.PubMedView ArticleGoogle Scholar
- Castiglione F, Piccoli B: Cancer immunotherapy, mathematical modeling and optimal control. J Theor Biol. 2007, 247 (4): 723-732. 10.1016/j.jtbi.2007.04.003.PubMedView ArticleGoogle Scholar
- Kirschner D, Marino S: Mycobacterium tuberculosis as viewed through a computer. Trends Microbiol. 2005, 13 (5): 206-211. 10.1016/j.tim.2005.03.005.PubMedView ArticleGoogle Scholar
- Kirschner DE: Toward a multiscale model of antigen presentation in immunity. Immunol Rev. 2007, 216 (1): 93-PubMedView ArticleGoogle Scholar
- Kim PS, Lee PP, Levy D: Modeling regulation mechanisms in the immune system. J Theor Biol. 2007, 246 (1): 33-69. 10.1016/j.jtbi.2006.12.012.PubMedView ArticleGoogle Scholar
- Agoram BM, Martin SW, van der Graaf PH: The role of mechanism-based pharmacokinetic-pharmacodynamic (PK-PD) modelling in translational research of biologics. Drug Discov Today. 2007, 12 (23–24): 1018-1024.PubMedView ArticleGoogle Scholar
- Mager DE, Jusko WJ: Development of Translational Pharmacokinetic-Pharmacodynamic Models. Clin Pharmacol Ther. 2008, 83 (6): 909-912. 10.1038/clpt.2008.52.PubMedPubMed CentralView ArticleGoogle Scholar
- Melder RJ: Systemic distribution and tumor localization of adoptively transferred lymphocytes in mice: comparison with physiologically based pharmacokinetic model. Neoplasia. 2002, 4 (1): 3-8. 10.1038/sj.neo.7900209.PubMedPubMed CentralView ArticleGoogle Scholar
- Garg A, Balthasar J: Physiologically-based pharmacokinetic (PBPK) model to predict IgG tissue kinetics in wild-type and FcRn-knockout mice. J Pharmacokinet Pharmacodyn. 2007, 34 (5): 687-709. 10.1007/s10928-007-9065-1.PubMedView ArticleGoogle Scholar
- Mager DE, Wyska E, Jusko WJ: Diversity of Mechanism-Based Pharmacodynamic Models. Drug Metab Dispos. 2003, 31 (5): 510-518. 10.1124/dmd.31.5.510.PubMedView ArticleGoogle Scholar
- Vicini P: Multiscale Modeling in Drug Discovery and Development: Future Opportunities and Present Challenges. Clin Pharmacol Ther. 2010, 88 (1): 126-129. 10.1038/clpt.2010.87.PubMedView ArticleGoogle Scholar
- Sorger PK: An NIH White Paper by the QSP Workshop Group. Quantitative and Systems Pharmacology in the Post-genomic Era: New Approaches to Discovering Drugs and Understanding Therapeutic Mechanisms. 2011,http://www.nigms.nih.gov/nr/rdonlyres/8ecb1f7c-be3b-431f-89e6-a43411811ab1/0/systemspharmawpsorger2011.pdf,Google Scholar
- Webb K, White T: UML as a cell and biochemistry modeling language. Biosystems. 2005, 80 (3): 283-302. 10.1016/j.biosystems.2004.12.003.PubMedView ArticleGoogle Scholar
- Ayyadurai V, Dewey C: CytoSolve: A Scalable Computational Method for Dynamic Integration of Multiple Molecular Pathway Models. Cell Mol Bioeng. 2011, 4 (1): 28-45. 10.1007/s12195-010-0143-x.PubMedPubMed CentralView ArticleGoogle Scholar
- Palsson B: Systems biology: properties of reconstructed networks. 2006, Cambridge: Cambridge Univ PrView ArticleGoogle Scholar
- Jamshidi N, Palsson BØ: Mass Action Stoichiometric Simulation Models: Incorporating Kinetics and Regulation into Stoichiometric Models. Biophys J. 2010, 98 (2): 175-185. 10.1016/j.bpj.2009.09.064.PubMedPubMed CentralView ArticleGoogle Scholar
- Weinberg RA: The biology of cancer. 2007, New York: Garland ScienceGoogle Scholar
- Colombo MP, Piconese S: Regulatory-T-cell inhibition versus depletion: the right choice in cancer immunotherapy. Nat Rev Cancer. 2007, 7 (11): 880-887. 10.1038/nrc2250.PubMedView ArticleGoogle Scholar
- Park J: Natural immunosurveillance against spontaneous, autochthonous breast cancers revealed and enhanced by blockade of IL-13-mediated negative regulation. Cancer Immunol Immunother. 2008, 57 (6): 907-912. 10.1007/s00262-007-0414-0.PubMedView ArticleGoogle Scholar
- Edwards JS, Palsson BO: Systems Properties of the Haemophilus influenzaeRd Metabolic Genotype. J Biol Chem. 1999, 274 (25): 17410-17416. 10.1074/jbc.274.25.17410.PubMedView ArticleGoogle Scholar
- Feist AM, Palsson BO: The growing scope of applications of genome-scale metabolic reconstructions using Escherichia coli. Nat Biotech. 2008, 26 (6): 659-667. 10.1038/nbt1401.View ArticleGoogle Scholar
- Duarte NC: Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proc Natl Acad Sci. 2007, 104 (6): 1777-1782. 10.1073/pnas.0610772104.PubMedPubMed CentralView ArticleGoogle Scholar
- Gianchandani EP: Functional States of the Genome-Scale Escherichia Coli Transcriptional Regulatory System. PLoS Comput Biol. 2009, 5 (6): e1000403-10.1371/journal.pcbi.1000403.PubMedPubMed CentralView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.