Molecular BioSystems View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 537
View Journal | View Issue
In silico identification of potential therapeutic targets in the TGF-b signal transduction pathway† Daniel Nicklas and Leonor Saiz* The transforming growth factor-b (TGF-b) superfamily of cytokines controls fundamental cellular processes, such as proliferation, motility, differentiation, and apoptosis. This fundamental role is emphasized by the widespread presence of mutations of the core components of the TGF-b signal transduction pathway in a number of human diseases. Therefore, there is an increasing interest in the development of therapies to specifically target this pathway. Here we develop a computational approach
Received 30th June 2013, Accepted 10th December 2013
to identify potential intervention points that are capable of restoring the normal signaling dynamics to
DOI: 10.1039/c3mb70259f
this approach explicitly to the TGF-b pathway to study the signaling dynamics of mutated and normal
the mutated system while maintaining the behavior of normal cells substantially unperturbed. We apply cells treated with inhibitory drugs and identify the processes in the pathway that are most susceptible to
www.rsc.org/molecularbiosystems
therapeutic intervention.
Introduction The transforming growth factor-b (TGF-b) signal transduction pathway controls multiple cellular processes, including proliferation, motility, differentiation, and apoptosis, that are essential for regulating embryonic development, tissue homeostasis, and wound healing.1 Dysregulation of the pathway through mutations of its core components is associated with a number of human diseases, such as cancer, fibrotic conditions of kidneys, liver, and lungs, and vascular and autoimmune disorders.2–4 Because of its involvement in many aspects of human health, a significant effort of clinical research focuses on developing therapies targeting the activity of the TGF-b signaling pathway.5,6 TGF-b-superfamily ligands signal through a network of interacting molecular components with complex patterns of coupled signaling and feedback regulation. At the plasma membrane, the ligands initiate signaling by binding to two types of serine–threonine kinases, known as type I and type II receptors. The active ligand–receptor complex propagates the signal through the intracellular mediator Smad proteins upon internalization into the endosome, where it may recruit and phosphorylate a receptor-regulated Smad (R-Smad). Phosphorylated R-Smads bind to Smad4 and translocate into the nucleus, where the resulting macromolecular complex regulates the
Modeling of Biological Networks Laboratory, Department of Biomedical Engineering, University of California, 451 East Health Sciences Drive, Davis, CA 95616, USA. E-mail: [email protected]; Fax: +1-530-754-5739; Tel: +1-530-752-6700 † Electronic supplementary information (ESI) available. See DOI: 10.1039/ c3mb70259f
This journal is © The Royal Society of Chemistry 2014
expression of several hundred genes, including Smad7. Smad7, an inhibitory Smad, negatively regulates signaling by binding to active ligand–receptor complexes at the plasma membrane and inducing their degradation. TGF-b-superfamily ligands use two parallel R-Smad channels, grouped as Smad1/5/8 and Smad2/3, to transmit the signal. Bone morphogenetic proteins (BMPs) and nodal/activin ligands induce phosphorylation of the Smad1/5/8 and Smad2/3 groups, respectively, while TGF-b induces phosphorylation of both R-Smad groups.7 A number of mutations affecting the TGF-b pathway signaling components have been identified in a variety of human cancers.2,6 While TGF-b is typically a growth suppressor, its dysregulated activation may promote invasion and growth of tumors in cancer cells.8 Therefore, a number of strategies have been developed to abrogate this effect by inhibiting TGF-b activity, such as the use of antisense oligonucleotides and monoclonal antibodies targeting TGF-b mRNA and protein, respectively, and small molecule kinase inhibitors that target the receptors, blocking R-Smad phosphorylation.6,9 Quantitative and predictive models have been successfully used to study key modules in the TGF-b signal transduction pathway, including receptor trafficking,10–12 Smad nucleocytoplasmic shuttling,13,14 and transcriptional feedback loops.15–20 Several models have combined these components into comprehensive mathematical representations of the signaling pathway.19,21–23 In addition to providing insight into the underlying mechanisms, computational models of signaling pathways have been used to study the effects of mutations in disease.22,24–26 In the TGF-b signaling pathway, for instance, Chung et al. investigated the effects of cancerous mutations on Smad-dependent signaling. Specifically, they used their model to predict how different
Mol. BioSyst., 2014, 10, 537--548 | 537
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
signaling components are affected by cancerous mutations by introducing a mutation resulting in downregulation of TGF-b receptors through decreased initial concentration and synthesis rate of the affected species.22 Furthermore, it has been shown that computational modeling can be used as a fundamental step in the development of novel pharmaceuticals.27,28 A prime example of the applicability of computational approaches is the epidermal growth factor receptor (EGFR) signaling pathway, which has been extensively studied using a number of computational models to investigate both receptor trafficking and downstream signaling.29 In particular, models of the EGF-induced extracellular signal regulated kinase (ERK) signal transduction pathway have been successfully used to study the effects of oncogenic mutations,25 to assess the role of mutation in the sensitivity of cells to an existing therapy,24 and to identify therapeutic targets and biomarkers.30 Specifically, Orton and coworkers developed a model of the ERK signal transduction pathway and compared the signaling dynamics of normal and cancer cells as modeled through knockouts and mutations of signaling proteins. From this, the authors hypothesized which components of the signaling pathway were susceptible to therapeutic intervention.25 Using similar strategies, Hendriks et al. introduced perturbations into their model to represent the effects of mutant receptors and tested potential mechanisms by which the mutated signal dynamics observed experimentally differed from those of normal cells.24 Using sensitivity analysis, Lebedeva and coworkers identified a number of potential therapeutic targets and biomarkers in the pathway.30 In addition to investigating the effects of mutations, computational modeling has been used to analyze how therapies interact with mutated signaling pathways and to propose novel targets for therapeutic intervention in a number of systems and diseases.31–36 Specifically, Sung and Simon compared the effects of novel inhibitor drugs on the signaling dynamics of the nuclear factor-kB (NF-kB) pathway.32 Yan et al. and Araujo et al. extended this analysis in the NF-kB and EGFR pathways, respectively, to simulations of drug combinations.31,33 In particular, Araujo and coworkers evaluated the effects of kinase inhibitors targeting one or two nodes in the pathway, revealing that simultaneous administration allows for a reduction in the required dose necessary to significantly affect the signal response.31 Mathematical modeling of transferrin binding and receptor trafficking was used to successfully propose a novel mechanism to improve its drug delivery capabilities,34 which was validated with in vitro and in vivo experiments.34–36 Here we apply a comprehensive computational model of the TGF-b signal transduction pathway19 to study the signaling dynamics of mutated and normal cells treated with inhibitory drugs or treatments. Using this model, we build upon strategies applied to the TGF-b signaling pathway and other systems to develop an effective in silico approach for identifying targetable processes and potential intervention points that are capable of restoring the normal signaling dynamics to the mutated system while maintaining the behavior of normal cells substantially unperturbed.
538 | Mol. BioSyst., 2014, 10, 537--548
Molecular BioSystems
We have chosen the TGF-b signaling system because of its prominent role in cancer.8 Cancer is a multistage process that involves the successive failure of many control points, such as cancer cells stimulating their own growth and resisting inhibitory signals that might otherwise stop their growth as well as the evasion of programmed cell death. Therefore, dysregulation of the TGF-b pathway in cancer cell lines is accompanied by the dysregulation of additional pathways, such as those arising from K-ras oncogene mutations.37,38 There are many therapies under development that target the components of the TGF-b pathway.6 One of the most sought features of a therapy is the selective targeting of the mutated cells while keeping the normal cells largely unperturbed. So far there has not been a systematic approach to identify those therapies. The in silico approach we present here is a first step in this direction, which can be applied to a number of other signaling systems as a guiding tool in the therapy development process. Subsequent aspects that would need to be taken into account in the full therapy development include the differences between cell lines and tumor samples as well as the cross-talking between different signaling pathways.
Methods The TGF-b pathway model To study the dynamics of the pathway upon stimulation with TGF-b, we use the comprehensive mathematical model developed in ref. 19, which has been shown to accurately reproduce the behavior of human keratinocytes (HaCaT), bovine aortic endothelial cells (BAEC), and the mouse mesenchymal C2C12 cell line. The model considers explicitly receptor trafficking, downstream signaling by intracellular Smad proteins, and negative feedback through Smad7. Signaling is initiated when a TGF-b ligand binds to its type II receptor, denoted as RII. This complex then recruits a type I receptor, one of RI1 or RI2, forming the ligand–receptor complex C1 or C2, respectively. The ligand–receptor complexes propagate the signal upon internalization into the endosome, where C1 may bind with cytosolic Smad1 (S1c) and C2 may bind with cytosolic Smad2 (S2c) and phosphorylate the R-Smad. We use the subscripts c and n to denote cytosolic and nuclear species, respectively, and apply a prefix p to phosphorylated species. pS1c and pS2c bind to Smad4 (S4c), forming the pS1S4c and pS2S4c complexes, respectively, which may translocate into the nucleus where pS1S4n and pS2S4n then may promote the expression of Smad7 (S7). This initiates a negative feedback loop in which S7 binds to the ligand–receptor complexes (C1 or C2) at the plasma membrane, preventing their association and phosphorylation of R-Smad proteins, and targeting them for degradation. A schematic representation of the model is shown in Fig. 1, where arrows indicate modeled reactions in the pathway. Reactions are mathematically represented using mass-action kinetics, which are then combined to form the system of ordinary differential equations (ODEs) to track the rate of change of each modeled species (Table S1, ESI†). Parameter values for the three cell types are listed in Table S2 (ESI†).10,19,22,56–64
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Paper
Fig. 1 Schematic illustration of the Smad-dependent TGF-b signaling pathway model upon stimulation with TGF-b. Arrows denote reaction steps in the pathway and are labeled with the rate constant for the reaction. Synthesis and degradation reactions are provided in the table (inset) with the notation for the different molecular species given in the text. Expression of Smad7 is under the control of a single gene, activated by nuclear phosphorylated Smad1–Smad4 and/or Smad2–Smad4 complexes.
Prior to ligand stimulation, we determine the steady state solution for the system of ODEs (Table S1, ESI†) by setting each time-derivative to zero and solving the resulting linear system of equations with the concentration of the ligand set to zero using the ‘linalg.solve’ method in Numpy 1.6.2 (http://numpy. scipy.org). Upon ligand stimulation, we numerically solve the system of ODEs using the CVODE method in the SUNDIALS 2.5.0 package.39 These correspond to the typical experimental conditions that measure the response of the system to a sudden change of the TGF-b ligand concentration from zero to a saturating value that is kept constant afterwards.
Mutation and molecular intervention simulations In order to capture the effects of mutations in the pathway, we use a strategy similar to that employed in ref. 22 and 31 to introduce altered protein expression or functional mutation. Focusing on the mutations outlined in Table 1, we define a mutation factor (MF) within the range [0,1] with which to reduce the parameter values associated with each mutation. For example, to introduce the ‘inhibited Smad4 nuclear import’ mutation with a 0.1 mutation factor, we multiply k7imp and k14imp by 0.1. Similarly, we simulate the effects of a molecular intervention by introducing an intervention factor (IF) within
This journal is © The Royal Society of Chemistry 2014
the range [0,1] that is multiplied with the rate constants affected by a particular molecular intervention. Sensitivity analysis We use a global sensitivity analysis to investigate how the model responds to perturbations of individual parameters in the parameter space.40 Concretely, we compute the derivativebased global sensitivity measures developed in ref. 41, which sample the effects of local parameter perturbation in a larger parameter space. To compute the effects of local parameter perturbation, we use the scaled sensitivity coefficients42 given by eki; j ¼
ki; j @Cj ; Cj @ki; j
(1)
where ki,j defines the value of the ith parameter of sample j and Cj is the time integral of the model output over 24 hours after stimulation with the ligand for the parameters of the jth sample. As model output, we examine the signaling dynamics of both nuclear phosphorylated Smad1–Smad4 (pS1S4n) and nuclear phosphorylated Smad2–Smad4 (pS2S4n) complexes since they control the long-term effects of ligand stimulation through regulation of gene expression.7 To approximate the partial derivative, we simulate the model with 0.5 percent
Mol. BioSyst., 2014, 10, 537--548 | 539
View Article Online
Paper
Molecular BioSystems
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Table 1 Potential targetable processes and their associated rate constants analyzed in the optimization of the molecular interventions. The parameters associated with each process are multiplied by an intervention factor in the range [0,1] to simulate the process inhibition
Inhibited process
Parameters
Type II receptor synthesis Type I receptor synthesis R-Smad synthesis Smad4 synthesis Smad7 synthesis Association of TGF-b with type II receptor Association of TGF-b–type II receptor complex with type I receptor Ligand–receptor complex–R-Smad association R-Smad–Smad4 association R-Smad nuclear import R-Smad nuclear export Smad4 nuclear import Smad4 nuclear export Ligand–receptor complex–Smad7 association Association of phospho-R-Smad–Smad4 complex with DNA
ksyn,RII ksyn,RI ksyn,RS ksyn,S4 ksyn,S7, klip,1, and klip,2 k1a k2a k4a k6a and k10a k7imp, k12imp, and k13imp k12exp k7imp and k14imp k14exp k20a,1 and k20a,2 KA,1 and KA,2
perturbations of parameter ki about its value at sample j and calculate the finite central difference.43 For each parameter ki, we calculate the local sensitivity coefficient eki,j for 4000 parameter values, i.e. j = 1,. . .,4000, by randomly sampling around its original value (Table S2, ESI†) in the parameter space. Specifically, we generate a random variable x with a uniform distribution within the range [1,1] and multiply the original parameter value by 10x. The scaled sensitivity coefficients are then used to compute the global sensitivity coefficient41 defined by Gki ¼
N 1X ek 2 ; N j¼1 i; j
(2)
where N is the number of samples, equal to 4000 here. We implemented the sensitivity analysis using Python 2.7.3 (http:// www.python.org), Numpy 1.6.2 (http://numpy.scipy.org), and Scipy 0.10.0.44 Molecular intervention optimization We optimize the intervention factor IF, defined as a value within the range [0,1], by fitting the simulation result of a mutated cell with the specific molecular intervention to that of a normal cell. Additionally, we fit the simulation result of the normal cell with the molecular intervention to that of a normal cell in order to minimize the potential risk of side effects from exposing normal cells to the molecular intervention since the goal of our approach is to identify those molecular interventions that lead to the recovery of normal behavior by mutated cells and that minimally affect the behavior of normal cells. Therefore, we define the objective function for minimization as the least-squares error: E¼
n X
ðSðti Þ M ðti ÞÞ2 þðSðti Þ S ðti ÞÞ2 ;
i¼1
540 | Mol. BioSyst., 2014, 10, 537--548
(3)
where S(ti) and M(ti) are the normal and mutated cell dynamics, respectively, at the time point ti of n total time points and the asterisk represents the corresponding dynamics with the molecular intervention. To determine the optimized intervention factor (OIF) we first evaluate the objective function with intervention factors within the range [0.1,0.9] with 0.1 increments. The intervention factor with the minimum objective function value becomes the starting point for a downhill simplex algorithm with the ‘optimize.fmin’ function in Scipy 0.10.0.44 In order to impose the bounds of [0,1] for the intervention factor, we use the transformations described in ref. 45. Our computational implementation uses Python 2.7.3 (http://www.python.org), Numpy 1.6.2 (http://numpy.scipy.org), and Scipy 0.10.0.44
Results and discussion Model dynamics To examine the ability of the model to reproduce the observed behavior in mutated cell lines, we compare the results of the simulation with experimental data in HaCaT cells and two pancreatic cancer cell lines, namely Colo-357 and PT45 (Fig. 2). The TGF-b signaling in these cell lines has been studied in detail in ref. 46. The results show that Colo-357 cells exhibit an epithelial phenotype, forming well-defined epithelial sheets and expressing E-cadherin at their adherens junctions. Colo357 cells respond to TGF-b in a way very similar to HaCaT cells, strongly inducing the expression of p21, JunB, c-Jun, Smad7, and PAI-1 mRNA within 1 hour of ligand stimulation. PT45 cells, in contrast, do not present the epithelial phenotype and the induction of p21 and JunB mRNA with TGF-b is only transient. Colo-357 cells, as HaCaT cells, display rapid nuclear accumulation of phosphorylated Smad2 that is still present, but reduced 6 hours after activation. In contrast, nuclear accumulation of phosphorylated Smad2 in PT45 cells disappears abruptly after 1 hour of stimulation. The kinetics of nuclear accumulation of Smad4 was observed to be very similar to that of phosphorylated Smad2 in both Colo-357 and PT45 cells. Using the parameters for HaCaT cells (Table S2, ESI†), the model reproduces the normal signaling dynamics of both nuclear phosphorylated Smad2 (Fig. 2A) and nuclear Smad4 (Fig. 2D) molecular species. To simulate the dynamics of pancreatic cancer cells Colo-357 and PT45, we introduce the ‘inhibited Smad4 nuclear import’ mutation found in pancreatic cancer cell lines (Table 2) to the normal cases with a mutation factor of MF = 0.1. As Colo-357 cells exhibit an epithelial phenotype,46 we use the model for HaCaT cells to reproduce their signaling dynamics. By introducing the mutation, the model results change only minimally for both nuclear phosphorylated Smad2 (Fig. 2B and Fig. S1, ESI†) and nuclear Smad4 (Fig. 2E and Fig. S1, ESI†), consistently with the experimental data. The differences between the model and experimental results observed in the early stages of signaling in Fig. 2B and E are attributable, as discussed in ref. 11, to the detailed trafficking dynamics of the receptors. We use the parameters for C2C12 cells to reproduce the PT45
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Paper
Fig. 2 Simulation of the effects of pancreatic cancer mutations on TGF-b signaling dynamics. We simulate the dynamics of the model using the parameters for HaCaT (black solid lines) and C2C12 cells (blue dashed lines) and show the results for nuclear phosphorylated Smad2 (top row; A–C) and nuclear Smad4 (bottom row; D–F) species upon stimulation with 2 ng mL1 TGF-b at t = 0 hours. Experimental data (symbols) was quantified using ImageJ 10.252 from Western blot experiments in ref. 46 for HaCaT cells (left column; A and D) and two pancreatic cancer cell lines: Colo-357 (middle column; B and E) and PT45 (right column; C and F). To simulate the dynamics of the pancreatic cancer cell lines, we use the ‘inhibited Smad4 nuclear import’ mutation (Table 2) with a mutation factor MF = 0.1. The root-mean-square deviation (RMSD) between the model results and the experimental data is 0.13, 0.18, 0.16, 0.22, 0.26, and 0.18 for panels A, B, C, D, E, and F, respectively.
Table 2 Effects of cancerous TGF-b pathway mutations and their implementation in the computational model. For each mutation, the affected parameters are multiplied by a mutation factor to capture the decreased synthesis or association rates
Mutation TGF-b receptor downregulation Smad4 downregulation Smad7 downregulation Inhibited Smad4 nuclear import Inhibited phospho-RSmad–Smad4 association
Disease
Ref.
Model alteration 2
Bladder, breast, esophageal, ovarian, and prostate cancer Breast cancer Colon cancer Pancreatic cancer
Levy & Hill
Decrease receptor synthesis rate (ksyn,RII and ksyn,RI)
Levy & Hill2 Levy & Hill2 ´n et al.50 More
Pancreatic cancer
Hata et al.51
Decrease Smad4 synthesis rate (ksyn,S4) Decrease Smad7 synthesis rate (ksyn,S7, klip,1, and klip,2) Decrease nuclear import rate constant of monomeric S mad4 and phospho-R-Smad–Smad4 complexes (k7imp and k14imp) Decrease association rate constant for cytosolic and nuclear phospho-R-Smad with Smad4 (k6a and k10a)
dynamics as the latter also exhibits a mesenchymal phenotype.46 The model captures the more transient signal in both nuclear phosphorylated Smad2 (Fig. 2C) and nuclear Smad4 (Fig. 2F) species in PT45 cells as compared to C2C12 dynamics without the mutation (Fig. 2A and D), as observed experimentally. To demonstrate that the model reproduces the observed experimental behavior in the presence of inhibitory drugs and treatments, we simulate the model upon stimulation with TGF-b and treatment with endocytosis inhibitors (Fig. 3 and Fig. S2, ESI†) and with a phosphatase inhibitor (Fig. 4). In Fig. 3, we show the results of simulating the effect of potassium depletion (–KCl) to inhibit clathrin-mediated endocytosis and of Nystatin to inhibit lipid raft-caveolar endocytosis on the dynamics of the ratio of phosphorylated Smad2 to total Smad2 species in Mv1Lu cells. In this case, we use the parameters for HaCaT cells in the model as both cells exhibit an epithelial phenotype. To simulate the effects of potassium depletion, we decreased the internalization and constitutive degradation rate constants for all receptor species
This journal is © The Royal Society of Chemistry 2014
(kdeg,RII, kdeg,RI, k3int, k16deg, k18int, and k19int) with an intervention factor of IF = 0.15. This results in a decreased phosphorylated Smad2 response upon TGF-b stimulation as it minimizes the pool of internalized ligand–receptor complexes to bind to Smad2 (Fig. 3B). To simulate the effects of Nystatin, we decrease the rate constants for Smad7–ligand–receptor complex association (k20a,1 and k20a,2) with an intervention factor of IF = 0.85. In contrast to the strong effects observed for potassium depletion, Nystatin treatment results just in a slight increase in the phosphorylated Smad2 concentration due to the inhibition of the Smad7-mediated negative feedback loop (Fig. 3C). The effects of phosphatase inhibition through treatment with sodium orthovanadate in BAEC cells are simulated by decreasing the dephosphorylation rate constants (k8dp and k11dp) with an intervention factor of IF = 0.3. The model reproduces a slightly more sustained phosphorylated Smad1 response to TGF-b in the presence of sodium orthovanadate as compared to normal signaling dynamics in this cell line (Fig. 4) as it eliminates one mode of negative regulation in the pathway.
Mol. BioSyst., 2014, 10, 537--548 | 541
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Molecular BioSystems
Fig. 3 Simulation of the effects of endocytosis inhibitors on TGF-b signaling dynamics. Using the parameters for HaCaT cells, we simulate the dynamics (lines) of the system and show the results for the ratio of phosphorylated Smad2 to total Smad2 upon stimulation with 0.5 nM (12.8 ng mL1) TGF-b at t = 0 hours. Experimental Western blot data (symbols) and quantification are from ref. 53 with Mv1Lu cells. The experimental data and simulation results are normalized to the maximum from all three conditions. (A) The control case with no inhibitory treatment. (B) Treatment with potassium depletion (–KCl). To simulate the effects of potassium depletion, which inhibits clathrin-mediated endocytosis, we decrease internalization and constitutive degradation rate constants for all receptor species (kdeg,RII, kdeg,RI, k3int, k16deg, k18int, and k19int) by an intervention factor IF = 0.15. (C) Treatment with Nystatin. To simulate the effects of Nystatin, which inhibits lipid raft-caveolar endocytosis, we decrease the rate constants for Smad7–ligand–receptor complex association (k20a,1 and k20a,2) with an intervention factor IF = 0.85. The root-mean-square deviation between the model results and the experimental data is 0.19, 0.06, and 0.06 for panels A, B, and C, respectively.
Fig. 4 Simulation of the effects of phosphatase inhibitors on TGF-b signaling dynamics. Using the parameters for BAEC cells, we simulate the dynamics (lines) of the system and show the results for phosphorylated Smad1 upon stimulation with 1 ng mL1 TGF-b at t = 0 hours. Experimental Western blot data (symbols) are from BAEC cells. (A) The control case with no inhibitory effects. Experimental data are from ref. 15 and 54 and quantified in ref. 15. (B) Treatment with sodium orthovanadate. To simulate the effects of sodium orthovanadate, a phosphatase inhibitor, we decrease the dephosphorylation rate constants (k8dp and k11dp) with an intervention factor IF = 0.3. Experimental data are from ref. 54 and quantified using ImageJ 10.2.52 The root-mean-square deviation between the model results and the experimental data is 0.18 and 0.07 for panels A and B, respectively.
Computational approach for identification of potential therapeutic targets As the TGF-b pathway model is capable of reproducing the effects of both mutations and treatment with inhibitory drugs observed experimentally, we use it as the foundation with which to demonstrate the applicability of the novel in silico approach to identify potential therapeutic targets in intracellular networks. The different steps of the computational approach, summarized in Fig. 5, are discussed in the following sections.
Step 1: identification of targetable network processes. Several strategies have been developed to specifically target biochemical processes that may be used here. Antisense oligonucleotides are used to specifically inhibit protein synthesis by binding to mRNA and preventing translation.47 In the modeled pathway, these could be used to inhibit synthesis of type II receptors, type I receptors, R-Smads, Smad4, and Smad7. Peptide aptamers are engineered to inhibit protein–protein or protein–DNA association by binding to the associated region of
Fig. 5 Flowchart of the approach used to identify potential therapeutic targets. See text (Results section) for additional details on the four steps of the approach.
542 | Mol. BioSyst., 2014, 10, 537--548
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
the protein, competing for the site with the normal binding partner.48 We envision that peptide aptamers could also be used to inhibit nuclear translocation by masking the sites on the Smad proteins responsible for interacting with the nuclear pore complex or accessory proteins mediating nuclear translocation.49 Finally, kinase inhibitors have already been developed to target the ligand–receptor complexes and inhibit phosphorylation of R-Smads.9 From these strategies, we have identified a set of potential targetable processes in the TGF-b pathway and their associated rate constants that would be decreased upon administration of the inhibitory treatment (Table 1). Step 2: sensitivity analysis. We use a global sensitivity analysis (see Methods section) in order to determine how the model responds to perturbation of its parameters. Processes with high sensitivity coefficients will be the most responsive to molecular intervention; therefore, we use sensitivity analysis as a guiding tool for identifying targets in the pathway. We perform the sensitivity analysis with each of the three cell lines mutated with one of the five mutations identified in cancer cells (Table 2), establishing 15 cell-type-mutation combinations, in addition to the 3 non-mutated cases. Each mutation is simulated with a mutation factor of MF = 0.1. In Fig. 6, we show a sample of these results, concretely those for the ‘receptor downregulation’ mutation. The results for the global sensitivity coefficients of the normal cells are included as a control to compare with the effects of the mutation. In general, the mutation decreases the sensitivity coefficients of each parameter as compared to those for the normal cell types, suggesting an increased difficulty in perturbing processes in cells with this mutation. Step 3: sensitive target identification. In order to identify targetable processes that may have the most significant effect on the signaling dynamics, we define as filtering criteria a lower bound of 0.01 for the global sensitivity coefficients Gki. From this, each of the 15 cell-type-mutation combinations will have a set of sensitive parameters with values satisfying the filtering criteria (Gki Z 0.01). To illustrate these results, we color the parameters in Fig. 6 in three groups. Black parameters were not identified as specifically targetable (i.e. not found in Table 1) and are not considered beyond this step. Those in red represent parameters with sensitivity coefficients satisfying the filtering criteria for at least one studied cell type for the considered mutation (H-, B-, or C-; denoting mutated HaCaT, BAEC, or C2C12 cells, respectively). It is important to note that the sensitivity analysis assesses single parameter perturbations, but several of the potential molecular interventions considered in Table 1 affect multiple parameters. Therefore, if at least one parameter from the set of parameters associated with a molecular intervention in Table 1 is selected at this step, we consider the molecular intervention as potentially effective and it passes the filtering procedure to the next step in the approach. Parameters are colored red if within a set of parameters for a particular intervention at least one satisfies the filtering criteria. Those parameters in green correspond to the rest of cases. This is illustrated in the parameter coloring of Fig. S3 (ESI†).
This journal is © The Royal Society of Chemistry 2014
Paper
Step 4: molecular intervention optimization. Using the set of molecular interventions filtered from the sensitivity analysis for each cell-type-mutation combination, we optimize the intervention factor for each case in order to restore normal signaling dynamics to the mutated system, while maintaining the signaling dynamics of the normal system in the presence of the molecular interventions. Fig. 7 displays the results of the optimization routine for all the considered cell types, mutations, and molecular interventions. The ‘No intervention’ case is included as a control in Fig. 7, where we evaluate the objective function (eqn (3)) without any treatment added to the normal or mutated systems. Successful molecular interventions will minimize the objective function value as compared to the ‘No intervention’ case for a given cell type and mutation. Mutations resulting in low values of the objective function signify that the signaling dynamics are robust to the perturbations introduced by that mutation in the specific cell type. In these cases, there is no need of targeting the system with an inhibitor because the mutated cell exhibits similar signaling dynamics as the corresponding normal cell. Therefore, any treatment that reduces the objective function will only result in minimal change from the mutated state as its dynamics already closely mirrors that of the normal state, as illustrated in Fig. 8A for the case of HaCaT cells. By using this optimization routine, rather than relying on the sensitivity analysis alone, we are able to ensure that targeting a particular process is capable of restoring normal signaling dynamics. For example, with the BAEC case affected by the ‘Smad7 downregulation’ mutation, the R-Smad–Smad4 association inhibition was selected as a potential molecular intervention from the sensitivity analysis results (Fig. 7B). However, this mutation results in a more sustained signaling response as compared to that of the normal system and inhibiting R-Smad–Smad4 association increases the signal duration further, moving in the opposite direction than desired (Fig. 8B). Similarly, when the C2C12 case with the ‘inhibited Smad4 nuclear import’ mutation is treated with an R-Smad synthesis inhibitor, the dynamics are unaffected as compared to those with the mutation alone (Fig. 8C). In Fig. 9, we show a selection of the best optimization results for each cell type, which minimize the objective function (eqn (3)) in cases where the mutation has a substantial effect on the normal signaling dynamics. Furthermore, these results minimally impact the signaling dynamics of normal cells in the presence of the molecular intervention, reducing the risk of potential side effects. In HaCaT cells, we show the effects of ‘Smad7 downregulation’ identified in colon cancer (Table 2), which results in a more sustained response to TGF-b for the nuclear phosphorylated Smad2–Smad4 species than without the mutation (Fig. 9A). By inhibiting type II receptor synthesis with an optimized intervention factor of OIF = 0.063, the model dynamics of the mutated case is restored to that of the normal case. For the BAEC case, the ‘inhibited phospho-R-Smad–Smad4 association’ mutation identified in pancreatic cancer (Table 2) results in a more permanent response to ligand stimulation for the nuclear phosphorylated Smad1–Smad4 species (Fig. 9B). By treating the system with Smad4 nuclear export inhibition at the
Mol. BioSyst., 2014, 10, 537--548 | 543
View Article Online
Molecular BioSystems
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Fig. 6 Sensitivity analysis for the ‘receptor downregulation’ mutation. We plot the results obtained for the global sensitivity coefficient Gki (eqn (2)) for each model parameter ki upon stimulation with 45 000 molecules (1.9 ng mL1) of TGF-b. Each column represents the result for one species (pS1S4n or pS2S4n) and one cell type. Results for normal HaCaT, BAEC, and C2C12 cells are denoted as H, B, and C, respectively. The mutated cells, with a mutation factor of MF = 0.1, are abbreviated by appending a dash symbol (–) to the cell type character. Rows represent the results for each parameter, grouped by the module in the pathway it represents. The heat map, generated with Matplotlib 1.1.0,55 represents the value of log10(Gki ) as the color intensities in the color key (righthand side of the heat map); empty spaces indicate a sensitivity coefficient equal to zero or a parameter value set to zero from Table S2 (ESI†). Parameter names in black are not specifically targetable, while those colored red or green represent processes that may be specifically targeted (Table 1). Those in red represent parameters satisfying the filtering criteria described in the Results section, while those in green do not satisfy the filtering criteria.
optimized intervention factor of OIF = 0.037, the dynamics is restored to its normal transient response, while minimally affecting the normal cell case. For the C2C12 case, we mutate the system with the ‘Smad7 downregulation’ mutation, which
544 | Mol. BioSyst., 2014, 10, 537--548
increases the duration of the response to TGF-b stimulation for the nuclear phosphorylated Smad1–Smad4 species (Fig. 9C). By inhibiting R-Smad–Smad4 association with an optimized intervention factor of OIF = 0.132, the mutated-cell dynamics returns
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Fig. 7 Results of the optimization of the molecular interventions. Bars represent the objective function value (eqn (3)) using one species, mutation, and molecular intervention for (A) HaCaT, (B) BAEC, and (C) C2C12 cells. Bars are grouped by the type of mutation and color coded with the molecular intervention used. Dots beneath an empty bar indicate that the molecular intervention was not selected as potentially effective from the sensitivity analysis step.
to the transient response from the normal case, while minimally impacting the normal C2C12 signaling dynamics. Finally, we filter the results of the molecular interventions optimization routine to produce a list of potential treatments for each cell-type-mutation combination (Table 3). For each cell-type-species combination, we first determine the set of mutations with objective functions for the ‘No intervention’
This journal is © The Royal Society of Chemistry 2014
case greater than or equal to 500 for each cell type. With this, we then determine the set of molecular interventions that minimize the objective function below 0.3 times the ‘No intervention’ value from the set of mutations for each cell-type-species combination. Although the HaCaT dynamics is only substantially responsive to the ‘Smad7 downregulation’ mutation, this method identified two potential molecular interventions to restore normal
Mol. BioSyst., 2014, 10, 537--548 | 545
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Molecular BioSystems
Fig. 8 Simulation results of mutations and molecular interventions with minimal or negative effects on normal signaling dynamics. Solid blue lines represent simulations of normal cells and solid red lines represent simulations of mutated cells with a mutation factor of MF = 0.1 upon stimulation with 45 000 molecules (1.9 ng mL1) of TGF-b at t = 0 hours. Dashed lines indicate simulations with the molecular interventions. (A) HaCaT cells with the ‘Smad4 downregulation’ mutation. (B) BAEC cells with the ‘Smad7 downregulation’ mutation and treated with an R-Smad–Smad4 association inhibitor using an intervention factor of IF = 0.1. (C) C2C12 cells with the ‘inhibited Smad4 nuclear import’ mutation and treated with an R-Smad synthesis inhibitor using an intervention factor of IF = 0.1.
Fig. 9 Simulation results for optimal molecular interventions. Same conditions as in Fig. 8 for combinations of species, mutation, and molecular interventions selected from the optimization routine for each cell type. (A) HaCaT cells with the ‘Smad7 downregulation’ mutation and treated with a type II receptor synthesis inhibitor using the optimized intervention factor OIF = 0.063. (B) BAEC cells with the ‘inhibited phospho-R-Smad–Smad4 association’ mutation and treated with a Smad4 nuclear export inhibitor using the optimized intervention factor OIF = 0.037. (C) C2C12 cells with the ‘Smad7 downregulation’ mutation and treated with an R-Smad–Smad4 association inhibitor using the optimized intervention factor OIF = 0.132.
Table 3 Potential targetable processes identified by optimization of the molecular interventions. For each cell type and mutations that significantly affect the normal signaling dynamics, we show the corresponding molecular interventions that restore normal signaling dynamics to the mutated system. The entries in boldface are those used in Fig. 9
HaCaT Mutation
Targeted process
Smad4 downregulation
Smad7 downregulation
BAEC Species Targeted process Type II receptor synthesis Type I receptor synthesis Smad4 nuclear export
Type II receptor synthesis pS2S4n Type I receptor synthesis R-Smad–Smad4 pS2S4n Smad4 nuclear export association
C2C12 Species Targeted process pS1S4n pS1S4n pS1S4n pS1S4n Type I receptor synthesis
pS1S4n
pS1S4n Smad4 nuclear import
pS1S4n
R-Smad–Smad4 association Inhibited Smad4 nuclear import
Inhibited R-Smad–Smad4 association
546 | Mol. BioSyst., 2014, 10, 537--548
Type I receptor synthesis R-Smad nuclear export Smad4 nuclear export
Species
pS1S4n
pS1S4n pS1S4n pS1S4n
Type I receptor synthesis pS1S4n R-Smad nuclear export pS1S4n Smad4 nuclear export pS1S4n Receptor–R-Smad association
pS2S4n
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
signaling dynamics to the mutated system: inhibited type II receptor synthesis and inhibited R-Smad–Smad4 association. For the BAEC case, this method identified a number of molecular interventions for four of the five mutations considered here, targeting receptor synthesis and R-Smad/Smad4 nucleocytoplasmic shuttling. The results for the C2C12 case reveal four distinct molecular interventions that may be used to target the system for two of the mutations we used in this analysis. Finally, our approach did not identify any molecular intervention for the ‘TGF-b receptor downregulation’ mutation in any of the considered cases.
Conclusions Effectively targeting the TGF-b signal transduction pathway requires a detailed dissection of the network interactions beyond insight gained from intuition alone. Here we have developed a strategy based on computational modeling to quantitatively analyze the signaling pathway in order to determine potential therapeutic targets. We have considered a model that captures the precise dynamics of the system, mutations that affect the system parameters, and a collection of potentially targetable components of the pathway, such as inhibition of association or synthesis of proteins. The key element of the approach is the identification of molecular interventions in the pathway that, in addition to bringing the mutated system back to its normal dynamics, do not substantially affect the non-mutated system. To do so, our approach uses global sensitivity analysis to identify parameters with high sensitivity coefficients and passes the molecular interventions containing these sensitive parameters to the optimization routine that selects the optimal strength of the molecular intervention. We have shown the capabilities of the computational approach for different cell types by reproducing experimental TGF-b signaling dynamics when the system is mutated in pancreatic cancer cells and in the presence of inhibitory drugs targeting distinct processes in the pathway. From the results of the molecular intervention optimization routine, we have identified a set of mutations that significantly affect the signaling dynamics for each cell type and a number of molecular interventions that may be used to effectively target the effects of these mutations. Interestingly, our approach identified multiple types of molecular interventions that can be used for a given mutation. This type of approach can generally be applied to a wide variety of signaling systems with characterized biochemical interactions and mutations affecting their components to identify novel targets for therapeutic intervention. Therefore, our methodology provides an efficient starting point to select potential targets for further experimental and clinical investigation.
Acknowledgements This work was supported by the University of California, Davis (to LS).
This journal is © The Royal Society of Chemistry 2014
Paper
References ´, Annu. Rev. Biochem., 1998, 67, 753–791. 1 J. Massague 2 L. Levy and C. S. Hill, Cytokine Growth Factor Rev., 2006, 17, 41–58. 3 P. ten Dijke and H. M. Arthur, Nat. Rev. Mol. Cell Biol., 2007, 8, 857–869. 4 G. C. Blobe, W. P. Schiemann and H. F. Lodish, N. Engl. J. Med., 2000, 342, 1350–1358. 5 R. Derynck and K. Miyazono, The TGF-b Family, Cold Spring Harbor Laboratory Press, Cold Spring Harbor, N.Y., 2008. 6 Y. Drabsch and P. ten Dijke, Cancer Metastasis Rev., 2012, 31, 553–568. ´, Cell, 2003, 113, 685–700. 7 Y. Shi and J. Massague ´, Cell, 2008, 134, 215–230. 8 J. Massague 9 J. M. Yingling, K. L. Blanchard and J. S. Sawyer, Nat. Rev. Drug Discovery, 2004, 3, 1011–1022. 10 J. M. G. Vilar, R. Jansen and C. Sander, PLoS Comput. Biol., 2006, 2, e3. 11 J. M. G. Vilar and L. Saiz, Biophys. J., 2011, 101, 2315–2323. 12 V. Becker, J. Timmer and U. Klingmuller, Adv. Exp. Med. Biol., 2012, 736, 313–323. 13 D. C. Clarke, M. D. Betterton and X. Liu, IEE Proc.: Syst. Biol., 2006, 153, 412–424. 14 B. Schmierer, A. L. Tournier, P. A. Bates and C. S. Hill, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 6608–6613. ¨nsson, E. Pardali, P. ten Dijke and 15 P. Melke, H. Jo C. Peterson, Biophys. J., 2006, 91, 4368–4380. 16 M. Paulsen, S. Legewie, R. Eils, E. Karaulanov and C. Niehrs, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10202–10207. 17 K. Wegner, A. Bachmann, J. U. Schad, P. Lucarelli, S. Sahle, ¨ller, S. Dooley and P. Nickel, C. Meyer, U. Klingmu U. Kummer, Biophys. Chem., 2012, 162, 22–34. 18 G. Celliere, G. Fengos, M. Herve and D. Iber, BMC Syst. Biol., 2011, 5, 184. 19 D. Nicklas and L. Saiz, J. R. Soc. Interface, 2013, 10, 20130363. 20 J. Ho and L. Saiz, Biophys. J., 2011, 100, 164a. 21 Z. Zi and E. Klipp, PLoS One, 2007, 2, e936. 22 S. W. Chung, F. L. Miles, R. A. Sikes, C. R. Cooper, M. C. Farach-Carson and B. A. Ogunnaike, Biophys. J., 2009, 96, 1733–1750. 23 Z. Zi, Z. Feng, D. A. Chapnick, M. Dahl, D. Deng, E. Klipp, A. Moustakas and X. Liu, Mol. Syst. Biol., 2011, 7, 492. 24 B. S. Hendriks, G. J. Griffiths, R. Benson, D. Kenyon, M. Lazzara, J. Swinton, S. Beck, M. Hickinson, J. M. Beusmans, D. Lauffenburger and D. de Graaf, IEE Proc.: Syst. Biol., 2006, 153, 457–466. 25 R. J. Orton, M. E. Adriaens, A. Gormand, O. E. Sturm, W. Kolch and D. R. Gilbert, BMC Syst. Biol., 2009, 3, 100. 26 S. W. Chung, C. R. Cooper, M. C. Farach-Carson and B. A. Ogunnaike, J. R. Soc. Interface, 2012, 9, 1389–1397. 27 N. Kumar, B. S. Hendriks, K. A. Janes, D. de Graaf and D. A. Lauffenburger, Drug Discovery Today, 2006, 11, 806–811. 28 L. H. Abbott and F. Michor, Br. J. Cancer, 2006, 95, 1136–1141.
Mol. BioSyst., 2014, 10, 537--548 | 547
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
29 K. Oda, Y. Matsuoka, A. Funahashi and H. Kitano, Mol. Syst. Biol., 2005, 1, 20050010. 30 G. Lebedeva, A. Sorokin, D. Faratian, P. Mullen, A. Goltsov, S. P. Langdon, D. J. Harrison and I. Goryanin, Eur. J. Pharm. Sci., 2012, 46, 244–258. 31 R. P. Araujo, E. F. Petricoin and L. A. Liotta, BioSystems, 2005, 80, 57–69. 32 M. H. Sung and R. Simon, Mol. Pharmacol., 2004, 66, 70–75. 33 H. Yan, B. Zhang, S. Li and Q. Zhao, BMC Syst. Biol., 2010, 4, 50. 34 B. J. Lao, W. L. Tsai, F. Mashayekhi, E. A. Pham, A. B. Mason and D. T. Kamei, J. Controlled Release, 2007, 117, 403–412. 35 D. J. Yoon, D. S. Chu, C. W. Ng, E. A. Pham, A. B. Mason, D. M. Hudson, V. C. Smith, R. T. MacGillivray and D. T. Kamei, J. Controlled Release, 2009, 133, 178–184. 36 D. J. Yoon, B. H. Kwan, F. C. Chao, T. P. Nicolaides, J. J. Phillips, G. Y. Lam, A. B. Mason, W. A. Weiss and D. T. Kamei, Cancer Res., 2010, 70, 4520–4527. 37 D. Hanahan and R. A. Weinberg, Cell, 2000, 100, 57–70. 38 D. Hanahan and R. A. Weinberg, Cell, 2011, 144, 646–674. 39 A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, ACM Trans. Math. Softw., 2005, 31, 363–396. 40 Z. Zi, IET Syst. Biol., 2011, 5, 336–346. 41 S. Kucherenko, M. Rodriguez-Fernandez, C. Pantelides and N. Shah, Reliab. Eng. Syst. Saf., 2009, 94, 1135–1148. 42 A. Varma, M. Morbidelli and H. Wu, Parametric sensitivity in chemical systems, Cambridge University Press, Cambridge, UK; New York, NY, 1999. 43 W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical recipes: the art of scientific computing, Cambridge University Press, New York, 2007. 44 E. Jones, T. Oliphant, P. Peterson and others, SciPy: Open source scientific tools for Python (http://www.scipy.org), 2001–. 45 F. James and M. Winkler, MINUIT user’s guide, CERN, Geneva, 2004. ´s and C. S. Hill, Oncogene, 2003, 22, 3698–3711. 46 F. J. Nicola
548 | Mol. BioSyst., 2014, 10, 537--548
Molecular BioSystems
47 C. F. Bennett and E. E. Swayze, Annu. Rev. Pharmacol. Toxicol., 2010, 50, 259–293. 48 A. Conidi, V. van den Berghe and D. Huylebroeck, Int. J. Mol. Sci., 2013, 14, 6690–6719. 49 C. S. Hill, Cell Res., 2009, 19, 36–46. ´n, S. Itoh, A. Moustakas, P. t. Dijke and 50 A. More C.-H. Heldin, Oncogene, 2000, 19, 4396–4404. 51 A. Hata, R. S. Lo, D. Wotton, G. Lagna and J. Massague, Nature, 1997, 388, 82–87. 52 M. D. Abramoff, P. J. Magalhaes and S. J. Ram, Biophoton. Int., 2004, 11, 36–42. 53 G. M. Di Guglielmo, C. Le Roy, A. F. Goodfellow and J. L. Wrana, Nat. Cell Biol., 2003, 5, 410–421. 54 G. Valdimarsdottir, M.-J. Goumans, F. Itoh, S. Itoh, C.-H. Heldin and P. ten Dijke, BMC Cell Biol., 2006, 7, 16. 55 J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95. ´, Nat. Cell Biol., 1999, 1, 472–478. 56 R. S. Lo and J. Massague 57 M. Wan, X. Cao, Y. Wu, S. Bai, L. Wu, X. Shi, N. Wang and X. Cao, EMBO Rep., 2002, 3, 171–176. 58 W. Liu, H. Rui, J. Wang, S. Lin, Y. He, M. Chen, Q. Li, Z. Ye, S. Zhang, S. C. Chan, Y.-G. Chen, J. Han and S.-C. Lin, EMBO J., 2006, 25, 1646–1658. 59 P. R. Segarini, D. M. Rosen and S. M. Seyedin, Mol. Endocrinol., 1989, 3, 261–272. 60 J.-F. Goetschy, O. Letourneur, N. Cerletti and M. A. Horisberger, Eur. J. Biochem., 1996, 241, 355–362. 61 S.-B. Peng, L. Yan, X. Xia, S. A. Watkins, H. B. Brooks, D. Beight, D. K. Herron, M. L. Jones, J. W. Lampe, W. T. McMillen, N. Mort, J. S. Sawyer and J. M. Yingling, Biochemistry, 2005, 44, 2293–2304. 62 M. Funaba and L. S. Mathews, Mol. Endocrinol., 2000, 14, 1583–1591. 63 X. Lin, X. Duan, Y.-Y. Liang, Y. Su, K. H. Wrighton, J. Long, M. Hu, C. M. Davis, J. Wang, F. C. Brunicardi, Y. Shi, Y.-G. Chen, A. Meng and X.-H. Feng, Cell, 2006, 125, 915–928. 64 B. Schmierer and C. S. Hill, Mol. Cell. Biol., 2005, 25, 9845–9858.
This journal is © The Royal Society of Chemistry 2014
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 537
View Journal | View Issue
In silico identification of potential therapeutic targets in the TGF-b signal transduction pathway† Daniel Nicklas and Leonor Saiz* The transforming growth factor-b (TGF-b) superfamily of cytokines controls fundamental cellular processes, such as proliferation, motility, differentiation, and apoptosis. This fundamental role is emphasized by the widespread presence of mutations of the core components of the TGF-b signal transduction pathway in a number of human diseases. Therefore, there is an increasing interest in the development of therapies to specifically target this pathway. Here we develop a computational approach
Received 30th June 2013, Accepted 10th December 2013
to identify potential intervention points that are capable of restoring the normal signaling dynamics to
DOI: 10.1039/c3mb70259f
this approach explicitly to the TGF-b pathway to study the signaling dynamics of mutated and normal
the mutated system while maintaining the behavior of normal cells substantially unperturbed. We apply cells treated with inhibitory drugs and identify the processes in the pathway that are most susceptible to
www.rsc.org/molecularbiosystems
therapeutic intervention.
Introduction The transforming growth factor-b (TGF-b) signal transduction pathway controls multiple cellular processes, including proliferation, motility, differentiation, and apoptosis, that are essential for regulating embryonic development, tissue homeostasis, and wound healing.1 Dysregulation of the pathway through mutations of its core components is associated with a number of human diseases, such as cancer, fibrotic conditions of kidneys, liver, and lungs, and vascular and autoimmune disorders.2–4 Because of its involvement in many aspects of human health, a significant effort of clinical research focuses on developing therapies targeting the activity of the TGF-b signaling pathway.5,6 TGF-b-superfamily ligands signal through a network of interacting molecular components with complex patterns of coupled signaling and feedback regulation. At the plasma membrane, the ligands initiate signaling by binding to two types of serine–threonine kinases, known as type I and type II receptors. The active ligand–receptor complex propagates the signal through the intracellular mediator Smad proteins upon internalization into the endosome, where it may recruit and phosphorylate a receptor-regulated Smad (R-Smad). Phosphorylated R-Smads bind to Smad4 and translocate into the nucleus, where the resulting macromolecular complex regulates the
Modeling of Biological Networks Laboratory, Department of Biomedical Engineering, University of California, 451 East Health Sciences Drive, Davis, CA 95616, USA. E-mail: [email protected]; Fax: +1-530-754-5739; Tel: +1-530-752-6700 † Electronic supplementary information (ESI) available. See DOI: 10.1039/ c3mb70259f
This journal is © The Royal Society of Chemistry 2014
expression of several hundred genes, including Smad7. Smad7, an inhibitory Smad, negatively regulates signaling by binding to active ligand–receptor complexes at the plasma membrane and inducing their degradation. TGF-b-superfamily ligands use two parallel R-Smad channels, grouped as Smad1/5/8 and Smad2/3, to transmit the signal. Bone morphogenetic proteins (BMPs) and nodal/activin ligands induce phosphorylation of the Smad1/5/8 and Smad2/3 groups, respectively, while TGF-b induces phosphorylation of both R-Smad groups.7 A number of mutations affecting the TGF-b pathway signaling components have been identified in a variety of human cancers.2,6 While TGF-b is typically a growth suppressor, its dysregulated activation may promote invasion and growth of tumors in cancer cells.8 Therefore, a number of strategies have been developed to abrogate this effect by inhibiting TGF-b activity, such as the use of antisense oligonucleotides and monoclonal antibodies targeting TGF-b mRNA and protein, respectively, and small molecule kinase inhibitors that target the receptors, blocking R-Smad phosphorylation.6,9 Quantitative and predictive models have been successfully used to study key modules in the TGF-b signal transduction pathway, including receptor trafficking,10–12 Smad nucleocytoplasmic shuttling,13,14 and transcriptional feedback loops.15–20 Several models have combined these components into comprehensive mathematical representations of the signaling pathway.19,21–23 In addition to providing insight into the underlying mechanisms, computational models of signaling pathways have been used to study the effects of mutations in disease.22,24–26 In the TGF-b signaling pathway, for instance, Chung et al. investigated the effects of cancerous mutations on Smad-dependent signaling. Specifically, they used their model to predict how different
Mol. BioSyst., 2014, 10, 537--548 | 537
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
signaling components are affected by cancerous mutations by introducing a mutation resulting in downregulation of TGF-b receptors through decreased initial concentration and synthesis rate of the affected species.22 Furthermore, it has been shown that computational modeling can be used as a fundamental step in the development of novel pharmaceuticals.27,28 A prime example of the applicability of computational approaches is the epidermal growth factor receptor (EGFR) signaling pathway, which has been extensively studied using a number of computational models to investigate both receptor trafficking and downstream signaling.29 In particular, models of the EGF-induced extracellular signal regulated kinase (ERK) signal transduction pathway have been successfully used to study the effects of oncogenic mutations,25 to assess the role of mutation in the sensitivity of cells to an existing therapy,24 and to identify therapeutic targets and biomarkers.30 Specifically, Orton and coworkers developed a model of the ERK signal transduction pathway and compared the signaling dynamics of normal and cancer cells as modeled through knockouts and mutations of signaling proteins. From this, the authors hypothesized which components of the signaling pathway were susceptible to therapeutic intervention.25 Using similar strategies, Hendriks et al. introduced perturbations into their model to represent the effects of mutant receptors and tested potential mechanisms by which the mutated signal dynamics observed experimentally differed from those of normal cells.24 Using sensitivity analysis, Lebedeva and coworkers identified a number of potential therapeutic targets and biomarkers in the pathway.30 In addition to investigating the effects of mutations, computational modeling has been used to analyze how therapies interact with mutated signaling pathways and to propose novel targets for therapeutic intervention in a number of systems and diseases.31–36 Specifically, Sung and Simon compared the effects of novel inhibitor drugs on the signaling dynamics of the nuclear factor-kB (NF-kB) pathway.32 Yan et al. and Araujo et al. extended this analysis in the NF-kB and EGFR pathways, respectively, to simulations of drug combinations.31,33 In particular, Araujo and coworkers evaluated the effects of kinase inhibitors targeting one or two nodes in the pathway, revealing that simultaneous administration allows for a reduction in the required dose necessary to significantly affect the signal response.31 Mathematical modeling of transferrin binding and receptor trafficking was used to successfully propose a novel mechanism to improve its drug delivery capabilities,34 which was validated with in vitro and in vivo experiments.34–36 Here we apply a comprehensive computational model of the TGF-b signal transduction pathway19 to study the signaling dynamics of mutated and normal cells treated with inhibitory drugs or treatments. Using this model, we build upon strategies applied to the TGF-b signaling pathway and other systems to develop an effective in silico approach for identifying targetable processes and potential intervention points that are capable of restoring the normal signaling dynamics to the mutated system while maintaining the behavior of normal cells substantially unperturbed.
538 | Mol. BioSyst., 2014, 10, 537--548
Molecular BioSystems
We have chosen the TGF-b signaling system because of its prominent role in cancer.8 Cancer is a multistage process that involves the successive failure of many control points, such as cancer cells stimulating their own growth and resisting inhibitory signals that might otherwise stop their growth as well as the evasion of programmed cell death. Therefore, dysregulation of the TGF-b pathway in cancer cell lines is accompanied by the dysregulation of additional pathways, such as those arising from K-ras oncogene mutations.37,38 There are many therapies under development that target the components of the TGF-b pathway.6 One of the most sought features of a therapy is the selective targeting of the mutated cells while keeping the normal cells largely unperturbed. So far there has not been a systematic approach to identify those therapies. The in silico approach we present here is a first step in this direction, which can be applied to a number of other signaling systems as a guiding tool in the therapy development process. Subsequent aspects that would need to be taken into account in the full therapy development include the differences between cell lines and tumor samples as well as the cross-talking between different signaling pathways.
Methods The TGF-b pathway model To study the dynamics of the pathway upon stimulation with TGF-b, we use the comprehensive mathematical model developed in ref. 19, which has been shown to accurately reproduce the behavior of human keratinocytes (HaCaT), bovine aortic endothelial cells (BAEC), and the mouse mesenchymal C2C12 cell line. The model considers explicitly receptor trafficking, downstream signaling by intracellular Smad proteins, and negative feedback through Smad7. Signaling is initiated when a TGF-b ligand binds to its type II receptor, denoted as RII. This complex then recruits a type I receptor, one of RI1 or RI2, forming the ligand–receptor complex C1 or C2, respectively. The ligand–receptor complexes propagate the signal upon internalization into the endosome, where C1 may bind with cytosolic Smad1 (S1c) and C2 may bind with cytosolic Smad2 (S2c) and phosphorylate the R-Smad. We use the subscripts c and n to denote cytosolic and nuclear species, respectively, and apply a prefix p to phosphorylated species. pS1c and pS2c bind to Smad4 (S4c), forming the pS1S4c and pS2S4c complexes, respectively, which may translocate into the nucleus where pS1S4n and pS2S4n then may promote the expression of Smad7 (S7). This initiates a negative feedback loop in which S7 binds to the ligand–receptor complexes (C1 or C2) at the plasma membrane, preventing their association and phosphorylation of R-Smad proteins, and targeting them for degradation. A schematic representation of the model is shown in Fig. 1, where arrows indicate modeled reactions in the pathway. Reactions are mathematically represented using mass-action kinetics, which are then combined to form the system of ordinary differential equations (ODEs) to track the rate of change of each modeled species (Table S1, ESI†). Parameter values for the three cell types are listed in Table S2 (ESI†).10,19,22,56–64
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Paper
Fig. 1 Schematic illustration of the Smad-dependent TGF-b signaling pathway model upon stimulation with TGF-b. Arrows denote reaction steps in the pathway and are labeled with the rate constant for the reaction. Synthesis and degradation reactions are provided in the table (inset) with the notation for the different molecular species given in the text. Expression of Smad7 is under the control of a single gene, activated by nuclear phosphorylated Smad1–Smad4 and/or Smad2–Smad4 complexes.
Prior to ligand stimulation, we determine the steady state solution for the system of ODEs (Table S1, ESI†) by setting each time-derivative to zero and solving the resulting linear system of equations with the concentration of the ligand set to zero using the ‘linalg.solve’ method in Numpy 1.6.2 (http://numpy. scipy.org). Upon ligand stimulation, we numerically solve the system of ODEs using the CVODE method in the SUNDIALS 2.5.0 package.39 These correspond to the typical experimental conditions that measure the response of the system to a sudden change of the TGF-b ligand concentration from zero to a saturating value that is kept constant afterwards.
Mutation and molecular intervention simulations In order to capture the effects of mutations in the pathway, we use a strategy similar to that employed in ref. 22 and 31 to introduce altered protein expression or functional mutation. Focusing on the mutations outlined in Table 1, we define a mutation factor (MF) within the range [0,1] with which to reduce the parameter values associated with each mutation. For example, to introduce the ‘inhibited Smad4 nuclear import’ mutation with a 0.1 mutation factor, we multiply k7imp and k14imp by 0.1. Similarly, we simulate the effects of a molecular intervention by introducing an intervention factor (IF) within
This journal is © The Royal Society of Chemistry 2014
the range [0,1] that is multiplied with the rate constants affected by a particular molecular intervention. Sensitivity analysis We use a global sensitivity analysis to investigate how the model responds to perturbations of individual parameters in the parameter space.40 Concretely, we compute the derivativebased global sensitivity measures developed in ref. 41, which sample the effects of local parameter perturbation in a larger parameter space. To compute the effects of local parameter perturbation, we use the scaled sensitivity coefficients42 given by eki; j ¼
ki; j @Cj ; Cj @ki; j
(1)
where ki,j defines the value of the ith parameter of sample j and Cj is the time integral of the model output over 24 hours after stimulation with the ligand for the parameters of the jth sample. As model output, we examine the signaling dynamics of both nuclear phosphorylated Smad1–Smad4 (pS1S4n) and nuclear phosphorylated Smad2–Smad4 (pS2S4n) complexes since they control the long-term effects of ligand stimulation through regulation of gene expression.7 To approximate the partial derivative, we simulate the model with 0.5 percent
Mol. BioSyst., 2014, 10, 537--548 | 539
View Article Online
Paper
Molecular BioSystems
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Table 1 Potential targetable processes and their associated rate constants analyzed in the optimization of the molecular interventions. The parameters associated with each process are multiplied by an intervention factor in the range [0,1] to simulate the process inhibition
Inhibited process
Parameters
Type II receptor synthesis Type I receptor synthesis R-Smad synthesis Smad4 synthesis Smad7 synthesis Association of TGF-b with type II receptor Association of TGF-b–type II receptor complex with type I receptor Ligand–receptor complex–R-Smad association R-Smad–Smad4 association R-Smad nuclear import R-Smad nuclear export Smad4 nuclear import Smad4 nuclear export Ligand–receptor complex–Smad7 association Association of phospho-R-Smad–Smad4 complex with DNA
ksyn,RII ksyn,RI ksyn,RS ksyn,S4 ksyn,S7, klip,1, and klip,2 k1a k2a k4a k6a and k10a k7imp, k12imp, and k13imp k12exp k7imp and k14imp k14exp k20a,1 and k20a,2 KA,1 and KA,2
perturbations of parameter ki about its value at sample j and calculate the finite central difference.43 For each parameter ki, we calculate the local sensitivity coefficient eki,j for 4000 parameter values, i.e. j = 1,. . .,4000, by randomly sampling around its original value (Table S2, ESI†) in the parameter space. Specifically, we generate a random variable x with a uniform distribution within the range [1,1] and multiply the original parameter value by 10x. The scaled sensitivity coefficients are then used to compute the global sensitivity coefficient41 defined by Gki ¼
N 1X ek 2 ; N j¼1 i; j
(2)
where N is the number of samples, equal to 4000 here. We implemented the sensitivity analysis using Python 2.7.3 (http:// www.python.org), Numpy 1.6.2 (http://numpy.scipy.org), and Scipy 0.10.0.44 Molecular intervention optimization We optimize the intervention factor IF, defined as a value within the range [0,1], by fitting the simulation result of a mutated cell with the specific molecular intervention to that of a normal cell. Additionally, we fit the simulation result of the normal cell with the molecular intervention to that of a normal cell in order to minimize the potential risk of side effects from exposing normal cells to the molecular intervention since the goal of our approach is to identify those molecular interventions that lead to the recovery of normal behavior by mutated cells and that minimally affect the behavior of normal cells. Therefore, we define the objective function for minimization as the least-squares error: E¼
n X
ðSðti Þ M ðti ÞÞ2 þðSðti Þ S ðti ÞÞ2 ;
i¼1
540 | Mol. BioSyst., 2014, 10, 537--548
(3)
where S(ti) and M(ti) are the normal and mutated cell dynamics, respectively, at the time point ti of n total time points and the asterisk represents the corresponding dynamics with the molecular intervention. To determine the optimized intervention factor (OIF) we first evaluate the objective function with intervention factors within the range [0.1,0.9] with 0.1 increments. The intervention factor with the minimum objective function value becomes the starting point for a downhill simplex algorithm with the ‘optimize.fmin’ function in Scipy 0.10.0.44 In order to impose the bounds of [0,1] for the intervention factor, we use the transformations described in ref. 45. Our computational implementation uses Python 2.7.3 (http://www.python.org), Numpy 1.6.2 (http://numpy.scipy.org), and Scipy 0.10.0.44
Results and discussion Model dynamics To examine the ability of the model to reproduce the observed behavior in mutated cell lines, we compare the results of the simulation with experimental data in HaCaT cells and two pancreatic cancer cell lines, namely Colo-357 and PT45 (Fig. 2). The TGF-b signaling in these cell lines has been studied in detail in ref. 46. The results show that Colo-357 cells exhibit an epithelial phenotype, forming well-defined epithelial sheets and expressing E-cadherin at their adherens junctions. Colo357 cells respond to TGF-b in a way very similar to HaCaT cells, strongly inducing the expression of p21, JunB, c-Jun, Smad7, and PAI-1 mRNA within 1 hour of ligand stimulation. PT45 cells, in contrast, do not present the epithelial phenotype and the induction of p21 and JunB mRNA with TGF-b is only transient. Colo-357 cells, as HaCaT cells, display rapid nuclear accumulation of phosphorylated Smad2 that is still present, but reduced 6 hours after activation. In contrast, nuclear accumulation of phosphorylated Smad2 in PT45 cells disappears abruptly after 1 hour of stimulation. The kinetics of nuclear accumulation of Smad4 was observed to be very similar to that of phosphorylated Smad2 in both Colo-357 and PT45 cells. Using the parameters for HaCaT cells (Table S2, ESI†), the model reproduces the normal signaling dynamics of both nuclear phosphorylated Smad2 (Fig. 2A) and nuclear Smad4 (Fig. 2D) molecular species. To simulate the dynamics of pancreatic cancer cells Colo-357 and PT45, we introduce the ‘inhibited Smad4 nuclear import’ mutation found in pancreatic cancer cell lines (Table 2) to the normal cases with a mutation factor of MF = 0.1. As Colo-357 cells exhibit an epithelial phenotype,46 we use the model for HaCaT cells to reproduce their signaling dynamics. By introducing the mutation, the model results change only minimally for both nuclear phosphorylated Smad2 (Fig. 2B and Fig. S1, ESI†) and nuclear Smad4 (Fig. 2E and Fig. S1, ESI†), consistently with the experimental data. The differences between the model and experimental results observed in the early stages of signaling in Fig. 2B and E are attributable, as discussed in ref. 11, to the detailed trafficking dynamics of the receptors. We use the parameters for C2C12 cells to reproduce the PT45
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Paper
Fig. 2 Simulation of the effects of pancreatic cancer mutations on TGF-b signaling dynamics. We simulate the dynamics of the model using the parameters for HaCaT (black solid lines) and C2C12 cells (blue dashed lines) and show the results for nuclear phosphorylated Smad2 (top row; A–C) and nuclear Smad4 (bottom row; D–F) species upon stimulation with 2 ng mL1 TGF-b at t = 0 hours. Experimental data (symbols) was quantified using ImageJ 10.252 from Western blot experiments in ref. 46 for HaCaT cells (left column; A and D) and two pancreatic cancer cell lines: Colo-357 (middle column; B and E) and PT45 (right column; C and F). To simulate the dynamics of the pancreatic cancer cell lines, we use the ‘inhibited Smad4 nuclear import’ mutation (Table 2) with a mutation factor MF = 0.1. The root-mean-square deviation (RMSD) between the model results and the experimental data is 0.13, 0.18, 0.16, 0.22, 0.26, and 0.18 for panels A, B, C, D, E, and F, respectively.
Table 2 Effects of cancerous TGF-b pathway mutations and their implementation in the computational model. For each mutation, the affected parameters are multiplied by a mutation factor to capture the decreased synthesis or association rates
Mutation TGF-b receptor downregulation Smad4 downregulation Smad7 downregulation Inhibited Smad4 nuclear import Inhibited phospho-RSmad–Smad4 association
Disease
Ref.
Model alteration 2
Bladder, breast, esophageal, ovarian, and prostate cancer Breast cancer Colon cancer Pancreatic cancer
Levy & Hill
Decrease receptor synthesis rate (ksyn,RII and ksyn,RI)
Levy & Hill2 Levy & Hill2 ´n et al.50 More
Pancreatic cancer
Hata et al.51
Decrease Smad4 synthesis rate (ksyn,S4) Decrease Smad7 synthesis rate (ksyn,S7, klip,1, and klip,2) Decrease nuclear import rate constant of monomeric S mad4 and phospho-R-Smad–Smad4 complexes (k7imp and k14imp) Decrease association rate constant for cytosolic and nuclear phospho-R-Smad with Smad4 (k6a and k10a)
dynamics as the latter also exhibits a mesenchymal phenotype.46 The model captures the more transient signal in both nuclear phosphorylated Smad2 (Fig. 2C) and nuclear Smad4 (Fig. 2F) species in PT45 cells as compared to C2C12 dynamics without the mutation (Fig. 2A and D), as observed experimentally. To demonstrate that the model reproduces the observed experimental behavior in the presence of inhibitory drugs and treatments, we simulate the model upon stimulation with TGF-b and treatment with endocytosis inhibitors (Fig. 3 and Fig. S2, ESI†) and with a phosphatase inhibitor (Fig. 4). In Fig. 3, we show the results of simulating the effect of potassium depletion (–KCl) to inhibit clathrin-mediated endocytosis and of Nystatin to inhibit lipid raft-caveolar endocytosis on the dynamics of the ratio of phosphorylated Smad2 to total Smad2 species in Mv1Lu cells. In this case, we use the parameters for HaCaT cells in the model as both cells exhibit an epithelial phenotype. To simulate the effects of potassium depletion, we decreased the internalization and constitutive degradation rate constants for all receptor species
This journal is © The Royal Society of Chemistry 2014
(kdeg,RII, kdeg,RI, k3int, k16deg, k18int, and k19int) with an intervention factor of IF = 0.15. This results in a decreased phosphorylated Smad2 response upon TGF-b stimulation as it minimizes the pool of internalized ligand–receptor complexes to bind to Smad2 (Fig. 3B). To simulate the effects of Nystatin, we decrease the rate constants for Smad7–ligand–receptor complex association (k20a,1 and k20a,2) with an intervention factor of IF = 0.85. In contrast to the strong effects observed for potassium depletion, Nystatin treatment results just in a slight increase in the phosphorylated Smad2 concentration due to the inhibition of the Smad7-mediated negative feedback loop (Fig. 3C). The effects of phosphatase inhibition through treatment with sodium orthovanadate in BAEC cells are simulated by decreasing the dephosphorylation rate constants (k8dp and k11dp) with an intervention factor of IF = 0.3. The model reproduces a slightly more sustained phosphorylated Smad1 response to TGF-b in the presence of sodium orthovanadate as compared to normal signaling dynamics in this cell line (Fig. 4) as it eliminates one mode of negative regulation in the pathway.
Mol. BioSyst., 2014, 10, 537--548 | 541
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Molecular BioSystems
Fig. 3 Simulation of the effects of endocytosis inhibitors on TGF-b signaling dynamics. Using the parameters for HaCaT cells, we simulate the dynamics (lines) of the system and show the results for the ratio of phosphorylated Smad2 to total Smad2 upon stimulation with 0.5 nM (12.8 ng mL1) TGF-b at t = 0 hours. Experimental Western blot data (symbols) and quantification are from ref. 53 with Mv1Lu cells. The experimental data and simulation results are normalized to the maximum from all three conditions. (A) The control case with no inhibitory treatment. (B) Treatment with potassium depletion (–KCl). To simulate the effects of potassium depletion, which inhibits clathrin-mediated endocytosis, we decrease internalization and constitutive degradation rate constants for all receptor species (kdeg,RII, kdeg,RI, k3int, k16deg, k18int, and k19int) by an intervention factor IF = 0.15. (C) Treatment with Nystatin. To simulate the effects of Nystatin, which inhibits lipid raft-caveolar endocytosis, we decrease the rate constants for Smad7–ligand–receptor complex association (k20a,1 and k20a,2) with an intervention factor IF = 0.85. The root-mean-square deviation between the model results and the experimental data is 0.19, 0.06, and 0.06 for panels A, B, and C, respectively.
Fig. 4 Simulation of the effects of phosphatase inhibitors on TGF-b signaling dynamics. Using the parameters for BAEC cells, we simulate the dynamics (lines) of the system and show the results for phosphorylated Smad1 upon stimulation with 1 ng mL1 TGF-b at t = 0 hours. Experimental Western blot data (symbols) are from BAEC cells. (A) The control case with no inhibitory effects. Experimental data are from ref. 15 and 54 and quantified in ref. 15. (B) Treatment with sodium orthovanadate. To simulate the effects of sodium orthovanadate, a phosphatase inhibitor, we decrease the dephosphorylation rate constants (k8dp and k11dp) with an intervention factor IF = 0.3. Experimental data are from ref. 54 and quantified using ImageJ 10.2.52 The root-mean-square deviation between the model results and the experimental data is 0.18 and 0.07 for panels A and B, respectively.
Computational approach for identification of potential therapeutic targets As the TGF-b pathway model is capable of reproducing the effects of both mutations and treatment with inhibitory drugs observed experimentally, we use it as the foundation with which to demonstrate the applicability of the novel in silico approach to identify potential therapeutic targets in intracellular networks. The different steps of the computational approach, summarized in Fig. 5, are discussed in the following sections.
Step 1: identification of targetable network processes. Several strategies have been developed to specifically target biochemical processes that may be used here. Antisense oligonucleotides are used to specifically inhibit protein synthesis by binding to mRNA and preventing translation.47 In the modeled pathway, these could be used to inhibit synthesis of type II receptors, type I receptors, R-Smads, Smad4, and Smad7. Peptide aptamers are engineered to inhibit protein–protein or protein–DNA association by binding to the associated region of
Fig. 5 Flowchart of the approach used to identify potential therapeutic targets. See text (Results section) for additional details on the four steps of the approach.
542 | Mol. BioSyst., 2014, 10, 537--548
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
the protein, competing for the site with the normal binding partner.48 We envision that peptide aptamers could also be used to inhibit nuclear translocation by masking the sites on the Smad proteins responsible for interacting with the nuclear pore complex or accessory proteins mediating nuclear translocation.49 Finally, kinase inhibitors have already been developed to target the ligand–receptor complexes and inhibit phosphorylation of R-Smads.9 From these strategies, we have identified a set of potential targetable processes in the TGF-b pathway and their associated rate constants that would be decreased upon administration of the inhibitory treatment (Table 1). Step 2: sensitivity analysis. We use a global sensitivity analysis (see Methods section) in order to determine how the model responds to perturbation of its parameters. Processes with high sensitivity coefficients will be the most responsive to molecular intervention; therefore, we use sensitivity analysis as a guiding tool for identifying targets in the pathway. We perform the sensitivity analysis with each of the three cell lines mutated with one of the five mutations identified in cancer cells (Table 2), establishing 15 cell-type-mutation combinations, in addition to the 3 non-mutated cases. Each mutation is simulated with a mutation factor of MF = 0.1. In Fig. 6, we show a sample of these results, concretely those for the ‘receptor downregulation’ mutation. The results for the global sensitivity coefficients of the normal cells are included as a control to compare with the effects of the mutation. In general, the mutation decreases the sensitivity coefficients of each parameter as compared to those for the normal cell types, suggesting an increased difficulty in perturbing processes in cells with this mutation. Step 3: sensitive target identification. In order to identify targetable processes that may have the most significant effect on the signaling dynamics, we define as filtering criteria a lower bound of 0.01 for the global sensitivity coefficients Gki. From this, each of the 15 cell-type-mutation combinations will have a set of sensitive parameters with values satisfying the filtering criteria (Gki Z 0.01). To illustrate these results, we color the parameters in Fig. 6 in three groups. Black parameters were not identified as specifically targetable (i.e. not found in Table 1) and are not considered beyond this step. Those in red represent parameters with sensitivity coefficients satisfying the filtering criteria for at least one studied cell type for the considered mutation (H-, B-, or C-; denoting mutated HaCaT, BAEC, or C2C12 cells, respectively). It is important to note that the sensitivity analysis assesses single parameter perturbations, but several of the potential molecular interventions considered in Table 1 affect multiple parameters. Therefore, if at least one parameter from the set of parameters associated with a molecular intervention in Table 1 is selected at this step, we consider the molecular intervention as potentially effective and it passes the filtering procedure to the next step in the approach. Parameters are colored red if within a set of parameters for a particular intervention at least one satisfies the filtering criteria. Those parameters in green correspond to the rest of cases. This is illustrated in the parameter coloring of Fig. S3 (ESI†).
This journal is © The Royal Society of Chemistry 2014
Paper
Step 4: molecular intervention optimization. Using the set of molecular interventions filtered from the sensitivity analysis for each cell-type-mutation combination, we optimize the intervention factor for each case in order to restore normal signaling dynamics to the mutated system, while maintaining the signaling dynamics of the normal system in the presence of the molecular interventions. Fig. 7 displays the results of the optimization routine for all the considered cell types, mutations, and molecular interventions. The ‘No intervention’ case is included as a control in Fig. 7, where we evaluate the objective function (eqn (3)) without any treatment added to the normal or mutated systems. Successful molecular interventions will minimize the objective function value as compared to the ‘No intervention’ case for a given cell type and mutation. Mutations resulting in low values of the objective function signify that the signaling dynamics are robust to the perturbations introduced by that mutation in the specific cell type. In these cases, there is no need of targeting the system with an inhibitor because the mutated cell exhibits similar signaling dynamics as the corresponding normal cell. Therefore, any treatment that reduces the objective function will only result in minimal change from the mutated state as its dynamics already closely mirrors that of the normal state, as illustrated in Fig. 8A for the case of HaCaT cells. By using this optimization routine, rather than relying on the sensitivity analysis alone, we are able to ensure that targeting a particular process is capable of restoring normal signaling dynamics. For example, with the BAEC case affected by the ‘Smad7 downregulation’ mutation, the R-Smad–Smad4 association inhibition was selected as a potential molecular intervention from the sensitivity analysis results (Fig. 7B). However, this mutation results in a more sustained signaling response as compared to that of the normal system and inhibiting R-Smad–Smad4 association increases the signal duration further, moving in the opposite direction than desired (Fig. 8B). Similarly, when the C2C12 case with the ‘inhibited Smad4 nuclear import’ mutation is treated with an R-Smad synthesis inhibitor, the dynamics are unaffected as compared to those with the mutation alone (Fig. 8C). In Fig. 9, we show a selection of the best optimization results for each cell type, which minimize the objective function (eqn (3)) in cases where the mutation has a substantial effect on the normal signaling dynamics. Furthermore, these results minimally impact the signaling dynamics of normal cells in the presence of the molecular intervention, reducing the risk of potential side effects. In HaCaT cells, we show the effects of ‘Smad7 downregulation’ identified in colon cancer (Table 2), which results in a more sustained response to TGF-b for the nuclear phosphorylated Smad2–Smad4 species than without the mutation (Fig. 9A). By inhibiting type II receptor synthesis with an optimized intervention factor of OIF = 0.063, the model dynamics of the mutated case is restored to that of the normal case. For the BAEC case, the ‘inhibited phospho-R-Smad–Smad4 association’ mutation identified in pancreatic cancer (Table 2) results in a more permanent response to ligand stimulation for the nuclear phosphorylated Smad1–Smad4 species (Fig. 9B). By treating the system with Smad4 nuclear export inhibition at the
Mol. BioSyst., 2014, 10, 537--548 | 543
View Article Online
Molecular BioSystems
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Fig. 6 Sensitivity analysis for the ‘receptor downregulation’ mutation. We plot the results obtained for the global sensitivity coefficient Gki (eqn (2)) for each model parameter ki upon stimulation with 45 000 molecules (1.9 ng mL1) of TGF-b. Each column represents the result for one species (pS1S4n or pS2S4n) and one cell type. Results for normal HaCaT, BAEC, and C2C12 cells are denoted as H, B, and C, respectively. The mutated cells, with a mutation factor of MF = 0.1, are abbreviated by appending a dash symbol (–) to the cell type character. Rows represent the results for each parameter, grouped by the module in the pathway it represents. The heat map, generated with Matplotlib 1.1.0,55 represents the value of log10(Gki ) as the color intensities in the color key (righthand side of the heat map); empty spaces indicate a sensitivity coefficient equal to zero or a parameter value set to zero from Table S2 (ESI†). Parameter names in black are not specifically targetable, while those colored red or green represent processes that may be specifically targeted (Table 1). Those in red represent parameters satisfying the filtering criteria described in the Results section, while those in green do not satisfy the filtering criteria.
optimized intervention factor of OIF = 0.037, the dynamics is restored to its normal transient response, while minimally affecting the normal cell case. For the C2C12 case, we mutate the system with the ‘Smad7 downregulation’ mutation, which
544 | Mol. BioSyst., 2014, 10, 537--548
increases the duration of the response to TGF-b stimulation for the nuclear phosphorylated Smad1–Smad4 species (Fig. 9C). By inhibiting R-Smad–Smad4 association with an optimized intervention factor of OIF = 0.132, the mutated-cell dynamics returns
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
Fig. 7 Results of the optimization of the molecular interventions. Bars represent the objective function value (eqn (3)) using one species, mutation, and molecular intervention for (A) HaCaT, (B) BAEC, and (C) C2C12 cells. Bars are grouped by the type of mutation and color coded with the molecular intervention used. Dots beneath an empty bar indicate that the molecular intervention was not selected as potentially effective from the sensitivity analysis step.
to the transient response from the normal case, while minimally impacting the normal C2C12 signaling dynamics. Finally, we filter the results of the molecular interventions optimization routine to produce a list of potential treatments for each cell-type-mutation combination (Table 3). For each cell-type-species combination, we first determine the set of mutations with objective functions for the ‘No intervention’
This journal is © The Royal Society of Chemistry 2014
case greater than or equal to 500 for each cell type. With this, we then determine the set of molecular interventions that minimize the objective function below 0.3 times the ‘No intervention’ value from the set of mutations for each cell-type-species combination. Although the HaCaT dynamics is only substantially responsive to the ‘Smad7 downregulation’ mutation, this method identified two potential molecular interventions to restore normal
Mol. BioSyst., 2014, 10, 537--548 | 545
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
Molecular BioSystems
Fig. 8 Simulation results of mutations and molecular interventions with minimal or negative effects on normal signaling dynamics. Solid blue lines represent simulations of normal cells and solid red lines represent simulations of mutated cells with a mutation factor of MF = 0.1 upon stimulation with 45 000 molecules (1.9 ng mL1) of TGF-b at t = 0 hours. Dashed lines indicate simulations with the molecular interventions. (A) HaCaT cells with the ‘Smad4 downregulation’ mutation. (B) BAEC cells with the ‘Smad7 downregulation’ mutation and treated with an R-Smad–Smad4 association inhibitor using an intervention factor of IF = 0.1. (C) C2C12 cells with the ‘inhibited Smad4 nuclear import’ mutation and treated with an R-Smad synthesis inhibitor using an intervention factor of IF = 0.1.
Fig. 9 Simulation results for optimal molecular interventions. Same conditions as in Fig. 8 for combinations of species, mutation, and molecular interventions selected from the optimization routine for each cell type. (A) HaCaT cells with the ‘Smad7 downregulation’ mutation and treated with a type II receptor synthesis inhibitor using the optimized intervention factor OIF = 0.063. (B) BAEC cells with the ‘inhibited phospho-R-Smad–Smad4 association’ mutation and treated with a Smad4 nuclear export inhibitor using the optimized intervention factor OIF = 0.037. (C) C2C12 cells with the ‘Smad7 downregulation’ mutation and treated with an R-Smad–Smad4 association inhibitor using the optimized intervention factor OIF = 0.132.
Table 3 Potential targetable processes identified by optimization of the molecular interventions. For each cell type and mutations that significantly affect the normal signaling dynamics, we show the corresponding molecular interventions that restore normal signaling dynamics to the mutated system. The entries in boldface are those used in Fig. 9
HaCaT Mutation
Targeted process
Smad4 downregulation
Smad7 downregulation
BAEC Species Targeted process Type II receptor synthesis Type I receptor synthesis Smad4 nuclear export
Type II receptor synthesis pS2S4n Type I receptor synthesis R-Smad–Smad4 pS2S4n Smad4 nuclear export association
C2C12 Species Targeted process pS1S4n pS1S4n pS1S4n pS1S4n Type I receptor synthesis
pS1S4n
pS1S4n Smad4 nuclear import
pS1S4n
R-Smad–Smad4 association Inhibited Smad4 nuclear import
Inhibited R-Smad–Smad4 association
546 | Mol. BioSyst., 2014, 10, 537--548
Type I receptor synthesis R-Smad nuclear export Smad4 nuclear export
Species
pS1S4n
pS1S4n pS1S4n pS1S4n
Type I receptor synthesis pS1S4n R-Smad nuclear export pS1S4n Smad4 nuclear export pS1S4n Receptor–R-Smad association
pS2S4n
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Molecular BioSystems
signaling dynamics to the mutated system: inhibited type II receptor synthesis and inhibited R-Smad–Smad4 association. For the BAEC case, this method identified a number of molecular interventions for four of the five mutations considered here, targeting receptor synthesis and R-Smad/Smad4 nucleocytoplasmic shuttling. The results for the C2C12 case reveal four distinct molecular interventions that may be used to target the system for two of the mutations we used in this analysis. Finally, our approach did not identify any molecular intervention for the ‘TGF-b receptor downregulation’ mutation in any of the considered cases.
Conclusions Effectively targeting the TGF-b signal transduction pathway requires a detailed dissection of the network interactions beyond insight gained from intuition alone. Here we have developed a strategy based on computational modeling to quantitatively analyze the signaling pathway in order to determine potential therapeutic targets. We have considered a model that captures the precise dynamics of the system, mutations that affect the system parameters, and a collection of potentially targetable components of the pathway, such as inhibition of association or synthesis of proteins. The key element of the approach is the identification of molecular interventions in the pathway that, in addition to bringing the mutated system back to its normal dynamics, do not substantially affect the non-mutated system. To do so, our approach uses global sensitivity analysis to identify parameters with high sensitivity coefficients and passes the molecular interventions containing these sensitive parameters to the optimization routine that selects the optimal strength of the molecular intervention. We have shown the capabilities of the computational approach for different cell types by reproducing experimental TGF-b signaling dynamics when the system is mutated in pancreatic cancer cells and in the presence of inhibitory drugs targeting distinct processes in the pathway. From the results of the molecular intervention optimization routine, we have identified a set of mutations that significantly affect the signaling dynamics for each cell type and a number of molecular interventions that may be used to effectively target the effects of these mutations. Interestingly, our approach identified multiple types of molecular interventions that can be used for a given mutation. This type of approach can generally be applied to a wide variety of signaling systems with characterized biochemical interactions and mutations affecting their components to identify novel targets for therapeutic intervention. Therefore, our methodology provides an efficient starting point to select potential targets for further experimental and clinical investigation.
Acknowledgements This work was supported by the University of California, Davis (to LS).
This journal is © The Royal Society of Chemistry 2014
Paper
References ´, Annu. Rev. Biochem., 1998, 67, 753–791. 1 J. Massague 2 L. Levy and C. S. Hill, Cytokine Growth Factor Rev., 2006, 17, 41–58. 3 P. ten Dijke and H. M. Arthur, Nat. Rev. Mol. Cell Biol., 2007, 8, 857–869. 4 G. C. Blobe, W. P. Schiemann and H. F. Lodish, N. Engl. J. Med., 2000, 342, 1350–1358. 5 R. Derynck and K. Miyazono, The TGF-b Family, Cold Spring Harbor Laboratory Press, Cold Spring Harbor, N.Y., 2008. 6 Y. Drabsch and P. ten Dijke, Cancer Metastasis Rev., 2012, 31, 553–568. ´, Cell, 2003, 113, 685–700. 7 Y. Shi and J. Massague ´, Cell, 2008, 134, 215–230. 8 J. Massague 9 J. M. Yingling, K. L. Blanchard and J. S. Sawyer, Nat. Rev. Drug Discovery, 2004, 3, 1011–1022. 10 J. M. G. Vilar, R. Jansen and C. Sander, PLoS Comput. Biol., 2006, 2, e3. 11 J. M. G. Vilar and L. Saiz, Biophys. J., 2011, 101, 2315–2323. 12 V. Becker, J. Timmer and U. Klingmuller, Adv. Exp. Med. Biol., 2012, 736, 313–323. 13 D. C. Clarke, M. D. Betterton and X. Liu, IEE Proc.: Syst. Biol., 2006, 153, 412–424. 14 B. Schmierer, A. L. Tournier, P. A. Bates and C. S. Hill, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 6608–6613. ¨nsson, E. Pardali, P. ten Dijke and 15 P. Melke, H. Jo C. Peterson, Biophys. J., 2006, 91, 4368–4380. 16 M. Paulsen, S. Legewie, R. Eils, E. Karaulanov and C. Niehrs, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10202–10207. 17 K. Wegner, A. Bachmann, J. U. Schad, P. Lucarelli, S. Sahle, ¨ller, S. Dooley and P. Nickel, C. Meyer, U. Klingmu U. Kummer, Biophys. Chem., 2012, 162, 22–34. 18 G. Celliere, G. Fengos, M. Herve and D. Iber, BMC Syst. Biol., 2011, 5, 184. 19 D. Nicklas and L. Saiz, J. R. Soc. Interface, 2013, 10, 20130363. 20 J. Ho and L. Saiz, Biophys. J., 2011, 100, 164a. 21 Z. Zi and E. Klipp, PLoS One, 2007, 2, e936. 22 S. W. Chung, F. L. Miles, R. A. Sikes, C. R. Cooper, M. C. Farach-Carson and B. A. Ogunnaike, Biophys. J., 2009, 96, 1733–1750. 23 Z. Zi, Z. Feng, D. A. Chapnick, M. Dahl, D. Deng, E. Klipp, A. Moustakas and X. Liu, Mol. Syst. Biol., 2011, 7, 492. 24 B. S. Hendriks, G. J. Griffiths, R. Benson, D. Kenyon, M. Lazzara, J. Swinton, S. Beck, M. Hickinson, J. M. Beusmans, D. Lauffenburger and D. de Graaf, IEE Proc.: Syst. Biol., 2006, 153, 457–466. 25 R. J. Orton, M. E. Adriaens, A. Gormand, O. E. Sturm, W. Kolch and D. R. Gilbert, BMC Syst. Biol., 2009, 3, 100. 26 S. W. Chung, C. R. Cooper, M. C. Farach-Carson and B. A. Ogunnaike, J. R. Soc. Interface, 2012, 9, 1389–1397. 27 N. Kumar, B. S. Hendriks, K. A. Janes, D. de Graaf and D. A. Lauffenburger, Drug Discovery Today, 2006, 11, 806–811. 28 L. H. Abbott and F. Michor, Br. J. Cancer, 2006, 95, 1136–1141.
Mol. BioSyst., 2014, 10, 537--548 | 547
View Article Online
Published on 10 December 2013. Downloaded by University of Connecticut on 29/10/2014 23:56:06.
Paper
29 K. Oda, Y. Matsuoka, A. Funahashi and H. Kitano, Mol. Syst. Biol., 2005, 1, 20050010. 30 G. Lebedeva, A. Sorokin, D. Faratian, P. Mullen, A. Goltsov, S. P. Langdon, D. J. Harrison and I. Goryanin, Eur. J. Pharm. Sci., 2012, 46, 244–258. 31 R. P. Araujo, E. F. Petricoin and L. A. Liotta, BioSystems, 2005, 80, 57–69. 32 M. H. Sung and R. Simon, Mol. Pharmacol., 2004, 66, 70–75. 33 H. Yan, B. Zhang, S. Li and Q. Zhao, BMC Syst. Biol., 2010, 4, 50. 34 B. J. Lao, W. L. Tsai, F. Mashayekhi, E. A. Pham, A. B. Mason and D. T. Kamei, J. Controlled Release, 2007, 117, 403–412. 35 D. J. Yoon, D. S. Chu, C. W. Ng, E. A. Pham, A. B. Mason, D. M. Hudson, V. C. Smith, R. T. MacGillivray and D. T. Kamei, J. Controlled Release, 2009, 133, 178–184. 36 D. J. Yoon, B. H. Kwan, F. C. Chao, T. P. Nicolaides, J. J. Phillips, G. Y. Lam, A. B. Mason, W. A. Weiss and D. T. Kamei, Cancer Res., 2010, 70, 4520–4527. 37 D. Hanahan and R. A. Weinberg, Cell, 2000, 100, 57–70. 38 D. Hanahan and R. A. Weinberg, Cell, 2011, 144, 646–674. 39 A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, ACM Trans. Math. Softw., 2005, 31, 363–396. 40 Z. Zi, IET Syst. Biol., 2011, 5, 336–346. 41 S. Kucherenko, M. Rodriguez-Fernandez, C. Pantelides and N. Shah, Reliab. Eng. Syst. Saf., 2009, 94, 1135–1148. 42 A. Varma, M. Morbidelli and H. Wu, Parametric sensitivity in chemical systems, Cambridge University Press, Cambridge, UK; New York, NY, 1999. 43 W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical recipes: the art of scientific computing, Cambridge University Press, New York, 2007. 44 E. Jones, T. Oliphant, P. Peterson and others, SciPy: Open source scientific tools for Python (http://www.scipy.org), 2001–. 45 F. James and M. Winkler, MINUIT user’s guide, CERN, Geneva, 2004. ´s and C. S. Hill, Oncogene, 2003, 22, 3698–3711. 46 F. J. Nicola
548 | Mol. BioSyst., 2014, 10, 537--548
Molecular BioSystems
47 C. F. Bennett and E. E. Swayze, Annu. Rev. Pharmacol. Toxicol., 2010, 50, 259–293. 48 A. Conidi, V. van den Berghe and D. Huylebroeck, Int. J. Mol. Sci., 2013, 14, 6690–6719. 49 C. S. Hill, Cell Res., 2009, 19, 36–46. ´n, S. Itoh, A. Moustakas, P. t. Dijke and 50 A. More C.-H. Heldin, Oncogene, 2000, 19, 4396–4404. 51 A. Hata, R. S. Lo, D. Wotton, G. Lagna and J. Massague, Nature, 1997, 388, 82–87. 52 M. D. Abramoff, P. J. Magalhaes and S. J. Ram, Biophoton. Int., 2004, 11, 36–42. 53 G. M. Di Guglielmo, C. Le Roy, A. F. Goodfellow and J. L. Wrana, Nat. Cell Biol., 2003, 5, 410–421. 54 G. Valdimarsdottir, M.-J. Goumans, F. Itoh, S. Itoh, C.-H. Heldin and P. ten Dijke, BMC Cell Biol., 2006, 7, 16. 55 J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95. ´, Nat. Cell Biol., 1999, 1, 472–478. 56 R. S. Lo and J. Massague 57 M. Wan, X. Cao, Y. Wu, S. Bai, L. Wu, X. Shi, N. Wang and X. Cao, EMBO Rep., 2002, 3, 171–176. 58 W. Liu, H. Rui, J. Wang, S. Lin, Y. He, M. Chen, Q. Li, Z. Ye, S. Zhang, S. C. Chan, Y.-G. Chen, J. Han and S.-C. Lin, EMBO J., 2006, 25, 1646–1658. 59 P. R. Segarini, D. M. Rosen and S. M. Seyedin, Mol. Endocrinol., 1989, 3, 261–272. 60 J.-F. Goetschy, O. Letourneur, N. Cerletti and M. A. Horisberger, Eur. J. Biochem., 1996, 241, 355–362. 61 S.-B. Peng, L. Yan, X. Xia, S. A. Watkins, H. B. Brooks, D. Beight, D. K. Herron, M. L. Jones, J. W. Lampe, W. T. McMillen, N. Mort, J. S. Sawyer and J. M. Yingling, Biochemistry, 2005, 44, 2293–2304. 62 M. Funaba and L. S. Mathews, Mol. Endocrinol., 2000, 14, 1583–1591. 63 X. Lin, X. Duan, Y.-Y. Liang, Y. Su, K. H. Wrighton, J. Long, M. Hu, C. M. Davis, J. Wang, F. C. Brunicardi, Y. Shi, Y.-G. Chen, A. Meng and X.-H. Feng, Cell, 2006, 125, 915–928. 64 B. Schmierer and C. S. Hill, Mol. Cell. Biol., 2005, 25, 9845–9858.
This journal is © The Royal Society of Chemistry 2014
