Computational modeling of drug response with applications to neuroscience.

The development of novel high-throughput technologies has opened up the opportunity to deeply characterize patient tissues at various molecular levels and has given rise to a paradigm shift in medicine towards personalized therapies. Computational analysis plays a pivotal role in integrating the various genome data and understanding the cellular response to a drug. Based on that data, molecular models can be constructed that incorporate the known downstream effects of drug-targeted receptor molecules and that predict optimal therapy decisions. In this article, we describe the different steps in the conceptual framework of computational modeling. We review resources that hold information on molecular pathways that build the basis for constructing the model interaction maps, highlight network analysis concepts that have been helpful in identifying predictive disease patterns, and introduce the basic concepts of kinetic modeling. Finally, we illustrate this framework with selected studies related to the modeling of important target pathways affected by drugs.

Introduction

The relatively poor efficacy of current drug treatments (illustrated, for example, with overall short survival times of cancer patients) along with the tremendous efforts in time and cost in drug development has given rise to a change in paradigm from blockbuster medicine to more personalized approaches.[1,2] While some neuropsychiatric diseases may remain barely treatable, individualized approaches have the potential to improve therapy in many other domains. Predicting individual molecular responses to a drug, for example, those caused by a certain mutation that foils the beneficial effects of the drug, is already addressed with gene tests prior to medical treatment, known as “precision medicine.”[3]

This paradigm shift has been enabled by new developments in biotechnology, in particular novel high-throughput sequencing technologies, along with quantitative proteomics. These have tremendously changed our abilities to implement personalized medicine.[4] On the other hand, these new developments impose huge challenges on the interpretation of these data in the clinical context. Large worldwide consortia, for example the 1000 Genomes Project[5] and International Cancer Genome Consortia,[6] generate huge amounts of sequencing information on human genetic variation that allow deep characterization of a patient's genome and contrasting of this information against the entire, “normal,” population. In contrast to the cancer field, neuroscience has been fairly reluctant to adopt genomics concepts, presumably because of the heterogeneity of neural material, the complexity of neural pathways, and, in the human situation, the lack of direct access to brain tissue, necessitating work with postmortem tissue.[7],[8] However, recent technology developments allow working with an ever-decreasing amount of material in more and more pure cell types, which makes genomics technology attractive for neuroscience. This has opened up the new research field of neurogenomics, which includes the genome-wide study of the nervous system and gene expression atlases of the brain, as well as network science.[9] Applications comprise DNA (re-)sequencing, the analysis of the transcriptome (RNA sequencing), the methylome (methyl sequencing) or the targeted enrichment of specific transcription factors using appropriate antibodies (chromatin immunoprecipitation [ChIP] sequencing) so that, in principle, dysfunctional molecular mechanisms in human diseases can be studied at a high level of detail (see Box for definitions of technical terms).[10] Additionally, proteomics technologies, such as the quantitative analysis of complex protein mixtures either with stable isotope labeling by amino acids in cell culture (SILAC)[11] or with isotope-coded affinity tags (ICAT),[12] or the measurement of protein expression with antibodies through reverse phase protein arrays (RPPA)[13] deliver quantitative measurements of proteins and thus information on functional changes in signaling pathways due to the disease under study.

While the measurement of the molecular landscape of a genome has quickly evolved during recent years, the bottleneck of molecular personalized medicine is clearly the computational analysis of the data along with its integration and interpretation in the disease context. Here, the sheer amount of data, comprising millions of data points, along with the complexity of the underlying data types that are often poorly correlated, is a challenge. Thus, new computational research fields such as “systems biology” or “systems medicine” have emerged that essentially aim to interpret genomic data at the molecular network level.[14]

Molecular networks are the key drivers of biological function. For example, a targeted drug exerts its effects by the inhibition or activation of drug targets that then activate or deactivate signaling cascades inside the cell, which in turn activate transcription factors that alter the gene expression response and, ultimately, lead to transient and rarely definite physiological or metabolic changes in the phenotype. Thus, the description of these networks is a key component in understanding the mode of action of a drug.

Knowledge about molecular interactions is spread through more than 500 dedicated data resources.[15] These databases are mostly curated, ie, their content is supervised by experts and annotators, and present detailed knowledge about specific organisms, specific interaction types (eg, protein-protein interactions, metabolic or signaling reactions) or specific disease domains (eg, cancer). On the other hand, there are ongoing attempts to try to integrate as much of these resources into meta-databases in order to derive more complete interaction networks.[16]

Molecular interaction networks have been used to infer function from high-throughput data and to draw hypotheses on the effects of drugs. There are essentially two strategies to explore drug action in the light of experimental data. The first approach maps experimental data onto large interaction networks and utilizes emerging properties of these networks through the computation of substructures or topological features. Results of this analysis are heavy-weighted subnetworks that can be used for qualitatively judging consequences of drug effects. The second approach uses kinetic modeling by describing the interaction network in mathematical terms, eg, with ordinary differential equation (ODE) systems, and by “translating” experimental data to kinetic model parameters (see Then, the dynamics of the system is simulated over time and results of key model parameters are used to judge the effect of the drug response.

In this article, we discuss the different steps in the course of modeling drug responses from experimental data, network building, as well as parameter fitting and model analysis. We review key aspects of the bioinformatics analyses and illustrate the conceptual framework with published examples from the neuroscience domain.

Resources for molecular pathway information

A mathematical model starts with the description of its species such as genes, proteins, protein complexes, or metabolites found to be dysregulated upon drug delivery. It is beyond the scope of this article to de scribe the different experimental omics and targeted approaches as well as primary data analysis methods that are used to quantify and identify the relevant molecules and, thus, we refer the reader to recent reviews on sequencing technology,[17] and proteomics applications,[18,19] as well as analysis methods, for example for transcriptome[20] and mass spectrometry-based proteome data (see Box).[21]

The next step in modeling is the definition of the interactions between the molecular components of interest. Interaction specification involves identifying the molecular drug targets along with their cellular signaling network and establishing relevant links to the metabolism and the gene regulatory system. A helpful concept in this direction is the molecular pathway concept. Pathway information, depicting the cellular information and mass flow from nuclear receptors to transcriptional regulators, is widely available. According to the recent pathway resource list, there are more than 500 databases that hold information on molecular interactions.[15] These databases comprise experimental results from large-scale studies on protein-protein interactions such as yeast-2-hybrid, experiments on protein-DNA interactions such as ChlP-sequencing, drug-target interactions or cellular signaling pathways annotated from functional experiments, literature, or expert knowledge. It should be noted, though, that the databases are manually curated so that variations in pathway annotation might occur. For example, for recent annotations of the well-studied “apoptosis” signaling pathway in three widely accessed databases, 76% of the annotated proteins are specific for a single database compared with only 24% of annotated proteins that are shared by at least two databases ( The reason for this is that pathway boundaries are not clearly defined, and that expert opinion on the extent of crosstalk with other pathways is highly variable. Additionally, pathway annotations are commonly focused on specific substructures or specific cellular context (eg, tissues, diseases, organisms), which may result in variations of the assembled protein lists.

This high variation in pathway annotation has given rise to the development of meta-databases, for example ConsensusPathDB[22] or STRING,[23] that integrate many different resources along with many different types of interactions, such as protein-protein, drug-target, signaling, and metabolic interactions, and, thus, allow a more complete view on pathways. Figure 1B exemplifies such an integrated view of drug action with fluoxetine (Prozac®). Fluoxetine is used as an antidepressant and belongs to the class of selective serotonin reuptake inhibitors (SSRIs). These compounds increase the extracellular level of the neurotransmitter serotonin in the brain. The compound has multiple interactions with metabolizing enzymes from the cytochrome P450 family as well as with the family of serotonin receptors and transporters.[24] By combining information on different interaction types, an integrated view can be generated in which drugs are linked to molecular targets and the molecular consequences of drug action can be monitored, building the basis for the generation of molecular models.

Pathway concepts can be directly used for the analysis of drug response with high-throughput data. For example, Yildirimman et al[25] investigated the carcinogenic hazard of different classes of chemicals in hepatocyte-like cells derived from human embryonic stem cells and found that transcriptome analysis based on pathways yields more robust classifiers as compared with singlegene expression analysis. Here, the expression changes for entire pathways were computed with a numerical score derived from all genes assigned to the pathway. Similar methodology has been carried out, eg, with T-profiler[26] or gene set enrichment analysis.[27]

Analyzing interaction networks

Besides the analysis of data in the light of preannotated pathways, researchers try to extract functional information from large interaction networks in an unsupervised way, ie, without prior knowledge on specific gene sets or pathways. In these approaches, molecular interactions are modeled with mathematical networks, ie, proteins correspond to network nodes and an interaction between two proteins is represented by an edge connecting the corresponding nodes. The human protein-protein interaction network is fairly large, comprising hundreds of thousands of interactions. Network theory aims at the computation of local and global properties of interaction networks and the deduction of emerging properties that might explain cellular function.[28,29] Network theory has a long tradition, and started in the 1950s with the seminal work of Erdösz and Rényi,[30] who invented the framework of random networks, ie, networks where nodes are connected by edges in a random way. Random networks are not what we observe in biological systems; however, local and global topological measures can be used to distinguish them from real-world networks. One such local measure is, eg, the clustering coefficient, C, of a node i, which measures its local connectedness. The measure evaluates the nodes that are connected with node i (the neighbors) by dividing the number of existing edges between the neighbors of node i by the number of all possible edges between the neighbors ( Another important local feature is the node degree, D., of node i, ie, the total number of edges of this node with other nodes. A global feature important for real-world graphs is the characteristic path length, L. This feature measures the average shortest path length between all possible pairs of nodes in the network. These characteristics were used in the seminal papers of Watts and Strogatz[31] and Barabasi and Albert[32] to distinguish biological networks from random ones. Watts and Strogatz found that socalled real-world networks, in contrast to random ones, have a high clustering coefficient, ie, a relatively high local connectivity while the average path length stays small. In addition, Barabasi and Albert introduced the notion of scale-free networks that show a small number of highly connected hubs and many nodes with only few connections. Typically, the node distribution of such networks follows a power law in contrast to the node distribution of random graphs that follow a Gaussian law. These characteristics explain the relative robustness of biological systems. Such scale-free networks are robust against random failures (eg, mutations) since it is very likely that the mutation occurs in a node with a low number of connections. On the other hand, the networks are very sensitive against failures of hub proteins since they are heavily connected with other proteins. As an example, we illustrate this with the human protein-protein interaction network composed of high quality interactions derived from the ConsensusPathDB. This network comprises 9 533 proteins with 80 422 interactions. It has an average cluster coefficient of 0.23 and characteristic path length of 3. The node degree follows the postulated power law (Figures 2B and 2C; see Box for explanation of terms used).

Large interaction networks have been applied in the analysis of high-throughput data. Ilofree et al,[33] for example, showed that protein-protein interaction networks can be used for classifying subtypes of diseases. Additionally, weighted gene coexpression network analysis (WGCNA)[34] has been introduced (see Box). This approach judges pairwise correlation of genes measured on their expression profiles from a set of experiments and organizes the gene pairs in a network structure. Next, multidimensional clustering is applied to generate clusters of local coexpression, which then can be associated with experimental phenotypes. Using this approach, several studies on human and animal brain generated stable network patterns with functional implications for neuroscience.[35,36,37] For example, in the latter work, Winden et al[37] used the WGCNA approach to investigate transcriptional profiles derived from microarray analysis of twelve populations of neurons in the adult mouse forebrain (cingulate cortex, pyramidal neurons, GABAergic interneurons, hippocampus, and thalamus, among others). They generated clusters of genes that could be associated with neuronal function and exert an underlying network structure of the neuronal transcriptome (see Box). They identified 13 modules with characteristic gene expression patterns that could be related to biological function, for example a module that corresponded to the subset of interneurons derived from the subpallium. As a highly connected gene in these samples they identified galectin-1, a previously unknown marker for this class of cells.

Basic concepts of dynamic modeling

While the approaches described so far generate static pictures of the cells under study in form of interaction networks or modules (subnetworks) thereof, dynamic computer modeling requires: (i) the construction of a directed network of model components (ie, proteins, complexes, metabolites), and (ii) the kinetic description of the interactions (see Box). Such models, when fitted to experimental data, can be used for predictive purposes, eg, predicting whether a perturbation introduced by a certain drug has specific effects on the model components, identifying model components that are highly sensitive with respect to the treatment, or predicting more general model features such as robustness and stability.[38,39]

The generic framework of such dynamic modeling is illustrated in Model design takes into account the experiments that interrogate the model along with a specification of the model components. Next, the interaction map has to be specified containing all interactions that interconnect the different species. The interaction map is typically retrieved from literature, interaction databases, or from the network approaches described in the previous section.[40] Prominent tools that assist in the automated construction of such interactions maps are, for example, CellDesigner,[41] Cytoscape,[42] and PathVisio.[43] Using these tools, large interaction maps of biochemical pathways have been designed such as the EGFR (epidermal growth factor receptor)[44] and the mTOR (mammalian target of rapamycin)[45] signaling pathways. Additionally, interaction maps for specific neurodegenerative diseases have been constructed, eg, for Alzheimer's disease (AD).[46] In the neuroscience field, there is specific interest in the modeling of neurons. Several software packages exist that allow construction and modeling of neurons such as NEURON[47] or neuroConstruct.[48]

In the deterministic kinetic modeling approach, the mathematical model is generated from the interaction map by defining reactions and the according reaction kinetics. Kinetic laws are mathematical terms that involve the vector of model components,, and kinetic parameters,. For example, the law of mass action for a unimolecular reaction has the kinetic law:

v = k S

Here, S indicates the total number of molecules of a species per cell, v denotes the rate and k the rate constant. Kinetic laws can be of any specificity and complexity, and are defined in order to reflect the dynamics and properties of reactions in networks. Additionally, positive and negative feedback loops can be introduced to realize specific influences on single reactions or oscillatory behavior.[49,50]

Figure 3B illustrates the general schema of a molecular model for drug response. Basically, three layers have to be implemented. The first layer describes the drug-target interactions and determines the effect of the drug on the target proteins; target proteins can be nuclear receptors, enzymes, ion channels, transporter proteins, among others, and are predominantly inhibited by the drug, however, there are several drugs that activate their targets. It should be noted that the space of cellular drug targets is rather limited; for example, early estimations by Drews and Ryser[51] suggested 487 potential targets, while after completion of the human genome sequence and the development of small-molecule drugs the estimated number of targets was increased to several thousands.[52] Current estimations range again in the order of hundreds of targets for a wide range of diseases.[53] The second layer describes the relevant cellular downstream effects such as signaling cascades, feedback loops, and crosstalk between different pathways. It has been discovered that signaling pathways are composed of a couple of typical building blocks such as Ras proteins, G protein cylces, phosphorelay systems, or MAPK (mitogen-activated protein kinase) cascades. The latter have been implicated in the pathogenesis of a variety of human disorders including cancer and neurodegenerative diseases such as AD, Parkinson's disease (PD), and amyotrophic lateral sclerosis (ALS). The third layer describes the essential readouts of the model that constitute the metabolic or physiological changes implied by the drug action. Responses of these readouts, eg, simulated changes in specific protein concentrations upon perturbation of the model by the drug, are used to predict drug effects as well as to judge and adapt the model parameters.

Generic kinetic models of basic molecular pathways are stored in model databases, for example BioModels,[54] and are represented with specific file formats that track information on all components, interactions, and mathematical descriptions of the model in order to ensure exchange, retrieval, and import to model analysis tools. Accepted community standards for model storage and representation are SBML[55] and BioPAX,[56] among others.

The critical bottleneck of molecular modeling is the uncertainty about the kinetic parameters, for example the rate constant in the above formula. These parameters have to be determined either by literature search or by computational analysis. Literature search involves databases that store already measured parameters from similar experiments and/or cellular models. BRENDA,[57] for example, is a database that includes comprehensive enzyme kinetic information. Parameters for metabolic models can be retrieved for many enzymes and, for example, Michaelis-Menten reaction rates. Similarly, SABIO-RK[58] holds information on reaction kinetics from various experimental conditions. Computational analysis uses selected model readouts and calibrates the model parameters to fit experimental data on these components, for example, derived from protein quantification (RPPA, Western blots), by minimizing a risk function. The set (or sets) of parameters that minimizes the predefined risk function is chosen as model parameters. Several mathematical approaches can be used, such as global optimization methods, either deterministic[59] or stochastic.[60] In particular, evolutionary approaches were shown to perform well in parameter optimization.[61] This class of algorithms uses evolutionary concepts such as mutations, selection, and recombination, in an iterative way in order to find the best set of parameters. Further stochastic approaches encompass, for example, Monte-Carlo sampling[62] or simulated annealing.[63]

Once the model is fitted to experimental data, simulation experiments can be carried out that assist in the identification of biomarkers (for example using sensitivity analysis[64]), in the testing of new hypotheses and the prediction of time and dosage effects of drugs and other external stressors.

Drug response modeling: selected examples

Here, we will describe recently published work that demonstrates the potentials and guide to possible further applications of drug response modeling.

Chang et al[65] investigated the cell signaling mechanisms induced by hallucinogens in human kidney cells. Using computational and experimental data they modeled the downstream effects of the activation of serotonin (5-hydroxytryptamine/5-HT) receptors 5-HT1AR and 5-HT2AR, both alone and in combination, and the resulting effects on the extracellular signal-regulated kinase (ERK) pathway that is a known target of hallucinogens. Dysfunction of serotonin receptors is a common observation in schizophrenia and depression, among others. Both receptors are coupled to different G proteins and induce different downstream signaling that results in a huge variability of drug response patterns dependent on which receptor is activated by the drug.[66] The authors build a computational model containing two branches of signaling cascades that result in MAPK activation by the different subtypes of serotonin receptors: while 5-HT1AR is coupled to Gi/G0 proteins and stimulates the MAPK cascade by phosphatidylinositol 3'-kinase (PI3K), 5-HT2AR is Gq/11 coupled and activates the MAPK cascade via protein kinase C (PKC)/Raf-1. Authors gave kinetic descriptions for transformation, translocation, and complex building of proteins and used mass action law as well as Michaelis-Menten kinetics to describe the dynamics. As the central model readout they used the amount of phosphorylated ERK1/2 at different time points. Simulation experiments were carried out with MATLAB. Kinetic parameters were partly taken from the literature and partly determined by parameter optimization using experimental data from Western blots measured on the pERK1/2 time courses as calibration. Using the model, the authors showed both experimentally and with model simulations that there are distinct responses dependent on the activation of the different receptors, in particular that 5-HT1AR agonists induced rather transient ERK activation by phosphorylation while 5-HT2AR agonists induced more sustained responses.

Wierling et al[67] examined the response of the epidermial growth factor receptor (EGFR)-mediated signaling network on targeted cancer therapies such as cetuximab (erbituxTM) and erlotinib (tarcevaTM) with a mere in silico approach. They generated an interaction map from literature[68] and described the kinetic interactions with a small set of different reaction types such as synthesis reactions, product formation, and degradation, which followed mass action law equations. In order to cope with the uncertainty in kinetic parameters, they used a massively parallel Monte-Carlo sampling strategy: kinetic parameters were sampled according to a log-normal probability distribution and simulations were carried out in the perturbed (drug action) and the nonperturbed state (no drug action) of the model. Resulting differences in model components were tracked. Next, this procedure was repeated multiple times with different parameter sets. Finally, an overall statistic was computed combining all repeated simulation runs. Simulation experiments were carried out with the custom modeling software PyBioS developed by the authors. Interestingly, they were able to show that the resulting variance in response of the model components was fairly low and, secondly, they were able to predict changes in drug response with respect to common knowledge on the effect of somatic mutations that are typically found in cancer patients, eg, a gain of function mutation in the B-Raf protein which leads to continuous proliferation even in the presence of an EGFR drug inhibitor. This paper exemplifies how computational modeling can be used to mimic different, and patient-specific, responses to targeted cancer therapy.

Chen et al[69] built and investigated a large computational model for the ErbB signaling pathway. The model contained different receptors (EGFR, HER2/neu2, ErbB3, and ErbB4) and immediate-early downstream responses involving MAPK and PI3K/AKT cascades. These receptors and activated signaling cascades control important cellular functions such as proliferation, survival, and motility, and dysregulation of the pathway is key in various cancer types[70] as well as neurodegenerative diseases.[71] The model contained kinetic descriptions for dimerization/phosphorylation reactions, complex binding, and receptor binding, and follows mass action kinetics. Parameter values were determined with simulated annealing starting with a literature-defined parameter set. Model readouts were phosphorylated ERK and AKT and simulated time courses were fitted to experimental data measured in different cancer cell lines. The computational analysis was done with MATLAB. In addition, the authors used the model to predict dose response curves from targeted inhibitors gefitinib and lapatinib.

The above selective examples demonstrate how the model-building and model-fitting process is used to predict effects of therapies and to answer complex questions such as the different response in activation of subtypes of receptors through different agonists. However, it should be noted that besides these potentials there are major pitfalls of such modeling approaches, in particular the lack of quantitative data, the low number of time points measured for the individual model components and their restriction to particular cells (or subcellular fractions), to specific time points, and to a relatively small number of model components compared with the full genome active in a “real-world” patient. Nonetheless, computational models are helpful complements to experimental medicine and once a model is constructed and fitted to experimental data, quantitative behavior of readouts can aid in the generation of new hypotheses and the testing of new therapies. Furthermore, in silico experiments can be conducted in cases that are hardly, or even not, measurable with experiments, eg, multitherapy approaches. In cases where clinical testing is possible, for example with tumor-carrying xenograft mice, such in silico experiments could be used to prioritize therapy combinations, which in turn can be validated experimentally.

Conclusion and outlook

In this article we present the basic concepts and steps in the workflow of molecular modeling. We focused on time-dependent modeling inside the cell, in particular receptor-mediated signaling. It should be noted though that computational modeling is by no means restricted to this situation and comprises, for example, concepts for modeling spatial processes,[72] pharmacokinetic modeling,[73] as well as whole -body[74] and multiscale approaches[75] that have not been discussed in this work.

Systems biology-driven modeling of drug response will gain more importance with the ever-increasing amounts of genome data becoming available for individual patients. The further decrease in sequencing costs along with the generation of purer materials, such as single-cell omics analyses[76,77] and material adapted closer to the individual patient, such as induced pluripotent stem (iPS) cells,[78,79] enable completely new chances for personalized therapies. Here, in particular the incorporation of mutations and larger human variations, eg, in neuro-oncology,[80,81] and allele-specific protein expression[82,83] will gain further importance. Additionally, the development of new imaging technology and the combination of different approaches will drive further insights.[84] Previously intractable to the approaches of human genetics, disorders of the brain are seeing accelerated gene discovery that was, until now, restricted to other branches of medicine. Genome analysis and genetic determinants in the research of neurodegenerative diseases will further give rise to linkage of molecular data with therapy decisions, for example with AD.[85] Neurogenomics as an emerging focus of neuroscience will further evolve with the results of genome-wide association studies and the knowledge of how these might be translated to the understanding of neurodegenerative diseases[86] and psychiatric disorders.[87] As brain research has attracted immense funding opportunities in recent years, and large initiatives are on the way such as the European Union's Human Brain Project or the US Brain Activity Map Project,[88] it can be anticipated that further data and modeling methodology will be generated and will increase the demand for in silico tools in neuroscience as a useful and efficient complement to the experiments.