Systems biology in immunology: a computational modeling perspective.

Systems biology is an emerging discipline that combines high-content, multiplexed measurements with informatic and computational modeling methods to better understand biological function at various scales. Here we present a detailed review of the methods used to create computational models and to conduct simulations of immune function. We provide descriptions of the key data-gathering techniques employed to generate the quantitative and qualitative data required for such modeling and simulation and summarize the progress to date in applying these tools and techniques to questions of immunological interest, including infectious disease. We include comments on what insights modeling can provide that complement information obtained from the more familiar experimental discovery methods used by most investigators and the reasons why quantitative methods are needed to eventually produce a better understanding of immune system operation in health and disease.

Migrating Myeloid Cells Sense Temporal Dynamics of Chemoattractant Concentrations.

Chemoattractant-mediated recruitment of hematopoietic cells to sites of pathogen growth or tissue damage is critical to host defense and organ homeostasis. Chemotaxis is typically considered to rely on spatial sensing, with cells following concentration gradients as long as these are present. Utilizing a microfluidic approach, we found that stable gradients of intermediate chemokines (CCL19 and CXCL12) failed to promote persistent directional migration of dendritic cells or neutrophils. Instead, rising chemokine concentrations were needed, implying that temporal sensing mechanisms controlled prolonged responses to these ligands. This behavior was found to depend on G-coupled receptor kinase-mediated negative regulation of receptor signaling and contrasted with responses to an end agonist chemoattractant (C5a), for which a stable gradient led to persistent migration. These findings identify temporal sensing as a key requirement for long-range myeloid cell migration to intermediate chemokines and provide insights into the mechanisms controlling immune cell motility in complex tissue environments.

Inflammation-induced interstitial migration of effector CD4⁺ T cells is dependent on integrin αV.

Leukocytes must traverse inflamed tissues to effectively control local infection. Although motility in dense tissues seems to be integrin independent and based on actomyosin-mediated protrusion and contraction, during inflammation, changes to the extracellular matrix (ECM) may necessitate distinct motility requirements. Indeed, we found that the interstitial motility of T cells was critically dependent on Arg-Gly-Asp (RGD)-binding integrins in the inflamed dermis. Inflammation-induced deposition of fibronectin was functionally linked to higher expression of integrin αV on effector CD4⁺ T cells. By intravital multiphoton imaging, we found that the motility of CD4⁺ T cells was dependent on αV expression. Selective blockade or knockdown of αV arrested T helper type 1 (TH1) cells in the inflamed tissue and attenuated local effector function. Our data demonstrate context-dependent specificity of lymphocyte movement in inflamed tissues that is essential for protective immunity.

An integrated approach to the characterization of immune repertoires using AIMS: An Automated Immune Molecule Separator.

The adaptive immune system employs an array of receptors designed to respond with high specificity to pathogens or molecular aberrations faced by the host organism. Binding of these receptors to molecular fragments-collectively referred to as antigens-initiates immune responses. These antigenic targets are recognized in their native state on the surfaces of pathogens by antibodies, whereas T cell receptors (TCR) recognize processed antigens as short peptides, presented on major histocompatibility complex (MHC) molecules. Recent research has led to a wealth of immune repertoire data that are key to interrogating the nature of these molecular interactions. However, existing tools for the analysis of these large datasets typically focus on molecular sets of a single type, forcing researchers to separately analyze strongly coupled sequences of interacting molecules. Here, we introduce a software package for the integrated analysis of immune repertoire data, capable of identifying distinct biophysical differences in isolated TCR, MHC, peptide, antibody, and antigen sequence data. This integrated analytical approach allows for direct comparisons across immune repertoire subsets and provides a starting point for the identification of key interaction hotspots in complementary receptor-antigen pairs. The software (AIMS-Automated Immune Molecule Separator) is freely available as an open access package in GUI or command-line form.

SBML Level 3: an extensible format for the exchange and reuse of biological models.

Systems biology has experienced dramatic growth in the number, size, and complexity of computational models. To reproduce simulation results and reuse models, researchers must exchange unambiguous model descriptions. We review the latest edition of the Systems Biology Markup Language (SBML), a format designed for this purpose. A community of modelers and software authors developed SBML Level 3 over the past decade. Its modular form consists of a core suited to representing reaction-based models and packages that extend the core with features suited to other model types including constraint-based models, reaction-diffusion models, logical network models, and rule-based models. The format leverages two decades of SBML and a rich software ecosystem that transformed how systems biologists build and interact with models. More recently, the rise of multiscale models of whole cells and organs, and new data sources such as single-cell measurements and live imaging, has precipitated new ways of integrating data with models. We provide our perspectives on the challenges presented by these developments and how SBML Level 3 provides the foundation needed to support this evolution.

Space-time histories approach to fast stochastic simulation of bimolecular reactions.

Computational models of reaction-diffusion systems involving low copy numbers or strongly heterogeneous molecular spatial distributions, such as those frequently found in cellular signaling pathways, require approaches that account for the stochastic dynamics of individual particles, as opposed to approaches representing them through their average concentrations. Efforts to remedy the high computational cost associated with particle-based stochastic approaches by taking advantage of Green's functions are hampered by the need to draw random numbers from complicated, and therefore costly, non-standard probability distributions to update particle positions. Here, we introduce an approach that permits the reconstruction of entire molecular trajectories, including bimolecular encounters, in retrospect, after a simulated time step, while avoiding inefficient draws from non-standard distributions. This means that highly accurate stochastic simulations can be performed for system sizes that would be prohibitively costly to simulate with conventional Green's function based methods. The algorithm applies equally well to one, two, and three dimensional systems and can be readily extended to include deterministic forces specified by an interaction potential, such as the Coulomb potential.

Targeted Proteomics-Driven Computational Modeling of Macrophage S1P Chemosensing.

Osteoclasts are monocyte-derived multinuclear cells that directly attach to and resorb bone. Sphingosine-1-phosphate (S1P)(1) regulates bone resorption by functioning as both a chemoattractant and chemorepellent of osteoclast precursors through two G-protein coupled receptors that antagonize each other in an S1P-concentration-dependent manner. To quantitatively explore the behavior of this chemosensing pathway, we applied targeted proteomics, transcriptomics, and rule-based pathway modeling using the Simmune toolset. RAW264.7 cells (a mouse monocyte/macrophage cell line) were used as model osteoclast precursors, RNA-seq was used to identify expressed target proteins, and selected reaction monitoring (SRM) mass spectrometry using internal peptide standards was used to perform absolute abundance measurements of pathway proteins. The resulting transcript and protein abundance values were strongly correlated. Measured protein abundance values, used as simulation input parameters, led to in silico pathway behavior matching in vitro measurements. Moreover, once model parameters were established, even simulated responses toward stimuli that were not used for parameterization were consistent with experimental findings. These findings demonstrate the feasibility and value of combining targeted mass spectrometry with pathway modeling for advancing biological insight.

Computational modeling of cellular signaling processes embedded into dynamic spatial contexts.

Cellular signaling processes depend on spatiotemporal distributions of molecular components. Multicolor, high-resolution microscopy permits detailed assessment of such distributions, providing input for fine-grained computational models that explore mechanisms governing dynamic assembly of multimolecular complexes and their role in shaping cellular behavior. However, it is challenging to incorporate into such models both complex molecular reaction cascades and the spatial localization of signaling components in dynamic cellular morphologies. Here we introduce an approach to address these challenges by automatically generating computational representations of complex reaction networks based on simple bimolecular interaction rules embedded into detailed, adaptive models of cellular morphology. Using examples of receptor-mediated cellular adhesion and signal-induced localized mitogen-activated protein kinase (MAPK) activation in yeast, we illustrate the capacity of this simulation technique to provide insights into cell biological processes. The modeling algorithms, implemented in a new version of the Simmune toolset, are accessible through intuitive graphical interfaces and programming libraries.

Sphingosine-1-phosphate mobilizes osteoclast precursors and regulates bone homeostasis.

Osteoclasts are the only somatic cells with bone-resorbing capacity and, as such, they have a critical role not only in normal bone homeostasis (called 'bone remodelling') but also in the pathogenesis of bone destructive disorders such as rheumatoid arthritis and osteoporosis. A major focus of research in the field has been on gene regulation by osteoclastogenic cytokines such as receptor activator of NF-kappaB-ligand (RANKL, also known as TNFSF11) and TNF-alpha, both of which have been well documented to contribute to osteoclast terminal differentiation. A crucial process that has been less well studied is the trafficking of osteoclast precursors to and from the bone surface, where they undergo cell fusion to form the fully differentiated multinucleated cells that mediate bone resorption. Here we report that sphingosine-1-phosphate (S1P), a lipid mediator enriched in blood, induces chemotaxis and regulates the migration of osteoclast precursors not only in culture but also in vivo, contributing to the dynamic control of bone mineral homeostasis. Cells with the properties of osteoclast precursors express functional S1P(1) receptors and exhibit positive chemotaxis along an S1P gradient in vitro. Intravital two-photon imaging of bone tissues showed that a potent S1P(1) agonist, SEW2871, stimulated motility of osteoclast precursor-containing monocytoid populations in vivo. Osteoclast/monocyte (CD11b, also known as ITGAM) lineage-specific conditional S1P(1) knockout mice showed osteoporotic changes due to increased osteoclast attachment to the bone surface. Furthermore, treatment with the S1P(1) agonist FTY720 relieved ovariectomy-induced osteoporosis in mice by reducing the number of mature osteoclasts attached to the bone surface. Together, these data provide evidence that S1P controls the migratory behaviour of osteoclast precursors, dynamically regulating bone mineral homeostasis, and identifies a critical control point in osteoclastogenesis that may have potential as a therapeutic target.

Opposing roles for RhoH GTPase during T-cell migration and activation.

T cells spend the majority of their time perusing lymphoid organs in search of cognate antigen presented by antigen presenting cells (APCs) and then quickly recirculate through the bloodstream to another lymph node. Therefore, regulation of a T-cell response is dependent upon the ability of cells to arrive in the correct location following chemokine gradients ("go" signal) as well as to receive appropriate T-cell receptor (TCR) activation signals upon cognate antigen recognition ("stop" signal). However, the mechanisms by which T cells regulate these go and stop signals remain unclear. We found that overexpression of the hematopoietic-specific RhoH protein in the presence of chemokine signals resulted in decreased Rap1-GTP and LFA-1 adhesiveness to ICAM-1, thus impairing T-cell chemotaxis; while in the presence of TCR signals, there were enhanced and sustained Rap1-GTP and LFA-1 activation as well as prolonged T:APC conjugates. RT-PCR analyses of activated CD4(+) T cells and live images of T-cell migration and immunological synapse (IS) formation revealed that functions of RhoH took place primarily at the levels of transcription and intracellular distribution. Thus, we conclude that RhoH expression provides a key molecular determinant that allows T cells to switch between sensing chemokine-mediated go signals and TCR-dependent stop signals.

High production rates sustain in vivo levels of PD-1high simian immunodeficiency virus-specific CD8 T cells in the face of rapid clearance.

Programmed Death 1 (PD-1) expression by human/simian immunodeficiency virus (HIV/SIV)-specific CD8 T cells has been associated with defective cytokine production and reduced in vitro proliferation capacity. However, the cellular mechanisms that sustain PD-1(high) virus-specific CD8 T cell responses during chronic infection are unknown. Here, we show that the PD-1(high) phenotype is associated with accelerated in vivo CD8 T cell turnover in SIV-infected rhesus macaques, especially within the SIV-specific CD8 T cell pool. Mathematical modeling of 5-bromo-2' deoxyuridine (BrdU) labeling dynamics demonstrated a significantly increased generation rate of PD-1(high) compared to PD-1(low) CD8 T cells in all memory compartments. Simultaneous analysis of Ki67 and BrdU kinetics revealed a complex in vivo turnover profile whereby only a small fraction of PD-1(high) cells, but virtually all PD-1(low) cells, returned to rest after activation. Similar kinetics operated in both chronic and acute SIV infection. Our data suggest that the persistence of PD-1(high) SIV-specific CD8 T cells in chronic infection is maintained in vivo by a mechanism involving high production coupled with a high disappearance rate.

The Simmune Modeler visual interface for creating signaling networks based on bi-molecular interactions.

MOTIVATION: Biochemical modeling efforts now frequently take advantage of the possibility to automatically create reaction networks based on the specification of pairwise molecular interactions. Even though a variety of tools exist to visualize the resulting networks, defining the rules for the molecular interactions typically requires writing scripts, which impacts the non-specialist accessibility of those approaches. We introduce the Simmune Modeler that allows users to specify molecular complexes and their interactions as well as the reaction-induced modifications of the molecules through a flexible visual interface. It can take into account the positions of the components of trans-membrane complexes relative to the embedding membranes as well as symmetry aspects affecting the reactions of multimeric molecular structures. Models created with this tool can be simulated using the Simmune Simulator or be exported as SBML code or as files describing the reaction networks as systems of ODEs for import into Matlab. AVAILABILITY: The Simmune Modeler and the associated simulators as well as extensive additional documentation and tutorials are freely available for Linux, Mac and Windows: http://go.usa.gov/QeH (Note shortened case-sensitive URL!).

Redirecting cell-type specific cytokine responses with engineered interleukin-4 superkines.

Cytokines dimerize their receptors, with the binding of the 'second chain' triggering signaling. In the interleukin (IL)-4 and IL-13 system, different cell types express varying numbers of alternative second receptor chains (γc or IL-13Rα1), forming functionally distinct type I or type II complexes. We manipulated the affinity and specificity of second chain recruitment by human IL-4. A type I receptor-selective IL-4 'superkine' with 3,700-fold higher affinity for γc was three- to ten-fold more potent than wild-type IL-4. Conversely, a variant with high affinity for IL-13Rα1 more potently activated cells expressing the type II receptor and induced differentiation of dendritic cells from monocytes, implicating the type II receptor in this process. Superkines showed signaling advantages on cells with lower second chain numbers. Comparative transcriptional analysis reveals that the superkines induce largely redundant gene expression profiles. Variable second chain numbers can be exploited to redirect cytokines toward distinct cell subsets and elicit new actions, potentially improving the selectivity of cytokine therapy.

Rate coefficients, binding probabilities, and related quantities for area reactivity models.

We further develop the general theory of the area reactivity model that describes the diffusion-influenced reaction of an isolated receptor-ligand pair in terms of a generalized Feynman-Kac equation and that provides an alternative to the classical contact reactivity model. Analyzing both the irreversible and reversible reaction, we derive the equation of motion of the survival probability as well as several relationships between single pair quantities and the reactive flux at the encounter distance. Building on these relationships, we derive the equation of motion of the many-particle survival probability for irreversible pseudo-first-order reactions. Moreover, we show that the usual definition of the rate coefficient as the reactive flux is deficient in the area reactivity model. Numerical tests for our findings are provided through Brownian Dynamics simulations. We calculate exact and approximate expressions for the irreversible rate coefficient and show that this quantity behaves differently from its classical counterpart. Furthermore, we derive approximate expressions for the binding probability as well as the average lifetime of the bound state and discuss on- and off-rates in this context. Throughout our approach, we point out similarities and differences between the area reactivity model and its classical counterpart, the contact reactivity model. The presented analysis and obtained results provide a theoretical framework that will facilitate the comparison of experiment and model predictions.

The area reactivity model of geminate recombination.

We investigate the reversible diffusion-influenced reaction of an isolated pair in the context of the area reactivity model that describes the reversible binding of a single molecule in the presence of a binding site in terms of a generalized version of the Feynman-Kac equation in two dimensions. We compute the corresponding exact Green's function in the Laplace domain for both the initially unbound and bound molecule. We discuss convolution relations that facilitate the calculation of the binding and survival probabilities. Furthermore, we calculate an exact analytical expression for the Green's function in the time domain by inverting the Laplace transform via the Bromwich contour integral and derive expressions for the binding and survival probability in the time domain as well. We numerically confirm the accuracy of the obtained expressions by propagating the generalized Feynman-Kac equation in the time domain. Our results should be useful for comparing the area reactivity model with the contact reactivity model.

Feedback regulation of proliferation vs. differentiation rates explains the dependence of CD4 T-cell expansion on precursor number.

The mechanisms regulating clonal expansion and contraction of T cells in response to immunization remain to be identified. A recent study established that there was a log-linear relation between CD4 T-cell precursor number (PN) and factor of expansion (FE), with a slope of ∼-0.5 over a range of 3-30,000 precursors per mouse. The results suggested inhibition of precursor expansion either by competition for specific antigen-presenting cells or by the action of other antigen-specific cells in the same microenvironment as the most likely explanation. Several molecular mechanisms potentially accounting for such inhibition were examined and rejected. Here we adopt a previously proposed concept, "feedback-regulated balance of growth and differentiation," and show that it can explain the observed findings. We assume that the most differentiated effectors (or memory cells) limit the growth of less differentiated effectors, locally, by increasing the rate of differentiation of the latter cells in a dose-dependent manner. Consequently, expansion is blocked and reversed after a delay that depends on initial PN, accounting for the dependence of the peak of the response on that number. We present a parsimonious mathematical model capable of reproducing immunization response kinetics. Model definition is achieved in part by requiring consistency with available BrdU-labeling and carboxyfluorescein diacetate succinimidyl ester (CFSE)-dilution data. The calibrated model correctly predicts FE as a function of PN. We conclude that feedback-regulated balance of growth and differentiation, although awaiting definite experimental characterization of the hypothetical cells and molecules involved in regulation, can explain the kinetics of CD4 T-cell responses to antigenic stimulation.

CountASAP: A Lightweight, Easy to Use Python Package for Processing ASAPseq Data.

Declining sequencing costs coupled with the increasing availability of easy-to-use kits for the isolation of DNA and RNA transcripts from single cells have driven a rapid proliferation of studies centered around genomic and transcriptomic data. Simultaneously, a wealth of new techniques have been developed that utilize single cell technologies to interrogate a broad range of cell-biological processes. One recently developed technique, transposase-accessible chromatin with sequencing (ATAC) with select antigen profiling by sequencing (ASAPseq), provides a combination of chromatin accessibility assessments with measurements of cell-surface marker expression levels. While software exists for the characterization of these datasets, there currently exists no tool explicitly designed to reformat ASAP surface marker FASTQ data into a count matrix which can then be used for these downstream analyses. To address this, we created CountASAP, an easy-to-use Python package purposefully designed to transform FASTQ files from ASAP experiments into count matrices compatible with commonly-used downstream bioinformatic analysis packages. CountASAP takes advantage of the independence of the relevant data structures to perform fully parallelized matches of each sequenced read to user-supplied input ASAP oligos and unique cell-identifier sequences.

Progressive CD4+ central memory T cell decline results in CD4+ effector memory insufficiency and overt disease in chronic SIV infection.

Primary simian immunodeficiency virus (SIV) infections of rhesus macaques result in the dramatic depletion of CD4(+) CCR5(+) effector-memory T (T(EM)) cells from extra-lymphoid effector sites, but in most infections, an increased rate of CD4(+) memory T cell proliferation appears to prevent collapse of effector site CD4(+) T(EM) cell populations and acute-phase AIDS. Eventually, persistent SIV replication results in chronic-phase AIDS, but the responsible mechanisms remain controversial. Here, we demonstrate that in the chronic phase of progressive SIV infection, effector site CD4(+) T(EM) cell populations manifest a slow, continuous decline, and that the degree of this depletion remains a highly significant correlate of late-onset AIDS. We further show that due to persistent immune activation, effector site CD4(+) T(EM) cells are predominantly short-lived, and that their homeostasis is strikingly dependent on the production of new CD4(+) T(EM) cells from central-memory T (T(CM)) cell precursors. The instability of effector site CD4(+) T(EM) cell populations over time was not explained by increasing destruction of these cells, but rather was attributable to progressive reduction in their production, secondary to decreasing numbers of CCR5(-) CD4(+) T(CM) cells. These data suggest that although CD4(+) T(EM) cell depletion is a proximate mechanism of immunodeficiency, the tempo of this depletion and the timing of disease onset are largely determined by destruction, failing production, and gradual decline of CD4(+) T(CM) cells.

Pathogenesis of HIV infection: what the virus spares is as important as what it destroys.

Upon transmission to a new host, HIV targets CCR5+ CD4+ effector memory T cells, resulting in acute, massive depletion of these cells from mucosal effector sites. This depletion does not initially compromise the regenerative capacity of the immune system because naive and most central memory T cells are spared. Here, we discuss evidence suggesting that frequent activation of these spared cells during the chronic phase of HIV infection supplies mucosal tissues with short-lived CCR5+ CD4+ effector cells that prevent life-threatening infections. This immune activation also facilitates continued viral replication, but infection and killing of target T cells by HIV are selective and the impact on effector-cell lifespan is limited. We propose, however, that persistent activation progressively disrupts the functional organization of the immune system, reducing its regenerative capacity and facilitating viral evolution that leads to loss of the exquisite target cell-sparing selectivity of viral replication, ultimately resulting in AIDS.

Theory of reversible diffusion-influenced reactions with non-Markovian dissociation in two space dimensions.

We investigate the reversible diffusion-influenced reaction of an isolated pair in the presence of a non-Markovian generalization of the backreaction boundary condition in two space dimensions. Following earlier work by Agmon and Weiss, we consider residence time probability densities that decay slower than an exponential and that are characterized by a single parameter 0 < σ ≤ 1. We calculate an exact expression for a Green's function of the two-dimensional diffusion equation subject to a non-Markovian backreaction boundary condition that is valid for arbitrary σ and for all times. We use the obtained expression to derive the survival probability for the initially unbound pair and we calculate an exact expression for the probability S(t[line]*) that the initially bound particle is unbound. Finally, we obtain an approximate solution for long times. In particular, we show that the ultimate fate of the bound state is complete dissociation, as in the Markovian case. However, the limiting value is approached quite differently: Instead of a ~t(-1) decay, we obtain 1 - S(t[line]*) ~ t(-σ)ln t. The derived expressions should be relevant for a better understanding of reversible membrane-bound reactions in cell biology.