The first mathematical model published in the journal Cell described the binding of EGF to EGFR. That paper ushered in an era of multidisciplinary research which brought together scientists and engineers to build mathematical models of intracellular signaling pathways. In fact, the field of Systems Biology, often thought of as the use of computational models to study multiple signaling pathways across variable time and length scales, grew directly from these first collaborations.
There are many approaches to computational modeling, each aiming to understand how information travels through a signaling network to ultimately produce a cell phenotype. We will discuss some of the various modeling techniques in lecture. In lab you will design and perform a high-throughput experiment that alters cell proliferation, a cell phenotype that is often dysregulated in cancer. There are several distinct signaling pathways downstream of EGFR activation that influence cell proliferation and viability. The best understood pathways are the Ras/Erk pathway, PI3K/Akt pathway, and the STAT transcriptional regulation pathway. One important tool that we can take advantage of as biological engineers is the wealth of mathematical models describing the dynamics of these pathways upon EGF stimulation. It is good practice to gather as much data as possible before designing an experiment -- therefore, we will use computational simulation to better understand how chemical inhibition will alter information flow through the Ras/Erk, PI3K/Akt, or STAT pathways.
In 2012, Bidkhori et al. employed an ordinary differential equation (ODE)-based mathematical model to study the effect of a common mutant EGFR protein on the pathways that control cell proliferation/viability. Today we will utilize two web-based programs, CellDesigner and COPASI, that facilitate building and visualizing cell signaling models to interogate the Bidkhori model. In fact, the computational study of cell signaling has become such a useful tool that a standardized programing language SBML has been developed. Concomitant with the development of SBML, a number of user friendly GUIs have popped up to make simulation of models accessible to everyone -- including the two that we are using today.
Simulation of the EGFR network
First, you will find the Bidkhori model on the BioModels database. Next, you will visualize the network structure using CellDesigner. Finally, you'll simulate the network to get a feel for the information flux through the Ras/Erk, PI3K/Akt, and STAT3 pathways using COPASI.
Import the SBML code
The BioModels database is a very large compendium of biochemical models that are for public use. The curators of the database accept SBML code from model builders and then validate that the original figures published with the model can be recapitulated using the submitted code. The Systems Biology community has tried very hard to make their computational models available to everyone, and you'll find that many engineers will start with a published model and then build from there. This is not only an acceptable practice, but an encouraged one -- with two important rules: 1. You must understand the output and design of the model that you start with and 2. The authors of the original model are always cited in your work.
- Go to the BioModels database.
- In the search box enter 'EGFR'.
- You should find a list of 24 'curated' models.
- A curated model is one that has been validated by the people who maintain the BioModels website.
- Scroll all the way down to the bottom of that list and choose BIOMD0000000452 Bidkhori2012 - normal EGFR signaling
- Note: take care not to scroll through the models that have not been curated (or validated).
- Click the disk icon to save the models SBML code to your computer.
- Now click on the BioModel database number associated with the model.
- In this window you will find information about the paper the model was described in, other models that the current one is derived from, and its creation date.
- Notice how computational modeling of biological systems is happening all around the world!
- Navigate to the Math tab.
- Here you will find all the reactions that are described by this mathematical model.
- For example: [EGF] + [EGFR] <--> [EGF-EGFR] is the equation that describes the binding of Epidermal Growth Factor to its receptor Epidermal Growth Factor Receptor. This reaction is the start of the signaling cascade!
- Navigate to the Curation tab.
- Here is the information provided by the BioModels curators as to how the model was validated. You'll see that the author of this tab also used COPASI to run the model.
- If you maximize the image of the plots, you can get an idea for what type of output you can expect from your simulation.
- Double check that you've downloaded the SBML file to your computer. If you have, you're good to move on to the next section.
Open the model in CellDesigner
Figure 3: The start of the EGFR signaling cascade as coded in the Bidkhori et al. model and visualized using CellDesigner.
- Open CellDesigner from your Applications menu.
- Use the Open command under File to open the SBML code that you downloaded form the BioModels database.
- Next, clean up the CellDesigner workspace a bit to make things easier to view:
- Maximize the CellDesigner window.
- Under View:
- unselect Show Reaction Id
- under Change Toolbar Visible, select Hide all
- deselect the Tree view
- Click the Proteins tab near the bottom of the workspace. How many proteins and protein-complexes are described by this model?
- Click Reactions. How many different ODEs are used to simulate this biochemical system?
- Now, minimize the bottom panel using the small, gray arrows in the upper left hand corner.
- The components of the EGFR network and their complexed states (for example EGF + EGFR --> EGF-EGFR) are illustrated in the large window in the middle of the screen. The number of model components makes it very difficult to view the structure of the model. To see a quick view of the system, under View choose Zoom Fit.
- The 'hairball' view of the model includes arrows (edges) and boxes (nodes) that represent the different biochemical reactions of the system. To better organize the model for ease of visualizing the information flow, go to Layout and choose Hierarchic Layout.
- Feel free to use the default settings.
- With this view it is a bit easier to picture the information flow through the network.
- Use the zoom-in button or View --> zoom-in to expand the upper right hand corner of the model.
- You should see something similar to Figure 3 on this wiki page.
- This diagram represents the first events in the EGFR signaling cascade. Specifically, the binding of EGF to EGFR and subsequent EGFR dimerization (as represented by EGFR2). The long, nonsense looking terms are default names for rate parameters given in the model. You don't have to worry about these, but unfortunately, sometimes we can't get rid of them using this program.
- Zoom around the model diagram to become familiar with what proteins (both phosphorylated/active -- denoted by pXXX and total -- denoted by XXX) are involved and when.
- Locate ppERK. This is the terminal component of the Ras/ERK pathway in this model. Doubly phosphorylated ERK serves two purposes: 1. It translocates to the nucleus to promote gene transcription and 2. It promotes a positive feedback loop driving amplification of the signaling pathway (see Figure 2).
- Locate the nuclear STAT3 phosphorylated dimer, pSTAT3n-pSTAT3n. Translocation of the STAT3 dimer to the nucleus also drives proliferation.
- Finally, locate pAkt. The phosphorylated from of Akt promotes cell survival by suppressing cell death, but does not necessarily promote cell proliferation.
- Use the handout that lists all of the equations that are included in the model as a guide.
Now that you are familiar with the topology of the model (or the order in which the proteins communicate with one another), it is time to run a simulation to see it in action!
Simulate the model using COPASI
Figure 4: View of COPASI start-up page with Bidkhori et al. model loaded into program.
Although CellDesigner includes code to perform simulations, there are elements of the software that are not as user-friendly as some other available programs. In fact, this type of deficiency is currently a challenge when working with computational models -- many software packages exist that do a few things very well, but it is difficult to locate a program that can do everything you need. Perhaps if computational biology is your area of interest, you can fix this problem for future biological engineers!
To simulate and interrogate the Bidkhori et al. model, we will utilize COPASI. This is another freeware that can be downloaded and used on lab or personal computers. You may find that too many windows open is distracting; if that is a problem, feel free to close CellDesigner.
Note: We are using COPASI v4.8 (Build 35) as it seems to be the most stable on our lab computers.
- Open COPASI.
- Under the File menu choose Import SBML... and open the SBML file that you downloaded from the BioModels database.
- On the COPASI GUI, change 'Quantity Unit' to μmol. The initial values for the model are all in units of μM. The time simulation occurs in seconds.
- Now that the model components are loaded into the program, you should explore how COPASI displays the model. There is no network map here, but you can look at the mathematical equations by opening the left-hand sidebar menu Model --> Mathematical --> Differential Equations. Here is another, perhaps more familiar view, of the ODEs underlying the model.
- Ensure that the local parameters are set to 'display numerical value' to get an idea for the order of magnitude of the rate constants used in this model.
- Alternatively, you may also browse the equations, parameters, and initial values for the reaction components under the Biochemical menu.
- Note that the initial concentration of EGF (look under the Species tab) -- you can think of it as adding EGF to a cell culture without any growth factor -- is set at 0.0081967 μM = 8.2 nM ~ 50 ng/mL EGF, a very common stimulating concentration in the scientific literature.
Before we run our simulation, we will prepare some output graphs to facilitate easy comparison once we start exploring the biology and parameter space.
Figure 5: Example of how to create a new plot in COPASI.
- Under Output Specifications menu, highlight Plots (0).
- Create a new plot by double clicking the 'New Plot' line (see Figure 5 as a guide).
- Title this plot 'pEGF-EGFR2'
- Click on New curve...
- For the x-axis, choose Time --> Model Time
- For the y-axis, choose Species --> Transient Concentrations --> [pEGF-EGFR2](t)
- Click OK
- Leave the remaining default settings
- Click Commit
- Next, create three more plots:
- Plot 2: Ras/ERK pathway
- Include pShc, Ras-GTP, and ppERK versus time in this plot. Hint: you can click on New curve... to add lines to the same plot.
- Plot 3: PI3K/Akt pathway
- Include pPI3K, pAkt versus time in this plot.
- Plot 4:STAT pathway
- Include pSTAT3c-pSTAT3c and pSTAT3n-pSTAT3n versus time in this plot.
- Once you have added all four plots, hit the Commit button at the bottom of the page in order to continue.
- Next, under Tasks --> Time Course change the duration of the simulation to 1000s and the interval size to 1s.
- Confirm that the Integration Interval and the Output Interval read 0 to 1000.
- Click on Run.
Figure 6: Simulation of the EGFR signaling pathway in response to 8.2 nM EGF in normal cells. This is a screenshot of Shannon's Powerpoint file
You should now have four plot windows open on the desktop. Take a minute to first save each image and paste them into a Powerpoint file. Use Figure 6 as a reference layout.
Use COPASI to simulate EGFR signaling in Cancer
One way that tumor cells can escape therapeutic targeting is by modulating their response to growth factor. This can be done in a variety of ways; signaling through a mutant form of the receptor, 're-wiring' the signaling cascade to rely on another pathway, or altering the way that the receptor is presented on the cell surface.
Today, we will study what happens when there are defects in how often EGFR is turned over from the surface of the cells. In normal cells, the binding of EGF to EGFR promotes receptor internalization and eventually degradation. This is one way that cells can control the level of response to EGF -- by causing EGFR to disappear from the cell surface, the cell can keep the intracellular signaling pathways in check. However, in some cancer cells, the level of EGFR on the surface is altered and the time taken for the receptor to internalize is significantly increased. Can you predict what effect this will have on the levels of phosphorylated and active signaling proteins within the cell?
Cbl is a protein that interacts with EGFR and modulates its internalization. Bidkhori et al. optimized the kinetic rate parameters describing Cbl interaction with EGFR and EGFR-containing protein complexes in order to 'fit' data in the literature describing the internalization of EGFR in some non-small cell lung cancer cell lines. In other words, the parameters were adjusted in a systematic manner until the model output matched the experimental data. Using the handout provided at the beginning of class (and linked here), change the model parameters to see what happens to our favorite Ras/Erk, PI3K/Akt, and STAT pathways when the level of the receptor on the cell surface is increased due to changes in Cbl-EGFR interaction.
- Go to 'Parameter Overview' under the 'Biochemical' menu.
- Scroll down to find the equations of interest.
- Unfortunately, those crazy long strings of letters and numbers appear again -- but trust us that the top parameter is the forward rate constant and the bottom one is the reverse rate constant.
- Change the values accordingly.
Generate the same plots that you made in the last section and add them to your PowerPoint file. You can print that file to take with you.
Finally, if you have a little extra time, find Equation r123. Set K123 = 0. This action simulates the loss of PTEN, and important brake on the PI3K/Akt pathway -- this is interesting to us because many cells have a very low level of PTEN, and deregulation of PTEN is one mechanism of acquired resistance. Run the simulation and compare the four graphs you obtained to the normal cells and the cancer cells with increased EGFR levels.
Choose your inhibitors for a high-throughput drug screen
Carefully examine the plots that you generated using the normal and cancer EGFR signaling model in COPASI:
- How have the dynamics of the Ras/Erk, PI3K/Akt, and STAT pathways changed?
- Is the steady state behavior of the pathways different? (To adequately answer this question, you might try increasing the time duration of the simulation.)
- Which pathway appears to be most affected by a change in Cbl-EGFR interaction (i.e. a change in EGFR trafficking)?
While it is unknown if the Cbl-EGFR interaction is altered in all NSLC cells, it is known that many cells have an increased level of the EGFR protein. At this point, we hope that you have proven to yourself that a change in EGFR levels has a large effect on intracellular signaling. To make matters even more complex, some tumor cells contain a mutant form of the EGFR protein -- these activating mutations turn ON the signaling cascade in the absence of EGF and alter the response of the pathway to drugs meant to shut it down.
Today you will use your newly gathered knowledge of the signaling pathway to choose one chemical inhibitor that mimics a clinically used (or close to the clinic!) therapeutic agent for cutting off the signaling through one arm of the EGFR pathway that we've studied today.
Your choices are:
- U0126 MEK Inhibitor
- LY924002 PI3K Inhibitor
- Stattic STAT3 Inhibitor
Consider what pathway you find most interesting -- at this point we have no data using our cell line of interest (SKOV3) regarding cell viability in the presence of dual EGFR and MEK, PI3K, or STAT inhibition. Which inhibitor do you predict will have the greatest impact on cell viability? What leads you to believe that?