Wolfram Computation Meets Knowledge

Modeling in the Search for New Drugs—G-Protein-Coupled Receptors (GPCRs)

Explore the contents of this article with a free Wolfram SystemModeler trial. Yesterday, the Nobel Prize in Chemistry was awarded to Robert J. Lefkowitz and Brian K. Kobilka for having mapped how a family of cellular receptors called G-protein-coupled receptors (GPCRs) work. The Nobel Prize winners’ research has proven to be very important in the development of novel therapeutic drugs—about 40–50% of all therapeutic drugs in use today are centered on GPCRs. The real beauty of GPCR-based response systems is that they include components that are used over and over again for the response to external signals in many kinds of cellular functions throughout our bodies. Sight, smell, and the adrenaline response are examples of these GPCR-mediated responses with physiologically important functions.

Identifying new targets for therapeutic drug intervention includes analysis of the complex webs of signaling pathways and feedback systems in our cells, extending beyond the first event of a signal connecting with the GPCR on the cell surface, which is non-trivial. Lately the cost-effective practice of using mathematical models as an initial step for finding those elusive new targets, and also as a tool for understanding how other reactions of a cell might be affected by a new drug, has been growing. In this blog post we are going to use modeling and simulation in order to illustrate how the GPCR-based cellular response to an external signal can be modified. And by performing this analysis, I thought we should also see how we can find promising targets for therapeutic drug design, which are then aimed at either increasing or decreasing the response. Since the first two steps in the pathways are identical in most of the GPCR-based signal responses in a cell, we can freely choose a representative model. One such well understood signal response pathway that uses GPCR is the mating pheromone response in yeast, which we are here going to explore using Mathematica and Wolfram SystemModeler.

Wolfram SystemModeler and Mathematica

The yeast-mating pheromone model

The model we’re going to use has been developed by Tau-Mu Yi et al. It can be retrieved from the BioModels Database, and by using the SBML import function and the BioChem library, we can quickly get the model into Wolfram SystemModeler. The imported model is shown in the figure below, where concentrations of a molecule are represented by the green circles and specific chemical reactions are shown as red arrows (indicating direction).

Model built in Wolfram SystemModeler

As can be seen in the figure above, the model describes two sub-cycles: the activation and deactivation of the G-protein-coupled receptor (GPCR) and the G-protein. The receptor’s activation is triggered when a signal substance (also known as an agonist) binds to the receptor. In general, the receptor consists of seven interconnected “rods” inserted into the membrane of the cell, and when the signal reaches the receptor, the arrangement of the rods is changed, which transmits the signal to the inside of the cell, as illustrated in the figure below. On the inside of the cell, the now activated receptor mediates the signal further by activating the G-protein (as seen on the right in the model), which causes the G-protein to divide into two subunits: Gα and Gβγ. Typically the Gα subunit then transmits the signal further into the cell, ultimately generating the desired physiological or cellular response.

Changing the arrangement of the rods

A first glance at the model

To familiarize ourselves with the model, we will start by simulating it using the input, initial conditions, and parameter values as they are presented in the original article. These values correspond to how a normal yeast cell responds to the external signal of a mating pheromone.

The simulations can be run in either SystemModeler or Mathematica, but here we are going to use the Mathematica environment. We will start by loading WSMLink:


And by using the WSMSimulate function, we simulate the model for 600 s.

Simulating the model for 600 s

As a measure of the response in the cell, we are going to calculate the fraction of activated G-proteins. Since this fraction will be computed multiple times, it is handy to define it as a small function, were we focus on the α subunit of the activated G-protein. Using the same notation as in the model, this fraction is defined as follows:

Fraction of activated G-proteins

With our response function sorted out, let’s see what the cellular response looks like.

Calculating the cellular response

Plotting the cellular response

Plot of the cellular response

Here we find that there is a rapid response, which is slowly being diminished; the exact timing and absolute height of the peak response is equal to:

Finding the exact timing and absolute height of the peak response

Exact timing and absolute height of the peak response

So, typically 45% of all the G-proteins are in an activated state 23 s after the signal is initially received.

Modifying the cellular response

Now, let’s start analyzing what happens to the cellular response if we alter specific functions in the cell. In the paper by Yi et al., mutated yeast cells with reduced ability to end the response were experimentally studied; this mutation can be implemented in the model by changing the value of the parameter reaction_6.k1 from 0.11 s-1 to 0.004 (s-1).

Parameter values in a model can be modified by using the WSMParameterValues option to WSMSimulate:

Modifying the cellular response

Modified cellular response

Plotting the modified cellular response

Plot of the modified cellular response

And once again we characterize the response quantitatively.

Finding the exact timing and absolute height of the peak response

Exact timing and absolute height of the peak response

From these results we can conclude that a cell without the proper ability to regulate its response has a much stronger response compared to the normal cell; the peak response is almost doubled. Also, there is a prolonged response that is very slowly being reduced.

How does the concentration of signal molecules affect the cell?

A standard analysis method in cell biological systems is the dose-response, that is, characterization of the response for different concentrations of a signal molecule. Let’s exemplify this method by evaluating how the yeast mutant responds to different concentrations of the signal molecules.

In this particular model the concentrations of the signal molecules are defined as molecules per cell (in contrast to moles per liter, M), and therefore the signal concentration has to be defined in numbers of molecules. So in order to have doses corresponding to an experimental setup, (e.g. a dose of 1μM), we define our input sweep, using Avogadro’s constant (which is used in the definition of the unit mole), as:

Defining the input sweep using Avogadro's constant

We now run the different doses for the mutant yeast cell, by including the WSMInitialValues option.

Running the different doses for the mutant yeast cell

Since we want to see the response as a function of the signal concentration, we extract the values at the peak response, which we established at 68 s for the mutant.

Extracting the values at the peak response

Defining the activeGDoseMutant function

We are now ready to explore the cellular response as a function of the signal strength. Typically there is a need for a dose that yields a response slightly higher than when the first signs of response are visible in order to ultimately induce any physiological effect. Furthermore, in the case of therapeutic drugs, at higher doses undesired side effects appear and grow stronger as the dose increases. Another reason for not administering a large dose is that after a certain dose, the response does not increase, which in the case of expensive drugs is a waste of money.

In the interactive diagram below you can play around with the cell response and find the optimal dose with respect to the safety margins and desired effects. Please note that the limits are just for illustration and are arbitrarily chosen.

Cell response

This is an example of how delicate the dosage of a drug can be; in our model there are basically yeast cells floating in a jar, but imagine that you need to get this just right within the human body. The dose and response can of course change in diseased cases, and it can also be the case that we want to use a drug in order to reduce an unwanted response to some signal, which we are going to focus the remainder of this blog post on.

Identifying potential drug targets

So far we have been focused on how a particular change in the cell affects the cellular response. Let’s now move on to an efficient method of analyzing a larger portion of the system simultaneously.

In order to quickly assess in what way a system reacts to changes in the reaction speed, we will use WSMSimulateSensitivity, which locally analyzes the sensitivity of the states in system (in this case, the various concentrations) due to shifts in the parameters.

Let’s test how the normal yeast cell responds to changes in the following reactions:

Analyzing a larger portion of the system simultaneously

Normal yeast cell responding to changes

Before we can fully evaluate the sensitivities of response, we need to specify the derivative of the response function with respect to the parameter that is being analyzed:

D[cellResponse[x[p], y[p], z[p]], p]

Derivative of the response function with respect to the parameter that is being analyzed

x, y, and z correspond to gαGTP, gαGDP, and gProt respectively; the renaming is just for simplification. We are now ready to define the derivative of the cell response as a new function:

Defining the derivative of the cell response as a new function

We are almost done. The final step is to get out the appropriate data from the simulations so that we might asses the sensitivities. These lists correspond to state values from the initial parameter simulations, state value derivations with respect to the parameter value, and the initial parameter values.

Getting out the appropriate data from the simulations

With everything gathered, we can now see how the cellular response is affected by a 50% change in the parameter values around the local linearization. Note that should we only want to look at a single state variable, instead of our response function, we can use WSMPlot directly on Msens, with just a few additional inputs (see the documentation for more information).

Building multiple plots

Multiple plots of data

From this analysis we find two promising reactions that we can modify in order to change the cellular response: G-protein activation and the receptor and signal degradation. The red lines (at the edge of the red area) correspond to an increase in parameter value. Altering these two parameters yields distinctly different cellular responses. For instance, an increase in how fast the G-protein is activated yields a rapid and stronger response, whereas an increase in how fast the receptor and signal degradation occurs only affects the “long-term” dynamics by lowering the response.


Once a part in the cell is identified for use as a drug target, it is wise to evaluate it more in-depth and find out how the rest of the system is affected if multiple reaction rates are lowered or increased. Since there will be many possible combinations, I’m going to reduce the complexity a little in this example. We are going use our previous results when choosing the three most important reactions and include one we so far have ignored, which is how fast the signal molecule interacts with the receptor. As another limit, only three different values for each reaction will be used, yielding a total of 81 combinations.

Running more simulations

And we calculate the response as usual.

Calculating the response

Now we can start exploring. In the upper panel of the figure is a diagram of the model, with popup menus superimposed. By changing the values, you can see how the response (in the lower panel) changes as a consequence of the different combinations of reaction speeds. I also added the maximal response in red; can you identify the solution that yields the strongest results?

Response changes

This is a rather simple example, and we only looked at a select few components in the signal response pathway, but it would be near impossible to figure out all the interplay without using mathematical modeling and simulation. As you can see in this example, if we want to get as small a response as possible, there are multiple targets to aim our research toward.

I mentioned earlier that most therapeutic drugs target GPCR pathways, and most of these work by competing with the signal for binding to the receptor (this effect can be seen if you change the setting in the top-left corner in the model above to “low”). Some specific examples of these drugs are caffeine (to reduce sleepiness), β-blockers (in heart medicine), and antihistamines (to stop allergic reactions), which all bind to the receptor so that the cell’s responses to the human body’s internal signals are dampened or completely hindered.

Some concluding thoughts

We have now done some in silica experiments on how one might alter the cellular response in signaling pathways that use the G-protein coupled receptor, for which the in-depth characterization was awarded this year’s Nobel Prize in Chemistry. And by performing these in silica experiments, we now have glimpsed how easy it is to evaluate many different modifications of a cellular signaling response compared to the tedious work needed to carry out such experiments in a laboratory. Thanks to the many years of diligent biochemical and cell biological research, there exist many building blocks for the design of mathematical models of signaling pathways. By merging the knowledge in SystemModeler and performing analyses in Mathematica, there is vast potential for identifying and characterizing novel drug targets, as well as defining appropriate doses, efficiently and cost-effectively.

Feel free to download this notebook and model and play around for yourself. I didn’t include any downstream events, since the intent was to focus on G-proteins and the GPCR cycles. However for the more interested reader, there is a model developed by Bente Kofahl and Edda Klipp that has an extensive downstream expansion. This model can rather easily be implemented in SystemModeler using the standard components in the BioChem library.

Download this ZIP file that includes the post as a CDF file and the model.


Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.


  1. I would like to see a physiological based pharmacokinetic model for humans built up in modeler.