|Home | About | Journals | Submit | Contact Us | Français|
Through a multidisciplinary approach involving experimental and computational studies, we address quantitative aspects of signaling mechanisms triggered in the cell by the receptor targets of hallucinogenic drugs, the serotonin 5-HT2A receptors. To reveal the properties of the signaling pathways, and the way in which responses elicited through these receptors alone and in combination with other serotonin receptors subtypes (the 5-HT1AR), we developed a detailed mathematical model of receptor-mediated ERK1/2 activation in cells expressing the 5-HT1A and 5-HT2A subtypes individually, and together. In parallel, we measured experimentally the activation of ERK1/2 by the action of selective agonists on these receptors expressed in HEK293 cells. We show here that the 5-HT1AR agonist Xaliproden HCl elicited transient activation of ERK1/2 by phosphorylation, whereas 5-HT2AR activation by TCB-2 led to higher, and more sustained responses. The 5-HT2AR response dominated the MAPK signaling pathway when co-expressed with 5-HT1AR, and diminution of the response by the 5-HT2AR antagonist Ketanserin could not be rescued by the 5-HT1AR agonist. Computational simulations produced qualitative results in good agreement with these experimental data, and parameter optimization made this agreement quantitative. In silico simulation experiments suggest that the deletion of the positive regulators PKC in the 5-HT2AR pathway, or PLA2 in the combined 5-HT1A/2AR model greatly decreased the basal level of active ERK1/2. Deletion of negative regulators of MKP and PP2A in 5-HT1AR and 5-HT2AR models was found to have even stronger effects. Under various parameter sets, simulation results implied that the extent of constitutive activity in a particular tissue and the specific drug efficacy properties may determine the distinct dynamics of the 5-HT receptor-mediated ERK1/2 activation pathways. Thus, the mathematical models are useful exploratory tools in the ongoing efforts to establish a mechanistic understanding and an experimentally testable representation of hallucinogen-specific signaling in the cellular machinery, and can be refined with quantitative, function-related information.
Through a collaborative multidisciplinary approach involving experimental and computational studies at various scales – from molecular to cellular and organismal – we are engaged in an effort to reveal the mechanisms that engender the complex effects of hallucinogenic drugs of abuse (for some reviews see (Aghajanian and Marek 1999; Gresch, Strickland et al. 2002; Nichols 2004)). Our computational work addresses quantitative and structural aspects of the mechanisms of hallucinogenic drugs in various chemical classes. Using methods of molecular biophysics and computational biology, we simulate the dynamic properties of the ligand-bound receptor systems for hallucinogens compared to non-hallucinogenic congeners (Weinstein 2006). Proceeding further up in the size and time scale of the relevant processes, we show here how mathematical models of receptor-mediated signaling properties can be used to connect to experimentally determined signaling. We focus on molecular complexes and interactions of these compounds with serotonin receptors in the G protein coupled receptors (GPCRs) family, and follow the mechanisms through ensuing interaction processes with other components of the signaling cascades such as membrane components (e.g., PIP2) and PDZ domains (Madsen, Beuming et al. 2005; Khelashvili, Weinstein et al. 2008).
The rigorous quantitative approach we apply is made possible by recent advances in many aspects of experimental and computational biology, from the molecular to the integrative level of cell signaling systems. These advances have improved our understanding of GPCR activation as an allosteric mechanism, triggered by ligand binding, that results in conformational changes transduced by the transmembrane domain (Han, Wang et al. 2008) (Urban, Clarke et al. 2007) (Weinstein 2006; Deupi and Kobilka 2007) (Kobilka and Deupi 2007). The advances have also allowed us to characterize the steps of intra-receptor activation mechanism by combing computational modeling and experimentation (Guo, Shi et al. 2005). Experimental evidence points to multiple conformations related to the activation of the receptor (Niv, Skrabanek et al. 2006), (Filizola, Wang et al. 2006; Han, Wang et al. 2008). Different ligands binding to the same receptor may induce different conformational states, which in turn can result in coupling to different signaling pathways (specifically for the hallucinogens, see Gonzales-Maeso, Weisstaub et al. 2007), and functional hetero-oligomarization (Gonzalez-Maeso, Ang et al. 2008). We have recently reviewed (Weinstein 2006) some key aspects of functional understanding achievable from computational modeling of hallucinogen mechanisms at the molecular and cellular level, emphasizing not only the structural context of the mechanisms of the receptor molecules and their interactions, but also the importance of bioinformatics and mathematical modeling tools in revealing the specific consequences of hallucinogen binding to GPCRs. The findings leading to this newly gained understanding include key mechanistic components such as (i)-modes of receptor response (conformational rearrangements and stabilization of “activated state(s)”) responsible for protein-protein interactions ranging from homo- and hetero-oligomerization to interactions with scaffold proteins (e.g., PDZ domains), and (ii)-the role of conformational rearrangements of the receptor due to hallucinogen binding in association/dissociation of specific protein-protein interactions and selective signaling. These developments show why models of the activated forms of GPCRs have become increasingly necessary for the development of a clear understanding of signal propagation into the cell (Niv, Skrabanek et al. 2006) (Filizola, Wang et al. 2006; Han, Wang et al. 2008).
Here we summarize briefly the recent progress along these lines, by presenting a topological network and a mathematical model that offer a detailed visual, quantitative and dynamic illustration of the 5-HT receptors-mediated ERK pathways, known to be targeted by hallucinogens in their actions (5-HT2A, and relations to 5-HT1A). The current understanding of detailed signaling mechanisms is still incomplete, and the determinants for the function of these GPCRs in cellular signaling triggered by hallucinogens are only partially delineated. Therefore, in order to gain quantitative presentations of this complicated network, we derive mathematical models by focusing on particular signaling processes activated by hallucinogens through 5-HT receptors. Quantitative understanding of the actions of hallucinogens must include the signaling pathways activated by these ligands after binding to the target GPCRs (for a discussion see (Niv, Skrabanek et al. 2006), (Weinstein 2006), and (Kholodenko 2006; Palsson 2006)). Here we illustrate the use of computational modeling in the quantitative interpretation of currently available data of serotonin receptors-mediated MAPK cascade, and collect all the pieces into function-related, time-dependent information. While network representations shown in Figures 1 and and22 may be incomplete, and the values of the parameters may carry significant uncertainties, these are likely to be remedied by results from continuing research. However, the integrative representation and quantitative summary of current literature offered by these models are the focus of interest. At the very least, the simulations of these models can aid ongoing efforts to construct an increasingly comprehensive mechanistic understanding by validating or eliminating specific assumptions, answering particular mechanistic questions (see below), and guide the necessary experimental effort by producing experimentally testable hypotheses based on the representation of hallucinogen-specific signaling in the cellular machinery.
We have focused on the activation of the eukaryotic MAPK cascade via serotonin receptors 5-HT1A and 5-HT2A as it is ubiquitously expressed in diverse biological processes. MAPK signal transduction pathways mediate short-term effects such as modulation of potassium channel (Yuan, Adams et al. 2002) and glutamate receptor function (Endo and Launey 2003) as well as long-term effects such as cell differentiation, long-term potentiation (LTP), learning and memory (Adams and Sweatt 2002). Signaling through MAPK pathways is known to positively regulate immediate early genes (IEGs). In addition, MAPK cascades in a variety of cells are tightly regulated by multiple feedback loops. Although the basic structure of all MAPK cascades is the same, differences in feedback control enable them to generate a plethora of biological responses, including oscillations, gradual and ultrasensitive responses(Huang and Ferrell 1996; Chang and Karin 2001) (Huang and Ferrell 1996) (Bhalla, Ram et al. 2002) (Heinrich, Neel et al. 2002) (Kholodenko, Hoek et al. 1997) (Kholodenko 2000).
The action of hallucinogens on 5-HT receptors is well documented (Nichols 2004). While drug discrimination experiments have singled out the 5-HT2AR subtype as the important target of hallucinogens, both the 5-HT1A and 5-HT2A receptors have been suggested to be involved. The serotonin receptor 5-HT1A is coupled to Gi/Go proteins, and stimulates the MAPK growth-signaling pathway possibly through G protein βγ complex-mediated phoshatidylinositol 3′-kinase (PI-3K) or phospholipase (PLC) β pathways (Della Rocca, Mukhin et al. 1999) (Adayev, El-Sherif et al. 1999) (Della Rocca, van Biesen et al. 1997).
The 5-HT2A serotonin receptor is Gq/11-coupled and has diverse roles in both the central nervous system (CNS) and peripheral vasculature, and is known to trigger MAPK activation via PKC/Raf-1 pathway (Hershenson, Chao et al. 1995) (Watts 1996), and also to stimulation of PLA2 to produce the second messenger arachidonic acid (AA) (Felder, Kanterman et al. 1990) (Tournois, Mutel et al. 1998). While the 5-HT2A receptors have been implicated as direct targets of hallucinogens, the balance of signaling activities that produce the hallucinogenic effect remains unknown. In particular, an essential issue in signal transduction is how the activated receptors are integrated into signaling pathways and how specific conformations of the activated receptor may establish the distinct patterns of signal transduction observed when they bind different ligands; notably, hallucinogens produced entirely different transcriptome fingureprints compared to their non-hallucinogenic congeners (Gonzalez-Maeso, Yuen et al. 2003). Some key questions then become: Which of the reactions in these complex networks are most important? Where are the cross-talk points that are regulated by upstream or downstream components? To address such questions it becomes essential to acquire qualitative information on which interactions take place, and quantitative data on their strength. Mathematical modeling can be then be used to integrate such information, to identify the key signaling molecules and predict the system behavior of the specified pathway.
To initiate such a study, we created the detailed topological representation of the 5-HT1A and 5-HT2A receptors-medicated ERK activation described here, based on known reactions and assumptions derived from canonical pathways. This representation of interactions was then translated into mathematical reactions that describe the network topology. We then used computational simulations to solve the mathematical equations, which then yielded predictions of species concentration profiles that vary with respect to time upon ligand stimulation. In order to validate the dynamics predicted from the simulations we carried out in parallel experimental analyses in Human Embryonic Kidney (HEK) 293 cells stably expressing serotonin receptors. Because the presence of Gi- and Gq-coupled receptors responsive to a common hormone/neurotransmitter may synergize, we hypothesized that co-expression of 5-HT1A and 5-HT2A receptors may result in enhanced activity of MAPK ERK in HEK 293 cell line. Here we report that HEK 293 cells recombinantly expressing 5-HT1A and 5-HT2A receptors produce dynamics of ERK activation distinct from receptors that are expressed alone. We also report that when the two receptors are expressed in similar ratios, the 5-HT2A receptors seem to dominate the ERK signal in HEK 293 cells, both in duration and magnitude. We show that treatment with the 5-HT2A receptor antagonist Ketanserin produces a switch in the ERK activation pattern from sustained to transient, suggesting that in vivo, the levels of the 5- HT2A receptor expression may play an important role in the ERK activity phenotype. Comparing experimental results with simulation data we demonstrate here that the individual activation pathway models produce results in qualitatively good agreement with the experimental data, and that following parameter optimization, the computational analysis agrees quantitatively with the experimental results.
The important parameters and intrinsic behaviors of the system were further revealed by sensitivity analysis. Moreover, we report simulation results from in silico experiments, including parameter perturbation and deletion of key regulators along MAPK feedback pathways. Knockout of the positive regulators PKC and PLA2 in 5-HT2AR and combined 5-HT1A/5-HT2AR models, respectively resulted in a great reduction of basal levels of active ERK1/2 (Supplementary figure 4). Compared with PKC and PLA2, the negative regulators PP2A and MKP are found to produce even stronger effects on cells in all three models.
In summary, the models suggest that constitutive activity in a particular tissue, combined with specific drug efficacy may determine distinct dynamics of a 5-HT receptors-mediated ERK1/2 pathway, and therefore affecting the receptor activation phenotypes. The in silico experiments provide insights into the underlying mechanisms of ERK pathways via 5-HT receptors, in which can be further validated by inhibitors or activators, siRNA or transfections to influence the activity and expression of target genes.
The full length cDNA clones for the human serotonin receptor gene subtypes were either purchased from the UMR cDNA Resource (5HT1AR) or were received as a generous gift (5H2AR) from Dr. Stuart Sealfon. Each cDNA was subcloned into a pcDNA3.1(+) plasmid vector (Invitrogen) containing drug resistance genes for either G418 (neomycin phosphotransferase) or hygromycin B (Hygromycin B phosphotransferase). The drugs were either purchased from Sigma/RBI (Natick, MA) or Tocris (Ellisville, MO). Analytical grade chemicals were purchased from Sigma-Aldrich (St. Louis, MO), and cell culture supplies were purchased from either Sigma-Aldrich or Hyclone Laboratories (Logan, UT). The radiolabeled antagonists [3H]methylspiperone (MSP) (NET-856, 80 Ci/mmol) and 4-(2′-Methoxy)-phenyl-1-[2′-(N-2″-pyridinyl)-p-fluorobenzamido]ethyl-piperazine (Methyl-[3H]-MPPF) (NET-1109, 79.8 Ci/mmol) were purchased from PerkinElmer Life and Analytical Sciences (Boston, MA). The antibiotics G418 and HygroGold were purchased from InvivoGen (San Diego, CA).
Plasmid constructs containing the wild type human serotonin receptors were transfected into HEK293 cells using CaPO4 precipitation (Invitrogen, Carlsbad, CA). Specifically, 20 μg of plasmid DNA was mixed with a final volume of 1 ml of CaPO4/HEPES solution, and the resulting precipitate was added drop-wise to 20 to 30% confluent cells attached to a 150 cm2 plate in a total volume of 20 ml Dulbecco’s modified Eagle’s media supplement with 10% bovine calf serum (BCS) and antibiotics. The following day, the media was removed by aspiration and replaced with fresh drug selection media: either G418 (2 mg/mL) or Hygromycin B (100 μg/mL) About two weeks later, single colonies were selected, expanded and checked for the presence of the serotonin receptors via the whole cell binding assay. Those clonal cell lines expressing high levels of receptor were passaged many times under conditions of constant selective pressure (4–6 weeks) until stable clones lines were achieved. Saturation isotherm binding assay was then used to examine receptor expression levels in stable clonal cell lines. For cells co-expressing two different serotonin receptors, 5-HT1AR gene was subcloned into pcDNA3.1/Hygro (+) plasmid and transfected into the 5-HT2AR (G418 resistant) stable cell line with similar conditions as described before.
HEK 293 cells expressing the desired receptor were dislodged by a 5 min incubation in Earle’s balanced saline solution lacking Ca2+ and Mg2+ and supplemented with 5 mM EDTA. After centrifugation, the cell pellet was lysed in lysis buffer (5 mM Tris and 5 mM MgCl2, pH 7.4) at 4 °C. The lysate was glass-glass homogenized (eight stroke), and the membranes were centrifuged at 35,000g for 30 min. The pellet was resuspended in ice-cold 50 mM Tris, pH 7.4, and centrifuged again. The washed membrane pellet was resuspended by light homogenization (three strokes) in binding buffer (see below) immediately before use.
Membranes containing wild type serotonin receptors were assayed for specific binding activity. [3H] MSP (0.3–10 nM) was used to estimate numbers of 5-HT2A receptors, while [3H] MPPF was used for measuring the expression of 5-HT1A receptors. The binding buffer consisted of 50 mM Tris, pH 7.4, at 25 °C. Non-specific bindings for 5-HT1AR and 5-HT2AR were defined in the presence of excess NAN-190 and mianserin, respectively. The incubation was carried out at 25°C for 1.5 h and terminated by rapid filtration through GF/C filters pretreated with 0.3% polyethylenimine. The filters were rinsed once using chilled washing buffer, consisting of ice-cold binding buffer (pH 7.4, 0 °C). Radioactivity bound to the filters was quantified using scintillation spectroscopy at a counting efficiency of 47%. Membrane protein concentrations were determined using the bicinchoninic acid protein reagent (Piere, IL) and a bovine serum albumin standard curve.
All points for each experiment were sampled in triplicates (unless specified otherwise). The geometric mean values of the data from three independent experiments are reported with their associated standard deviation. The equilibrium dissociation constant (KD) of the primary radioligand was measured by saturation isotherm analysis.
Prism 4 (GraphPad Software, Inc., San Diego, CA) was used to plot data and determine drug binding affinity values by analyzing the saturation isotherm curves. A 95% confidence interval was used for all curve-fitting procedures and to compare different curve-fitting models. The statistical measures of fit employed were the F-test, the run test, and a correlation coefficient.
5-HT and selective agonists, such as Xaliproden hydrochloride (HCl) and (4-Bromo-3,6- dimethoxybenzocyclobuten-1-yl)methylamine hydrobromide (TCB-2) for 5-HT1AR and 5-HT2AR, respectively, were used to study the activation of ERK. Stable cell lines were maintained in selective medium containing low concentrations (100 μg/ml) of G418 and/or HygroGold. Prior to the drug treatment, cells were starved by serum-free DMEM for at least 24 h. Various concentrations of stimuli (10−10–10−4 M) were tested on cells for 5 min; and optimal concentrations of agonists determined from dose-response curves were tested from 5 min up to 1 hour. Agonist stimulations were stopped by removing the media and rinsing once with cold PBS, solubilized with pre-chilled lysis buffer containing protease inhibitors. After brief sonication and centrifugation at 13,000g at 0°C for 45 min, cellular proteins were extracted and protein concentrations were measured with Bio-Rad protein assay. An equal amount of protein was analyzed on 10% SDS-PAGE and electroblotted onto nitrocellulose membranes (Whatman Inc., Dassel, Germany). Membranes were blocked in a solution containing TBS-Tween 20 (50 mM Tris Base, 0.9% NaCl, 0.1% Tween 20, pH 7.6) and nonfat dry milk (5% wt/vol) and then incubated in a solution containing TBS-Tween 20, 5% wt/vol nonfat dry milk and primary antibody (1:1,000) overnight at 4 °C. Membranes were then washed with TBS-Tween 20 for three times for 5 min each time, and incubated with secondary antibody at room temperature for 1 h. Three washes with TBS-Tween 20 were repeated as before. Note that experiments were carried out with two blots running parallely in each same amount of identical samples were loaded, each blot were then incubated with different primary antibody (phospho and non-phospho).
A developing solution implementing Enzymatic Chemiluminescence (ECL) was used for visualization of proteins. Phosphorylated forms of ERK1/2 were detected by immunoblotting using a rabbit phospho- p42/p44 (ERK1/2) polyclonal primary antibody (Cell Signaling, Inc., Boston, MA), and a secondary anti-rabbit IgG antibody was used for amplifying the signal; while non-phosphorylated ERK1/2 as basal level control was validated with a rabbit p42/p44 ERK polyclonal primary antibody after the immunoblots. The data from both the phospho and non-phospho- ERK immunoblots were quantified by densitometric analysis using KODAK 1D image software. First, all collective data from the same experiment/blot were subtracted from its background nose. Second, to quantify the fold of activation, all experimental data were divided by their respective control density. The data point where no treatment was given was set to one, therefore representing the basal level. Each experiment was repeated at least three times.
Based on the simplified cartoon diagram (Figure 1), the kinetic scheme (Figure 2) of MAP kinase signaling cascades for 5-HT1A and 5-HT2A receptors in HEK 293 cells was constructed using Visio (Microsoft, Inc.). This scheme describes the connectivity of the reaction network (Figure 2).
The temporal dynamics of the proposed reaction scheme (Figure 2) are described by a deterministic Ordinary Differential Equation (ODE) model. In the ODE model, the rate of change of a species concentration with time (time derivative) is equal to the sum of all reactions that produce this species minus all reactions that consume this species. The model inputs are kinetic parameters and initial concentrations. Enzyme kinetics are described by the standard Michaelis-Menten formulation. The cytosolic space is assumed to be well-mixed and diffusion reactions are ignored. The rate constants are either from the literature or estimated; the detailed rate equations and parameter values used in this study are given in the Supplementary Tables 1, 2, 3.
Three types of biochemical reactions were considered: translocation, transformation and binding. Translocation is the movement of a chemical species from one location to another. For instance, PKC can be recruited from cytosol to membrane by DAG activation; IP3, a second messenger generated from PIP2 hydrolysis, is released from plasma membrane to cytosolic space. Transformation is the conversion of one molecule to another. Examples include conversion of PIP2 to IP3 and DAG by PLCβ, phosphorylation of Raf by active form of PKC. Binding consists of two molecules combing, usually non-covalently, to form a single complex. Examples are agonist binding to the 5-HT2A receptor, and RGS binding to Gαq-GTP. Binding can be regarded as a special case of transformation since second order or higher order reactions can be treated as one-step processes. In our model, the transformation reactions were represented by the Michaelis-Menten formulation or two-step binding reactions, with k−2 = 0. Individual steps were translated into differential equations based on the proposed kinetic scheme.
We adopted the two-state mechanism for 5-HT receptors activation, i.e., prior to agonist stimulation receptors exist in equilibrium between an inactive (R) and an active (R*) state with a stability (equilibrium) constant J (Leff 1995). In the R state, receptors are uncoupled from G-proteins, whereas in the R* state, receptors can couple to and activate G-proteins. The conformational change in GPCRs associated with R- to R* isomerization enables GPCRs to promote the dissociation of GDP from G-proteins, which is the rate-limiting step in the G-protein (Gilman 1987). The binding affinities for L+R (lowaffinity site: , where M is the equilibrium binding constant) and L+RG (high affinity site: , where β is defined as ligand’s efficacy) were taken from the literature (Pou, Nenonene et al. 1997; Roth, Choudhary et al. 1997), and the β constant was set as β > 1 for an agonist which promotes the functional R* form (the higher the efficacy, the larger the value of β), and β < 1 for an inverse agonist; β= 1 characterizes a neutral agonist. Assuming that J, the ratio between R* and R =~0.1, i.e., that 10% of the GPCRs were constitutively active.
Sensitivity analysis was carried out to capture the essential characteristics of the model and to determine the critical parameters that control the peak phosphorylation of ERK1/2. These critical parameters include either the initial conditions of the species states or the kinetic rate constants. The time-dependent sensitivities of the phosphorylated ERK1/2 to the changes in these parameters were plotted and the sensitivity values at the time point of phosphorylated ERK1/2 peaks were identified. These values were fully normalized to enable their comparison with each other and are listed in the Supplementary Tables 4, 5. The details of sensitivity analysis calculations can be found in the Supplementary Information.
The mathematical model consists of a set of coupled ODEs and is parametric in reaction rate constants and initial concentrations. After the critical parameters that affect the phosphorylated peak ERK1/2 concentration were identified by sensitivity analysis (Supplementary Table 6), we investigated whether modification of these parameters led to a better fit of the model to the experimental phosphorylated ERK1/2 data using dynamic local optimization. This was done by minimizing the sum of squared errors between the experimentally observed and mathematically simulated phosphorylated ERK1/2 concentration at 4 different time points by manipulating these parameter values. The experimental error was assumed to be normally distributed with zero mean and the covariance matrix was assumed to be same in each time point and diagonal.
Note that local optimization was performed to identify a set of parameter combinations that lead to the observed phosphorylated ERK1/2 behavior. We did not employ a global optimization algorithm to identify the parameters that best fit the experimental data as these algorithms are computationally expensive and, furthermore, the best parameter combination is not unique; different sets of parameter combinations can lead to the same observed output, an inherent property of multiparameter nonlinear systems biology models (Gutenkunst, Waterfall et al. 2007).
Simulations were carried out using SimBiology Toolbox within MATLAB software (The Mathworks; Natick, MA). Differential equations were integrated using the function ode15s, which is a variable order solver, based on the numerical differentiation formulas (NDFs) and is designed for stiff systems. Local optimization was performed using the fmincon and lsqnonlin functions in SimBiology Toolbox of MATLAB. Both functions employ a subspace trust region method and are based on an interiorreflective Newtonian method. Each iteration involves the approximate solution of a large linear system using the method of preconditioned conjugate gradients.
Saturation isotherm binding assay showed that HEK 293 cells stably expressed h5-HT1AR, h5-HT2AR, and h5-HT1A/2AR (Figure 3.) The ratio of 5-HT1AR and 5-HT2AR is close to 1:1 in co-expressing cell line (Table 1).
Dose-response experiments were carried out in 5-HT1AR and 5-HT2AR cells, using Xaliproden HCl and TCB-2, respectively. Minimal concentrations of Xaliproden HCl 100 nM and TCB-2 1 nM were required to obtain observed ERK activation in individual cell lines, respectively (Supplementary figure 1C and 1D). We found that 1 μM of Xaliproden HCl and 10 nM of TCB-2 served best in the characterization of the 5-HT1AR and 5-HT2AR-mediated ERK activation, respectively, in both individual and co-expressed cell lines.
Concentration-inhibition curves were obtained for Ketanserin in the 5-HT2AR clone, and the IC50 values ranging from 9–50 nM (see Supplementary figure 2A). With 50 nM Ketanserin present, about half of the initial ERK activation was observed, corresponding to the measured IC50 value in this system.
Overall, 5-HT1AR displayed transient ERK activation whether expressed alone or co-expressed with 5-HT2AR; whereas 5-HT2AR exhibited sustained and intense response of ERK phosphorylation, either alone or co-expressed with 5-HT1AR. Although co-stimulation of 5-HT1AR and 5-HT2AR by both Xaliproden HCl and TCB-2 enhanced the ERK activation, it did not vary much from 5-HT2AR-mediated ERK induced by TCB-2 (Figure 4, ,55).
Interestingly, prolonged exposure (2 hours) of the preparation to the antagonist Ketanserin (50 nM) before treatment with TCB-2 alone, or with both TCB-2 and Xaliproden HCl, not only attenuated the response but also switched it from sustained to transient. Addition of Xaliproden HCl did not increase the response, indicating that the 5-HT2AR signaling pathway may dominate the ERK activation in the co-expressing cell line, which is in agreement with our result showing that co-stimulation of 5-HT1AR and 5-HT2AR by both Xaliproden HCl and TCB-2 is similar to 5-HT2AR-mediated ERK induced by TCB-2.
Literature review enabled the construction of the model structure shown in Figure 1 as a simplified cartoon diagram. The corresponding detailed kinetic scheme is given in Figure 2. Note that since the Adenylyl Cyclase (AC) inhibition induced by 5-HT1A agonists did not seem to be responsible for the mitogenic effects of 5-HT1A agonists, the effect of cAMP on MAPK activation was not considered. In addition, we did not include the PLC effect because PLC activation by 5-HT1AR is cell-type specific (Liu and Albert 1991; Boddeke, Fargin et al. 1992; Varrault, Bockaert et al. 1992; Ni, Panicker et al. 1997) and is not clear in HEK 293 cells (this is indicated by coloring the reactions involving PLC via 5-HT1AR in a light gray color in Figure 2). For the Gq-coupled 5-HT2A receptors, activation of PLC appears to be almost universal. Thus, the model suggests that 5-HT2A receptors activate PKC via Gαq/11, which in turn acts on Ras and Raf; and Gi and Gq converge downstream of Ras. The ERK regulation consists of positive and negative feedback loops, as shown in green and red thick lines (Figure 1), respectively. Elevation of the active form of MAPK (ERK) stimulates PLA2 (Lin, Chuang et al. 2003), which releases AA and furthermore, AA has a positive effect on PKC. These pair-wise connections create a potential positive feedback loop in the regulation of MAPK. On the other hand, the two key phosphatases in the systems are PP2A and MKP. PP2A is a broad-specificity Ser-Thr phosphatase that dephosphorylates both Raf and MEK. MKP is a dual-specificity phosphatase that dephosphorylates both Tyr and Thr residues on MAPKs, and whose expression is transcriptionally regulated by MAPK. MAPK also phosphorylates MKP, which reduces ubiquitination and degradation of MKP. The increased amount of MKP produces a negative feedback.
The temporal dynamics of the proposed reactions scheme (Figure 2) are described by a deterministic Ordinary Differential Equation (ODE) model. The ODE model consisted of 112 species, 228 parameters and 128 real reactions. Note that four reactions involving PI generation were included but not active in the model. PIP2 serves as a substrate for two apparently independent signaling mechanisms: one involving activation of PIP2-directed PLC and the other one involving activation of PIP2-directed PI3K. In most cells, the levels of PIP2 do no rise dramatically in response to agonist stimulation, possibly due to containing substantial concentration of PIP2 (Stephens, Jackson et al. 1993). The PIP2 level therefore was kept constant. Details of the rate equations, parameter values and parameter estimations can be found in the Supplementary information. The goals for this model are to reflect the experimental data measured in this study, and to provide insights into mechanisms that drive the observed phenomena. The reaction network diagram of the 5-HT1A and 2AR signaling pathways was shown in Figure 2.
Comparison of the simulation results with experimental data enhances understanding of the processes underlying the signaling pathways. Two distinct types of simulations were carried out using ordinary differential equations (ODE) (Schoeberl, Eichler-Jonsson et al. 2002): 1) steady state, in which the variables are determined when no protein concentration changes for the active states prior to the ligand stimulation; and 2) time course, in which the values of the variables are determined as a time series upon ligand treatment, using the steady-state concentrations as starting values.
5-HT1AR- and 5-HT2AR-mediated ERK activation was simulated for each receptor alone, and together. The initial normalized simulation results for 5-HT1AR and 5-HT2AR-mediated ERK phosphorylation qualitatively match the experimental data we obtained (Figure 6 A and B).
For 5-HT1AR, the ligand efficiency, β, and the R*/R ratio, J, parameter sets of (J, β) = (0.01, 41), and (0.25, 2) explained the experimental data the best; we chose the model with (J, β) = (0.01, 41) for analysis, in view of the high efficacy of the agonist used in the experiments.
TCB-2 was introduced as a high affinity 5-HT2AR agonist (Ki is 0.75 nM for human receptors) (Tocris Inc. (McLean, Parrish et al. 2006)). However, there is no substantial information regarding its efficacy, KL and KH. Experimentally we found that 1 nM of TCB-2 was able to elicit considerable activation of ERK in HEK293 cells stably expressing human 5-HT2AR, but had no effect on 5-HT1AR expressing cells (Supplementary figure 1B). Models with (J, β) = (0.6, 357) and (J, β) = (0.25, 357) matched the experimental data best and we chose (J, β) = (0.6, 357) for the analysis. Note that although the concentration of TCB-2 (10 nM) used in the experiment is 100 fold lower than that of the 5-HT1A agonist Xaliproden HCl (1 μM), the ERK activation response is much greater in amplitude, and longer in duration than the Xaliproden-triggered ERK activation (5-HT1AR specific).
To create the MAPK activation model through both the 5-HT1A and 5-HT2A mediated pathways, we combined the two models described above. The basal MAPK level is about 17 and 8 fold higher than individually simulated 5-HT1AR and 5-HT2AR alone, respectively, suggesting a synergistic effect. However, this was not confirmed in our experimental system, possibly because the basal activity is reduced by constitutive desensitization which was not taken into consideration in our mathematical model.
Differences in ERK activation observed between the two receptor-activated pathways were characterized in the simulations of the corresponding models. One example is the difference in time-course of ERK activation, i.e., sustained for the 5-HT2A-triggered response, and transient for the 5-HT1A response. Another example is the effect of the antagonist Ketanserin (at 10 nM (data not shown) and 50 nM) on the nature of the response elicited by 10 nM TCB-2 on the 5-HT2A receptor. The intensity of ERK activation was significantly attenuated at 50 nM of Ketanserin (Figure 7), which is in agreement with our experimental data, in which about half of the ERK activation was inhibited in the presence of 50 nM of Ketanserin. In addition, the signal became transient in duration when Ketanserin was at 50nM, which again corresponds to our experimental observation.
The effects of the positive regulators of the MAPK feedback loop, PKC and PLA2, were addressed as follows: deleting PKC in the 5-HT2A receptor model greatly reduced the activation of ERK. Similarly, deleting PLA2 in the combined 5-HT1A/2AR model caused a significant reduction in the magnitude of ERK activation (Supplementary figure 4). Notably, similar deletion experiments in silico for the negative regulators of MAPK, MKP and PP2A did not lead to a feasible simulation, indicating that the cells could not survive under such conditions.
Sensitivity analysis and parameter optimization were performed as described in the Methods section. The robustness properties of ERK1/2 phosphorylation were evaluated from the statistics of the sensitivity coefficient ‘population’: if the mean of the population is close to zero, then a small standard deviation indicates robustness (since most sensitive coefficients will then be close to zero). Illustrative results are summarized below for specific pathways.
Forward sensitivity analysis revealed that the population mean is close to zero (the absolute mean |<S>| = 0.079, and standard deviation σ = 0.151), suggesting that the model is generally robust and only very sensitive when some of the parameters are perturbed. 29 out of 228 parameters (~12.72%) had normalized absolute sensitivity, |S| > 0.3 (Supplementary Table 4). Parameters with the largest impact are k82, k63, and k51, which are involved in the dephosphorylatation of phosphorylated MKP, ERK and Rafstar, respectively, suggesting significant effects for the phosphatases MKP and PP2A. Local optimization of these 29 parameters led to an improved R-squared value (R2) from −1.2350 to 0.9825 (Figure 6C).
Time-dependent sensitivity analysis of phosphorylated ERK to small changes in initial concentrations of other species, identified 10 out of 112 (~8.93%) with |S| > 0.2 (|<S>| = 0.0564, and σ = 0.178) (Supplementary Table 5). These include MKP, ERK, PP2A, MEK, Raf, RasGAP, RasGDP, Src, Shc, and Grb_sos. MKP and PP2A are the most critical species with the largest impact on the model.
Forward sensitivity analysis (|<S>| = 0.045, σ =0.112) identified 18 of 228 kinetic parameters (~7.89%) with |S| > 0.3 (Supplementary Table 4). Local optimization of these 18 parameters improved the R2 from −0.29813 to 0.9919 (Figure 6D). Most are turnover rates involved in the activation and deactivation of ERK_PP by MEK_PP and MKP_PP, respectively. The most sensitive parameters are k82 and k63, with |S| > 0.6. These parameters participate in the MKP_PP and ERK_PP dephosphorylation respectively, suggesting that the ERK activation dynamics may be determined by the induction of MKP through activation of ERK and MKP_PP proteolysis (Supplementary Table 5). Notably, optimizing k82 alone dramatically improved the R2 from −0.29813 to 0.9924. As k82 is the turnover rate constant of MKP_PP degradation, this result suggests that the dephosphorylation of MKP_PP may play a major role in the sustained activation of MAPK. However, the experimental value of k82 is currently unavailable.
Sensitivity analysis of ERK1/2 phosphorylation to small changes in initial concentrations identified 9 out of 112 species (~8.04%) with |S| > 0.2 (|<S>| = 0.03, and σ = 0.11). These are cpxERK_MEK_PP, MKP, ERK, PP2A, cpxMEK_Rafstar, MEK, Rafstar, cpxRasGTP_Raf, and Raf. Again, PP2A and MKP are the relatively most influential species in this list (Supplementary Table 5).
Forward sensitivity analysis of the combined 5-HT1A/2AR-mediated MAPK pathway produced |<S>| = 0.0186 and σ = 0.064. These values are much lower for 5-HT1AR or 5-HT2AR alone, suggesting that there are fewer sensitive parameters in the co-stimulated pathway. Therefore we chose a higher threshold and 10 out of 228 parameters (~ 4.39%) had |S| > 0.1 (Supplementary Table 4). The parameters exhibiting the largest impact are k82 and k63.
Sensitivity analysis of phosphorylated ERK with respect to small changes in initial concentrations of other species identified 6 out of 112 (~5.36%) with |S| > 0.1 (|<S>| = 0.015 and σ = 0.049). These are cpxERK_PP_MKP_PP, ERK_PP, cpxERK_MEK_PP, MEK_PP, MKP and PP2A. MKP and ERK_PP, In particular, have the highest sensitivity values (Supplementary Table 5).
Overall, k63 and k63r are the most sensitive parameters in 5-HT1AR, 5-HT2AR, and 5-HT1A/2AR models, implying that the induction of MKP by active ERK is the key element in determining the dynamics of ERK phosphorylation. Furthermore, in all three models, the species initial concentrations with the largest impact are involved in phosphatase reactions (Supplementary Table 6).
Parameter perturbations were performed to describe how variations of initial concentrations of species identified from the previous sensitivity analyses affect levels of ERK activation. A range of values (for example, 0.1 ~ 1 μM, see Figure 8) for these species initial concentrations was scanned and iterative simulation results were summarized in one figure to compare their significance.
Parameter scan was done for ten species identified from the sensitivity analysis of the 5-HT1AR model (MKP, ERK, PP2A, MEK, Raf, RasGAP, RasGDP, Src, Shc and Grb2_sos). Variations of PP2A, RasGAP and MKP were found to yield opposite patterns of ERK phosphorylation, i.e. activation vs. inhibition. Most notably, about 5-fold increases of Raf and RasGDP changed the transient ERK activation to sustained response. These two species are involved in integrating Gi and Gq pathways. Our experimental data showed transient activation of ERK by 5-HT1AR (Figure 4B, ,5B),5B), thereby suggesting lower values of RasGDP and Raf were induced in 5-HT1AR signaling pathway.
Parameter scan was also carried out for the nine species identified from the sensitivity analysis of the 5-HT2AR model (cpxERK_MEK_PP, MKP, ERK, PP2A, cpxMEK_Rafstar, MEK, Rafstar, cpxRasGTP_Raf, and Raf). Interestingly, variations of cpxMEK_Rafstar, Rafstar and cpxRasGTP_Rafstar showed that initial values of those proteins cannot be zeros, suggesting an intrinsic activation of the 5-HT2AR–mediated ERK pathway. Increases of Raf and RasGDP greatly elevated ERK activation. Comparing the 5-HT1A and 5-HT2A models, Ras and Raf (at 0.1 μM each) produced a transient ERK activation by the 5-HT1AR, but a sustained one by 5-HT2AR. Note that 5-HT2AR has higher basal activation of ERK than 5-HT1AR, in addition to higher J (J=0.6 in 5-HT2AR and J=0.01 in 5-HT1AR). We note, however, that positive feedback regulation of PKC, activated by AA, is missing in the 5-HT1AR model. Finally, the amount variation of the species cpxERK_MEK_PP changed the active ERK to look like transient activation (see Figure 2, Reactions 24, 25).
For a more detailed look at how the change of positive regulators affect levels of ERK activity, parameter scan was carried out for PKC and PLA2 in 5-HT2AR and 5-HT1A/2AR models, although they were not among the most sensitive state variables (see Supplementary figure 4 and for the parameter scan of PKC in 5-HT1A/2AR, see Figure 8D). Moreover, knockout experiments showed that the ERK activity was significantly decreased when PKC and PLA2 were absent in 5-HT2AR and 5-HT1A/2AR cells, respectively, suggesting the importance of feedback loops enabled by these two players (Supplementary figure 4).
Representative results of parameter scan are shown in Figure 8. Note that the Y axes are in different scales for the different figures. The larger effect of changes in MKP and PP2A (Figure 8B, 8C) than in PKC and PLA2 (Figure 8D, Supplementary figure 4A and 4B) suggests that the variations of negative regulators have a higher impact on the ERK activation than those of positive regulators.
Little is known about the interactions between 5-HT receptor subtypes, although Andrade & Nicoll (Andrade and Nicoll 1987), for example, had demonstrated long ago that 5-HT has multiple actions on single neurons in the CA1 region of the hippocampus and that they are mediated through distinct mechanisms and signaled by different serotonin receptors. Such findings also appear in the more recent literature. Berg et al. (Berg, Maayani et al. 1996) have shown that Gq-coupled 5-HT2C receptors can reduce Gi-coupled 5-HT1B receptor-mediated inhibition of AC. Johnson-Farley et al. (Johnson-Farley, Kertesy et al. 2005) have found that Gq-coupled 5-HT2A and Gs-coupled 5-HT7A receptors positively interact to activate ERK and Akt in PC12 cells. But such individually characterized interactions are not usually integrated into a system-level scheme that describes the evolution of a complex phenotype such as ERK MAPK activity. The coexpression of 5-HT1A and 5-HT2A receptors in the HEK293 cells allows us to evaluate important elements in their interacting pathways in the activation of ERK. As it is known that the kinetics and amplitude of ERK signaling can regulate various cell fate decisions (Marshall 1995; Werlen, Hausmann et al. 2003), our present findings regarding the duration of ERK activation as a function of receptor population are especially interesting.
Differences in duration, magnitude and compartmentalization of ERK MAPK activity, as observed from simulations of our model, are known to provide signaling specificity in cell fate decisions. Regulators of ERK duration include PKC, Ras-GAP, and Sprouty (Bhalla, Ram et al. 2002) (Sasagawa, Ozaki et al. 2005) (Hanafusa, Torii et al. 2002). Regulators of ERK activation are considered to be mainly scaffold proteins, including KSR (kinase suppressor of Ras) (Therrien, Michaud et al. 1996) (Morrison 2001), E3 ubiquitin ligase IMP (impedes mitogenic signal propagation) (Matheny, Chen et al. 2004), and MP1 (MEK partner1) (Teis, Wunderlich et al. 2002). Spatial regulators of ERK activity, such as PEA-15 (Formstecher, Ramos et al. 2001) and Sef (similar expression to FGF genes) (Torii, Kusakabe et al. 2004), play roles in ERK nuclear localization, thereby influencing the IEGs expression. Using the models in the simulations, we were able to predict the duration of ERK activity by considering the drug efficacy and constitutive activity inside the cells (e.g. Figure 6). However, we did not consider spatial factors that influence the ERK activity, and did not include the scaffold proteins in our models (the approximations inherent in these simplifications were necessary to keep the modeling problem computationally tractable for this initial stage). To delineate some of the factors producing the distinct patterns of ERK activation, we used the simulations to test the effects of various parameters and state variables, and examined as well the effects of deletion of positive and negative regulators in the feedback loop of MAP kinase signaling pathway. Thus, upon deletion of PKC in silico the basal level of phosphorylated ERK started very low and instead of increasing, decreased in the early time period. In silico deletion of PLA2 in the combined 5-HT1A/2AR model, still showed activation of ERK, but with much less intensity compared with wild type, and it might be that PKC remained active in the model which contributed to the feedback.
Parameter perturbation of state variables based on sensitivity analysis, allowed us to further reveal intrinsic dynamics of the system. In 5-HT2AR, if initial amounts of proteins like Rafstar, cpxRasGTP_Raf, cpxMEK_Rafstar and cpxERK_MEK_PP are zero, instead of activation we observe a drop in the calculated values of ERK, suggesting the presence of a basal level of active ERK. In addition, converging points of Gi and Gq pathways such as RasGDP and Raf were shown to play significant roles in the ERK dynamics. Increases of RasGDP and Raf, particularly in 5-HT1AR model, switched the duration from transient to sustained. All three models (for 5-HT1AR, 5-HT2AR, and combined 5-HT1A/2AR triggered pathway) showed that levels of MKP and PP2A must be finely regulated or else will lead to opposite results of ERK activation, i.e. increase vs. decrease. Moreover, variations of the initial concentration of the species cpxERK_PP_MKP_PP in the 5-HT1A/2AR model, displayed a switch-like behavior of ERK activation, and our simulation indicated that the ratio between ERK_PP and MKP_PP determine the dynamics of ERK activation, which is consistent with experimental studies (Lin and Yang 2006).
Note that state variables, such as ligand-receptor (LR), or ligand-receptor-G protein (LRG) complexes, were not detected as sensitive variables. These state variables are located upstream in the pathway, producing instant changes as soon as the ligand was introduced, but ERK activation did not dramatically change after 300 seconds of ligand addition, so that the sensitive state variables and rate constants became those involved in the reactions around 300 seconds. Nevertheless, we note that the ligand-specific parameter most correlated with both measures of response generation was the effectiveness with which the ligand induces an active receptor conformation. This is an important finding in relation to the specific mechanisms of hallucinogens acting on the 5-HT2A receptors, because hallucinogens, such as LSD, differ in their modes of binding to the 5-HT2A receptors compared to non-hallucinogenic agonists, and we have already shown (Ebersole, Visiers et al. 2003) that the pharmacological efficacy at these receptors depends on specific modes of ligand-receptor interaction. Moreover, we found here that the nature of responses at different levels of efficacy is affected by several cell-specific parameters including receptor and G protein expression (as the ratio of G protein and receptor increases, so does the degree of precoupling - see Supplementary figure 3), as well as parameters involved in the G protein activation loop, the equilibrium ratio of active and inactive receptor and the rate and efficacy of receptor-G protein coupling.
It is undoubtedly true that in the current state of this art of cell signaling models, there are major pitfalls and concern about such modeling efforts. Chief among these are (i)-the gathering and curation of the qualitative (pathway) data required for the models, and the quantitative data required for the simulations; (ii)-the sheer paucity of such data; and (iii)-the conceptual and methodological problems related to modeling interactions of small numbers of molecules in spatial subdivisions of the cell, at relevant time points. We are continuing to address the first problem mentioned above with the development of SigPath (Campagne, Neves et al. 2004). The second issue, regarding the availability of the data, has been discussed and reviewed responsibly in the current literature in systems biology (e.g., see discussions and articles in (Bock 2002; Gutenkunst, Waterfall et al. 2007) and references therein; as well as in (Slepchenko, Schaff et al. 2002; Slepchenko, Schaff et al. 2003; Sauro and Kholodenko 2004; Wachman, Poage et al. 2004). Still, the value of modeling with the available levels of data for both hypothesis testing and experiment design have been considered to remain highly significant (in particular, see discussions in (Bock 2002; Palsson 2006) and as illustrated more recently for several specific signaling pathways (Hautaniemi, Kharait et al. 2005; Janes, Albeck et al. 2005; Hendriks, Cook et al. 2006; Kumar, Hendriks et al. 2006; Gutenkunst, Waterfall et al. 2007)). Our results presented here should add some weight to the argument in favor of application of currently available information to explore mechanistic hypotheses and suggest experimental probing of inferences that could not have emerged without the system simulation.
The methodological problems raised in point (iii) above remain the most formidable. These methodological problems are not only subjects of vigorous discussion in the literature (Ideker, Winslow et al. 2006; Kriete and Eils 2006; Palsson 2006) (Palsson 2006), but also of intense activity. Because of the many advantages that quantification and formal representation in mathematical models could offer, we are maintaining close attention on this topic in the context of mechanisms of drugs of abuse, as discussed in the recent review (Niv, Skrabanek et al. 2006). We plan to proceed carefully in tracking such approaches and contributing to their development, as they appear likely to be increasingly necessary for an understanding of how the complexity of hallucinogen mechanisms in which multiple functions of a cell are coordinated and regulated.
Our study is an illustration of current capabilities and shortcomings of a new, emerging field of study, providing the first detailed kinetic interaction map and quantitative representation of a fundamental component of cell signaling- the interplay of two receptors responding to the same neurotransmitter- with the objective of identifying the characteristics of that interaction. This hypothesis (regarding the interplay) is explored through simulations. It is possible that the proposed signaling pathways used by these receptors may be quite different in neurons than in HEK293 cells. On the other hand, the reproduction of the phosphorylation behavior of the target suggests that this is a plausible model of receptor interplay. Clearly, we have left out elements of the signaling pathway-both some of those that are known, and the currently unknown ones- so that if the parameters or some pathway components are actually in error, the agreement we find with the experimental results may mean either that there is a compensation that balances out the errors, or that the other details are not essential for the behavior simulated by this model. The knockdown experiments provide further experimentally testable hypotheses. In summary, our study offers new perspectives and attractive working hypotheses for the field. The ability to predict the correct response should help much in the design and evaluation of therapeutic strategies that target specific components of the pathways activated by hallucinogens. The value of integrating computational and experimental results illustrated by this work is clearly dependent on the ability of modelers and theorists to communicate with the experimentally derived knowledge and with the researchers who uncover it. We can only hope that our manuscript will help initiate and support this dialog in the neuropharmacology field.
We thank Dr. Xin-Yun Huang for access to experimental facilities and patient guidance of C-w C’s work in his lab, Dr. Ravi Iyengar and the Virtual Cell group (National Resource for Cell Analysis and Modeling (NRCAM)) for general consultation; and Angela Chao for generation of some data in early stages of the experimental effort. The work was supported in part by NIH grants P01 DA012923 (to HW) and R01 MH063162 and funds G67673 (to J.A.S). Computational resources of the HRH Prince Alwaleed Bin Talal Bin Abdulaziz Alsaud Institute for Computational Biomedicine are gratefully acknowledged.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errorsmaybe discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.