Introduction

Today we will begin our journey through Systems Biology and the careful consideration of experimental and computational design that is required for large-scale biological engineering studies. At the end of lab today, you will choose a chemical inhibitor that targets the Epidermal Growth Factor Receptor (EGFR) signaling network with the goal of overcoming a growing common problem in the treatment of disease -- drug resistance.

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 an 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.

Part 1: 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. After your brain is full of EGFR signaling, you'll choose a combination of two inhibitors in an attempt to curb signaling and kill SKOV3 cells, a human model of ovarian cancer.

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.

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 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 your 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 2: 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 2 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, 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 1).

Locate the nuclear STAT3 dimer, STAT3n-STAT3n. Translocation of the STAT3 dimer to the nucleus also drives proliferation.

Finally, locate pAkt. The phosphorylated from of Akt promotes cell survival by suppressing apoptosis, 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, it is time to run a simulation to see it in action!

Simulate the model using COPASI

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 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 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 Menu --> 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 number' 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 -- 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 4: 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 4 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

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 5: 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 5 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 (See Figure 5). 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.

Generate the same plots that you made in the last section and add them to your PowerPoint file. Attach this file to your EN notebook for easy reference.

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)?

Form a concise hypothesis about how the viability or proliferation of the simulated cancer cell might change. State your evidence. Write this hypothesis in bold type in your EN notebook.

While it is unknown in our SKOV3 cells if the Cbl-EGFR interaction is altered, it is known that SKOV3 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.

In our next lab class we will use PCR and DNA sequencing to determine if the EGFR protein expressed in SKOV3 cells contains a common EGFR mutation, but today you will use you newly gather 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 STAT 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?

As you discovered in Module 1, experiments are more powerful when performed in replicates. Therefore, a minimum of two groups will team up to perform biological replicates of each inhibition condition. Please sign up on the talk page to choose your inhibitor.

Part 2: Mock drug combination exercise

Note: We have a total of 6 multichannel pipettes to complete this exercise, therefore groups will need to share pipettes or stagger their start time.

Although we won't be preparing our high-throughput experiment until M2D6, today we will practice designing a 96-well plate to facilitate efficient mixing and transfer of our chemical inhibitors. On M2D6-7 we will perform a dose-response study using two inhibitors -- an EGFR inhibitor and an inhibitor of your choice based upon the simulations you performed today -- to quantify cell viability. Here are some important details to consider:

We will test 5 inhibitor concentrations + 1 control (no inhibitor).

Inhibitor combinations are of interest to determine any synergies.

Your experiment will be completed in duplicate on your plate and in at least duplicate across the class.

The total volume of media + inhibitor will be 100 μL per well.

One way set up this experiment is to duplicate a 6x6 grid in a 96-well plate -- this allows for all combinations of inhibitor to be tested. To facilitate this plan, we will first set up a dilution plate to prepare 10x serial dilutions of inhibitor that can easily be transferred to our cells using a multi-channel pipetteman. There are several reasons to set up the experiment this way:

Finally, let's prepare our experiment plate - a mock-up of what your high throughput cell viability experiment will look like on M2D6-7.

Using the multi-channel pipette, add 50 μL from column 12 from your dilution plate to columns 6 and 12 of a new 96-well plate (your experiment plate).

Note: We deliberately started at the lowest dye/drug concentration (zero!) because that will allow you to re-use the same pipette tips for the next steps.

Add 50 μL from column 11 of the dilution plate to columns 5 and 11 of the experiment plate.

Add 50 μL from column 10 of the dilution plate to columns 4 and 10 of the experiment plate.

Repeat these steps until all of the columns contain red dye.

Add 6 new pipette tips to your multi-channel pipetteman.

Now add the blue dye -- begin by pipetting the blue contents (50 μL) of row F from the dilution plate into row F of the experiment plate.

Note: You will need to perform two pipetting steps to fill all of row F (columns 1-6 and then 7-12).

Continue by adding 50 μL from row E of the dilution plate to row E of the experiment plate.

Fill the remaining rows with increasing blue dye.

Finally, check the contents of your experiment plate with the experiment plate on the front bench. If they don't match, simply dump the contents of your plate in the lab sink and repeat filling the experiment plate until they do.

The next time you perform these dilution and pipetting steps, they will be completed inside the TC hood, so take a moment to reflect on how you will minimize the contamination risk while setting up your M2D6 high throughput experiment.

For Next Time

Complete and submit your Module 1 Abstract and Data Summary assignment by 11am on the day of your next lab class. Submit this assignment to Stellar. Review the submission guidelines here.

The first Reflection of the semester is due on October 14th at 5pm. Read the prompt for Reflection #1 and write it as a separate note in your EN notebook. Email that note to Shannon (skalford at mit dot edu) or Agi (astachow at mit dot edu) by 5pm on the due date.

One of two major communication assignments in Module 2 is the presentation of a Journal Club. The two JC dates are M2D5 and M2D8. Please check the Schedule page to determine which date fits better with your busy schedule. After M2D2 you will have the opportunity to sign up for a day and a paper covering either High Throughput Screening technologies or Cancer Systems Biology.