 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Mathematical Modelling and Scientific Computing in the Biosciences

27 March 2007

Lecture 2: Overview

● Singular Perturbation:  Michaelis-Menten Kinetics
● Cooperative Phenomena, Hill-Function
● Motifs in Biological Circuits, Dynamical Properties
    - Negative Auto-Regulation
    - Feed-Forward Loop

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Catalytic reaction: enzyme (Enz) and substrate (S) forming complex EnzS, then giving rise to product (P) and enzyme:
Enz+S Underoverscript[⇔, k1b, arg3] EnzSOverscript[→, k2d]Enz+P

Get the corresponding ODEs using Cellerator:

In[1]:=

<<cellerator.m ;

EnzymeReactions = { {Enz + S⇔EnzS, k1f, k1b}, {EnzS→P + Enz, k2f}} ;

Cellerator™ 1.5.8 (1-July-2005) loaded 27-Mar-2007 09:32 using Mathematica Version 5.2 for Linux (June 20, 2005)

In[3]:=

LargeBlue[x_] := StyleForm[x, FontColor→ RGBColor[0, 0, 1], FontWeight→ "Bold", FontSize→ 15] ;

EnzymeODE = interpret[EnzymeReactions]//First//ColumnForm//LargeBlue

Out[4]//StyleForm=

Enz^′[t] == k1b EnzS[t] + k2f EnzS[t] - k1f Enz[t] S[t]
EnzS^′[t] == -k1b EnzS[t] - k2f EnzS[t] + k1f Enz[t] S[t]
P^′[t] == k2f EnzS[t]
S^′[t] == k1b EnzS[t] - k1f Enz[t] S[t]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Last time, we simplified the ODE system to:

In[5]:=

HighlightEpsRule = {ε→StyleForm[ε, FontColor→ RGBColor[1, 0, 0], FontWeight→ "Bold"] } ;

%/. HighlightEpsRule//ColumnForm//LargeBlue

Out[7]//StyleForm=

ε Overscript[EnzS,^]^′[Overscript[t,^]] == -(1 + (k1b + k2f)/(k1f STotal[0])) Overscript[EnzS,^][Overscript[t,^]] + Overscript[S,^][Overscript[t,^]]
Overscript[S,^]^′[Overscript[t,^]] == -Overscript[S,^][Overscript[t,^]] + Overscript[EnzS,^][Overscript[t,^]] (k1b/(k1f STotal[0]) + Overscript[S,^][Overscript[t,^]])

with the small (dimensionless) parameter ε ≡ EnzTotal[0]/STotal[0].

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

In[210]:=

[Graphics:HTMLFiles/Lecture3_40.gif]

 ◂ ▸  1 of 1

Singular Perturbation: Michaelis-Menten Kinetics

In[244]:=

EpsList = {0.1, 0.03, 0.01} ;

PlotOrig = SolveAndPlot/@ EpsList//Show ;

PlotBlownUp = SolveAndPlotBlownUp/@EpsList//Show ;

Show[ ({PlotOrig, PlotBlownUp}//GraphicsArray), DisplayFunction→ $DisplayFunction, ImageSize→ {500, 180}] ;

[Graphics:HTMLFiles/Lecture3_51.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Observation: it seems that in terms of the new (magnified) time variable τ ≡ Overscript[t,^]/ε, the inner (singular) solution is essentially independent of parameter ε.

Regular perturbation failed because the solution is not analytic in the perturbative parameter ε. Maybe the inner solution is an analytic function when in terms of the new variable τ. Try this time-transformation.

In[124]:=

Inner—ODE—IC = Join[%, { EnzS—I[0] == 0, S—I[0] == 1}] ;

(%/.{ Derivative[n_][Func_][t_] →  Derivative[n][Func][t]/ε^n})//Simplify[#, Assumptions→ {ε>0}] & ;

Inner—ODE—IC = %/.{Overscript[t,^] →τ} ;

%/.HighlightEpsRule//ColumnForm

Out[128]=

EnzS—I[τ] (1 + (k1b + k2f)/(k1f STotal[0])) + EnzS—I^′[τ] == S—I[τ]
EnzS—I[τ] ε (k1b/(k1f STotal[0]) + S—I[τ]) == ε S—I[τ] + S—I^′[τ]
EnzS—I[0]==0
S—I[0]==1

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Now do perturbative analysis: expand the solution in terms of the small parameter,  ε:

In[222]:=

MaxOrd = 1 ;   

That is, we have expanded the inner solution as:

In[224]:=

{S—Inner[t],   EnzS—Inner[t]}/.HighlightEpsRule

Out[224]=

{S—I[0][t] + ε S—I[1][t], EnzS—I[0][t] + ε EnzS—I[1][t]}

Therefore, the inner ODE system is:

In[225]:=

Expanded—ODE—IC = Inner—ODE—IC/.{ S—I→ S—Inner, EnzS—I→EnzS—Inner}//Simplify ;

%/.HighlightEpsRule//ColumnForm//StyleForm[#, FontSize-> 7] &

Out[226]//StyleForm=

EnzS—I[0][0]+ε EnzS—I[1][0]==0
S—I[0][0]+ε S—I[1][0]==1

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Now let's obtain ODEs for each order. Zeroth order equations for the inner (singular) solution:

In[259]:=

<<PerturbationODE.m

OrderList = PolyOrderList[Inner—ODE—IC, {S—I[τ], EnzS—I[τ]}, ε, 1] ;

OrderList[[1]]//ColumnForm

Out[261]=

EnzS—I[0][τ] + (k1b EnzS—I[0][τ])/(k1f STotal[0]) + (k2f EnzS—I[0][τ])/(k1f STotal[0]) + EnzS—I[0]^′[τ] == S—I[0][τ]
0 == S—I[0]^′[τ]
EnzS—I[0][0]==0
S—I[0][0]==1

In[262]:=

Order0Soln = DSolve[OrderList[[1]], { EnzS—I[0], S—I[0]}, τ]//Simplify//Flatten ;

%//ColumnForm

Out[263]=

EnzS—I[0] →Function[{τ}, -((-1 + ^(τ (-k1b - k2f - k1f STotal[0]))/(k1f STotal[0])) k1f STotal[0])/(k1b + k2f + k1f STotal[0])]
S—I[0]→Function[{τ},1]

First order equations for the inner (singular) solution:

In[264]:=

Order1Soln = DSolve[OrderList[[2]]/.Order0Soln, { EnzS—I[1], S—I[1]}, τ]//Simplify//Flatten ;

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

In[265]:=

PositivityAssumption = {STotal[0] >0, k1f> 0, k1b>0,   k2f>0} ;

Limit—EnzS—Inner—Order0 = Limit[EnzS—I[0]/.Order0Soln//Simplify//Last, τ→ Infinity, Assumptions→%]

Out[266]=

(k1f STotal[0])/(k1b + k2f + k1f STotal[0])

In[267]:=

Limit—S—Inner—Order0 = Limit[S—I[0]/.Order0Soln//Simplify//Last, τ→ Infinity]

Out[267]=

1

Inner/Outer Solutions: Matching Condition

In the limit of ε → 0, the width of the boundary layer tends to zero.        
Therefore, for the outer solution, its value at the boundary layer → value at t = 0.  
However, in the same limit of ε → 0, we have τ ≡ Overscript[t,^]/ε → ∞.
Therefore, matching condition from continuity of solution requires that:
     Underscript[lim, τ→∞] {S—Inner[τ],  EnzS—Inner[τ]} =   Underscript[lim, t→0]  {S—Outer[t],  EnzS—Outer[t]}

Thus, the inner solution provides the initial conditions (t=0) for the outer solution.

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

[Graphics:HTMLFiles/Lecture3_110.gif]

Inner/Outer Solutions: Matching Condition

In the limit of ε → 0, the width of the boundary layer tends to zero.        
Therefore, for the outer solution, its value at the boundary layer → value at t = 0.  
However, in the same limit of ε → 0, we have τ ≡ Overscript[t,^]/ε → ∞.
Therefore, matching condition from continuity of solution requires that:
     Underscript[lim, τ→∞] {S—Inner[τ],  EnzS—Inner[τ]} =   Underscript[lim, t→0]  {S—Outer[t],  EnzS—Outer[t]}

Thus, the inner solution provides the initial conditions (t=0) for the outer solution.

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Singular Perturbation: Michaelis-Menten Kinetics

Using singular perturbation, we have derived the correct IC for the outer solution:

In[44]:=

CorrectOuterIC = {S—Outer[0] == 1, EnzS—Outer[0] == (k1f STotal[0])/(k1b + k2f + k1f STotal[0])} ;

Compare this to the (inconsistent) IC obtained from regular perturbation (last lecture):

In[45]:=

WrongOuterIC = {S—Outer[0] == 1, EnzS—Outer[0] == 0} ;

In particular, look at the outer ODE, at zeroth order in ε:

In[46]:=

OuterOrder0ODE = SimpEnzymeODE/.{Overscript[EnzS,^] → EnzS—Outer, Overscript[S,^] → S—Outer }/.{ε→ 0}//Simplify ;

%//ColumnForm

Out[47]=

EnzS—Outer[Overscript[t,^]] (1 + (k1b + k2f)/(k1f STotal[0])) == S—Outer[Overscript[t,^]]
EnzS—Outer[Overscript[t,^]] (k1b/(k1f STotal[0]) + S—Outer[Overscript[t,^]]) == S—Outer[Overscript[t,^]] + S—Outer^′[Overscript[t,^]]

In[82]:=

Together//@(OuterOrder0ODE[[1]]/.{Overscript[t,^] → 0})

Out[82]=

(EnzS—Outer[0] (k1b + k2f + k1f STotal[0]))/(k1f STotal[0]) == S—Outer[0]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Michaelis-Menten Kinetics: Inner Solution

In[49]:=

EnzS—IOrder0[τ_] := Evaluate[EnzS—I[0]/.Order0Soln//Last] ;

HighlightRule[x_] := x→StyleForm[x, FontColor→ RGBColor[0., 0., 1], FontSize→ 12] ;

EnzS—IOrder0[τ]/.HighlightRule[τ]

Out[51]=

-((-1 + ^((-k1b - k2f - k1f STotal[0]) τ)/(k1f STotal[0])) k1f STotal[0])/(k1b + k2f + k1f STotal[0])

In[52]:=

UnitConstantsRule = {STotal[0] → 1, k1b→ 1, k1f→ 1, k2f→ 1} ;

EnzS—IOrder0Plot = Plot[EnzS—IOrder0[τ]/.%, {τ, 0, 1}, PlotStyle→ {Hue[1], Thickness[0.008]}, DisplayFunction→ Identity] ;

Show[PlotBlownUp, EnzS—IOrder0Plot, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/Lecture3_139.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Michaelis-Menten Kinetics: Outer Solution

In[55]:=

EnzS—Order0OuterRule = Solve[OuterOrder0ODE[[1]], EnzS—Outer[Overscript[t,^]]]//Flatten ;

S—OuterODE = OuterOrder0ODE[[2]]/.EnzS—Order0OuterRule//Simplify ;

S—OuterSoln = NDSolve[ Join[{%}, {S—Outer[0] == 1}]/.UnitConstantsRule, S—Outer,    {Overscript[t,^], 0, 1}]//Flatten ;

EnzS—OuterSoln = EnzS—Order0OuterRule/.S—OuterSoln/.UnitConstantsRule ;

EnzS—OuterOrder0Plot = Plot[EnzS—Outer[Overscript[t,^]]/.%, {Overscript[t,^], 0, 1}, PlotStyle→ {Hue[1], Thickness[0.008]}, DisplayFunction→ Identity] ;

Show[PlotOrig, %, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/Lecture3_151.gif]

Conclusion: the outer perturbative solution captures the long-time dynamics well, especially for small values of ε; the initial behavior is well-captured by the inner solution.

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Cooperative Phenomenon

● Substrate binds to the enzyme on the binding site.
● There are some enzymes with more than 1 binding site.
● Reaction between enzyme and substrate is called cooperative, if when the enzyme is bond to a substrate molecule, other substrate molecule can bind to the remaining binding sites.
● Indirect interaction between distinct binding sites: allosteric effect.

        S+EUnderoverscript[⇔, k1 -, arg3]C_1Overscript[→, k2] E+P

               S+C_1Underoverscript[⇔, k3 -, arg3]C_2Overscript[→, k4]C_1+P


● Again, assume ε ≡ EnzymeTotal/STotal≪ 1
● From singular perturbation analysis, obtain expression for the substrate flux/reaction velocity (see book by J. Murray, Intro. Math. Biol.):

        |ds/dt| =  (EnzTotal  STotal   (k2  K_m^'    + k4  STotal))/(K_m K_m^' + K_m^' STotal   + STotal^2)

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Cooperative Phenomenon: Hill-Function

● When cooperativity is suspected, usually assume reaction velocity to take the following Hill-function form:
    |ds/dt| =  (Q * STotal^n )/(K_m + STotal^n);
where n is usually taken to be an integer. Typically, 2 ≤ n ≤ 4.
● n takes the role number of molecules taking part in the reaction.
   In the case n = 1, there is no cooperativity; reduce to Michaelis-Menten, as analyzed.
   Larger n, higher non-linearity.

In[78]:=

Hill—Function[x_, n_] := Q x^n/(Km + x^n) ; ListOfN = {1, 2, 4, 8} ;

Hill—Plots = Map[ Plot[Hill—Function[x, #]/.{Q→ 1, Km→ 1}, {x, 0, 2}, DisplayFunction→ Identity, PlotStyle→ {Hue[#/5.], Thickness[0.015]}] &, ListOfN] ;

[Graphics:HTMLFiles/Lecture3_178.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Cooperative Phenomenon: Hill-Function

Besides activating Hill-functions, repressing Hill-functions are also commonly used in biological modelling.

In[81]:=

Repressing—Hill—Function[x_, n_] := Q /(Km + x^n) ;

[Graphics:HTMLFiles/Lecture3_187.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Regulation

In[85]:=

Show[GraphicsArray[{ActivatingHill, RepressingHill}], ImageSize→ { 500, 120}] ;

[Graphics:HTMLFiles/Lecture3_194.gif]

● Such activators (positive control) and repressors (negative control) are common in gene regulation
● Cell: thousands of interacting genes, acting together to respond to diverse signals
    - external: temperature, level of nutrients, harmful chemicals
    - internal: level of metabolites, damage to proteins, DNA, etc
● Information-processing capability: biological circuitry
● Environment ⇔ transcription factors  (proteins transiting between active and inactive molecular states)
                ⇔ genes
● Set of interactions: transcription network

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Regulation

● Schematic of interactions between the environment, transcription factors, genes (U. Alon, Intro. Systems Biology):

In[86]:=

Import["~/Teaching/NoteBooks/Lecture3/GeneNetwork.jpg", "JPG"]//Show[#, ImageSize→ {500, 250}] & ;

[Graphics:HTMLFiles/Lecture3_201.gif]

● How do transcription factors influence the on/off states of genes?

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Regulation

In[87]:=

Show[Import["~/Teaching/NoteBooks/Lecture3/State_Active.jpg", "JPG"], ImageSize→ {500, 320}] ;

[Graphics:HTMLFiles/Lecture3_208.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Regulation

In[88]:=

Show[Import["~/Teaching/NoteBooks/Lecture3/State_Inactive.jpg", "JPG"], ImageSize→ {500, 320}] ;

[Graphics:HTMLFiles/Lecture3_215.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Transcription Network

In[89]:=

Import["~/Teaching/NoteBooks/Lecture3/GeneNetwork.jpg", "JPG"] ;

Show[%, ImageSize→ {500, 280}] ;

[Graphics:HTMLFiles/Lecture3_223.gif]

● What kind of patterns/motifs are observed in gene transcription networks?

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Networks

We can generate a list of random, connected, 3-node, activator/repressor graphs:

In[112]:=

<<DiscreteMath`Combinatorica`

PlotNum = 80 ; Prob = 0.3 ; NumVert = 3 ;

RandomGraphList = Table[ RandomGraph[NumVert, Prob, Type→ Directed], {PlotNum} ] ;

DistinctConnectedGraphs = Pick[%, Map[(ConnectedQ[#] == True) &, %]]//Union ;

Show[GraphicsArray[Partition[%, 7]], ImageSize→ {500, 170}, PlotRange-> All] ;

[Graphics:HTMLFiles/Lecture3_236.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Networks

A larger list of 4-node graphs:

In[119]:=

NumVert = 4 ;

RandomGraphList = Table[ RandomGraph[NumVert, Prob, Type→ Directed], {PlotNum} ] ;

DistinctConnectedGraphs = Pick[%, Map[(ConnectedQ[#] == True) &, %]]//Union ;

Show[GraphicsArray[Partition[%, 10]], ImageSize→ {500, 210}, PlotRange-> All] ;

[Graphics:HTMLFiles/Lecture3_247.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Gene Networks

● It has been found that gene transcription networks are not random graphs:
    - Motifs occur much more abundantly
    - Anti-Motifs occur much more rarely

Relative abundance of 8 Feed-Forward-Loop (FFL) types in the transcription network of yeast and E. coli. (from Mangan et. al, 2006):

In[78]:=

Import["~/Teaching/NoteBooks/Lecture3/FeedForwardMotifOccurence.jpg", "JPG"] ;

Show[%, ImageSize→ {500, 220}] ;

[Graphics:HTMLFiles/Lecture3_255.gif]

 ◂ ▸  FormBox[StyleBox[RowBox[{CounterBox[SlideShowNavigationBar], of , CounterBox[SlideShowNavigationBar, {None, SlideShowHeader, -1}]}], SR], StandardForm]

Network Motifs and Dynamical Implications

● Given that mutations occur randomly, it is probable that the abundance of observed motifs confer some dynamical properties that for its selection:
    - robustness to noise (environmental fluctuations)
    - information processing capabilities


● Negative Autoregulation:
    - insensitive to protein production rate (a noisy process)

● Feed-Forward Loop:
    - sign-sensitive delay
    - persistence detector/filtering of brief fluctuating signals
    

Created by Mathematica  (March 27, 2007) Valid XHTML 1.1!