Impacts of Vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 Variants Alpha and Delta on Coronavirus Disease 2019 Transmission Dynamics in Four Metropolitan Areas of the United States

To characterize Coronavirus Disease 2019 (COVID-19) transmission dynamics in each of the metropolitan statistical areas (MSAs) surrounding Dallas, Houston, New York City, and Phoenix in 2020 and 2021, we extended a previously reported compartmental model accounting for effects of multiple distinct periods of non-pharmaceutical interventions by adding consideration of vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) variants Alpha (lineage B.1.1.7) and Delta (lineage B.1.617.2). For each MSA, we found region-specific parameterizations of the model using daily reports of new COVID-19 cases available from January 21, 2020 to October 31, 2021. In the process, we obtained estimates of the relative infectiousness of Alpha and Delta as well as their takeoff times in each MSA (the times at which sustained transmission began). The estimated infectiousness of Alpha ranged from 1.1x to 1.4x that of viral strains circulating in 2020 and early 2021. The estimated relative infectiousness of Delta was higher in all cases, ranging from 1.6x  to 2.1x. The estimated Alpha takeoff times ranged from February 1 to February 28, 2021. The estimated Delta takeoff times ranged from June 2 to June 26, 2021. Estimated takeoff times are consistent with genomic surveillance data.

Impacts of vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 variants Alpha and Delta on Coronavirus Disease 2019 transmission dynamics in four metropolitan areas of the United States

To characterize Coronavirus Disease 2019 (COVID-19) transmission dynamics in each of the metropolitan statistical areas (MSAs) surrounding Dallas, Houston, New York City, and Phoenix in 2020 and 2021, we extended a previously reported compartmental model accounting for effects of multiple distinct periods of non-pharmaceutical interventions by adding consideration of vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) variants Alpha (lineage B.1.1.7) and Delta (lineage B.1.617.2). For each MSA, we found region-specific parameterizations of the model using daily reports of new COVID-19 cases available from January 21, 2020 to October 31, 2021. In the process, we obtained estimates of the relative infectiousness of Alpha and Delta as well as their takeoff times in each MSA (the times at which sustained transmission began). The estimated infectiousness of Alpha ranged from 1.1x to 1.4x that of viral strains circulating in 2020 and early 2021. The estimated relative infectiousness of Delta was higher in all cases, ranging from 1.6x to 2.1x. The estimated Alpha takeoff times ranged from February 1 to February 28, 2021. The estimated Delta takeoff times ranged from June 2 to June 26, 2021. Estimated takeoff times are consistent with genomic surveillance data.

One-Sentence Summary

Using a compartmental model parameterized to reproduce available reports of new Coronavirus Disease 2019 (COVID-19) cases, we quantified the impacts of vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) variants Alpha (lineage B.1.1.7) and Delta (lineage B.1.617.2) on regional epidemics in the metropolitan statistical areas (MSAs) surrounding Dallas, Houston, New York City, and Phoenix.

Evaluation of FluSight influenza forecasting in the 2021-22 and 2022-23 seasons with a new target laboratory-confirmed influenza hospitalizations

Accurate forecasts can enable more effective public health responses during seasonal influenza epidemics. Forecasting teams were asked to provide national and jurisdiction-specific probabilistic predictions of weekly confirmed influenza hospital admissions for one through four weeks ahead for the 2021-22 and 2022-23 influenza seasons. Across both seasons, 26 teams submitted forecasts, with the submitting teams varying between seasons. Forecast skill was evaluated using the Weighted Interval Score (WIS), relative WIS, and coverage. Six out of 23 models outperformed the baseline model across forecast weeks and locations in 2021-22 and 12 out of 18 models in 2022-23. Averaging across all forecast targets, the FluSight ensemble was the 2nd most accurate model measured by WIS in 2021-22 and the 5th most accurate in the 2022-23 season. Forecast skill and 95% coverage for the FluSight ensemble and most component models degraded over longer forecast horizons and during periods of rapid change. Current influenza forecasting efforts help inform situational awareness, but research is needed to address limitations, including decreased performance during periods of changing epidemic dynamics.

Quantification of early nonpharmaceutical interventions aimed at slowing transmission of Coronavirus Disease 2019 in the Navajo Nation and surrounding states (Arizona, Colorado, New Mexico, and Utah)

During an early period of the Coronavirus Disease 2019 (COVID-19) pandemic, the Navajo Nation, much like New York City, experienced a relatively high rate of disease transmission. Yet, between January and October 2020, it experienced only a single period of growth in new COVID-19 cases, which ended when cases peaked in May 2020. The daily number of new cases slowly decayed in the summer of 2020 until late September 2020. In contrast, the surrounding states of Arizona, Colorado, New Mexico, and Utah all experienced at least two periods of growth in the same time frame, with second surges beginning in late May to early June. Here, we investigated these differences in disease transmission dynamics with the objective of quantifying the contributions of non-pharmaceutical interventions (NPIs) (e.g., behaviors that limit disease transmission). We considered a compartmental model accounting for distinct periods of NPIs to analyze the epidemic in each of the five regions. We used Bayesian inference to estimate region-specific model parameters from regional surveillance data (daily reports of new COVID-19 cases) and to quantify uncertainty in parameter estimates and model predictions. Our results suggest that NPIs in the Navajo Nation were sustained over the period of interest, whereas in the surrounding states, NPIs were relaxed, which allowed for subsequent surges in cases. Our region-specific model parameterizations allow us to quantify the impacts of NPIs on disease incidence in the regions of interest.

Quantification of early nonpharmaceutical interventions aimed at slowing transmission of Coronavirus Disease 2019 in the Navajo Nation and surrounding states (Arizona, Colorado, New Mexico, and Utah)

During an early period of the Coronavirus Disease 2019 (COVID-19) pandemic, the Navajo Nation, much like New York City, experienced a relatively high rate of disease transmission. Yet, between January and October 2020, it experienced only a single period of growth in new COVID-19 cases, which ended when cases peaked in May 2020. The daily number of new cases slowly decayed in the summer of 2020 until late September 2020. In contrast, the surrounding states of Arizona, Colorado, New Mexico, and Utah all experienced at least two periods of growth in the same time frame, with second surges beginning in late May to early June. To investigate the causes of this difference, we used a compartmental model accounting for distinct periods of non-pharmaceutical interventions (NPIs ) ( e.g., behaviors that limit disease transmission) to analyze the epidemic in each of the five regions. We used Bayesian inference to estimate region-specific model parameters from regional surveillance data (daily reports of new COVID-19 cases) and to quantify uncertainty in parameter estimates and model predictions. Our results suggest that NPIs in the Navajo Nation were sustained over the period of interest, whereas in the surrounding states, NPIs were relaxed, which allowed for subsequent surges in cases. Our region-specific model parameterizations allow us to quantify the impacts of NPIs on disease incidence in the regions of interest.

Implementation of a practical Markov chain Monte Carlo sampling algorithm in PyBioNetFit

SUMMARY

Bayesian inference in biological modeling commonly relies on Markov chain Monte Carlo (MCMC) sampling of a multidimensional and non-Gaussian posterior distribution that is not analytically tractable. Here, we present the implementation of a practical MCMC method in the open-source software package PyBioNetFit (PyBNF), which is designed to support parameterization of mathematical models for biological systems. The new MCMC method, am, incorporates an adaptive move proposal distribution. For warm starts, sampling can be initiated at a specified location in parameter space and with a multivariate Gaussian proposal distribution defined initially by a specified covariance matrix. Multiple chains can be generated in parallel using a computer cluster. We demonstrate that am can be used to successfully solve real-world Bayesian inference problems, including forecasting of new Coronavirus Disease 2019 case detection with Bayesian quantification of forecast uncertainty.

AVAILABILITY AND IMPLEMENTATION

PyBNF version 1.1.9, the first stable release with am, is available at PyPI and can be installed using the pip package-management system on platforms that have a working installation of Python 3. PyBNF relies on libRoadRunner and BioNetGen for simulations (e.g. numerical integration of ordinary differential equations defined in SBML or BNGL files) and Dask.Distributed for task scheduling on Linux computer clusters. The Python source code can be freely downloaded/cloned from GitHub and used and modified under terms of the BSD-3 license (https://github.com/lanl/pybnf). Online documentation covering installation/usage is available (https://pybnf.readthedocs.io/en/latest/). A tutorial video is available on YouTube (https://www.youtube.com/watch?v=2aRqpqFOiS4&t=63s).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Bayesian Inference of State-Level COVID-19 Basic Reproduction Numbers across the United States

Although many persons in the United States have acquired immunity to COVID-19, either through vaccination or infection with SARS-CoV-2, COVID-19 will pose an ongoing threat to non-immune persons so long as disease transmission continues. We can estimate when sustained disease transmission will end in a population by calculating the population-specific basic reproduction number ℛ0, the expected number of secondary cases generated by an infected person in the absence of any interventions. The value of ℛ0 relates to a herd immunity threshold (HIT), which is given by 1-1/ℛ0. When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely (barring mutations allowing SARS-CoV-2 to escape immunity). Here, we report state-level ℛ0 estimates obtained using Bayesian inference. Maximum a posteriori estimates range from 7.1 for New Jersey to 2.3 for Wyoming, indicating that disease transmission varies considerably across states and that reaching herd immunity will be more difficult in some states than others. ℛ0 estimates were obtained from compartmental models via the next-generation matrix approach after each model was parameterized using regional daily confirmed case reports of COVID-19 from 21 January 2020 to 21 June 2020. Our ℛ0 estimates characterize the infectiousness of ancestral strains, but they can be used to determine HITs for a distinct, currently dominant circulating strain, such as SARS-CoV-2 variant Delta (lineage B.1.617.2), if the relative infectiousness of the strain can be ascertained. On the basis of Delta-adjusted HITs, vaccination data, and seroprevalence survey data, we found that no state had achieved herd immunity as of 20 September 2021.

Bayesian Inference of State-Level COVID-19 Basic Reproduction Numbers across the United States

Although many persons in the United States have acquired immunity to COVID-19, either through vaccination or infection with SARS-CoV-2, COVID-19 will pose an ongoing threat to non-immune persons so long as disease transmission continues. We can estimate when sustained disease transmission will end in a population by calculating the population-specific basic reproduction number ℛ ₀ , the expected number of secondary cases generated by an infected person in the absence of any interventions. The value of ℛ ₀ relates to a herd immunity threshold (HIT), which is given by 1 - 1/ℛ ₀ . When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely (barring mutations allowing SARS-CoV-2 to escape immunity). Here, we report state-level ℛ ₀ estimates obtained using Bayesian inference. Maximum a posteriori estimates range from 7.1 for New Jersey to 2.3 for Wyoming, indicating that disease transmission varies considerably across states and that reaching herd immunity will be more difficult in some states than others. ℛ ₀ estimates were obtained from compartmental models via the next-generation matrix approach after each model was parameterized using regional daily confirmed case reports of COVID-19 from 21-January-2020 to 21-June-2020. Our ℛ ₀ estimates characterize infectiousness of ancestral strains, but they can be used to determine HITs for a distinct, currently dominant circulating strain, such as SARS-CoV-2 variant Delta (lineage B.1.617.2), if the relative infectiousness of the strain can be ascertained. On the basis of Delta-adjusted HITs, vaccination data, and seroprevalence survey data, we find that no state has achieved herd immunity as of 20-September-2021.

Significance Statement

COVID-19 will continue to threaten non-immune persons in the presence of ongoing disease transmission. We can estimate when sustained disease transmission will end by calculating the population-specific basic reproduction number ℛ ₀ , which relates to a herd immunity threshold (HIT), given by 1 - 1/ℛ ₀ . When the immune fraction of a population exceeds this threshold, sustained disease transmission becomes exponentially unlikely. Here, we report state-level ℛ ₀ estimates indicating that disease transmission varies considerably across states. Our ℛ ₀ estimates can also be used to determine HITs for the Delta variant of COVID-19. On the basis of Delta-adjusted HITs, vaccination data, and serological survey results, we find that no state has yet achieved herd immunity.

Daily Forecasting of Regional Epidemics of Coronavirus Disease with Bayesian Uncertainty Quantification, United States

To increase situational awareness and support evidence-based policymaking, we formulated a mathematical model for coronavirus disease transmission within a regional population. This compartmental model accounts for quarantine, self-isolation, social distancing, a nonexponentially distributed incubation period, asymptomatic persons, and mild and severe forms of symptomatic disease. We used Bayesian inference to calibrate region-specific models for consistency with daily reports of confirmed cases in the 15 most populous metropolitan statistical areas in the United States. We also quantified uncertainty in parameter estimates and forecasts. This online learning approach enables early identification of new trends despite considerable variability in case reporting.

Daily Forecasting of New Cases for Regional Epidemics of Coronavirus Disease 2019 with Bayesian Uncertainty Quantification

To increase situational awareness and support evidence-based policy-making, we formulated a mathematical model for COVID-19 transmission within a regional population. This compartmental model accounts for quarantine, self-isolation, social distancing, a non-exponentially distributed incubation period, asymptomatic individuals, and mild and severe forms of symptomatic disease. Using Bayesian inference, we have been calibrating region-specific models daily for consistency with new reports of confirmed cases from the 15 most populous metropolitan statistical areas in the United States and quantifying uncertainty in parameter estimates and predictions of future case reports. This online learning approach allows for early identification of new trends despite considerable variability in case reporting.

ARTICLE SUMMARY LINE

We report models for regional COVID-19 epidemics and use of Bayesian inference to quantify uncertainty in daily predictions of expected reporting of new cases, enabling identification of new trends in surveillance data.

Daily Forecasting of New Cases for Regional Epidemics of Coronavirus Disease 2019 with Bayesian Uncertainty Quantification

To increase situational awareness and support evidence-based policy-making, we formulated two types of mathematical models for COVID-19 transmission within a regional population. One is a fitting function that can be calibrated to reproduce an epidemic curve with two timescales (e.g., fast growth and slow decay). The other is a compartmental model that accounts for quarantine, self-isolation, social distancing, a non-exponentially distributed incubation period, asymptomatic individuals, and mild and severe forms of symptomatic disease. Using Bayesian inference, we have been calibrating our models daily for consistency with new reports of confirmed cases from the 15 most populous metropolitan statistical areas in the United States and quantifying uncertainty in parameter estimates and predictions of future case reports. This online learning approach allows for early identification of new trends despite considerable variability in case reporting. We infer new significant upward trends for five of the metropolitan areas starting between 19-April-2020 and 12-June-2020.

PyBioNetFit and the Biological Property Specification Language

In systems biology modeling, important steps include model parameterization, uncertainty quantification, and evaluation of agreement with experimental observations. To help modelers perform these steps, we developed the software PyBioNetFit, which in addition supports checking models against known system properties and solving design problems. PyBioNetFit introduces Biological Property Specification Language (BPSL) for the formal declaration of system properties. BPSL allows qualitative data to be used alone or in combination with quantitative data. PyBioNetFit performs parameterization with parallelized metaheuristic optimization algorithms that work directly with existing model definition standards: BioNetGen Language (BNGL) and Systems Biology Markup Language (SBML). We demonstrate PyBioNetFit's capabilities by solving various example problems, including the challenging problem of parameterizing a 153-parameter model of cell cycle control in yeast based on both quantitative and qualitative data. We demonstrate the model checking and design applications of PyBioNetFit and BPSL by analyzing a model of targeted drug interventions in autophagy signaling.

Modeling cell line-specific recruitment of signaling proteins to the insulin-like growth factor 1 receptor

Receptor tyrosine kinases (RTKs) typically contain multiple autophosphorylation sites in their cytoplasmic domains. Once activated, these autophosphorylation sites can recruit downstream signaling proteins containing Src homology 2 (SH2) and phosphotyrosine-binding (PTB) domains, which recognize phosphotyrosine-containing short linear motifs (SLiMs). These domains and SLiMs have polyspecific or promiscuous binding activities. Thus, multiple signaling proteins may compete for binding to a common SLiM and vice versa. To investigate the effects of competition on RTK signaling, we used a rule-based modeling approach to develop and analyze models for ligand-induced recruitment of SH2/PTB domain-containing proteins to autophosphorylation sites in the insulin-like growth factor 1 (IGF1) receptor (IGF1R). Models were parameterized using published datasets reporting protein copy numbers and site-specific binding affinities. Simulations were facilitated by a novel application of model restructuration, to reduce redundancy in rule-derived equations. We compare predictions obtained via numerical simulation of the model to those obtained through simple prediction methods, such as through an analytical approximation, or ranking by copy number and/or K_D value, and find that the simple methods are unable to recapitulate the predictions of numerical simulations. We created 45 cell line-specific models that demonstrate how early events in IGF1R signaling depend on the protein abundance profile of a cell. Simulations, facilitated by model restructuration, identified pairs of IGF1R binding partners that are recruited in anti-correlated and correlated fashions, despite no inclusion of cooperativity in our models. This work shows that the outcome of competition depends on the physicochemical parameters that characterize pairwise interactions, as well as network properties, including network connectivity and the relative abundances of competitors.

Cells rely on networks of interacting biomolecules to sense and respond to environmental perturbations and signals. However, it is unclear how information is processed to generate appropriate and specific responses to signals, especially given that these networks tend to share many components. For example, receptors that detect distinct ligands and regulate distinct cellular activities commonly interact with overlapping sets of downstream signaling proteins. Here, to investigate the downstream signaling of a well-studied receptor tyrosine kinase (RTK), the insulin-like growth factor 1 (IGF1) receptor (IGF1R), we formulated and analyzed 45 cell line-specific mathematical models, which account for recruitment of 18 different binding partners to six sites of receptor autophosphorylation in IGF1R. The models were parameterized using available protein copy number and site-specific affinity measurements, and restructured to allow for network generation. We find that recruitment is influenced by the protein abundance profile of a cell, with different patterns of recruitment in different cell lines. Furthermore, in a given cell line, we find that pairs of IGF1R binding partners may be recruited in a correlated or anti-correlated fashion. We demonstrate that the simulations of the model have greater predictive power than protein copy number and/or binding affinity data, and that even a simple analytical model cannot reproduce the predicted recruitment ranking obtained via simulations. These findings represent testable predictions and indicate that the outputs of IGF1R signaling depend on cell line-specific properties in addition to the properties that are intrinsic to the biomolecules involved.

A Step-by-Step Guide to Using BioNetFit

BioNetFit is a software tool designed for solving parameter identification problems that arise in the development of rule-based models. It solves these problems through curve fitting (i.e., nonlinear regression). BioNetFit is compatible with deterministic and stochastic simulators that accept BioNetGen language (BNGL)-formatted files as inputs, such as those available within the BioNetGen framework. BioNetFit can be used on a laptop or stand-alone multicore workstation as well as on many Linux clusters, such as those that use the Slurm Workload Manager to schedule jobs. BioNetFit implements a metaheuristic population-based global optimization procedure, an evolutionary algorithm (EA), to minimize a user-defined objective function, such as a residual sum of squares (RSS) function. BioNetFit also implements a bootstrapping procedure for determining confidence intervals for parameter estimates. Here, we provide step-by-step instructions for using BioNetFit to estimate the values of parameters of a BNGL-encoded model and to define bootstrap confidence intervals. The process entails the use of several plain-text files, which are processed by BioNetFit and BioNetGen. In general, these files include (1) one or more EXP files, which each contains (experimental) data to be used in parameter identification/bootstrapping; (2) a BNGL file containing a model section, which defines a (rule-based) model, and an actions section, which defines simulation protocols that generate GDAT and/or SCAN files with model predictions corresponding to the data in the EXP file(s); and (3) a CONF file that configures the fitting/bootstrapping job and that defines algorithmic parameter settings.

New insights into RAS biology reinvigorate interest in mathematical modeling of RAS signaling

RAS is the most frequently mutated gene across human cancers, but developing inhibitors of mutant RAS has proven to be challenging. Given the difficulties of targeting RAS directly, drugs that impact the other components of pathways where mutant RAS operates may potentially be effective. However, the system-level features, including different localizations of RAS isoforms, competition between downstream effectors, and interlocking feedback and feed-forward loops, must be understood to fully grasp the opportunities and limitations of inhibiting specific targets. Mathematical modeling can help us discern the system-level impacts of these features in normal and cancer cells. New technologies enable the acquisition of experimental data that will facilitate development of realistic models of oncogenic RAS behavior. In light of the wealth of empirical data accumulated over decades of study and the advancement of experimental methods for gathering new data, modelers now have the opportunity to advance progress toward realization of targeted treatment for mutant RAS-driven cancers.

A phase 2 trial of personalized cytotoxic therapy based on tumor immunohistochemistry in previously treated metastatic pancreatic cancer patients

Background

The choice of a regimen in metastatic pancreatic cancer patients following progression on 1st line therapy is empiric and outcomes are unsatisfactory. This phase II study was performed to evaluate the efficacy of therapy selected by immunohistochemistry (IHC) in these patients following progression after one or more therapies.

Methods

Eligible patients underwent a percutaneous biopsy of a metastatic lesion and treatment selection was determined by IHC. The study required 35 evaluable patients (power of 86%) for detecting a true 1-year survival rate of >20%.

Results

A tumor biopsy was performed in 48 of 49 accrued patients. Study therapy was not given (n=13) either due to insufficient tumor on biopsy (n=8) or due to worsening cancer related symptoms after biopsy (n=5). The demographics of evaluable patients (n=35) are male/female (59%/41%), with age range 34-78 years (median 63 years). Patients had 1-6 prior regimens (median of 2). The most common IHC targets were topoisomerase 1 or 2, thymidylate synthase, excision repair cross-complementation group 1 protein (ERCC1), and osteonectin secreted protein acidic and rich in cysteine (SPARC). Commercially available treatment regimens prescribed included FOLFIRI, FOLFOX, irinotecan, and doxorubicin. The response (RECIST) was 9%, the median survival was 5.6 months (94% CI, 3.8-8.2), and the 1-year survival was 20% (95% CI, 7-33%).

Conclusions

In all patients, IHC assays resulted in identification of at least two targets for therapy and a non-cross resistant regimen could be prescribed for therapy with evidence of some benefit. An IHC based treatment strategy is feasible and needs validation in larger studies.

AutoDock VinaXB: implementation of XBSF, new empirical halogen bond scoring function, into AutoDock Vina

Background

Halogen bonding has recently come to play as a target for lead optimization in rational drug design. However, most docking program don’t account for halogen bonding in their scoring functions and are not able to utilize this new approach. In this study a new and improved halogen bonding scoring function (XBSF) is presented along with its implementation in the AutoDock Vina molecular docking software. This new improved program is termed as AutoDock VinaXB, where XB stands for the halogen bonding parameters that were added.

Results

XBSF scoring function is derived based on the X···A distance and C–X···A angle of interacting atoms. The distance term was further corrected to account for the polar flattening effect of halogens. A total of 106 protein-halogenated ligand complexes were tested and compared in terms of binding affinity and docking poses using Vina and VinaXB. VinaXB performed superior to Vina in the majority of instances. VinaXB was closer to native pose both above and below 2 Å deviation categories almost twice as frequently as Vina.

Conclusions

Implementation of XBSF into AutoDock Vina has been shown to improve the accuracy of the docking result with regards to halogenated ligands. AutoDock VinaXB addresses the issues of halogen bonds that were previously being scored unfavorably due to repulsion factors, thus effectively lowering the output RMSD values.

Electronic supplementary material

The online version of this article (doi:10.1186/s13321-016-0139-1) contains supplementary material, which is available to authorized users.

BioNetFit: a fitting tool compatible with BioNetGen, NFsim and distributed computing environments

Rule-based models are analyzed with specialized simulators, such as those provided by the BioNetGen and NFsim open-source software packages. Here, we present BioNetFit, a general-purpose fitting tool that is compatible with BioNetGen and NFsim. BioNetFit is designed to take advantage of distributed computing resources. This feature facilitates fitting (i.e. optimization of parameter values for consistency with data) when simulations are computationally expensive.

AVAILABILITY AND IMPLEMENTATION

BioNetFit can be used on stand-alone Mac, Windows/Cygwin, and Linux platforms and on Linux-based clusters running SLURM, Torque/PBS, or SGE. The BioNetFit source code (Perl) is freely available (http://bionetfit.nau.edu).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

CONTACT

bionetgen.help@gmail.com.

Use of mechanistic models to integrate and analyze multiple proteomic datasets

Proteins in cell signaling networks tend to interact promiscuously through low-affinity interactions. Consequently, evaluating the physiological importance of mapped interactions can be difficult. Attempts to do so have tended to focus on single, measurable physicochemical factors, such as affinity or abundance. For example, interaction importance has been assessed on the basis of the relative affinities of binding partners for a protein of interest, such as a receptor. However, multiple factors can be expected to simultaneously influence the recruitment of proteins to a receptor (and the potential of these proteins to contribute to receptor signaling), including affinity, abundance, and competition, which is a network property. Here, we demonstrate that measurements of protein copy numbers and binding affinities can be integrated within the framework of a mechanistic, computational model that accounts for mass action and competition. We use cell line-specific models to rank the relative importance of protein-protein interactions in the epidermal growth factor receptor (EGFR) signaling network for 11 different cell lines. Each model accounts for experimentally characterized interactions of six autophosphorylation sites in EGFR with proteins containing a Src homology 2 and/or phosphotyrosine-binding domain. We measure importance as the predicted maximal extent of recruitment of a protein to EGFR following ligand-stimulated activation of EGFR signaling. We find that interactions ranked highly by this metric include experimentally detected interactions. Proteins with high importance rank in multiple cell lines include proteins with recognized, well-characterized roles in EGFR signaling, such as GRB2 and SHC1, as well as a protein with a less well-defined role, YES1. Our results reveal potential cell line-specific differences in recruitment.

Abstract 1911: Analysis of the pancreatic cancer (PC) genome by immunohistochemistry (IHC) and comparative genome hybridization (CGH) in patients with long term survival with metastatic disease (mPC) treated on a phase II molecular profiling trial

Pancreatic cancer (PC) is a heterogeneous disease. Several mutations have been reported to develop along the path of PC evolution, and it is important to identify targets for potential therapy. CGH and IHC analysis of pancreatic tumor samples may reveal prognostic and predictive markers.

Methods: We evaluated tumor samples of patients (n=35) who were accrued to a phase II molecular profiling study (Ramanathan RK et al.12-LB-8953-AACR 2012). All subjects, with metastatic PC, with progression after 1st line of therapy had a core biopsy for CGH, IHC, and microarray analysis. Samples interrogated by CGH were flow sorted prior to analysis. Treatment recommendations were based on the IHC panel (Target Now®, Caris Life Sciences, Irving, TX).

Results: In 30 patients, adequate tumor was available for CGH. We identified 9 long term survivors, living > 18 months from first diagnosis of mPC. Of the 9 patients, 5 had available CGH data.

(1) had two distinct aneuploid tumor cell populations. There was evidence suggestive of increased WNT signaling based on amplicon targeting FZD3;. IHC identified TS and SPARC as potential targets. Recommended treatment: Capecitabine and Nab-Paclitaxel

(2) had a focal amplicon targeting WNT6 and WNT10A suggestive of increased WNT signaling, and homozygous deletions targeting RB1 and PARD3B. An intragenic break within NOTCH2 was also observed. IHC identified TS and TOPO1 as potential targets. Recommended treatment: FOLFIRI

(3) was identified (and confirmed to be) KRAS wild-type, but had a homozygous deletion of RASA1, the gene that codes for p120-RasGAP a known negative regulator of KRAS signaling. Commonly seen deletions of CDKN2A and SMAD4 were not observed. IHC identified TS and ERCC1 as potential targets. Recommended treatment: FOLFOX

(4) had an amplification of a PARP gene (PARP8) and of a base excision enzyme (NEIL2), chromosomal copy number changes to genes involved in chromatin regulation (LCMT, EZH2, KLF13) and EMT with dedifferentiation (ZNF703). IHC identified MGMT as a potential target. Recommended treatment: Temazolamide.

(5) had amplifications of PIK3CA and of KRAS, and a homozygous deletion of CDKN2A. IHC identified TOPO2 as a potential target. Recommended treatment: Doxorubicin

Conclusion: Due to the heterogeneity of the genome, even in a select group with long term survival, each patient may require an individualized treatment plan. Pathway analysis of all patients treated in the study is ongoing. Treatment strategies to focus on the stroma and of methylation of the genome are being explored to circumvent the genetic heterogeneity. (Supported by a SU2C pancreatic dream team grant)

Citation Format: Safi Shahda, Michael T. Barrett, Elizabeth Lenkiewicz, Lisa Evers, Richard G. Posner, Meraj Aziz, Glen J. Weiss, Gayle Jameson, Daniel D. Von Hoff, Ramesh K. Ramanathan. Analysis of the pancreatic cancer (PC) genome by immunohistochemistry (IHC) and comparative genome hybridization (CGH) in patients with long term survival with metastatic disease (mPC) treated on a phase II molecular profiling trial. [abstract]. In: Proceedings of the 104th Annual Meeting of the American Association for Cancer Research; 2013 Apr 6-10; Washington, DC. Philadelphia (PA): AACR; Cancer Res 2013;73(8 Suppl):Abstract nr 1911. doi:10.1158/1538-7445.AM2013-1911

Computational model for autophagic vesicle dynamics in single cells

Macroautophagy (autophagy) is a cellular recycling program essential for homeostasis and survival during cytotoxic stress. This process, which has an emerging role in disease etiology and treatment, is executed in four stages through the coordinated action of more than 30 proteins. An effective strategy for studying complicated cellular processes, such as autophagy, involves the construction and analysis of mathematical or computational models. When developed and refined from experimental knowledge, these models can be used to interrogate signaling pathways, formulate novel hypotheses about systems, and make predictions about cell signaling changes induced by specific interventions. Here, we present the development of a computational model describing autophagic vesicle dynamics in a mammalian system. We used time-resolved, live-cell microscopy to measure the synthesis and turnover of autophagic vesicles in single cells. The stochastically simulated model was consistent with data acquired during conditions of both basal and chemically-induced autophagy. The model was tested by genetic modulation of autophagic machinery and found to accurately predict vesicle dynamics observed experimentally. Furthermore, the model generated an unforeseen prediction about vesicle size that is consistent with both published findings and our experimental observations. Taken together, this model is accurate and useful and can serve as the foundation for future efforts aimed at quantitative characterization of autophagy.

Tumor-stromal interactions in pancreatic cancer

The tumor associated stroma has been described in recent years as being complicit in tumor growth in pancreatic cancer. The stroma hosts a variety of components of both cellular and molecular makeup. In normal tissues, the stroma provides nutrients and regulatory signals for proper cellular polarity and function. However, following oncogenic transformation, the stromal compartment is conscripted to provide stimulatory signals and protection to tumor cells. It is these tumor-stromal interactions that are currently of great therapeutic interest. Several key reports have suggested that therapeutic targeting of the tumor-stromal interactions in pancreatic cancer has the potential to offer survival benefit. In this review, we will discuss the tumor-stromal interactions that contribute to tumor growth and progression, and ways in which we might counter these interactions.

Targeting the tumor microenvironment in cancer: why hyaluronidase deserves a second look

Increased extracellular matrix (ECM) deposition is a characteristic observed in many solid tumors. Increased levels of one ECM component-namely, hyaluronan (HA)-leads to reduced elasticity of tumor tissue and increased interstitial fluid pressure. Multiple initial reports showed that the addition of hyaluronidase (HYAL) to chemotherapeutic regimens could greatly improve efficacy. Unfortunately, the bovine HYAL used in those studies was limited therapeutically by immunologic responses to treatment. Newly developed recombinant human HYAL has recently been introduced into clinical trials. In this article, we describe the role of HA in cancer, methods of targeting HA, and clinical studies performed to date, and we propose that targeting HA could now be an effective treatment option for patients with many different types of solid tumors.

Specification, annotation, visualization and simulation of a large rule-based model for ERBB receptor signaling

BACKGROUND

Mathematical/computational models are needed to understand cell signaling networks, which are complex. Signaling proteins contain multiple functional components and multiple sites of post-translational modification. The multiplicity of components and sites of modification ensures that interactions among signaling proteins have the potential to generate myriad protein complexes and post-translational modification states. As a result, the number of chemical species that can be populated in a cell signaling network, and hence the number of equations in an ordinary differential equation model required to capture the dynamics of these species, is prohibitively large. To overcome this problem, the rule-based modeling approach has been developed for representing interactions within signaling networks efficiently and compactly through coarse-graining of the chemical kinetics of molecular interactions.

RESULTS

Here, we provide a demonstration that the rule-based modeling approach can be used to specify and simulate a large model for ERBB receptor signaling that accounts for site-specific details of protein-protein interactions. The model is considered large because it corresponds to a reaction network containing more reactions than can be practically enumerated. The model encompasses activation of ERK and Akt, and it can be simulated using a network-free simulator, such as NFsim, to generate time courses of phosphorylation for 55 individual serine, threonine, and tyrosine residues. The model is annotated and visualized in the form of an extended contact map.

CONCLUSIONS

With the development of software that implements novel computational methods for calculating the dynamics of large-scale rule-based representations of cellular signaling networks, it is now possible to build and analyze models that include a significant fraction of the protein interactions that comprise a signaling network, with incorporation of the site-specific details of the interactions. Modeling at this level of detail is important for understanding cellular signaling.

Abstract 3697: Phase II study of therapy selected by molecular profiling in patients with previously treated metastatic pancreatic cancer - SU2C-001

We have completed a second line Phase II study of therapy to be selected by molecular profiling in 35 patients with previously treated metastatic pancreatic cancer. On this study all patients have a biopsy of a metastatic lesion (usually liver). Therapeutic targets are identified by profiling the metastatic biopsies with the combined use of immunohistochemistry (IHC), comparative genomic hybridization (CGH) on sorted tumor populations, and correlating gene expression microarray data to that of a panel of xenografts for which drug sensitivity has been previously determined. The IHC assays are performed in a CLIA certified laboratory and include the following therapeutic targets: MRP1, TOPO1, MGMT, PGP, Her2/Neu, PTEN, c-kit, RRM1, EGFR, BCRP, PDGFR, TOP2A, TS, ERCC1, SPARC, ER, PR, and AR. In addition each sample is screened by PCR for the presence of the most common K-RAS mutations. Gene expression microarray analysis is performed at John Hopkins University using Affymetrix U133 Plus 2.0 gene arrays. For CGH analysis flow sorted nuclei of diploid and aneuploid tumor cell populations are processed and hybridized to 400,000 feature CGH arrays. We have identified distinct high level focal amplicons targeting AKT2, LYN, RET, CDK6, K-RAS, MYCBP, BCL11A, RASA1 and WNT6, and gene specific homozygous deletions including CDKN2A, PTEN, MAP2K4, RASA1 and PARK2 in the sorted aneuploid tumor nuclei. In an effort to elucidate new targets and potential contexts of vulnerability based on the results of this trial, we have extracted combined gene level data for each patient and present these mapped to 32 pancreatic specific pathways and processes.

Citation Format: {Authors}. {Abstract title} [abstract]. In: Proceedings of the 103rd Annual Meeting of the American Association for Cancer Research; 2012 Mar 31-Apr 4; Chicago, IL. Philadelphia (PA): AACR; Cancer Res 2012;72(8 Suppl):Abstract nr 3697. doi:1538-7445.AM2012-3697

Guidelines for visualizing and annotating rule-based models

Rule-based modeling provides a means to represent cell signaling systems in a way that captures site-specific details of molecular interactions. For rule-based models to be more widely understood and (re)used, conventions for model visualization and annotation are needed. We have developed the concepts of an extended contact map and a model guide for illustrating and annotating rule-based models. An extended contact map represents the scope of a model by providing an illustration of each molecule, molecular component, direct physical interaction, post-translational modification, and enzyme-substrate relationship considered in a model. A map can also illustrate allosteric effects, structural relationships among molecular components, and compartmental locations of molecules. A model guide associates elements of a contact map with annotation and elements of an underlying model, which may be fully or partially specified. A guide can also serve to document the biological knowledge upon which a model is based. We provide examples of a map and guide for a published rule-based model that characterizes early events in IgE receptor (FcεRI) signaling. We also provide examples of how to visualize a variety of processes that are common in cell signaling systems but not considered in the example model, such as ubiquitination. An extended contact map and an associated guide can document knowledge of a cell signaling system in a form that is visual as well as executable. As a tool for model annotation, a map and guide can communicate the content of a model clearly and with precision, even for large models.

RuleMonkey: software for stochastic simulation of rule-based models

Background

The system-level dynamics of many molecular interactions, particularly protein-protein interactions, can be conveniently represented using reaction rules, which can be specified using model-specification languages, such as the BioNetGen language (BNGL). A set of rules implicitly defines a (bio)chemical reaction network. The reaction network implied by a set of rules is often very large, and as a result, generation of the network implied by rules tends to be computationally expensive. Moreover, the cost of many commonly used methods for simulating network dynamics is a function of network size. Together these factors have limited application of the rule-based modeling approach. Recently, several methods for simulating rule-based models have been developed that avoid the expensive step of network generation. The cost of these "network-free" simulation methods is independent of the number of reactions implied by rules. Software implementing such methods is now needed for the simulation and analysis of rule-based models of biochemical systems.

Results

Here, we present a software tool called RuleMonkey, which implements a network-free method for simulation of rule-based models that is similar to Gillespie's method. The method is suitable for rule-based models that can be encoded in BNGL, including models with rules that have global application conditions, such as rules for intramolecular association reactions. In addition, the method is rejection free, unlike other network-free methods that introduce null events, i.e., steps in the simulation procedure that do not change the state of the reaction system being simulated. We verify that RuleMonkey produces correct simulation results, and we compare its performance against DYNSTOC, another BNGL-compliant tool for network-free simulation of rule-based models. We also compare RuleMonkey against problem-specific codes implementing network-free simulation methods.

Conclusions

RuleMonkey enables the simulation of rule-based models for which the underlying reaction networks are large. It is typically faster than DYNSTOC for benchmark problems that we have examined. RuleMonkey is freely available as a stand-alone application http://public.tgen.org/rulemonkey. It is also available as a simulation engine within GetBonNie, a web-based environment for building, analyzing and sharing rule-based models.

Modeling multivalent ligand-receptor interactions with steric constraints on configurations of cell-surface receptor aggregates

We use flow cytometry to characterize equilibrium binding of a fluorophore-labeled trivalent model antigen to bivalent IgE-FcepsilonRI complexes on RBL cells. We find that flow cytometric measurements are consistent with an equilibrium model for ligand-receptor binding in which binding sites are assumed to be equivalent and ligand-induced receptor aggregates are assumed to be acyclic. However, this model predicts extensive receptor aggregation at antigen concentrations that yield strong cellular secretory responses, which is inconsistent with the expectation that large receptor aggregates should inhibit such responses. To investigate possible explanations for this discrepancy, we evaluate four rule-based models for interaction of a trivalent ligand with a bivalent cell-surface receptor that relax simplifying assumptions of the equilibrium model. These models are simulated using a rule-based kinetic Monte Carlo approach to investigate the kinetics of ligand-induced receptor aggregation and to study how the kinetics and equilibria of ligand-receptor interaction are affected by steric constraints on receptor aggregate configurations and by the formation of cyclic receptor aggregates. The results suggest that formation of linear chains of cyclic receptor dimers may be important for generating secretory signals. Steric effects that limit receptor aggregation and transient formation of small receptor aggregates may also be important.

Simulation of large-scale rule-based models

MOTIVATION

Interactions of molecules, such as signaling proteins, with multiple binding sites and/or multiple sites of post-translational covalent modification can be modeled using reaction rules. Rules comprehensively, but implicitly, define the individual chemical species and reactions that molecular interactions can potentially generate. Although rules can be automatically processed to define a biochemical reaction network, the network implied by a set of rules is often too large to generate completely or to simulate using conventional procedures. To address this problem, we present DYNSTOC, a general-purpose tool for simulating rule-based models.

RESULTS

DYNSTOC implements a null-event algorithm for simulating chemical reactions in a homogenous reaction compartment. The simulation method does not require that a reaction network be specified explicitly in advance, but rather takes advantage of the availability of the reaction rules in a rule-based specification of a network to determine if a randomly selected set of molecular components participates in a reaction during a time step. DYNSTOC reads reaction rules written in the BioNetGen language which is useful for modeling protein-protein interactions involved in signal transduction. The method of DYNSTOC is closely related to that of StochSim. DYNSTOC differs from StochSim by allowing for model specification in terms of BNGL, which extends the range of protein complexes that can be considered in a model. DYNSTOC enables the simulation of rule-based models that cannot be simulated by conventional methods. We demonstrate the ability of DYNSTOC to simulate models accounting for multisite phosphorylation and multivalent binding processes that are characterized by large numbers of reactions.

AVAILABILITY

DYNSTOC is free for non-commercial use. The C source code, supporting documentation and example input files are available at http://public.tgen.org/dynstoc/.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

GetBonNie for building, analyzing and sharing rule-based models

SUMMARY

GetBonNie is a web-based application for building, analyzing and sharing rule-based models encoded in the BioNetGen language (BNGL). Tools accessible within the GetBonNie environment include (i) an applet for drawing graphs that correspond to BNGL code; (ii) a network-generation engine for translating a set of rules into a chemical reaction network; (iii) simulation engines that implement generate-first, on-the-fly and network-free methods for simulating rule-based models; and (iv) a database for sharing models, parameter values, annotations, simulation tasks and results.

AVAILABILITY

GetBonNie is free at (http://getbonnie.org).

Predominance of activated EGFR higher-order oligomers on the cell surface

The epidermal growth factor receptor (EGFR) kinase is generally considered to be activated by either ligand-induced dimerisation or a ligand-induced conformational change within pre-formed dimers. We report the relationship between ligand-induced higher-order EGFR oligomerization and EGFR phosphorylation on the surface of intact cells. We have combined lifetime-detected Forster resonance energy transfer, as a probe of the receptor phosphorylation state and image correlation spectroscopy, to extract the relative association state of activated versus unactivated EGFR, to determine the ratio of the average number of receptors for active (phosphorylated) and inactive clusters. There are at least four times as many receptors in the ligand-induced active clusters than inactive clusters. Contrary to the prevailing view that the EGFR dimer is the predominant, active form, our data determine that higher-order EGFR oligomers are the dominant species associated with the ligand activated EGFR tyrosine kinase.

Trivalent antigens for degranulation of mast cells

Degranulation of basophils and mast cells plays a central role in allergic reactions. Degranulation is a response to cell surface receptor aggregation caused by association of receptors with antibodies bound to multivalent antigens. Tools used in studying this process have included small-molecule divalent antigens, but they suffer from weak signaling apparently due to small aggregate size. We have prepared trivalent antigens that allow formation of larger aggregates and potent responses from mast cells.

Rules for modeling signal-transduction systems

Formalized rules for protein-protein interactions have recently been introduced to represent the binding and enzymatic activities of proteins in cellular signaling. Rules encode an understanding of how a system works in terms of the biomolecules in the system and their possible states and interactions. A set of rules can be as easy to read as a diagrammatic interaction map, but unlike most such maps, rules have precise interpretations. Rules can be processed to automatically generate a mathematical or computational model for a system, which enables explanatory and predictive insights into the system's behavior. Rules are independent units of a model specification that facilitate model revision. Instead of changing a large number of equations or lines of code, as may be required in the case of a conventional mathematical model, a protein interaction can be introduced or modified simply by adding or changing a single rule that represents the interaction of interest. Rules can be defined and visualized by using graphs, so no specialized training in mathematics or computer science is necessary to create models or to take advantage of the representational precision of rules. Rules can be encoded in a machine-readable format to enable electronic storage and exchange of models, as well as basic knowledge about protein-protein interactions. Here, we review the motivation for rule-based modeling; applications of the approach; and issues that arise in model specification, simulation, and testing. We also discuss rule visualization and exchange and the software available for rule-based modeling.

Glutathione-S-transferase-green fluorescent protein fusion protein reveals slow dissociation from high site density beads and measures free GSH

BACKGROUND

Glutathione, a ubiquitous tripeptide, is an important cellular constituent, and measurement of reduced and oxidized glutathione is a measure of the redox state of cells. Glutathione-S-transferase (GST) fusion proteins bind naturally to beads derivatized with glutathione, and elution of such bead-bound fusion proteins with buffer containing millimolar glutathione is a commonly used method of protein purification. Many protein-protein interactions have been established by using GST fusion proteins and measuring binding of fusion protein binding partners by GST pulldown assays, usually monitored by Western blot methodology.

METHODS

Dextran beads suitable for flow cytometry were derivatized with glutathione. A fusion protein of GST and green fluorescent protein was used to define kinetic and equilibrium binding characteristics of GST fusion proteins to glutathione beads. Free glutathione competes with this binding, and this competition was used to measure free glutathione concentration.

RESULTS

A 10 microl assay can measure 5 microl of 20 microM glutathione (100 pmol glutathione) in 2 h by flow cytometry. This concentration is two orders of magnitude lower than cellular glutathione concentrations, and three orders of magnitude lower than affinity chromatography eluates. One important result is that by generating high site density, the GST fusion proteins can be constrained to the surface of one bead without hopping to the next bead in multiplex assays.

CONCLUSIONS

Glutathione in cellular lysates and GST-fusion protein affinity chromatography eluates can be measured by flow cytometry. Many interactions between GST fusion proteins and their fluorescent binding partners should be quantifiable by flow cytometry. Although a system may have the disadvantage that it has a low affinity and a correspondingly quick off-rate in solution, it may remain on beads if the site density can be increased to offer a slow apparent off rate.

CO2 transduction mechanisms in avian intrapulmonary chemoreceptors: experiments and models

Intrapulmonary chemoreceptors (IPC) are neurons that sense tonic and phasic CO2 stimuli in the lungs of birds and diapsid reptiles. IPC are different from most other vertebrate respiratory CO2 receptors because: (1) they are stimulated by low PCO2 and inhibited by high PCO2, (2) they have extremely rapid response characteristics, (3) their CO2 sensitivity is nearly abolished by intracellular inhibitors of carbonic anhydrase, and (4) their CO2 sensitivity is strongly depressed by inhibiting Na+/H+ antiport exchange. Experimental evidence suggests that IPC respond to intracellular pH, not CO2 directly, and that intracellular pH and IPC discharge are determined by a kinetic balance between CO2 hydration/dehydration rates, transmembrane acid/base exchange rates, and intracellular buffering. We review experimental evidence for and against various mechanisms of IPC CO2 chemotransduction, present a conceptual and mathematical model of the proposed mechanisms, and compare this model to CO2 transduction in other respiratory chemoreceptors.

Interaction of a monoclonal IgE-specific antibody with cell-surface IgE-Fc epsilon RI: characterization of equilibrium binding and secretory response

Aggregation of FcepsilonRI, the high-affinity cell-surface receptor for IgE antibody, is required for degranulation of basophils and mast cells, but not all receptor aggregates elicit this cellular response. The stereochemical constraints on clusters of FcepsilonRI that are able to signal cellular responses, such as degranulation, have yet to be fully defined. To improve our understanding of the properties of FcepsilonRI aggregates that influence receptor signaling, we have studied the interaction of 23G3, a rat IgG(1)(kappa) IgE-specific monoclonal antibody, with IgE-FcepsilonRI complexes on rat mucosal-type mast cells (RBL-2H3 line). We find that 23G3 is a potent secretagogue. This property and the structural features of 23G3 (two symmetrically arrayed IgE-specific binding sites) make 23G3 a potentially valuable reagent for investigating the relationship between FcepsilonRI clustering and FcepsilonRI-mediated signaling events. To develop a mathematical model of 23G3-induced aggregation of FcepsilonRI, we used fluorimetry and flow cytometry to quantitatively monitor equilibrium binding of FITC-labeled 23G3 intact Ab and its Fab' fragment to cell-surface IgE. The results indicate that IgE bound to FcepsilonRI expresses two epitopes for 23G3 binding; that 23G3 binds IgE resident on the cell surface with negative cooperativity; and that 23G3 appears to induce mostly but not exclusively noncyclic dimeric aggregates of FcepsilonRI. There is no simple relationship between receptor aggregation at equilibrium and the degranulation response. Further studies are needed to establish how 23G3-induced aggregation of IgE-FcepsilonRI correlates with cellular responses.

A quantitative approach for studying IgE-FcepsilonRI aggregation

Aggregation of cell surface receptors is a ubiquitous means of initiating signal transduction in many cellular systems. In this manuscript, we describe a combined theoretical and experimental approach based on multiparameter flow cytometry for measuring the time course of ligand induced aggregation of IgE-FcepsilonRI on RBL cells. By fluorescently labeling both the ligand and surface IgE (sIgE), we have developed an assay that permits us to simultaneously measure both occupancy of sIgE combining sites and association of antigen with the cell surface. This allows for a direct calculation of the degree of receptor aggregation present on the cell. By employing new mixing technologies developed for flow cytometry, we are able to look at aggregation in the sub second time domain. To extend our work, we have synthesized a new set of chemically well defined ligands (of valences 1-3) to use as probes in our studies. We show that the magnitude of the cellular response is dramatically increased as the valence of our ligand is raised from two to three.

Steric effects on multivalent ligand-receptor binding: exclusion of ligand sites by bound cell surface receptors

Steric effects can influence the binding of a cell surface receptor to a multivalent ligand. To account for steric effects arising from the size of a receptor and from the spacing of binding sites on a ligand, we extend a standard mathematical model for ligand-receptor interactions by introducing a steric hindrance factor. This factor gives the fraction of unbound ligand sites that are accessible to receptors, and thus available for binding, as a function of ligand site occupancy. We derive expressions for the steric hindrance factor for various cases in which the receptor covers a compact region on the ligand surface and the ligand expresses sites that are distributed regularly or randomly in one or two dimensions. These expressions are relevant for ligands such as linear polymers, proteins, and viruses. We also present numerical algorithms that can be used to calculate steric hindrance factors for other cases. These theoretical results allow us to quantify the effects of steric hindrance on ligand-receptor kinetics and equilibria.

Quantifying aggregation of IgE-FcepsilonRI by multivalent antigen

Aggregation of cell surface receptors by multivalent ligand can trigger a variety of cellular responses. A well-studied receptor that responds to aggregation is the high affinity receptor for IgE (FcepsilonRI), which is responsible for initiating allergic reactions. To quantify antigen-induced aggregation of IgE-FcepsilonRI complexes, we have developed a method based on multiparameter flow cytometry to monitor both occupancy of surface IgE combining sites and association of antigen with the cell surface. The number of bound IgE combining sites in excess of the number of bound antigens, the number of bridges between receptors, provides a quantitative measure of IgE-FcepsilonRI aggregation. We demonstrate our method by using it to study the equilibrium binding of a haptenated fluorescent protein, 2,4-dinitrophenol-coupled B-phycoerythrin (DNP25-PE), to fluorescein isothiocyanate-labeled anti-DNP IgE on the surface of rat basophilic leukemia cells. The results, which we analyze with the aid of a mathematical model, indicate how IgE-FcepsilonRI aggregation depends on the total concentrations of DNP25-PE and surface IgE. As expected, we find that maximal aggregation occurs at an optimal antigen concentration. We also find that aggregation varies qualitatively with the total concentration of surface IgE as predicted by an earlier theoretical analysis.

Kinetics of adhesion of IgE-sensitized rat basophilic leukemia cells to surface-immobilized antigen in Couette flow

Antigen-antibody systems provide the flexibility of varying the kinetics and affinity of molecular interaction and studying the resulting effect on adhesion. In a parallel-plate flow chamber, we measured the extent and rate of adhesion of rat basophilic leukemia cells preincubated with anti-dinitrophenyl IgE clones SPE-7 or H1 26. 82 to dinitrophenyl-coated polyacrylamide gel substrates in a linear shear field. Both of these IgEs bind dinitrophenyl, but H1 26.82 has a 10-fold greater on rate and a 30-fold greater affinity. Adhesion was found to be binary; cells either arrested irreversibly or continued at their unencumbered hydrodynamic velocity. Under identical conditions, more adhesion was seen with the higher affinity (higher on rate) IgE clone. At some shear rates, adhesion was robust with H1 26.82, but negligible with SPE-7. Reduction in receptor number or ligand density reduced the maximum level of adhesion seen at any shear rate, but did not decrease the shear rate at which adhesion was first observed. The spatial pattern of adhesion for both IgE clones is well represented by the first-order kinetic rate constant kad, and we have determined how kad depends on ligand and receptor densities and shear rate. The rate constant kad found with H1 26.82 was approximately fivefold greater than with SPE-7. The dependence of kad on site density and shear rate for SPE-7 is complex: kad increases linearly with antigen site density at low to moderate shear rates, but is insensitive to site density at high shear. kad increases with shear rate at low site density but decreases with shear at high site density. With H1 26.82, the functional dependence of kad with shear rate was similar. Although these data are consistent with the hypothesis that we have sampled both transport and reaction-limited adhesion regimes, they point out deficiencies in current theories describing cell attachment under flow.

A rapid mix flow cytometer with subsecond kinetic resolution

Kinetic approaches are valuable tools for mechanistic studies of cell function. Flow cytometry is well suited to make sensitive kinetic measurements, but the time required to deliver mixed samples to the point of measurement (10-20 s in a conventional cytometer) limits analysis of rapidly occurring events. To address this limitation, we adapted a syringe-based stopped-flow rapid mixing device to a modified commercial flow cytometer to achieve mixing and measurement of sample in under 1 s. Because such screw-driven mixers are designed to deliver fluid at rates of microliters per millisecond and cytometers accept samples at microliters per second, the syringe mixer was modified with a screw to allow sample delivery at rates as low as 1.8 microliters/s. A custom-made nozzle holder featuring a fast-acting three-way sample delivery valve and a 1.5- microliters dead volume was designed for a Becton Dickinson FACS stream-in-air flow nozzle. Syringe motors and valves are computer controlled, as is the start signal for an adjustable time ramp. A stable sample stream can be established within the sheath stream in less than 1 s, enabling fluorescence measurements of microspheres with coefficients of variation of approximately 5%. Light scatter gating to select particles in the center of the laser beam enables fluorescence measurements at times of under 300 ms. Efficient mixing of reagents is demonstrated by the iodide quenching of microspheres surface labeled with fluorescein isothiocyanate (FITC). The instrument is capable of quantitatively proportioning cells and reagent, thereby allowing precise control of reagent concentration and dilution. Rapid kinetic measurements of intact cells are demonstrated by FITC-formyl peptide binding to cell surface receptors.

Simultaneous cross-linking by two nontriggering bivalent ligands causes synergistic signaling of IgE Fc epsilon RI complexes

We have used two bivalent ligands that bind IgE to study the relationship between the aggregation of receptors with high affinity for IgE (Fc epsilon RI) and the responses (receptor immobilization, Ca2+ influx, and degranulation) of rat basophilic leukemia (RBL-2H3) cells. One of these is a symmetric bivalent ligand, N,N'-bis[[epsilon-[(2,4-dinitrophenyl)amino]caproyl]-L-tyrosyl]-L- cystine ((DCT)2-cys), which binds specifically to the combining sites of a mAb anti-DNP IgE and efficiently cross-links cell surface IgE, but does not trigger significant degranulation or increases in intracellular Ca2+. Several lines of evidence, including lateral mobility measurements, indicate that this ligand preferentially forms stable cyclic complexes containing two (DCT)2-cys and two IgE. The second ligand is a mAb anti-IgE, B1E3, which causes lateral mobility changes consistent with dimerized IgE-Fc epsilon RI and also does not trigger increases in intracellular Ca2+ or degranulation. The two ligands together trigger robust responses. In the presence of B1E3, (DCT)2-cys causes immobilization of IgE-Fc epsilon RI in a broad concentration range; in a more narrow concentration range, it is a potent stimulant of changes in both degranulation and Ca2+. We have compared the dose-response curves for cellular activation to simulated IgE aggregation curves, i.e., curves that predict the equilibrium IgE aggregate size distribution as a function of the (DCT)2-cys concentration. Our results indicate that maximal cellular activation occurs at a much higher (DCT)2-cys concentration than maximal IgE aggregation. When IgE aggregation is maximal, almost all aggregated IgE is in cyclic dimers. Thus, cyclic dimers appear to be functionally ineffective, even after they have been cross-linked by B1E3. Aggregated IgE-Fc epsilon RI that is effective in stimulating a cellular response may have particular structural or dynamic properties that allow critical interactions for initiating the signaling cascade.

Simultaneous cross-linking by two nontriggering bivalent ligands causes synergistic signaling of IgE Fc epsilon RI complexes

We have used two bivalent ligands that bind IgE to study the relationship between the aggregation of receptors with high affinity for IgE (Fc epsilon RI) and the responses (receptor immobilization, Ca2+ influx, and degranulation) of rat basophilic leukemia (RBL-2H3) cells. One of these is a symmetric bivalent ligand, N,N'-bis[[epsilon-[(2,4-dinitrophenyl)amino]caproyl]-L-tyrosyl]-L- cystine ((DCT)2-cys), which binds specifically to the combining sites of a mAb anti-DNP IgE and efficiently cross-links cell surface IgE, but does not trigger significant degranulation or increases in intracellular Ca2+. Several lines of evidence, including lateral mobility measurements, indicate that this ligand preferentially forms stable cyclic complexes containing two (DCT)2-cys and two IgE. The second ligand is a mAb anti-IgE, B1E3, which causes lateral mobility changes consistent with dimerized IgE-Fc epsilon RI and also does not trigger increases in intracellular Ca2+ or degranulation. The two ligands together trigger robust responses. In the presence of B1E3, (DCT)2-cys causes immobilization of IgE-Fc epsilon RI in a broad concentration range; in a more narrow concentration range, it is a potent stimulant of changes in both degranulation and Ca2+. We have compared the dose-response curves for cellular activation to simulated IgE aggregation curves, i.e., curves that predict the equilibrium IgE aggregate size distribution as a function of the (DCT)2-cys concentration. Our results indicate that maximal cellular activation occurs at a much higher (DCT)2-cys concentration than maximal IgE aggregation. When IgE aggregation is maximal, almost all aggregated IgE is in cyclic dimers. Thus, cyclic dimers appear to be functionally ineffective, even after they have been cross-linked by B1E3. Aggregated IgE-Fc epsilon RI that is effective in stimulating a cellular response may have particular structural or dynamic properties that allow critical interactions for initiating the signaling cascade.

The kinetics of bivalent ligand-bivalent receptor aggregation: ring formation and the breakdown of the equivalent site approximation

When bivalent ligands capable of bridging binding sites on two different receptors interact with bivalent receptors, aggregates form. The aggregates can be of two types: chains (open structures containing n receptors, n-1 doubly bound ligands and 0, 1, or 2 singly bound ligands) and rings (closed structures containing n receptors and n doubly bound ligands). Both types of aggregates have been detected experimentally. In general, to determine the time dependence of the concentration of any particular aggregate requires solving an infinite set of coupled ordinary differential equations (ODEs). Perelson and DeLisi [19] showed that great simplification results if all receptor binding sites are equivalent, i.e., the binding properties of a site on a receptor are independent of the size of the aggregate the receptor is in. If only chains form, the problem reduces to solving two coupled ODEs for the concentrations of singly and doubly bound ligands. From the solutions to these ODEs, the time dependence of the entire aggregate size distribution can be determined. We show that the equivalent site approximation as formulated by Perelson and DeLisi [19] is incompatible with ring formation. We then present a modified equivalent site approximation that is useful if chains of any size can form but rings above a certain size (k) cannot. We show how to reduce the resulting infinite set of coupled ODEs to a closed system of at most 4k + 2 ODEs for the ligand concentrations, the ring concentrations, and the concentrations of all chains up to size k. Although we can only predict the kinetics of aggregate formation for aggregates of size k or less, at equilibrium the modified equivalent site approximation yields the complete aggregate size distribution.

Flow cytometric analysis of transmembrane phospholipid movement in bull sperm

Fluorescent phospholipids are useful to investigate phospholipid dynamics in biological membranes. We used flow cytometry to investigate transbilayer phospholipid movement in live sperm cells. Acyl-labeled N-4-nitrobenzo-2-oxa-1,3-diazole (NBD) -phosphatidylcholine (-PC), -phosphatidylethanolamine (-PE), or -phosphatidylserine (-PS) were incorporated into sperm cells, and the transbilayer location was determined by extraction of probe from cell with excess bovine serum albumin (BSA) or by chemical destruction of probe by sodium dithionite. Using these methods, we have measured the head group specific outer leaflet to inner leaflet movement (flip) of the aminophospholipids NBD-PS and NBD-PE. The fluorescent phospholipids moved inward across the plasma membrane with half-times of 1.8, 2.5, and 11.2 min, for NBD-PS, NBD-PE, and NBD-PC and reached apparent equilibrium levels of 88%, 94%, and 32% inside, respectively. The inward movement of NBD-PE was inhibited by sulfhydryl reagents, elevated intracellular Ca2+, and depletion of cellular ATP. Analysis of the kinetics of NBD-PE and -PS extraction by BSA allows determination of the rates for outward movement (flop) across the plasma membrane. Half-times for flop were 4.7 and 4.5 min for NBD-PS and -PE, respectively. Based on these measurements, a simple model of NBD-phospholipid equilibria was developed and fit to the kinetic data. Computer-generated fits reflected major features of the experimental data and provide a potential tool for predicting the dynamics of endogenous lipids.

Binding of bivalent ligand to cell surface IgE: can one detect ring formation?

It is well established that aggregation of cell surface immunoglobulin is involved in signal transduction by cells of the immune system. It is less well understood what special properties of these cell surface aggregates are important in initiating the signal cascade. Several authors have proposed that cells respond to the size (Fewtrell and Metzger (1980) J. Immun. 125, 701-710) as well as the stereochemistry (Ortega et al. (1989) Eur. J. Immun. 19, 2251-2256) of receptor aggregates. One approach to arriving at data relevant to this question has been to construct simple bivalent ligands that can bind to surface immunoglobulin. Several authors have suggested that when these bivalent ligands interact with surface immunoglobulin the formation of small stable cyclic complexes is highly favored. In this paper we consider whether it is possible to completely determine the parameters that describe the binding of a bivalent ligand to a bivalent receptor with the available experimental technology. We show that with the appropriate analysis procedure, using a modified equivalent site model, these parameters can be reliably determined from only three experiments even when there is a large amount of ring formation.

Multiparameter flow cytometric analysis of a pH sensitive formyl peptide with application to receptor structure and processing kinetics

Environmentally sensitive molecules have many potential cellular applications. We have investigated the utility of a pH sensitive ligand for the formyl peptide receptor, CHO-Met-Leu-Phe-Phe-Lys (SNAFL)-OH (SNAFL-seminaphtho-fluorescein), because in previous studies (Fay et al.: Biochemistry 30:5066-5075, 1991) protonation has been used to explain the quenching when the fluoresceinated formyl pentapeptide ligand binds to this receptor. Moreover, acidification in intracellular compartments is a general mechanism occurring in cells during processing of ligand-receptor complexes. Because the protonated form of SNAFL is excited at 488 nm with emission at 530 nm and the unprotonated form is excited at 568 nm with emission at 650 nm, the ratio of protonated and unprotonated forms can be examined by multiparameter flow cytometry. We found that the receptor-bound ligand is sensitive to both the extracellular and intracellular pH. There is a small increase in the pKa of the ligand upon binding to the receptor consistent with protonation in the binding pocket. Once internalized, spectral changes in the probe consistent with acidification and ligand dissociation from the receptor are observed.

Continuous spectrofluorometric analysis of formyl peptide receptor ternary complex interactions

Fluorescent formyl peptides have made it possible to study ligand-receptor-G protein (ternary complex) dynamics in real-time, but limitations to sample mixing and delivery in flow cytometry have interfered with continuous observation. We have taken advantage of the quenching of a fluoresceinated N-formyl pentapeptide upon binding to its receptor on permeabilized neutrophils to extend the analysis of the ternary complex dynamics to the second time scale. The association and dissociation of ligand in the presence and absence of saturating concentrations of GTP[S] were examined continuously and the results were found to be in agreement with results predicted previously from flow cytometry. We observe comparable initial rates for the formation of ligand-receptor (LR) binary complexes and ligand-receptor guanine nucleotide binding protein (LRG) ternary complexes, dissociation rates differing by two orders of magnitude, and slow interconversions between LR and LRG in the absence of guanine nucleotide. When fit by the ternary complex model, at least three sides of the model are required and the fit is improved if a significant fraction of receptors (RG) are allowed to be precoupled to G protein. One of the limitations of the analysis is that data fits are insensitive to additional parameters in the calculation which would permit analysis of all four sides of the ternary complex model. Experiments performed with subsaturating GTP[S] identified coexisting classes of LR and LRG and allowed analysis of the altered distribution between coupled and uncoupled receptors. At saturating nucleotide levels, the binding of GTP[S] and the breakup of the ternary complex occur on a subsecond time frame. This result is consistent with the idea that inside a neutrophil where GTP levels are several hundred microM, once ternary complex forms, ternary complex decomposition is rapid. Taken together, the observed rapid assembly and disassembly of ternary complex account for subsecond cell responses to ligand.

Aggregation of IgE-receptor complexes on rat basophilic leukemia cells does not change the intrinsic affinity but can alter the kinetics of the ligand-IgE interaction

The aggregation of IgE anchored to high-affinity Fc epsilon receptors on rat basophilic leukemia (RBL) cells by multivalent antigens initiates transmembrane signaling and ultimately cellular degranulation. Previous studies have shown that the rate of dissociation of bivalent and multivalent DNP ligands from RBL cells sensitized with anti-DNP IgE decreases with increasing ligand incubation times. One mechanism proposed for this effect is that when IgE molecules are aggregated, a conformational change occurs that results in an increase in the intrinsic affinity of IgE for antigen. This possibility was tested by measuring the equilibrium constant for the binding of monovalent DNP-lysine to anti-DNP IgE under two conditions, where the cell-bound IgE is dispersed and where it has been aggregated into visible patches on the cell surface using anti-IgE and a secondary antibody. No difference in the equilibrium constant in these two cases was observed. We also measured the rate of dissociation of a monovalent ligand from cell surface IgE under these two conditions. Whereas the affinity for monovalent ligand is not altered by IgE aggregation, we observe that the rate of ligand dissociation from IgE in clusters is slower than the rate of ligand dissociation from unaggregated IgE. These results are discussed in terms of recent theoretical developments concerning effects of receptor density on ligand binding to cell surfaces.

Real-time analysis of the assembly of ligand, receptor, and G protein by quantitative fluorescence flow cytometry

We describe a general approach for the quantitative analysis of the interaction among fluorescent peptide ligands (L), receptors (R), and G proteins (G) using fluorescence flow cytometry. The scheme depends upon the use of commercially available fluorescent microbeads as standards to calibrate the concentration of fluorescent peptides in solution and the receptor number on cells in suspension. We have characterized a family of fluoresceinated formyl peptides and analyzed both steady-state and dynamic aspects of ligand formyl peptide-receptor interactions in digitonin-permeabilized human neutrophils. Detailed receptor-binding studies were performed with the pentapeptide N-formyl-Met-Leu-Phe-Phe-Lys-fluorescein. Equilibrium studies showed that GTP [S] caused a loss of binding affinity of approximately two orders of magnitude, from approximately 0.04 nM (LRG) to approximately 3 nM (LR), respectively. Kinetic studies revealed that this change in affinity was principally due to an increase in the dissociation rate constants from approximately 1 x 10(-3) s-1 (LRG) to approximately 1 x 10(-1) s-1 (LR). In contrast, the association rate constants in the presence and absence of guanine nucleotide (approximately 3 x 10(7) s-1 M-1) were statistically indistinguishable and close to the diffusion limit. In the presence of guanine nucleotide (LR), the kinetic data were adequately fit by a single-step reversible-binding model. In the absence of guanine nucleotides, not all receptors have rapid access to G to form the LRG ternary complex. Mathematically, those R that have rapid access to G are either precoupled to R or the association of G with R is fast compared to the association of L with R. The physiological consequences of coupling heterogeneity are discussed.