PathwayLab Mathematica Help Browser Contents

Introduction

PathwayLab is a software tool for modeling, analysis, and information management of biochemical pathways. The tool streamlines the pathway building process by its rich and flexible set of graphical building blocks for specifying biochemical entities, reaction, and control mechanisms. The pathway models are built using drag and drop from stencils into a drawing page. This makes it very easy to rapidly build models of biochemical reaction networks such as metabolic pathways, signaling and gene regulatory networks.

A pathway model defined in PathwayLab can be exported to a Mathematica m-file as a pathway model object. This file can be loaded into Mathematica and the PathwayLab application package provides functions to extract information from such pathway model objects. The built-in capabilities of Mathematica can the be utilized for computational analysis such as transient and steady state analysis, stability analysis, sensitivity analysis, and parameter estimation. Hence, PathwayLab and its integration with Mathematica provides a powerful environment for modeling and computational analysis of biochemical reaction networks.

Overview

The PathwayLab package provides functions to extract information from pathway model objects exported from the PathwayLab software.

This loads the package.

In[1]:=

In the Examples folder there is a file FirstExample.vsd that contains a small pathway model, see the figure below. There is also a file FirstExample.m that contains a Mathematica expression for the pathway model. In PathwayLab any pathway model can be exported to an m-file by the menu command Export to ... Mathematica.

        A small pathway model.

Load the model and store it in a variable. The standard print form of a pathway model object Pathway[ … ] contains information about the number of states, rates, and parameters of the model.

In[2]:=

Out[2]=

Functions for extracting equations from pwobjects.

The rate equations of the model.

In[3]:=

Out[3]=

The rate equations can also be presented with the parameters substituted by their default numerical values.

In[4]:=

Out[4]=

They can also be given as mass balances, i.e., formulated in terms of the rate variables.

In[5]:=

Out[5]=

Assuming constant compartment volume (V'[t]==0) gives even simplier expressions.

In[6]:=

Out[6]=

The initial conditions of the model.

In[7]:=

Out[7]=

Functions for extracting different types of symbol names from a pwobject.

The state variables of the model

In[8]:=

Out[8]=

The rate variables of the model

In[9]:=

Out[9]=

The model parameters. Note that parameter names also indicate to which rate expression the parameter belongs.

In[10]:=

Out[10]=

The model variables are all variables of the model except the state variables.

In[11]:=

Out[11]=

Functions returning lists of rules for parameters and model variables.

The parameters and their default values.

In[12]:=

Out[12]=

The model variables and their corresponding expression.

In[13]:=

Out[13]=

The rate expression for the first reaction.

In[14]:=

Out[14]=

NDSolve[pwobject, {t, 0, }]

solve numerically the rate equations for the state variables, with time in the range 0 to

Extension of NDSolve.

The solution of the rate equations for the first 20 seconds.

In[15]:=

Out[15]=

A plot of the concentration curves for each entity.

In[16]:=

BioCircuit

This example is a purely articficial biochemical network that illustrates different pathway modeling objects and how they are transformed into mathematical expressions. The model can be launched in PathwayLab from the Start Menu but its pathway diagram is also shown as plain graphics in the figure below. The file BioCircuit.m used below has been generated by Export to ...  Mathematica in Visio's PathwayLab menu.

    The BioCircuit pathway.

Load the package.

In[1]:=

Load the model and store it in a variable.

In[2]:=

Out[2]=

State variables in a pathway model correspond to concentrations of entities described by a set of ordinary differential equations - the rate equations. The full name of these variables are given by their location, see the figure above.

The state variables of the pathway.

In[3]:=

Out[3]=

A rate variable of a pathway model corresponds to the rate of a transformation (reaction or translocation) between two entities. The unit of a rate variable are either given in amount per time or in concentration per time depending on the settings in the Formula Settings dialog window for the transformation object in PathwayLab. In addition to the state and rate variables there are variables describing compartment volumes and auxiliary entities (described by explicit functions of time). The model variables are defined by expressions containing parameters and states.

The rate variables of the pathway.

In[4]:=

Out[4]=

The variables of the pathway model (all non-state variables).

In[5]:=

Out[5]=

The model variables that are not rates.

In[6]:=

Out[6]=

The expressions defining the model variables can be can be extracted as a list of rules.

In[7]:=

Out[7]=

The rate expression for rate r1.

In[8]:=

Out[8]=

The parameters of the model.

In[9]:=

Out[9]=

The values of the parameters is given by a list of rules.

In[10]:=

Out[10]=

The list of parameter rules using TableForm for nice formatting.

In[11]:=

Out[11]//TableForm=

The rate expression for rate r1 when the parameters have been substituted with their default numerical values.

In[12]:=

Out[12]=

The rate equations (or reaction rate equations) of a pathway model is a set of differential equations describing the time evolution of the entity concentrations.

The rate equations for the pathway model.

In[13]:=

Out[13]=

The rate equations when the parameters have been have been substituted with their numerical values.

In[14]:=

Out[14]=

Rate equations can also be presented in mass balance format relating rate of change of entity concentrations with the transformation rates. Since entity concentrations are used as state variables there is also a term in each differential equation corresponding to dilution for time varying volumes. Rates defined as absolute (in terms of amount per time) will be divided by the volume of the compartment or location containing the entity to obtain the corresponding concentration flow rate.

The rate equations in mass balance format.

In[15]:=

Out[15]=

The rate equations in mass balance format assuming constant volume and area of the loctions.

In[16]:=

Out[16]//TableForm=

The rate equations in symbolic form (given by the default usage of RateEquations) is obtained by substituting the expressions for the model variables into the rate equations in mass balance format.

In[17]:=

Out[17]=

Given the the rate equations and a set of initial conditions the model can be simulated using NDSolve.

Store the rate equations, states, and initial conditions in variables.

In[18]:=

Out[20]=

Simulate the system for 40 seconds.

In[21]:=

Out[21]=

A plot of state variables x4 and NucleusVolume`x4.

In[22]:=

A change of the parameter r5`Vm to 2.

In[24]:=

A new plot of x4 and NucleusVolume`x4.

In[26]:=

A comparison of x4 and NucleusVolume`x4 for the default and new value of r5`Vm.

In[27]:=

Since simulation of the rate equations with default parameter values and initial conditions is a common operation the PathwayLab package implements an extension of NDSolve for this particular task.

Simulate the rate equations for the model using default parameter values and initial conditions.

In[28]:=

Out[29]=

A plot of rate variables NucleusVolume`Membrane`r6[t] and r7[t]. Additional options defined by the Graphics`Legend` standard package are used to add a legend to the plot. Observe the usage of ModelVariableRules and ParameterRules to obtain the symbolic expressions for the rates. Note the use of ReplaceAllRepeated (//.) instead of ReplaceAll (/.) below which is required in general but not for this particular model.

In[30]:=

One way to understand how sensitive a pathway model is with respect to changes in a particular parameter is to make a number of simulations for different values of the parameter.

Replace the parameter r5`Vm with a dummy variable Vm and substitute the remaining parameters with their default values.

In[34]:=

Solve the rate equations for Vm between 0.5 and 1.5 in steps of 0.1.

In[35]:=

A plot of the solutions for x4.

In[36]:=

A MAPK Cascade

This model is taken from the paper Untangling the wires:  A strategy to trace functional interactions in signaling and gene networks by B. N. Kholodenko et al. (www.pnas.org/cgi/doi10.1073/pnas.192442699). The model can be launched in PathwayLab from the Start Menu but its pathway diagram is also shown as plain graphics in the figure below. The file MAPKCascade.m used below has been generated by Export to ...  Mathematica in Visio's PathwayLab menu.

    A  mitogen-activated protein kinase (MAPK) cascade.

Load the package.

In[1]:=

Load the model and store it in a variable.

In[2]:=

Out[2]=

The state variables of the pathway.

In[3]:=

Out[3]=

The rate variables of the pathway.

In[4]:=

Out[4]=

The mass balance equations.

In[5]:=

Out[6]//TraditionalForm=

Once the model equations are available in Mathematica it is easy to analyze them in various ways. Here we will compute the stoichiometric matrix.

The right hand side of the mass balance equations.

In[7]:=

Out[7]=

The stoichiometric matrix (in Mathematica version 5.1 this can be directly accomplished by D[rhs,{#[t]&/@rates}]).

In[8]:=

Out[9]//MatrixForm=

In matrix notation the mass balance equations become.

In[10]:=

Out[10]=

Linear dependent rows of the stoichiometric matrix implies that there are conserved moieties in the system. In this particular case it is easy to find these linear combinations by inspection. The corresponding linear relations between the state derivatives show that certain combinations of states are constant over time.

The linear combinations and a test that they map zero.

In[11]:=

Out[14]//MatrixForm=

The conservation equations.

In[15]:=

Out[15]=

A list of rules for the total amounts of MKKK, MKK, and MAPK.

In[16]:=

Out[16]=

Utilizing the existence of conserved moieties in the system the rate equations can be reduced.

The rate equation for all states except MKKK—PP, MKK—PP

In[17]:=

In Mathematica 5 NDSolve can handle some differential algebraic equations (NDSolve[mapk,{t,0,300}] gives the same result and can be used in older Mathematica versions).

In[18]:=

Out[18]=

A plot of MKKK—PP, MKK—PP, and MAPK—PP.

In[19]:=

A plot of MAPK, MAPK—P, MAPK—PP and their sum.

In[21]:=

An advantage of having the system in reduced form is that its steady state can be computed algebraically. For this example the rate expressions all are rational functions of the state. Hence, NSolve can be used to find all steady state solutions.

The equations for the steady state with functions that are explicitly dependent on time fixed to their values at t=0.

In[22]:=

Solutions of the steady state equations (this may take several seconds). There are quite many solutions!

In[23]:=

Out[24]=

The subset of solutions that are both real and positive.

In[25]:=

Out[25]=

The states approach these values close to t=100 in the transient simulation above.

In[26]:=

Out[26]=

In this case FindRoot also finds the steady state starting in the initial conditions for the rate equations (which is not a good approximation to the steady state).

In[27]:=

Out[27]=

Integral Feedback Control

The following model illustrates how robust perfect adaptation (integral feedback control) can be implemented by a small reaction network. It can be shown that only the blue and red interactions are crucial for this behavior and that the steady state value of Y only depends on parameters of reaction v3 and v4. The pathway model is taken from Robust perfect adaptation in bacterial chemotaxis through integral feedback control, Tau-Mu Yi, Yun Huang, Melvin I. Simon, and John Doyle. PNAS Vol. 97, Issue 9, 4649-4653, April 25, 2000, see www.pnas.org/cgi/content/abstract/97/9/4649. The model can be launched in PathwayLab from the Start Menu but its pathway diagram is also shown as plain graphics in the figure below. The file IntegralFeedbackControl.m used below has been generated by Export to ...  Mathematica in Visio's PathwayLab menu.

    A biochemical pathway implementing integral feedback control (robust perfect adaptation).

Load the package.

In[1]:=

Load the model and store it in a variable.

In[2]:=

Out[2]=

The state variables of the pathway.

In[3]:=

Out[3]=

The rate equations. Note the use of ReplaceAll (/.) and not ReplaceAllRepeated (//.) to prevent E1 and E2 to be substituted with the expressions containing UnitStep.

In[4]:=

Out[4]=

The steady state equations.

In[5]:=

Out[5]=

The steady state solution.

In[6]:=

Out[7]//TableForm=

Observe that the value of Y[∞] is independent of E1[∞] and E2[∞]. Hence, the steady state is independent of these signals. The condition v3`Vm-v4`Vm>0 has to be satisfied for the solution to be physically realizable.

Simulate the system.

In[8]:=

A plot of the solution for Y when E1 and E2 are step and pulse signals.

In[9]:=

To illustrate the robustness of the pathway behavior one can make a number of simulations for different values of some parameter.

Replace the parameter v1`Ki with a dummy variable Ki and substitute the remaining parameters with their default values.

In[13]:=

Solve the rate equations for Ki between 0.1 and 0.5 in steps of 0.1.

In[14]:=

A plot of the solutions for Y and A.

In[15]:=

Observe that the value of Y returns to the fixed level 0.125 for any setting of the parameter Ki and regardless of the changing flux profiles due to the changes in E1 and E2.

Cell Cycle Control in Xenopus Frog Eggs

This model is presented in Cell Cycle Control in Xenopus Frog Eggs (Novak and Tyson, Journal of Theoretical Biology, 1993) and illustrates oscillations in concentration of biochemical entities involved in cell cycle control. The model can be launched in PathwayLab from the Start Menu but its pathway diagram is also shown as plain graphics in the figure below. The file CellCycle.m used below has been generated by Export to ...  Mathematica in Visio's PathwayLab menu.

    A pathway model for cell cycle control in Xenopus frog eggs.

Load the package.

In[1]:=

Load the model and store it in a variable.

In[2]:=

Out[2]=

The state variables of the pathway.

In[3]:=

Out[3]=

In this example the given initial conditions correspond to a steady state of the system. The stability properties of the steady state can be investigated by linearizing the model around the steady state and study the stability properties of the obtained linear system.

The initial conditions.

In[4]:=

Out[4]=

Store the rate equations in a variable.

In[5]:=

The right hand side of the rate equations.

In[6]:=

The right hand side of the rate equations evaluated at t=0 for the given initial conditions and default parameter values.

In[7]:=

Out[7]=

The Jacobian matrix in symbolic form (in Mathematica version 5.1 this can be directly accomplished by D[rhs,{#[t]&/@states}]).

In[8]:=

The Jacobian matrix evaluated at the steady state gives the system matrix for the linearized system.

In[9]:=

Out[10]//MatrixForm=

The eigenvalues of the linearized system.

In[11]:=

Out[11]=

There is one pair of eigenvalues with positive real part which shows that the steady state is unstable. This can also clearly be seen in a simulation of the non-linear system. First the solution is very close to the initial conditions but then it starts to deviate and eventually enters a limit cycle.

Simulate the non-linear system for 1000 seconds.

In[12]:=

Out[12]=

A plot of MPF and APC.

In[13]:=

A phase diagram of MPF and APC.

In[15]:=

Store the numerical value of the steady state in a variable x0 and construct a symbolic vector representing deviations from the steady state.

In[16]:=

Out[17]=

The linearized system equations.

In[18]:=

Initial conditions for the linearized system (slightly off zero to excite the unstable behavior).

In[19]:=

Simulate the linear system for 1000 seconds.

In[20]:=

Out[20]=

A plot of MPF and APC for the first 400 seconds.

In[21]:=

The behavior of the linear system resembles the non-linear one for small deviations although the onset of visible oscillations are different (due to the arbitrary choice of small linear initial conditions).

References

InNetics - www.innetics.com.

Created by Mathematica  (August 17, 2005)