Mathematica 9 is now available
PathwayLab Research Edition Products
-----
 /
PathwayLab Research Edition
*Features
*Examples
<BioCircuit
*A MAPK Cascade
*Integral Feedback Control
*Cell Cycle Control in Xenopus Frog Eggs
*Buy Online
*Trial Version
*For More Information
*Life & Medical Sciences Solutions
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

Example: BioCircuit

This example is a purely artificial biochemical network that illustrates different pathway modeling objects and how they are transformed into mathematical expressions. The model can be launched in PathwayLab Research Edition from the Start menu, but its pathway diagram is also shown as plain graphics in the figure below. The file BioCircuit.nb* used below was generated by Export to > Mathematica in the PathwayLab menu.

    [Graphics:HTMLFiles/biocircuit_1.gif]

    The BioCircuit pathway

Load the package.

In[1]:=

Needs["PathwayLab`"]

Load the model and store it in a variable.

In[2]:=

biocircuit = Get[ToFileName[{"PathwayLab", "Examples"}, "BioCircuit.m"]]

Out[2]=

-Pathway[8 state variables, 10 rate variables, 19 parameters]-

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

This lists the state variables of the pathway.

In[3]:=

states = StateVariables[biocircuit]

Out[3]=

{x2, x1, x3, x5, x4, NucleusVolume`x4, MembraneArea`x1, MembraneArea`x1P}

A rate variable of a pathway model corresponds to the rate of a transformation (reaction or translocation) between two entities. The units 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 Research Edition. 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.

This lists the rate variables of the pathway.

In[4]:=

rates = RateVariables[biocircuit]

Out[4]=

{r1, r10, r3, r4, NucleusVolume`Membrane`r6, r2, r5, MembraneArea`r8, MembraneArea`r9, r7}

This lists all non-state variables of the pathway model.

In[5]:=

modelvariables = ModelVariables[biocircuit]

Out[5]=

{V, u1, r1, u2, x0, r10, r3, r4, NucleusVolume, NucleusVolume`Membrane, NucleusVolume`Membrane`r6, r2, r5, MembraneArea, MembraneArea`r8, MembraneArea`r9, r7}

This lists the model variables that are not rates.

In[6]:=

Complement[modelvariables, rates]

Out[6]=

{NucleusVolume`Membrane, MembraneArea, NucleusVolume, u1, u2, V, x0}

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

In[7]:=

ModelVariableRules[biocircuit]

Out[7]=

This gives the rate expression for rate r1.

In[8]:=

r1[t]/.ModelVariableRules[biocircuit]

Out[8]=

(r1`Vm x1[t])/(r1`Km (1 + x1[t]/r1`Km + x2[t]/r1`Ki))

Here are the parameters of the model.

In[9]:=

Parameters[biocircuit]

Out[9]=

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

In[10]:=

ParameterRules[biocircuit]

Out[10]=

This provides the list of parameter rules, using TableForm for nice formatting.

In[11]:=

ParameterRules[biocircuit]//TableForm

Out[11]//TableForm=

r1`Vm→1
r1`Km→1
r1`Ki→0.2
r10`Vm→1
r10`Km→1
r3`alpha→0.2
r4`consrateE→1
NucleusVolume`Membrane`r6`Dconst→0.2
r2`Vm→1
r2`Km→1
r5`Vm→1
r5`Km→1
MembraneArea`Area→0.01
MembraneArea`r8`Vm→0.1
MembraneArea`r8`Km→0.4
MembraneArea`r9`k0→0.2
MembraneArea`r9`Km→1
r7`kf→0.5
r7`kr→0.01

This provides the rate expression for rate r1 when the parameters have been replaced by their default numerical values.

In[12]:=

r1[t]/.ModelVariableRules[biocircuit]/.ParameterRules[biocircuit]

Out[12]=

x1[t]/(1 + x1[t] + 5. x2[t])

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

This gives the rate equations for the pathway model.

In[13]:=

RateEquations[biocircuit]

Out[13]=

Here are the rate equations when the parameters have been replaced by their numerical values.

In[14]:=

RateEquations[biocircuit, EquationFormat→Numeric]

Out[14]=

Rate equations can also be presented in mass balance format, relating rate of change of entity concentrations with the transformation rates. Because 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.

Here are the rate equations in mass balance format.

In[15]:=

RateEquations[biocircuit, EquationFormat→MassBalances]

Out[15]=

This provides the rate equations in mass balance format assuming constant volume and area of the locations.

In[16]:=

RateEquations[biocircuit, EquationFormat→MassBalances]/.{V '[t] →0, MembraneArea '[t] →0, NucleusVolume '[t] →0}//TableForm

Out[16]//TableForm=

x2^′[t] == r1[t] - r2[t]
x1^′[t] == -r1[t] + r10[t] - r7[t]
x3^′[t] == r2[t] - r5[t]
x5^′[t] == r3[t] - r4[t]
x4^′[t] == r5[t] - NucleusVolume`Membrane`r6[t]/V[t]
NucleusVolume`x4^′[t] == NucleusVolume`Membrane`r6[t]/NucleusVolume[t]
MembraneArea`x1^′[t] == -MembraneArea`r8[t] + MembraneArea`r9[t] + (r7[t] V[t])/MembraneArea[t]
MembraneArea`x1P^′[t] == MembraneArea`r8[t] - MembraneArea`r9[t]

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

In[17]:=

RateEquations[biocircuit] === (RateEquations[biocircuit, EquationFormat→MassBalances]//.ModelVariableRules[biocircuit])

Out[17]=

True

Given 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]:=

deq = RateEquations[biocircuit, EquationFormat→Numeric] ;

states = StateVariables[biocircuit] ;

ic = InitialConditions[biocircuit]

Out[20]=

{x2[0] == 0.1, x1[0] == 1, x3[0] == 1, x5[0] == 1, x4[0] == 0, NucleusVolume`x4[0] == 0, MembraneArea`x1[0] == 0.5, MembraneArea`x1P[0] == 0}

Simulate the system for 40 seconds.

In[21]:=

sol = NDSolve[Join[deq, ic], states, {t, 0, 40}]//First

Out[21]=

This plots state variables x4 and NucleusVolume`x4.

In[22]:=

Needs["Graphics`Colors`"] ;

pl1 = Plot[Evaluate[{x4[t], NucleusVolume`x4[t]}/.sol], {t, 0, 40}, PlotRange→All, PlotStyle→ {Red, Blue}] ;

[Graphics:HTMLFiles/biocircuit_49.gif]

This changes the parameter r5`Vm to 2.

In[24]:=

deq2 = RateEquations[biocircuit]/.r5`Vm→2/.ParameterRules[biocircuit] ;

sol2 = NDSolve[Join[deq2, ic], states, {t, 0, 40}]//First ;

Here is a new plot of x4 and NucleusVolume`x4.

In[26]:=

pl2 = Plot[Evaluate[{x4[t], NucleusVolume`x4[t]}/.sol2], {t, 0, 40}, PlotRange→All, PlotStyle→ {{Red, Dashing[{0.02, 0.02}]}, {Blue, Dashing[{0.02, 0.02}]}}] ;

[Graphics:HTMLFiles/biocircuit_53.gif]

This compares x4 and NucleusVolume`x4 for the default and new value of r5`Vm.

In[27]:=

Show[pl1, pl2] ;

[Graphics:HTMLFiles/biocircuit_55.gif]

Because 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]:=

sol3 = NDSolve[biocircuit, {t, 0, 40}]//First ;

sol === sol3

Out[29]=

True

This plots 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 below of ReplaceRepeated (//.) instead of ReplaceAll (/.), which is required in general but not for this particular model.

In[30]:=

Needs["Graphics`Legend`"] ;

mvRules = ModelVariableRules[biocircuit] ;

paramRules = ParameterRules[biocircuit] ;

[Graphics:HTMLFiles/biocircuit_63.gif]

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 replace the remaining parameters with their default values.

In[34]:=

deq3 = RateEquations[biocircuit]/.r5`Vm→Vm/.ParameterRules[biocircuit] ;

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

In[35]:=

multiplesol = Table[NDSolve[Join[deq3, ic], states, {t, 0, 40}]//First, {Vm, .5, 1.5, .1}] ;

Plot the solutions for x4.

In[36]:=

Plot[Evaluate[x4[t]/.multiplesol], {t, 0, 40}, PlotStyle→Red, PlotRange→All] ;

[Graphics:HTMLFiles/biocircuit_67.gif]

* Note: If you do not own a copy of Mathematica, you will need to download a trial version in order to view the notebook document.