Deciphering the mechanism behind Fibroblast Growth Factor (FGF) induced biphasic signal-response profiles
- Jitendra Kanodia†1Email author,
- Diana Chai†1,
- Jannik Vollmer1,
- Jaeyeon Kim1,
- Andreas Raue1,
- Greg Finn1 and
- Birgit Schoeberl1Email author
© Kanodia et al.; licensee BioMed Central Ltd. 2014
Received: 21 December 2013
Accepted: 28 April 2014
Published: 15 May 2014
The Fibroblast Growth Factor (FGF) pathway is driving various aspects of cellular responses in both normal and malignant cells. One interesting characteristic of this pathway is the biphasic nature of the cellular response to some FGF ligands like FGF2. Specifically, it has been shown that phenotypic behaviors controlled by FGF signaling, like migration and growth, reach maximal levels in response to intermediate concentrations, while high levels of FGF2 elicit weak responses. The mechanisms leading to the observed biphasic response remains unexplained.
A combination of experiments and computational modeling was used to understand the mechanism behind the observed biphasic signaling responses. FGF signaling involves a tertiary surface interaction that we captured with a computational model based on Ordinary Differential Equations (ODEs). It accounts for FGF2 binding to FGF receptors (FGFRs) and heparan sulfate glycosaminoglycans (HSGAGs), followed by receptor-phosphorylation, activation of the FRS2 adapter protein and the Ras-Raf signaling cascade. Quantitative protein assays were used to measure the dynamics of phosphorylated ERK (pERK) in response to a wide range of FGF2 ligand concentrations on a fine-grained time scale for the squamous cell lung cancer cell line H1703. We developed a novel approach combining Particle Swarm Optimization (PSO) and feature-based constraints in the objective function to calibrate the computational model to the experimental data. The model is validated using a series of extracellular and intracellular perturbation experiments. We demonstrate that in silico model predictions are in accordance with the observed in vitro results.
Using a combined approach of computational modeling and experiments we found that competition between binding of the ligand FGF2 to HSGAG and FGF receptor leads to the biphasic response. At low to intermediate concentrations of FGF2 there are sufficient free FGF receptors available for the FGF2-HSGAG complex to enable the formation of the trimeric signaling unit. At high ligand concentrations the ligand binding sites of the receptor become saturated and the trimeric signaling unit cannot be formed. This insight into the pathway is an important consideration for the pharmacological inhibition of this pathway.
Signaling pathways are arguably one of the most important components driving how biological systems respond to environmental cues [1, 2]. The ability of cells to perceive and respond to their microenvironment is the basis of development, tissue repair, and immunity as well as normal tissue homeostasis. Errors in cellular information processing contribute towards diseases such as cancer. By understanding the intricacies of cell signaling and processing, diseases may be treated more effectively.
The FGF pathway plays a pivotal role in the stimulation of cell proliferation, cell migration and differentiation of a large number of cell types including muscle, neurons, cartilage and bone cells [3, 4]. FGF ligand - receptor signaling is regulated both by primary sequence differences between the 18 FGF ligands, the 7 main FGF receptors (FGFR1b, FGFR1c, FGFR2b, FGFR2c, FGFR3b, FGFR3c and FGFR4), by temporal and spatial expression patterns of FGFs, FGFRs and HSGAGs and FGF Binding Proteins. Tissue-specific alternative splicing in the second half of Ig-like domain 3 (D3) of fibroblast growth factor receptors 1–3 generates epithelial FGFR1b-FGFR3b and mesenchymal FGFR1c-FGFR3c splice isoforms. This splicing event establishes a selectivity filter to restrict the ligand binding specificity of FGFRb and FGFRc isoforms to mesenchymally and epithelially derived fibroblast growth factors (FGFs), respectively . FGF Binding Proteins (FGF-BP) have been described to function as a chaperone molecule that can mobilize FGF locally and present it to the FGF receptor. The FGF-BPs have been described to enhance proliferation and signaling in NIH-3T3 cells .
HSGAG has been assigned multiple roles: it is known to serve as a co-receptor essential for signaling, as transport mediator to increase the local concentration of growth factors close to the cell surface or as a regulator to accelerate or attenuate signaling . HSGAG are known to be essential for FGF signaling  and typical HSGAG levels on the cell surface are with 105-106 molecules per cell  much higher than typical FGFR levels (<50,000 molecules per cell Merrimack unpublished data). Thus, the benefits of understanding this pathway and the role of HSGAG in regulating FGF signaling are several fold: reveal greater insights into the fundamental principles of signaling pathway regulation by HSGAG in the case of FGF and more generally to understand the effect of inhibitors targeting the FGF pathway that are currently in development [10, 11].
Interestingly, some cells respond to activation by FGF ligands atypically [3, 4, 19]. Instead of the typical saturable response (Figure 1A), these cells respond in a biphasic manner (Figure 1B). Specifically, cellular response increases from low to intermediate levels of FGF ligand but then decreases at high levels of the ligand. For instance, Williams et al. demonstrated that FGF2-induced neurite outgrowth reaches a maximum at intermediate concentrations of FGF2 and that the outgrowth decreases at higher levels of FGF2 . Garcia-Maya et al. demonstrated that FGF-induced proliferation of NIH3T3 cells also reaches a maximum at intermediate concentrations of FGF2 [4, 20]. Similar results have been shown for fibroblasts and osteoblasts as well . Notably, biphasic response to FGF ligands has most commonly been reported at the level of cellular phenotype, while the underlying molecular signaling events that lead to the biphasic response remain unexplored.
Recently, Zhu et al. measured protein levels using Western Blots and indicated that FRS2- (FGFR adaptor protein) and ERK-phosphorylation is biphasic in response to FGFR activation by FGF2 . Thus, they provide hints that FGF signaling response does not follow the typical Michaelis-Menten reaction kinetics. However, owing to the labor-intensive and low throughput nature of the experimental technique, protein levels were measured at a small number of time points and only for a small number of ligand concentrations. Moreover, they did not provide any mechanistic insight into how the FGFR pathway activation cascade drives this biphasic behavior.
Efforts to model the FGFR pathway previously have primarily focused on the interactions taking place on the surface of the cell while completely ignoring the intracellular cascade [21–24]. Their efforts were directed towards understanding the effect of adding heparin, a soluble source of HSGAGs, on FGFR pathway activation. However, the time-course and dose–response of pERK to activation by ligand alone remains unexplored. In contrast, Yamada et al. built a combined mathematical model for the extracellular and intracellular components of the FGFR pathway but completely ignored the biphasic nature of response upon activation by the ligand .
Here, we investigate the mechanistic rationale for biphasic response to FGFR pathway activation utilizing a combination of high-density signaling data and Ordinary-Differential-Equation (ODE)-based mathematical modeling. We report the dynamic pERK response to FGF2 stimulation using a fine-grained time grid over 120 minutes. The model was calibrated to the experimental results using a published particle swarm optimization (PSO) algorithm  combined with a novel feature-based constrained optimization. Finally, we apply Markov-chain Monte Carlo sampling to explore nonidentifiability and uncertainty in the model calibration . The computational model represents the extracellular and intracellular components as well as the feedback regulation of the pathway. To validate the model, multiple simulations were conducted to predict signaling response to extracellular and the intracellular perturbations of the pathway. The model predictions were validated by in vitro experiments. We demonstrate that without any fitting, the model explains each of these perturbation experiments. We demonstrate that the complex protein interactions at the cell surface are necessary to explain the observed biphasic dose–response while the negative feedback loop from pERK to FRS2 controls the magnitude of the biphasic response. At low to intermediate concentrations of FGF2 there are sufficient free FGF receptors available for the FGF2-HSGAG complex to enable the formation of the trimeric signaling unit. At high ligand concentrations the ligand binding sites on FGFR become saturated and the trimeric signaling unit cannot be formed because binding of FGF2-HSGAG is weak, thereby leading to a decrease in pERK response.
Results and discussion
These detailed measurements increase the identifiability of the system, and allows the signaling cascade to be modeled mechanistically using a set of Ordinary Differential Equations (ODEs) describing the biochemical reactions. We constructed a mathematical description of the FGFR pathway that includes all the ligand-receptor interactions on the cell surface as well as the intracellular MAP kinase signaling cascade.
Description of the computational model
It should be recognized that the reactions included in this scheme are not necessarily direct protein-protein interactions. Some of these reactions represent a combination of multiple interactions within the cascade. The underlying principle of such a reductionist approach is that for large non-identifiable systems like signaling pathways, most components remain unmeasured. Therefore, a reduced mathematical model with appropriate features can better describe the experimental results while allowing for easy interpretation as compared to a model that accounts for each protein-protein interaction individually. The applicability of such reduced models is tested by verifying whether it can quantitatively fit experimental results (Figure 2) with physically meaningful rate constants followed by validation of independent perturbation experiments without any fitting.
Most of the kinetic rate constants described in the mathematical model remain unmeasured to date. Even for those rate constants that have been measured, for instance the ligand-receptor binding constants, a wide-range of values have been reported in the literature depending upon the experimental system and technique used [22, 32, 34]. Therefore, model calibration against the experimental data with physically meaningful parameters is critical in order to gain a better understanding of the system.
Calibration of the computational model using particle swarm optimization and feature-based constraints
As is typical for signaling pathways, the mathematical model is parametrically non-identifiable with many more components than the number of measurements available. Accordingly, there exists no unique solution for the parameter values and purely deterministic parameter estimation techniques will fail to provide good estimates unless the initial parameter guess is in the vicinity of the global optimum itself. Therefore a global parameter estimation technique, Particle Swarm Optimization (PSO), was utilized to estimate the parameters of the model . The approach allows for a more exhaustive and efficient exploration of the objective function manifold to find good parameter fits. The version of PSO utilized in this report has been published previously and was particularly useful for fitting parameters to signaling models [26, 35]. The only inputs required for this approach are the physical constraints for the parameter values.
Given the network structure, a critical aspect of estimating the parameters is the choice for the objective function to be minimized. We noted that minimization of the most commonly used objective function, root mean square error (RMSE) between experimental data and simulation predictions was insufficient and the estimated parameters provide poor solutions. Therefore we developed a hybrid objective function that accounts for RMSE as well as differences between qualitative features of experimental and simulation data. Figure 2A summarizes the feature-based constraints used in the objective function. The objective function was penalized if the pERK level reached peak values in response to ligand stimulation too early or too late (detailed penalty function expressions in Materials and methods Section 3). Similarly, the objective function was penalized if the pERK dose response at the relevant time points was not qualitatively biphasic. Minimizing the hybrid objective function using PSO led to physically meaningful parameter estimates and was used to obtain the parameters used in this report. In contrast, using RMSE as the objective function provided parameter estimates that had low objective function values but provided biophysically meaningless solutions. A detailed comparison between the two objective functions used and the corresponding parameter estimates obtained is described in Materials and methods Section 3. Based on these results, we propose a more general hypothesis that optimization of all signaling pathway models might benefit from utilizing a combination of qualitative constraints and RMSE values as compared to the simple RMSE-based objective function. The applicability of this approach to other models remains to be tested and is the subject of future research.
The model recapitulates pERK activation data at all different FGF2 concentrations. The model captures the experimental observations that at low FGF2 concentrations (below 4 ng/ml), the peak level of ERK activation increases with FGF2 concentration; however, at higher concentrations (above 4 ng/ml), the peak level of ERK activation remains constant (Figure 2B). Careful observation reveals that at high FGF2 concentrations, the time to peak ERK activation decreases as FGF2 levels increase. Most importantly, the model also captures the fact that pERK de-activation varies with FGF2 concentration. Specifically, at intermediate levels of FGF2, pERK levels reach a peak after 5 minutes and then slowly decrease over the next hour. In contrast, at the high levels of FGF2, pERK levels peak before 5 minutes and decrease to low levels within 20-30 minutes. Therefore the dose response curve at time points later than 10 minutes shows a strong biphasic response (Figure 2C). Thus, the reduced model adequately captures all of the essential features of pERK dynamics in response to FGF2 stimulation.
As described previously, even with the fine-grained pERK response measurements, the model is highly unidentifiable and therefore multiple parameter sets that might correspond to different response mechanisms explain the data equally well (Additional file 1: Figure S1A). Accordingly, before utilizing the aforementioned model and parameter sets shown in Figure 2 for further investigation of the underlying mechanisms, we tested and validated the model against multiple perturbation experiments. The validation will help rule out parameter sets that fit the training data-set well but do not describe the FGFR signaling pathway. Specifically, we perturbed the extracellular and intracellular components of the model and tested model validity by comparing the predictions with in vitro experimental measurements of pERK. The objective was to validate the underlying mechanisms predicted by the model rather than exact quantitative numbers. Therefore, the model was not trained on any these new experimental results and all model predictions were qualitatively compared to the experimental results.
Addition of soluble heparin decreases the pERK response at low and intermediate FGF2 levels but increases the pERK response at high levels of FGF2
One of the essential components of the model is the binding of FGF2 to HSGAGs. It is known that FGF2 also binds to soluble heparin and this binding in solution can compete with the binding of FGF2 to HSGAG. Soluble heparin can also rescue signaling behavior in cells that have been stripped of cell surface HSGAGs. Therefore, to validate the ligand binding and signaling complex formation module of the model, we tested how the pERK response changes with addition of soluble heparin to media. The effect of heparin on the FGF pathway has been explored previously at the level of surface interactions and ligand-receptor complex formation [21, 24]. However, those results and interpretations have never been extended to include the effect on the intracellular signaling cascade. Using our model, we explore the effect of heparin on pERK dynamics at the same fine time-grid as discussed in the previous section (Materials and methods Section 4). Model simulation of heparin addition to extracellular medium provided some non-intuitive insights (Figure 4A, B solid curves). At low to intermediate levels of FGF2, addition of heparin was predicted to decrease the level of ERK activation at all the time points. This is in line with the role of heparin as a ligand trap. However, surprisingly, at high levels of FGF2, addition of heparin was predicted to increase the level of pERK response. Model prediction for change in time-course and dose–response of pERK upon addition of heparin were validated using in vitro experiments (Figure 4A, C symbols and dotted curves).
Without any fitting, the model accurately captures pERK response to FGF2 in the absence/presence of heparin (Figure 4). Note that the same qualitative predictions were made by all the different parameter sets (Additional file 1: Figure S1B). Therefore, the level of confidence in the explanation provided by the simulations is high: The model indicates that at low levels of FGF2, heparin acts as a ligand trap and reduces the level of FGF2 binding to HSGAG and subsequent formation of signaling complexes while at high levels of FGF2, there is enough excess FGF2 present in the extracellular medium to overcome the ligand trap and to form signaling complexes. Additionally, a small number of heparin-FGF2 complexes bind to FGFR and initiate signaling complex formation, leading to an increase in pERK response. Thus, without fitting any of the parameters for soluble heparin, the model captures pERK response to extracellular perturbations. It is noteworthy here that some experimental data points do not match perfectly with simulation results. We hypothesize that the differences could be due to non-pathway specific interactions or other higher order interactions that are indeed a part of the FGFR pathway but have a small influence on the overall response and are not captured by the current model. Explanation of these differences would require measurement of more proteins within the cascade combined with a more detailed model and will be the subject of future work.
Delayed inhibition of MEK post FGF2 stimulation amplifies the biphasic pERK response
To validate the intracellular ERK activation module of the pathway, MEK was inhibited using a small molecule inhibitor U0126 (henceforth referred to as 'MEKi’) with two different schedules. In the first experiment, MEKi was added to the media at the same time as the ligand itself (Figure 5A). Five different concentrations of MEKi were used in combination with seven different concentrations of FGF2. For each combination of MEKi-FGF2, the pERK response was measured in cells 10 minutes after FGF2 stimulation. As predicted by the computational model only the 0.37 μM and 1.11 μM concentrations of MEKi lead to a significant pERK inhibition across the entire range of FGF2 concentrations. It is noteworthy that the biphasic shape of the FGF2 dose–response curve remained the same for different MEKi levels. pERK levels reach a peak at 4-20 ng/ml of FGF2 and then decrease at higher FGF2 concentrations. These results were in line with the predictions made by the model. The ODEs used to describe MEKi interactions are outlined in Materials and methods Section 5 in the Appendix.
It is not surprising that co-incubation of MEKi with FGF2 was expected to show uniform inhibition of pERK. Therefore, another experiment was conducted where the MEK inhibitor was added 5 minutes post ligand stimulation and pERK response was measured 5 minutes post MEKi addition (Figure 5B). In this experiment, the cells respond to ligand stimulation for 5 minutes under MEKi-free conditions and then for 5 minutes in the presence of MEKi. Thus, this perturbation experiment is a more stringent test of the model since it is required to accurately predict the combination of control as well as perturbation conditions within the same simulation. Similar to the co-inhibition experiment in Figure 5A the model predicted reduced pERK levels – however to a much lesser extent than in the co-inhibtion experiment. Specifically, the model predicted that the peak ERK response at intermediate FGF concentrations is significantly inhibited while the signal at the low and high FGF concentrations are inhibited to a lesser extent. This prediction was validated as showed in Figure 5B in the delayed addition experiment. The fact that the model captures the partial inhibition of pERK in the delayed inhibition experiment indicates that the negative feedback loop from pERK to FRS2 is accurately captured by this parameter set. It should be noted here that there are differences between pERK response under control conditions (i.e. no MEKi) for the two experimental setups. This difference can be attributed to the super-sensitivity of pERK to fluid movement and addition rather than biological response to growth factors or inhibitors. Therefore, it is further emphasized in this case to focus on the qualitative prediction of model compared with experimental results and interpret the changes in pERK response with respect to control conditions rather than absolute values.
It is also noted that not all the parameter sets predicted the effect of MEKi accurately. One of the sets gave almost no response to MEKi while two others incorrectly predicted pERK response to stimulation by FGF2. As expected previously, matching model prediction with the more stringent experimental test of fitting pERK response to FGF2 and MEKi scheduled with a time lag helps identify physiologically relevant parameter sets for the FGFR pathway. Finally, only four out of the seven parameter sets that predicted both perturbation experiments correctly were used for subsequent model interpretation and analysis. Overall, one can be fairly confident that the model significantly captures the behavior of intracellular ERK activation by FGF pathway and that the observed biphasic pERK dose response can be explained by the complex interplay between FGF receptor, FGF2 ligand and HSGAG.
Biphasic phenotypic responses to FGF signaling have been observed in multiple cell types in varying contexts. However, an investigation of the mechanism behind biphasic response at the level of signaling components remains unexplored. Here, we utilize high-throughput measurements of ERK activation for measuring the dynamic signaling response to a wide range of FGF2 concentrations. Such detailed quantification allows us to tease apart subtle changes in ERK activation dynamics and thus provide a more complete picture of the system. The quantification also serves as an excellent training dataset for a mechanistic ODE-based model of the FGFR pathway. Given the lack of identifiability of the model, the global parameter estimation method PSO was utilized with a custom-built objective function to estimate the parameters. Finally, we validated that the identified sets of parameters are representative using an MCMC approach. We demonstrate that the model successfully captures the subtle changes in the timing of pERK peak levels as well as the change in pERK decay rates. To validate the assumptions made to build the pathway network structure and the estimated parameter values, the model was tested against multiple perturbation experiments. Without any fitting, the model accurately predicted changes in pERK response to both extracellular and intracellular perturbation conditions. Thus, we assert that the reduced model captures the essential features of the system appropriately and can be utilized for further investigation. It is noteworthy that the objective of our approach was to build a mathematical model based on commonly understood aspects of FGFR pathway biology and uncover a plausible mechanism for biphasic pERK response. Therefore, we did not exhaustively explore the space of model structures for alternative theoretical models that can also explain the data. Such an exploration can be instrumental for uncovering completely unknown or less well-understood biology and will be the subject of future research.
The computational model was explored in silico and then tested experimentally by using another cell-line. According to the model, a system with a large number of FGF receptors is predicted to show a sigmoidal pERK response (Figure 1A). At high FGF receptor levels, there is no competition between FGF2 binding to FGFRs and HSGAG and thus, as expected, pERK response saturates at high FGF2 levels rather than showing a biphasic response. This model interpretation was verified experimentally by using another cancer-cell line MDA-MB-134 (Figure 6B). This cell line has about 8× more FGFR1 mRNA and ~2× more HSGAG and syndecan-1, a cell-surface proteoglycan containing HSGAGs . Under high FGFR and high HSGAG conditions, the model predicts that the pERK response looks sigmoidal up to 500 ng/ml of FGF2. The model prediction was verified by measuring pERK in MDA-MB-134 cells 10 min post FGF2 stimulation. These results highlight the importance of HSGAG and soluble heparin in regulating FGFR pathway activation. They indicate that HSGAG is not just a passive scaffold that facilitates binding of FGF ligands to FGFRs but actively participates in the local regulation of the dynamics of intracellular signal activation . It is noteworthy at this point that currently there are no good methods for accurately quantifying levels of HSGAGs. HSGAGs are not homogeneous, but rather a mixture of HSGAG species of different lengths and sulfation patterns. The sugar sequence and level of sulfation can impact the level of binding and its role in the signaling complex . Therefore, in the absence of such measurements, the model provides an effective alternative approach for investigating the importance of HSGAGs in modulating signaling response to FGF ligands.
The FGFR pathway signaling model and estimated parameter fit presented in this report is specific for the FGF2 ligand interacting with the FGFR1 receptor. However, the quantitative approach, including measurements, model topology and estimation approaches are generalizable. HSGAGs affect the kinetics of ligand-binding and thus signaling of several other growth factors including HB-EGF, HGF, PDGF, MBPs, TGF-β or wingless (a member of the Wnt family) [40–45]. In the case of HB-EGF binding to the EGF receptor, binding is enhanced at relatively low levels of heparin but higher heparin concentrations lead to inhibition of receptor, as is the case in FGF system . In contrast, despite strong in vivo evidence for a stimulatory role for endogenous HSGAGs in bone-morphogenetic-protein (BMP) signaling, the results of in vitro studies have been inconsistent [46, 47]. Since FGF and BMPs differ in the way they form signaling competent dimers, the biphasic response that can be observed with FGF may not be observable with BMPs or other growth factors . Therefore, the role of HSGAGs in different signaling pathways will need to be evaluated individually.
In summary, this work increases the quantitative understanding of the FGF signaling pathway and its intricate regulation on the cell surface as well as downstream of the receptor, which is important to understand cellular decision making.
Materials and methods
Section 1. Experimental details of quantitative pERK measurements & perturbation experiments
NCI-H1703 cells were obtained from American Type Culture Collection (ATCC CRL-5889). Cells were maintained in RPMI 1640 medium (Lonza) supplemented with 10% fetal calf serum (Tissue Culture Biologicals, Lot # 107197), 2 mM L-Glutamine (Gibco), penicillin (100 U/mL) and streptomycin (100 μg/mL) (Gibco), and grown in a humidified atmosphere of 5% CO2 and 95% air at 37°C.
In vitro signaling studies
For in vitro signaling studies, cells were seeded in 96-well tissue culture plates, allowed to attach overnight, then switched to serum-free media supplemented with 0.5% bovine serum albumin (Sigma-Aldrich) for 16 to 20 hours. For dose-time matrices, serum-starved cells were stimulated with serial dilutions of FGF-2 (R&D Systems, 4114-TC-01 M) starting at 500 ng/mL for 1, 2, 3, 4, 5, 7, 10, 15, 30, 60, 120 minutes. After stimulation, cells were washed with ice cold phosphate-buffered saline (PBS), then lysed in cold 1× lysis buffer from the AlphaScreen SureFire p-Erk1/2 Assay Kit (Perkin-Elmer TGRES10K).
For perturbation with exogenous heparin addition, cells were seeded as described above, then switched to serum-free media supplemented with 0.5% bovine serum albumin with or without 500 μg/mL heparin sodium salt (Sigma-Aldrich H3149-250KU). Serum-starved cells were stimulated with serial dilutions of FGF-2 in the starvation media and lysed as before. For perturbation with MEK inhibition, cells were seeded and starved as described above. Serum-starved cells were treated with one of two protocols: 1) simultaneous addition of a dose matrix of serially diluted FGF-2 starting at 500 ng/mL and U0126 starting at 10 μM (EMD Chemicals 662005-1MG) for 10 minutes followed by PBS wash and lysis, 2) addition of serially diluted FGF-2 followed by addition of serially diluted U0126 after 5 minutes of incubation and finally PBS wash and lysis 10 minutes after ligand addition.
Measurement of phospho-Erk1/2 levels on all samples were performed with the AlphaScreen SureFire pErk1/2 assay kit (Perkin-Elmer TGRES10K), according to the manufacturer protocol in 384-well plate format (Alpha Plate, PerkinElmer 6008350) and read using the EnVision 2013 Multilabl Plate Reader (Perkin-Elmer).
FGF Receptor Expression by qPCR
RNA was isolated from 2×10^6 cells using RNeasy kit (Qiagen) following manufacturer protocols. Genomic DNA was eliminated using deoxyribonuclease (DNase) treatment using DNase I (Roche). 6 g RNA was reverse transcribed to cDNA in a 60 L reaction using a High Capacity cDNA Reverse Transcription kit (Life Technologies Cat#4368814), and samples were stored at -80°C until qPCR.
Primer sequences for FGFRs
qPCR data was analyzed by Viaa 7 Software v1.2.2 . Baseline values were set automatically, and threshold values were kept constant. Samples with Ct values of >35 were considered below the limit of detection. Expression levels were normalized to the average of the housekeeping genes using the delta Ct method, and normalized expression is calculated as 2^-deltaCt. These results are plotted below.
Results - The relative expression of each of the FGF receptor isoforms was measured in NCI-H1703 cells as part of a wider screen. Based on the expression levels, FGFR1c contributed more than 95% of the total while FGFR1b, FGFR2b, FGFR2c, FGFR3b and FGFR3c combined contribute less than 5% of the total. Consistent with previous reports, NCI-H1703 cells predominantly express FGFR1c.
Section 2. Equations and parameters (table of seven sets) of the mathematical model
List of species
Set of ordinary differential equations
Section 3. Comparing different forms of objective function for optimization
One of the most commonly used forms of objective function for minimization using local and global optimization methods is the root mean square error (RMSE). The error term serves as a scale-dependent estimation of the deviation between observed and predicted values. This functional form is the exact equation that needs to be minimized for fitting a linear relationship between input and output. However, for highly non-linear models, such as the ODEs described above, the RMSE-based objective function manifold in the multi-parameter space is rife with a large number of local minima and therefore not suitable for parameter estimation. An alternative to RMSE-based objective functions is to construct hybrid objective functions that combine RMSE error with feature-based constraints incorporated as cost functions. The details of the hybrid objective function used for fitting the pERK data (Figure 2B) to FGFR model (Materials and methods, Section 2) are as follows.
- 1.Total number of phosphorylated molecules of all intracellular molecules should exceed 100 at peak activation (This constraint needs to be satisfied so that the mean field approximation required for building ODE-based models is valid). If this was not fulfilled the following penalty was added
- 2.pERK should be a sizable fraction of total-ERK present in the system (This constraint needs to be satisfied to ensure that the parameters don’t fall into a regime where only a small fraction of pERK molecules get phosphorylated and control the system).
- 3.The number of phosphorylated molecules down the cascade should not decrease by more than 25% within a single step (This constraint is expected to be true since kinases act as catalyst for phosphorylation of the molecule down the cascade and therefore it is expected that the input signal amplifies down the cascade).
where m is equal to the number of FGF concentrations, n+1 is equal to the number of kinases in the system (FRS2, MEK and ERK) and max(kinase(i,j)) or max(substrate(i,j)) is maximum level of upstream or downstream kinase respectively for FGF dose i and kinase j.
- 4.The pERK time response curve should reach a peak between 4-10 min after ligand stimulation
- 5.The pERK time response curve should be smooth and not oscillate like an under-damped second order control system. If there were more than 3 peaks, the following penalty was defined:
- 6.pERK dose-response curve should be biphasic at certain time-points as observed in experiments.
where k is the number of time points, σ(i) is the mimimal concentration of pERK at high FGF2 concentrations at time point i, ζ(i) is the approximate factor by which pERK should decrease from its maximal to its minimal response (determined from the experimental data) and μ(i) is the maximal level of pERK for time point i. ζ was chosen to be equal to [1.0, 1.0, 1.0, 1.0, 1.0, 0.95, 0.85, 0.5, 0.5, 0.9, 0.9] for the different time points.
The final objective function value was calculated as the product of all penalty values and the RMSE error.
Although the final objective function value obtained was, in general, lower for the RMSE-based objective function than the hybrid objective function, parameter estimates from RMSE-based function provided physically meaningless solutions. For instance, 19 out of 20 parameter estimates predicted that less than 2% of total ERK gets phosphorylated by FGFR pathway. In contrast, 12 out of 20 parameter estimates obtained using hybrid objective function satisfied all the experimental and physicality constraints. Each of these fits can be used as good initial estimates for obtaining more refined parameter estimates.
Visual inspection of the fits showed that many of the parameter sets obtained using RMSE-based function showed only minimal to no biphasic dose response. Further investigation revealed that these parameter sets fit the low and intermediate doses of FGF2 extremely well but completely failed to fit the response to high doses. In addition, many fits failed to match the time point of peak phosphorylation.
Based on these observations, we hypothesize that the objective function manifold that uses a combination of RMSE and constraints gets rid of various local minima that correspond to physically meaningless parameter estimates and thus facilitates a much more efficient estimation of parameter values. It can be expected that these findings can be extended to parameter estimation problems for various other signaling pathway models and other systems of ODEs. However, a formal comparison of the RMSE-based and hybrid objective function for benchmarked problems is beyond the scope of this work and remains to be verified in future studies.
Section 4. Additional/modified equations and parameters for the heparin addition experiment (Table 4)
List of additional species and parameters in heparin perturbation model
Section 5: Additional equations and parameters for the MEKi addition experiment (Table 5)
List of additional species and parameters in MEKi perturbation model
Section 6: Analyzing parameter uncertainty using Markov-chain Monte Carlo sampling
Parameter uncertainty for the seven sets identified using PSO was analyzed using Markov-chain Monte Carlo (MCMC) sampling approach as described previously . Hug et. al. showed that a multi-chain sampling approach can efficiently sample the posterior probability distribution for high-dimensional and non-linear models such as the FGF signaling model presented in this paper. Initializing from each of the seven parameter sets, Markov chains of length 100,000 were run using a modified version of 'adaptive metropolis’ algorithm as described in Haario et. al. . Additional file 2: Figure S2 shows the results of MCMC sampling for each of the 39 parameters with the initial positions marked as red stars.
For some parameters like 'H’ or 'Kd1a’ , the estimated optimal value for all seven optimal parameter sets are close together and located at the boundary of the allowed parameter range. Accordingly, the MCMC samples show parameter distribution close to the boundary. Importantly, the distribution shows that other parameter sets with values away from the boundary are able to explain the experimental data equally well. Thus, this analysis gives a more realistic impression of the objective function manifold as a function of parameter values. For some other parameters like 'kf44’ , two distinct clusters of MCMC samples were observed. This indicates that the parameter sets belong to two distinct basins of attraction. In summary, these results indicate that the seven sets of parameters identified by PSO are representative for the distribution of good fits that can be identified by MCMC sampling and provide the argument for limiting the further model analysis to these representative sets. The model was further analyzed using this large set of parameter values. For each of the parameter set represented in Additional file 2: Figure S2, the model satisfies all the qualitative constraints as described in Materials and methods section 3, including the constraint for biphasic pERK response. Thus, the model predicts biphasic pERK response for a wide range of parameter values represented in Additional file 2: Figure S2.
The authors acknowledge Dan Kirouac, Yasmin Hashambhoy-Ramsay and Andrijana Radivojevic for discussions on the results and feedback on the manuscript. The authors also thank Jeff Kearns and Sergio Iadevaia for support with using particle swarm optimization algorithm. Finally, the authors would like to thank Marco Muda for giving feedback on the manuscript.
- Sporn MB, Roberts AB: Peptide growth factors are multifunctional. Nature. 1988, 332: 217-219. 10.1038/332217a0.PubMedView ArticleGoogle Scholar
- Sporn MB, Roberts AB, Wakefield LM, de Crombrugghe B: Some recent advances in the chemistry and biology of transforming growth factor-beta. J Cell Biol. 1987, 105: 1039-1045. 10.1083/jcb.105.3.1039.PubMedView ArticleGoogle Scholar
- Williams EJ, Furness J, Walsh FS, Doherty P: Characterisation of the second messenger pathway underlying neurite outgrowth stimulated by FGF. Development. 1994, 120: 1685-1693.PubMedGoogle Scholar
- Garcia-Maya M, Anderson AA, Kendal CE, Kenny AV, Edwards-Ingram LC, Holladay A, Saffell JL: Ligand concentration is a driver of divergent signaling and pleiotropic cellular responses to FGF. J Cell Physiol. 2006, 206: 386-393. 10.1002/jcp.20483.PubMedView ArticleGoogle Scholar
- Delehedde M, Seve M, Sergeant N, Wartelle I, Lyon M, Rudland PS, Fernig DG: Fibroblast growth factor-2 stimulation of p42/44MAPK phosphorylation and IkappaB degradation is regulated by heparan sulfate/heparin in rat mammary fibroblasts. J Biol Chem. 2000, 275: 33905-33910. 10.1074/jbc.M005949200.PubMedView ArticleGoogle Scholar
- Tassi E, Al-Attar A, Aigner A, Swift MR, McDonnell K, Karavanov A, Wellstein A: Enhancement of fibroblast growth factor (FGF) activity by an FGF-binding protein. J Biol Chem. 2001, 276: 40247-40253.PubMedView ArticleGoogle Scholar
- Lander AD: Proteoglycans: master regulators of molecular encounter?. Matrix Biol. 1998, 17: 465-472. 10.1016/S0945-053X(98)90093-2.PubMedView ArticleGoogle Scholar
- Burgess WH, Maciag T: The heparin-binding (fibroblast) growth factor family of proteins. Annu Rev Biochem. 1989, 58: 575-606. 10.1146/annurev.bi.58.070189.003043.PubMedView ArticleGoogle Scholar
- Höök M, Kjellén L, Johansson S: Cell-surface glycosaminoglycans. Annu Rev Biochem. 1984, 53: 847-869. 10.1146/annurev.bi.53.070184.004215.PubMedView ArticleGoogle Scholar
- Dieci MV, Arnedos M, Andre F, Soria JC: Fibroblast growth factor receptor inhibitors as a cancer treatment: from a biologic rationale to medical perspectives. Cancer Discov. 2013, 3: 264-279. 10.1158/2159-8290.CD-12-0362.PubMedView ArticleGoogle Scholar
- Nielsen UB, Schoeberl B: Using computational modeling to drive the development of targeted therapeutics. IDrugs. 2005, 8: 822-826.PubMedGoogle Scholar
- Beenken A, Mohammadi M: The FGF family: biology, pathophysiology and therapy. Nat Rev Drug Discov. 2009, 8: 235-253. 10.1038/nrd2792.PubMedPubMed CentralView ArticleGoogle Scholar
- Anandappa SY, Winstanley JH, Leinster S, Green B, Rudland PS, Barraclough R: Comparative expression of fibroblast growth factor mRNAs in benign and malignant breast disease. Br J Cancer. 1994, 69: 772-776. 10.1038/bjc.1994.146.PubMedPubMed CentralView ArticleGoogle Scholar
- Mori S, Wu C-Y, Yamaji S, Saegusa J, Shi B, Ma Z, Kuwabara Y, Lam KS, Isseroff RR, Takada YK, Takada Y: Direct binding of integrin alphavbeta3 to FGF1 plays a role in FGF1 signaling. J Biol Chem. 2008, 283: 18066-18075. 10.1074/jbc.M801213200.PubMedPubMed CentralView ArticleGoogle Scholar
- Kolitz SE, Lauffenburger D a: Measurement and modeling of signaling at the single-cell level. Biochemistry. 2012, 51: 7433-7443. 10.1021/bi300846p.PubMedView ArticleGoogle Scholar
- Poste G: New insights into receptor regulation. J Appl Physiol. 1984, 57: 1297-1305.PubMedGoogle Scholar
- Schoeberl B, Eichler-Jonsson C, Gilles ED, Müller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002, 20: 370-375. 10.1038/nbt0402-370.PubMedView ArticleGoogle Scholar
- Himpe E, Kooijman R: Insulin-like growth factor-I receptor signal transduction and the Janus Kinase/Signal Transducer and Activator of Transcription (JAK-STAT) pathway. Biofactors. 2008, 35: 76-81.View ArticleGoogle Scholar
- Zhu H, Duchesne L, Rudland PS, Fernig DG: The heparan sulfate co-receptor and the concentration of fibroblast growth factor-2 independently elicit different signalling patterns from the fibroblast growth factor receptor. Cell Commun Signal. 2010, 8: 14-10.1186/1478-811X-8-14.PubMedPubMed CentralView ArticleGoogle Scholar
- Dupree M a, Pollack SR, Levine EM, Laurencin CT: Fibroblast growth factor 2 induced proliferation in osteoblasts and bone marrow stromal cells: a whole cell model. Biophys J. 2006, 91: 3097-3112. 10.1529/biophysj.106.087098.PubMedPubMed CentralView ArticleGoogle Scholar
- Fannon M, Forsten KE, Nugent M a: Potentiation and inhibition of bFGF binding by heparin: a model for regulation of cellular response. Biochemistry. 2000, 39: 1434-1445. 10.1021/bi991895z.PubMedView ArticleGoogle Scholar
- Filion RJ, Popel AS: A reaction–diffusion model of basic fibroblast growth factor interactions with cell surface receptors. Ann Biomed Eng. 2004, 32: 645-663.PubMedView ArticleGoogle Scholar
- Gopalakrishnan M, Forsten-Williams K, Täuber UC: Ligand-induced coupling versus receptor pre-association: cellular automaton simulations of FGF-2 binding. J Theor Biol. 2004, 227: 239-251. 10.1016/j.jtbi.2003.11.004.PubMedView ArticleGoogle Scholar
- Forsten KE, Fannon M, Nugent M a: Potential mechanisms for the regulation of growth factor binding by heparin. J Theor Biol. 2000, 205: 215-230. 10.1006/jtbi.2000.2064.PubMedView ArticleGoogle Scholar
- Yamada S, Taketomi T, Yoshimura A: Model analysis of difference between EGF pathway and FGF pathway. Biochem Biophys Res Commun. 2004, 314: 1113-1120. 10.1016/j.bbrc.2004.01.009.PubMedView ArticleGoogle Scholar
- Iadevaia S, Lu Y, Morales FC, Mills GB, Ram PT: Identification of optimal drug combinations targeting cellular networks: integrating phospho-proteomics and computational network analysis. Cancer Res. 2010, 70: 6704-6714. 10.1158/0008-5472.CAN-10-0460.PubMedPubMed CentralView ArticleGoogle Scholar
- Hug S, Raue A, Hasenauer J, Bachmann J, Klingmüller U, Timmer J, Theis FJ: High-dimensional Bayesian parameter estimation: case study for a model of JAK2/STAT5 signaling. Math Biosci. 2013, 246: 293-304. 10.1016/j.mbs.2013.04.002.PubMedView ArticleGoogle Scholar
- Marek L, Ware KE, Fritzsche A, Hercule P, Helton WR, Smith JE, McDermott LA, Coldren CD, Nemenoff RA, Merrick DT, Helfrich BA, Bunn PA, Heasley LE: Fibroblast growth factor (FGF) and FGF receptor-mediated autocrine signaling in non-small-cell lung cancer cells. Mol Pharmacol. 2009, 75: 196-207. 10.1124/mol.108.049544.PubMedPubMed CentralView ArticleGoogle Scholar
- Schlessinger J, Plotnikov AN, Ibrahimi OA, Eliseenkova AV, Yeh BK, Yayon A, Linhardt RJ, Mohammadi M: Crystal structure of a ternary FGF-FGFR-heparin complex reveals a dual role for heparin in FGFR binding and dimerization. Mol Cell. 2000, 6: 743-750. 10.1016/S1097-2765(00)00073-3.PubMedView ArticleGoogle Scholar
- Pellegrini L, Burke DF, von Delft F, Mulloy B, Blundell TL: Crystal structure of fibroblast growth factor receptor ectodomain bound to ligand and heparin. Nature. 2000, 407: 1029-1034. 10.1038/35039551.PubMedView ArticleGoogle Scholar
- Nugent MA, Edelman ER: Kinetics of basic fibroblast growth factor binding to its receptor and heparan sulfate proteoglycan: a mechanism for cooperactivity. Biochemistry. 1992, 31: 8876-8883. 10.1021/bi00152a026.PubMedView ArticleGoogle Scholar
- Ibrahimi OA, Zhang F, Hrstka SCL, Mohammadi M, Linhardt RJ: Kinetic model for FGF, FGFR, and proteoglycan signal transduction complex assembly. Biochemistry. 2004, 43: 4724-4730. 10.1021/bi0352320.PubMedView ArticleGoogle Scholar
- Wu Y, Chen Z, Ullrich A: EGFR and FGFR signaling through FRS2 is subject to negative feedback control by ERK1/2. Biol Chem. 2003, 384: 1215-1226.PubMedView ArticleGoogle Scholar
- Forsten-Williams K, Chua CC, Nugent MA: The kinetics of FGF-2 binding to heparan sulfate proteoglycans and MAP kinase signaling. J Theor Biol. 2005, 233: 483-499. 10.1016/j.jtbi.2004.10.020.PubMedView ArticleGoogle Scholar
- Eberhart R, Kennedy J: A new optimizer using particle swarm theory. MHS’95 Proc Sixth Int Symp Micro Mach Hum Sci. 1995, 39-43.View ArticleGoogle Scholar
- Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. 2009, 25: 1923-1929. 10.1093/bioinformatics/btp358.PubMedView ArticleGoogle Scholar
- Brown A, Robinson CJ, Gallagher JT, Blundell TL: Cooperative heparin-mediated oligomerization of fibroblast growth factor-1 (FGF1) precedes recruitment of FGFR2 to ternary complexes. Biophys J. 2013, 104: 1720-1730. 10.1016/j.bpj.2013.02.051.PubMedPubMed CentralView ArticleGoogle Scholar
- Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S, Wilson CJ, Lehár J, Kryukov GV, Sonkin D, Reddy A, Liu M, Murray L, Berger MF, Monahan JE, Morais P, Meltzer J, Korejwa A, Jané-Valbuena J, Mapa FA, Thibault J, Bric-Furlong E, Raman P, Shipway A, Engels IH, Cheng J, Yu GK, Yu J, Aspesi P, De Silva M,et al: The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature. 2012, 483: 603-607. 10.1038/nature11003.PubMedPubMed CentralView ArticleGoogle Scholar
- Naimy H, Buczek-Thomas JA, Nugent M, Leymarie N, Zaia J: Highly sulfated nonreducing end-derived heparan sulfate domains bind fibroblast growth factor-2 with high affinity and are enriched in biologically active fractions. J Biol Chem. 2011, 286: 19311-19319. 10.1074/jbc.M110.204693.PubMedPubMed CentralView ArticleGoogle Scholar
- Rapraeger A, Krufka A, Olwin BB: Requirement of heparan sulfate for bFGF-mediated fibroblast growth and myoblast differentiation. Science. 1991, 252 ((80-)): 1705-1708.PubMedView ArticleGoogle Scholar
- Yayon A, Klagsbrun M, Esko JD, Leder P, Ornitz DM: Cell surface, heparin-like molecules are required for binding of basic fibroblast growth factor to its high affinity receptor. Cell. 1991, 64: 841-848. 10.1016/0092-8674(91)90512-W.PubMedView ArticleGoogle Scholar
- Aviezer D, Yayon A: Heparin-dependent binding and autophosphorylation of epidermal growth factor (EGF) receptor by heparin-binding EGF-like growth factor but not by EGF. Proc Natl Acad Sci U S A. 1994, 91: 12173-12177. 10.1073/pnas.91.25.12173.PubMedPubMed CentralView ArticleGoogle Scholar
- Tessler S, Rockwell P, Hicklin D, Cohen T, Levi BZ, Witte L, Lemischka IR, Neufeld G: Heparin modulates the interaction of VEGF165 with soluble and cell associated flk-1 receptors. J Biol Chem. 1994, 269: 12456-12461.PubMedGoogle Scholar
- Zioncheck TF, Richardson L, Liu J, Chang L, King KL, Bennett GL, Fügedi P, Chamow SM, Schwall RH, Stack RJ: Sulfated oligosaccharides promote hepatocyte growth factor association and govern its mitogenic activity. J Biol Chem. 1995, 270: 16871-16878. 10.1074/jbc.270.28.16871.PubMedView ArticleGoogle Scholar
- Reichsman F, Smith L, Cumberledge S: Glycosaminoglycans can modulate extracellular localization of the wingless protein and promote signal transduction. J Cell Biol. 1996, 135: 819-827. 10.1083/jcb.135.3.819.PubMedView ArticleGoogle Scholar
- Ruppert R, Hoffmann E, Sebald W: Human bone morphogenetic protein 2 contains a heparin-binding site which modifies its biological activity. Eur J Biochem. 1996, 237: 295-302. 10.1111/j.1432-1033.1996.0295n.x.PubMedView ArticleGoogle Scholar
- Takada T, Katagiri T, Ifuku M, Morimura N, Kobayashi M, Hasegawa K, Ogamo A, Kamijo R: Sulfated polysaccharides enhance the biological activities of bone morphogenetic proteins. J Biol Chem. 2003, 278: 43229-43235. 10.1074/jbc.M300937200.PubMedView ArticleGoogle Scholar
- Kuo W, Digman MA, Lander AD: Heparan sulfate acts as a bone morphogenetic protein coreceptor by facilitating ligand-induced receptor hetero-oligomerization. Mol Biol Cell. 2010, 21: 4028-4041. 10.1091/mbc.E10-04-0348.PubMedPubMed CentralView ArticleGoogle Scholar
- Fon Tacer K, Bookout AL, Ding X, Kurosu H, John GB, Wang L, Goetz R, Mohammadi M, Kuro-o M, Mangelsdorf DJ, Kliewer SA: Research resource: Comprehensive expression atlas of the fibroblast growth factor system in adult mouse. Mol Endocrinol. 2010, 24: 2050-2064. 10.1210/me.2010-0142.PubMedPubMed CentralView ArticleGoogle Scholar
- Haario H, Saksman E, Tamminen J: An Adaptive Metropolis Algorithm. Bernoulli. 2001, 7: 223-242. 10.2307/3318737.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.