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.

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

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.

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.

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:

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

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:

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:

And once again we characterize the response quantitatively.

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:

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

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.

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.

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:

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:

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:

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.

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 M_sens, with just a few additional inputs (see the documentation for more information).

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.

Interactions

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.

And we calculate the response as usual.

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?

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.