Introduction to Mathematica
Part II
May 14, 2009



Shirley Pepke, CACR

Last time: overview, basic syntax, examples

This time:
Advanced concepts: parallelization, packages, automatic code generation
Complex example: cell cycle regulation

Parallelization

•Multicore parallelism --> Local Kernels

Use one command: 'Parallelize' for operations on lists.
Parallel Kerenels Status Monitor has useful info.

AbsoluteTiming[ProteinData[#,"MoleculePlot"]&/@{"GGA2","PAH"}]

AbsoluteTiming[Parallelize[ProteinData[#,"MoleculePlot"]&/@{"GGA2","PAH"}]]

•Automatic Process Generation (like 'fork')

Idea is that every instance of a function f will be submitted as a separate process. Use 'ParallelSubmit' and 'WaitAll'.

AbsoluteTiming[Table[RandomReal[50,{5000}],{2}];]

AbsoluteTiming[WaitAll[Table[ParallelSubmit[RandomReal[50,{5000}]],{2}];]]

Example 2.

AbsoluteTiming[WaitAll[Sum[ParallelSubmit[{k},NIntegrate[x^(-1-1/k),{x,1,∞}]],{k,1,24}]]]

AbsoluteTiming[TableForm[Table[{k,Integrate[x^(-1-1/k),{x,1,∞}]},{k,1,24}],TableDepth->2]]

AbsoluteTiming[TableForm[WaitAll[Table[ParallelSubmit[{k},{k,Integrate[x^(-1-1/k),{x,1,∞}]}],{k,1,24}]],TableDepth->2]]

AbsoluteTiming[TableForm[ParallelTable[{k,Integrate[x^(-1-1/k),{x,1,∞}]},{k,1,24}],TableDepth->2]]

Much advanced support for process control as well as parallelization over networks (including heterogeneous clusters).

Automatic Code Generation

Consider a simple coupled ODE system.

sol=NDSolve[{y'[x]==y[x]Cos[z[x]],z'[x]==y[x],y[0]==1,z[0]==0},{y,z},{x,0,15}]

Plot the results using options defined in the Mathematica package 'PlotLegends'.

<<PlotLegends`
Plot[{y[x]/.sol,z[x]/.sol},{x,0,10},PlotRange->All,PlotStyle->Thick,PlotLegend->{"y","z"},LegendPosition->{1,-.4}]

Now suppose I want to add two more variables:

sol=NDSolve[{y'[x]==y[x]Cos[z[x]],z'[x]==y[x],q'[x]==q[x]Cos[2*z[x]],w'[x]==q[x]/2,y[0]==1,z[0]==0,q[0]==-1,w[0]==0},{y,z,q,w},{x,0,15}]
Plot[{y[x]/.sol,z[x]/.sol,q[x]/.sol,w[x]/.sol},{x,0,10},PlotRange->All,PlotStyle->Thick,PlotLegend->{"y","z","q","w"},LegendPosition->{1,-.4}]

Now suppose I want to add two more variables with the same pattern, and two more, and two more...I can use 'Table' within functions to automatically generate equations.

d1=a[1]'[x]==a[1][x]Cos[a[2][x]];
dn2[nn_]:=Table[a[n]'[x]==a[n-1][x]/(n/2),{n,2,nn,2}]
dn3[nn_]:=Table[a[n]'[x]==a[n][x]Cos[(n-1)*a[n-1][x]],{n,3,nn,2}]
initsn[nn_]:=Table[a[n][0]==0,{n,2,nn,2}]
inits2n[nn_]:=Table[a[n][0]==(-1)^(Floor[n/2]),{n,1,nn,2}]

nn=6;
eqns=Join[{d1},dn2[nn],dn3[nn],initsn[nn],inits2n[nn]];
vars=Table[a[n],{n,1,nn}];
sol=NDSolve[eqns,vars,{x,0,15}]
Plot[Table[a[n][x]/.sol,{n,1,nn}],{x,0,10},PlotRange->All,PlotStyle->Thick]

Packages

I can make my set of functions into a package to be used by others by declaring some cells as initialization cells.

Automatic Code Generation - Application to Chemical Kinetics

xCellerator (Shapiro, BE, Levchenko A, Meyerowitz EM, Wold BJ, Mjolsness ED. Cellerator: extending a computer algebra system to include biochemical arrows for signal transduction simulations. Bioinformatics 19(5):677-678 (2003) translates chemical reaction networks to systems of ordinary differential equations, solves them, and plots the results.
Online: http://xlr8r.info

Good for networks that have combinatoric reactions or for networks with arbitrary numbers of components.

<<xlr8r.m

Basic example: a single equilibrium reaction

Consider a substance A which transforms at a rate kf to substance B, i.e. A->B at rate kf and also its reverse reaction B->A at rate kr. You can find this reaction on the xCellerator Palette under "Mass Action".

The Law of Mass Action states that the rate of change of the concentration of a reactant due to a particular reaction is equal to the reaction rate times the concentration of the contributing components.

system={{A⇄B, kf, kr}};
interpret[system]

kf=2;
kr=3;
sim=run[interpret[system], {t,0,15}, initialConditions->{A[0]==1, B[0]==0}]

runPlot[sim, PlotStyle->Thickness[.01]]

A (much) more complex example: Regulation of Cell Division

Murray AW and Kirschner MW. Dominoes and Clocks: The Union of Two Views of the Cell Cycle, Science, Vol. 246(4930), 1989.
Gonze D and Goldbeter A. A model for a network of phosphorylation-dephosphorylation cycles displaying the dynamics of dominoes and clocks, Journal of Theoretical Biology, Vol. 210, 2000.

Graphic Depiction of Two Views of the Cell Cycle

SetDirectory["/Users/spepke"];
diagrams=Import["dc.jpg","JPEG"]

Graphic Depiction of Gonze and Goldbeter's Model

SetDirectory["/Users/spepke"];
diagrams=Import["gg.jpg","JPEG"]

Bruce Shapiro's xCellerator implementation: http://xlr8r.info/examples