Classic and contemporary approaches to modeling biochemical reactions.

Recent interest in modeling biochemical networks raises questions about the relationship between often complex mathematical models and familiar arithmetic concepts from classical enzymology, and also about connections between modeling and experimental data. This review addresses both topics by familiarizing readers with key concepts (and terminology) in the construction, validation, and application of deterministic biochemical models, with particular emphasis on a simple enzyme-catalyzed reaction. Networks of coupled ordinary differential equations (ODEs) are the natural language for describing enzyme kinetics in a mass action approximation. We illustrate this point by showing how the familiar Briggs-Haldane formulation of Michaelis-Menten kinetics derives from the outer (or quasi-steady-state) solution of a dynamical system of ODEs describing a simple reaction under special conditions. We discuss how parameters in the Michaelis-Menten approximation and in the underlying ODE network can be estimated from experimental data, with a special emphasis on the origins of uncertainty. Finally, we extrapolate from a simple reaction to complex models of multiprotein biochemical networks. The concepts described in this review, hitherto of interest primarily to practitioners, are likely to become important for a much broader community of cellular and molecular biologists attempting to understand the promise and challenges of "systems biology" as applied to biochemical mechanisms.

Properties of cell death models calibrated and compared using Bayesian approaches.

Using models to simulate and analyze biological networks requires principled approaches to parameter estimation and model discrimination. We use Bayesian and Monte Carlo methods to recover the full probability distributions of free parameters (initial protein concentrations and rate constants) for mass-action models of receptor-mediated cell death. The width of the individual parameter distributions is largely determined by non-identifiability but covariation among parameters, even those that are poorly determined, encodes essential information. Knowledge of joint parameter distributions makes it possible to compute the uncertainty of model-based predictions whereas ignoring it (e.g., by treating parameters as a simple list of values and variances) yields nonsensical predictions. Computing the Bayes factor from joint distributions yields the odds ratio (∼20-fold) for competing 'direct' and 'indirect' apoptosis models having different numbers of parameters. Our results illustrate how Bayesian approaches to model calibration and discrimination combined with single-cell data represent a generally useful and rigorous approach to discriminate between competing hypotheses in the face of parametric and topological uncertainty.

Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis.

Stochastic fluctuations in gene expression give rise to cell-to-cell variability in protein levels which can potentially cause variability in cellular phenotype. For TRAIL (TNF-related apoptosis-inducing ligand) variability manifests itself as dramatic differences in the time between ligand exposure and the sudden activation of the effector caspases that kill cells. However, the contribution of individual proteins to phenotypic variability has not been explored in detail. In this paper we use feature-based sensitivity analysis as a means to estimate the impact of variation in key apoptosis regulators on variability in the dynamics of cell death. We use Monte Carlo sampling from measured protein concentration distributions in combination with a previously validated ordinary differential equation model of apoptosis to simulate the dynamics of receptor-mediated apoptosis. We find that variation in the concentrations of some proteins matters much more than variation in others and that precisely which proteins matter depends both on the concentrations of other proteins and on whether correlations in protein levels are taken into account. A prediction from simulation that we confirm experimentally is that variability in fate is sensitive to even small increases in the levels of Bcl-2. We also show that sensitivity to Bcl-2 levels is itself sensitive to the levels of interacting proteins. The contextual dependency is implicit in the mathematical formulation of sensitivity, but our data show that it is also important for biologically relevant parameter values. Our work provides a conceptual and practical means to study and understand the impact of cell-to-cell variability in protein expression levels on cell fate using deterministic models and sampling from parameter distributions.

Input-output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data.

The ErbB signaling pathways, which regulate diverse physiological responses such as cell survival, proliferation and motility, have been subjected to extensive molecular analysis. Nonetheless, it remains poorly understood how different ligands induce different responses and how this is affected by oncogenic mutations. To quantify signal flow through ErbB-activated pathways we have constructed, trained and analyzed a mass action model of immediate-early signaling involving ErbB1-4 receptors (EGFR, HER2/Neu2, ErbB3 and ErbB4), and the MAPK and PI3K/Akt cascades. We find that parameter sensitivity is strongly dependent on the feature (e.g. ERK or Akt activation) or condition (e.g. EGF or heregulin stimulation) under examination and that this context dependence is informative with respect to mechanisms of signal propagation. Modeling predicts log-linear amplification so that significant ERK and Akt activation is observed at ligand concentrations far below the K(d) for receptor binding. However, MAPK and Akt modules isolated from the ErbB model continue to exhibit switch-like responses. Thus, key system-wide features of ErbB signaling arise from nonlinear interaction among signaling elements, the properties of which appear quite different in context and in isolation.

Energetics of protein-DNA interactions.

Protein-DNA interactions are vital for many processes in living cells, especially transcriptional regulation and DNA modification. To further our understanding of these important processes on the microscopic level, it is necessary that theoretical models describe the macromolecular interaction energetics accurately. While several methods have been proposed, there has not been a careful comparison of how well the different methods are able to predict biologically important quantities such as the correct DNA binding sequence, total binding free energy and free energy changes caused by DNA mutation. In addition to carrying out the comparison, we present two important theoretical models developed initially in protein folding that have not yet been tried on protein-DNA interactions. In the process, we find that the results of these knowledge-based potentials show a strong dependence on the interaction distance and the derivation method. Finally, we present a knowledge-based potential that gives comparable or superior results to the best of the other methods, including the molecular mechanics force field AMBER99.

All-atom ab initio folding of a diverse set of proteins.

Natural proteins fold to a unique, thermodynamically dominant state. Modeling of the folding process and prediction of the native fold of proteins are two major unsolved problems in biophysics. Here, we show successful all-atom ab initio folding of a representative diverse set of proteins by using a minimalist transferable-energy model that consists of two-body atom-atom interactions, hydrogen bonding, and a local sequence-energy term that models sequence-specific chain stiffness. Starting from a random coil, the native-like structure was observed during replica exchange Monte Carlo (REMC) simulation for most proteins regardless of their structural classes; the lowest energy structure was close to native-in the range of 2-6 A root-mean-square deviation (rmsd). Our results demonstrate that the successful folding of a protein chain to its native state is governed by only a few crucial energetic terms.

A Systems Toxicology Approach for the Prediction of Kidney Toxicity and Its Mechanisms In Vitro.

The failure to predict kidney toxicity of new chemical entities early in the development process before they reach humans remains a critical issue. Here, we used primary human kidney cells and applied a systems biology approach that combines multidimensional datasets and machine learning to identify biomarkers that not only predict nephrotoxic compounds but also provide hints toward their mechanism of toxicity. Gene expression and high-content imaging-derived phenotypical data from 46 diverse kidney toxicants were analyzed using Random Forest machine learning. Imaging features capturing changes in cell morphology and nucleus texture along with mRNA levels of HMOX1 and SQSTM1 were identified as the most powerful predictors of toxicity. These biomarkers were validated by their ability to accurately predict kidney toxicity of four out of six candidate therapeutics that exhibited toxicity only in late stage preclinical/clinical studies. Network analysis of similarities in toxic phenotypes was performed based on live-cell high-content image analysis at seven time points. Using compounds with known mechanism as reference, we could infer potential mechanisms of toxicity of candidate therapeutics. In summary, we report an approach to generate a multidimensional biomarker panel for mechanistic de-risking and prediction of kidney toxicity in in vitro for new therapeutic candidates and chemical entities.

A knowledge-based move set for protein folding.

The free energy landscape of protein folding is rugged, occasionally characterized by compact, intermediate states of low free energy. In computational folding, this landscape leads to trapped, compact states with incorrect secondary structure. We devised a residue-specific, protein backbone move set for efficient sampling of protein-like conformations in computational folding simulations. The move set is based on the selection of a small set of backbone dihedral angles, derived from clustering dihedral angles sampled from experimental structures. We show in both simulated annealing and replica exchange Monte Carlo (REMC) simulations that the knowledge-based move set, when compared with a conventional move set, shows statistically significant improved ability at overcoming kinetic barriers, reaching deeper energy minima, and achieving correspondingly lower RMSDs to native structures. The new move set is also more efficient, being able to reach low energy states considerably faster. Use of this move set in determining the energy minimum state and for calculating thermodynamic quantities is discussed.

Entropic stabilization of proteins and its proteomic consequences.

Evolutionary traces of thermophilic adaptation are manifest, on the whole-genome level, in compositional biases toward certain types of amino acids. However, it is sometimes difficult to discern their causes without a clear understanding of underlying physical mechanisms of thermal stabilization of proteins. For example, it is well-known that hyperthermophiles feature a greater proportion of charged residues, but, surprisingly, the excess of positively charged residues is almost entirely due to lysines but not arginines in the majority of hyperthermophilic genomes. All-atom simulations show that lysines have a much greater number of accessible rotamers than arginines of similar degree of burial in folded states of proteins. This finding suggests that lysines would preferentially entropically stabilize the native state. Indeed, we show in computational experiments that arginine-to-lysine amino acid substitutions result in noticeable stabilization of proteins. We then hypothesize that if evolution uses this physical mechanism as a complement to electrostatic stabilization in its strategies of thermophilic adaptation, then hyperthermostable organisms would have much greater content of lysines in their proteomes than comparably sized and similarly charged arginines. Consistent with that, high-throughput comparative analysis of complete proteomes shows extremely strong bias toward arginine-to-lysine replacement in hyperthermophilic organisms and overall much greater content of lysines than arginines in hyperthermophiles. This finding cannot be explained by genomic GC compositional biases or by the universal trend of amino acid gain and loss in protein evolution. We discovered here a novel entropic mechanism of protein thermostability due to residual dynamics of rotamer isomerization in native state and demonstrated its immediate proteomic implications. Our study provides an example of how analysis of a fundamental physical mechanism of thermostability helps to resolve a puzzle in comparative genomics as to why amino acid compositions of hyperthermophilic proteomes are significantly biased toward lysines but not similarly charged arginines.

Lessons from the design of a novel atomic potential for protein folding.

We investigate all-atom potentials of mean force for estimating free energies in protein folding and fold recognition. We search through the space potentials and design novel atomic potentials with a random mixing approximation and a contact-correlated Gaussian approximation of decoy states. We show that the two derived potentials are highly correlated, supporting the use of the random energy model as an accurate statistical description of protein conformational states. The novel atomic potentials perform well in a Z-score and fold decoy recognition test. Furthermore, the designed atomic potential performs slightly and significantly better than atomic potentials derived under a quasi-chemical assumption. While accounting for connectivity correlations between atom types does not improve the performance of the designed potential, we show these correlations lead to ambiguities in the distribution of energetic contributions for atoms on the same residue. Within the confines of the model then, many potentials may exist which stabilize all native folds in subtly different ways. Comparison of different protein conformations under the various atomic potentials reveals both a remarkable degree of correspondence in the estimated free energies and a remarkable degree of correspondence in the identity of the contacts types that make the dominant contributions to the estimated free energies. This consistency may be interpreted as a sign that the design procedure is extracting physically meaningful quantities.

Recognition of nucleic acid bases and base-pairs by hydrogen bonding to amino acid side-chains.

Sequence-specific protein-nucleic acid recognition is determined, in part, by hydrogen bonding interactions between amino acid side-chains and nucleotide bases. To examine the repertoire of possible interactions, we have calculated geometrically plausible arrangements in which amino acids hydrogen bond to unpaired bases, such as those found in RNA bulges and loops, or to the 53 possible RNA base-pairs. We find 32 possible interactions that involve two or more hydrogen bonds to the six unpaired bases (including protonated A and C), 17 of which have been observed. We find 186 "spanning" interactions to base-pairs in which the amino acid hydrogen bonds to both bases, in principle allowing particular base-pairs to be selectively targeted, and nine of these have been observed. Four calculated interactions span the Watson-Crick pairs and 15 span the G:U wobble pair, including two interesting arrangements with three hydrogen bonds to the Arg guanidinum group that have not yet been observed. The inherent donor-acceptor arrangements of the bases support many possible interactions to Asn (or Gln) and Ser (or Thr or Tyr), few interactions to Asp (or Glu) even though several already have been observed, and interactions to U (or T) only if the base is in an unpaired context, as also observed in several cases. This study highlights how complementary arrangements of donors and acceptors can contribute to base-specific recognition of RNA, predicts interactions not yet observed, and provides tools to analyze proposed contacts or design novel interactions.

MicroRNA-155 deficient mice experience heightened kidney toxicity when dosed with cisplatin.

The development of nephrotoxicity limits the maximum achievable dosage and treatment intervals for cisplatin chemotherapy. Therefore, identifying mechanisms that regulate this toxicity could offer novel methods to optimize cisplatin delivery. MicroRNAs are capable of regulating many different genes, and can influence diverse cellular processes, including cell death and apoptosis. We previously observed miR-155 to be highly increased following ischemic or toxic injury to the kidneys and, therefore, sought to determine whether mice deficient in miR-155 would respond differently to kidney injury. We treated C57BL/6 and miR-155(-/-) mice with 20 mg/kg of cisplatin and found a significantly higher level of kidney injury in the miR-155(-/-) mice. Genome-wide expression profiling and bioinformatic analysis indicated the activation of a number of canonical signaling pathways relating to apoptosis and oxidative stress over the course of the injury, and identified potential upstream regulators of these effects. One predicted upstream regulator was c-Fos, which has two confirmed miR-155 binding sites in its 3' UTR and, therefore, can be directly regulated by miR-155. We established that the miR-155(-/-) mice had significantly higher levels of c-Fos mRNA and protein than the C57BL/6 mice at 72 h after cisplatin exposure. These data indicate a role for miR-155 in the cisplatin response and suggest that targeting of c-Fos could be investigated to reduce cisplatin-induced nephrotoxicity.