Search

Making drug development a more efficient and cost-effective process will have a transformative effect on human health. A key, yet underutilized, tool to aid in this transformation is mechanistic computational modeling. By incorporating decades of hard-won prior knowledge of molecular interactions, cellular signaling, and cellular behavior, mechanistic models can achieve a level of predictiveness that is not feasible using solely empirical characterization of drug pharmacodynamics. These models can integrate diverse types of data from cell culture and animal experiments, including high-throughput systems biology experiments, and translate the results into the context of human disease. This provides a framework for identification of new drug targets, measurable biomarkers for drug action in target tissues, and patient populations for which a drug is likely to be effective or ineffective. Additionally, mechanistic models are valuable in virtual screening of new therapeutic strategies, such as gene or cell therapy and tissue regeneration, identifying the key requirements for these approaches to succeed in a heterogeneous patient population. These capabilities, which are distinct from and complementary to those of existing drug development strategies, demonstrate the opportunity to improve success rates in the drug development pipeline through the use of mechanistic computational models.

Vascular Endothelial Growth Factor (VEGF) signal transduction is central to angiogenesis in development and in pathological conditions such as cancer, retinopathy and ischemic diseases. We constructed and validated a computational model of VEGFR2 trafficking and signaling, to study the role of receptor trafficking kinetics in modulating ERK phosphorylation in VEGF-stimulated endothelial cells. Trafficking parameters were optimized and validated against four previously published in vitro experiments. Based on these parameters, model simulations demonstrated interesting behaviors that may be highly relevant to understanding VEGF signaling in endothelial cells. First, at moderate VEGF doses, VEGFR2 phosphorylation and ERK phosphorylation are related in a log-linear fashion, with a stable duration of ERK activation; but with higher VEGF stimulation, phosphoERK becomes saturated, and its duration increases. Second, a large endosomal fraction of VEGFR2 makes the ERK activation reaction network less sensitive to perturbations in VEGF dosage. Third, extracellular-matrix-bound VEGF binds and activates VEGFR2, but by internalizing at a slower rate, matrix-bound VEGF-induced intracellular ERK phosphorylation is predicted to be greater in magnitude and more sustained, in agreement with experimental evidence. Fourth, different endothelial cell types appear to have different trafficking rates, which result in different levels of endosomal receptor localization and different ERK response profiles.

The vascular endothelial growth factor (VEGF) family binds multiple endothelial cell surface receptors. Our goal is to build comprehensive models of these interactions for the purpose of simulating angiogenesis. In view of low concentrations of growth factors in vivo and in vitro, stochastic modeling of molecular interactions may be necessary. Here, we compare Monte Carlo simulations of the stochastic binding of VEGF and two of its major receptors on cells in vitro to equivalent deterministic simulations. In the range of typical VEGF concentrations, the stochastic and deterministic models are in agreement. However, we observe significant variability in receptor binding, which may be linked to biological stochastic events, e.g., blood vessel sprout initiation. We study patches of cell surface of varying sizes to investigate spatial integration of the signal by the cell, which impacts directly the variability of binding, and find significant variability up to the single-cell level. Dimerization of VEGF receptors does not significantly alter the variability in ligand binding. A ‘sliding window’ approach demonstrated no reduction in the variability of binding by temporal integration. The variability is expected to be more prominent in in vivo situations where the number of ligand molecules available for binding is less.

In this issue of the Biophysical Journal, Biedermann et al. (1) describe a new computational package that implements a split GPU-CPU algorithm for speeding up the calculations of particle dynamics and reaction kinetics in heterogeneous biomolecular mixtures. The movement and interactions of proteins and other biomolecules are central to understanding drug action and improving drug design. Computational models of these molecular movements and reactions are indispensable, and many different simulation scales are available (Fig. 1 A), from atom-level resolution using molecular-dynamics packages such as CHARMM (www.charmm.org) (2) to continuum reaction-diffusion models at the tissue level (3) and compartmental models at the whole-body level (4). Figure 1 Simulation of heterogeneous molecular mixtures. (A) As the level of particle detail and particle discrimination increases (to the right), the computational time increases, while the ability to simulate large volumes and long times decreases. (B) Simplified ... If we focus on cell-level and subcellular-level simulations, then low-protein-copy number, spatial crowding, compartmentalization, and heterogeneity can make stochasticity a key consideration. At this scale, discriminating between the individual molecules becomes more important, and continuum approaches cannot do this. Methods that can explicitly simulate individual molecules are computationally costly; one approach to reduce the cost is to couple different model types. For example, molecular-dynamics simulations over short times can be used to estimate binding and unbinding energies; these can then be used in continuum models to estimate the longer-term behavior of a population of molecules in response to a perturbation. However, in many cases, to understand the system we will want to track each of the molecules individually; for these purposes the appropriate simulation scale is not atomistic, but instead that of individual molecular particles, each trackable over time as it moves and becomes involved in intermolecular reactions. Several packages have been developed to facilitate this level of simulation of individual particle movement and reaction, including MCell by groups at Salk Institute for Biological Studies (La Jolla, CA) and the University of Pittsburgh (Pittsburgh, PA; mcell.org) (5); Smoldyn by the Andrews group at the Fred Hutchinson Cancer Research Center (Seattle, WA; www.smoldyn.org) (6); and ReaDDy by the Noe group at the Free University of Berlin (Berlin, Germany; www.readdy-project.org) (7). Several of these modules have been or are being integrated into the larger simulation environments that model cell signaling, e.g., Spatiocyte in E-Cell (www.e-cell.org) and Smoldyn in VCell (University of Connecticut Health Center, Farmington, CT; www.nrcam.uchc.edu). Packages that use a lattice-based movement system and do not discriminate between all particles are also available, such as the following: MesoRD by the Elf lab at Uppsala University (Uppsala, Sweden; http://mesord.sourceforge.net/); StochSim by the Bray group at The University of Cambridge (Cambridge, United Kingdom; http://sourceforge.net/projects/stochsim/). Spatiocyte by the Tomita group at the Institute for Advanced Biosciences, Keio University (Fujisawa, Japan; spatiocyte.org) (8) uses a lattice, but individual molecules are still represented. The key problem with this level of simulation—tracking individual particles—is the computational cost; the time required to simulate all proteins within a single cell can be prohibitive. This cost increases with particle density, with simulation volume, and with the complexity and heterogeneity of the spatial domain. For the most part, this computational cost is not due to the calculation of the reactions that the proteins undergo with each other, but instead to the particle dynamics. Methods to speed up the calculation of particle dynamics could greatly enhance the ability to simulate subcellular processes such as signaling. Biedermann et al. (1), in this issue, report a new implementation of the ReaDDy package, ReaDDyMM, that leverages OpenMM and splits the calculations—the particle dynamics (movement, e.g., Brownian dynamics) to be calculated by the GPU and the reactions to be calculated by the CPU (Fig. 1 B). Three conditions must be met for this to obtain the reported two-orders of magnitude speedup, as follows. 1. There must be considerable separation in timescales between the particle dynamics and the reaction kinetics. ReaDDyMM takes advantage of the separation in timescales by running many smaller timesteps of particle dynamics for each larger timestep of reaction kinetics (Fig. 1 C). If this were not the case, then communication between the processors would make the algorithm inefficient. Greater timescale separation—faster diffusion, slower reactions—allows increased n (the number of GPU dynamics steps per CPU kinetics step; Fig. 1 B). 2. The level of speedup increases with particle density. At lower densities, both the overall computational time and the speedup are lower. This methodology thus appears particularly suited to understanding highly crowded biological environments. 3. Both the diffusion timestep (τ) and the reaction timestep (nτ) must be less than their corresponding characteristic timescales; otherwise, accuracy can be severely affected. ReaDDyMM split CPU-GPU calculation is distinct from Smoldyn’s GPU implementation, in which both diffusion and reactions are calculated on the GPU (9), and thus may provide an improvement for complex reaction networks. Benchmarking for both accuracy and computation time in various scenarios may indicate application-specific benefits to one approach over the other. Are two orders of magnitude of speedup sufficient to make this form of molecular simulation feasible? The answer depends on the application. The volumes simulated in the presented examples (1) are small—an attoLiter, or one-thousandth of an Escherichia coli cell. Mammalian cells are larger still. The protein densities used are similar to the total protein densities in E. coli. Thus, simulating a whole cell using this method would take one or more weeks, but this may have the potential to be further improved through parallelization. The benefit is the level of detail and accuracy possible, and for subcellular components this method will be both quick and highly detailed. ReaDDyMM is a useful advance and increases the tools and feature sets available for simulation of heterogeneous multimolecular simulations. This is an opportune time for such simulation packages, as superresolution microscopy and other advances make possible the experimental tracking of individual molecules within a complex cellular context.

We have created an experimentally-based computational model describing spatial transport of vascular endothelial growth factor (VEGF) and its receptors to quantitatively understand how guidance cues may modulate blood vessel sprout growth. VEGF binds to endothelial cells and initiates angiogenesis. Both VEGF concentration and VEGF gradients may control sprout formation. Soluble VEGF receptor 1 (sFlt1) can bind and sequester VEGF. Based on observations in developing vasculature, we hypothesize that a local reduction in sFlt1 expression can increase locally available VEGF and thus control angiogenesis. However, the complex VEGF interaction network makes it difficult to isolate how individual proteins contribute to the spatial distribution of the growth factor using experiments alone. Our computational model represents the local environment of a single blood vessel and nearby tissue and directly incorporates the network of VEGF interactions. In the model, parenchymal cells secrete VEGF, which diffuses through interstitial space and binds extracellular matrix (ECM) and sFlt1. VEGF binds endothelial cells via membrane-bound receptors Flt1 and Flk1, and endothelial cells secrete sFlt1. Additionally, the model accounts for degradation of VEGF and sFlt1 as well as internalization of receptor-bound ligands. Using partial differential equations, we simulate this system, which is constrained by experimentally-derived parameters. Our simulations show that when a sprout-leading tip cell secretes less sFlt1 than neighboring cells, there is decreased local sFlt1 sequestration of VEGF, thus resulting in augmented VEGF-Flk1 levels on the surface of the low-sFlt1 secreting tip cell. This could lead to sprout generation. We also show how variations in sFlt1 secretion and tip cell configuration may affect the gradients of guidance cues and directionality of sprout growth.

Vascular endothelial growth factor (VEGF) is an important regulator of blood vessel growth. Matrix-binding and non-matrix-binding isoforms of VEGF stimulate production of blood vessel networks that are structurally and functionally different. Here, we develop a computational model of the binding of soluble and immobilized ligands to VEGF receptor 2 (VEGFR2), endosomal trafficking of VEGFR2, and site-specific VEGFR2 tyrosine phosphorylation to study differences in induced signaling between these VEGF isoforms. We build our model using experimental data in multiple cell lines, under different conditions, from several groups. In capturing essential features of VEGFR2 signaling and trafficking, our model suggests that VEGFR2 trafficking parameters are largely consistent across multiple endothelial cell lines. Simulations demonstrate distinct cellular localization of VEGFR2 phosphorylated on tyrosines 1175 and 1214. This is the first model to clearly show that differences in site-specific VEGFR2 activation when stimulated with immobilized VEGF compared to soluble VEGF can be accounted for by altered trafficking of VEGFR2 without an intrinsic difference in receptor activation. The model predicts that Neuropilin-1 can induce differences in the surface-to-internal distribution of VEGFR2, and that ligated VEGFR2 and phosphorylated VEGFR2 levels diverge over time following stimulation. We validated our model by successfully predicting changes in VEGFR2 phosphorylation resulting from perturbations of multiple phosphatases acting on VEGFR2. Using this model, we identify multiple key levers that alter how VEGF binding to VEGFR2 results in different coordinated patterns of multiple downstream signaling pathways. Specifically, simulations predict that perturbations of VEGFR2 trafficking (e.g. via VEGF immobilization), interactions with Neuropilin-1, and changes in expression or activity of phosphatases acting on VEGFR2 all affect the magnitude, duration, and relative strength of VEGFR2 phosphorylation on Y1175 and Y1214, and they do so predictably within our single consistent model framework.