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

Mathematical Modelling and Scientific Computing in the Biosciences

20 March 2007

Lecture 2: Overview

Topics


● Mathematica Preliminaries
● Mass-Action Kinetics
● Michaelis-Menten Kinetics
● Perturbation Analysis:
    -regular perturbation
    -singular perturbation

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

Mathematica Preliminaries

List: collection of general expressions        { ... }

In[8]:=

Table[f[i], {i, 1, 3}]

Out[8]=

{f[1], f[2], f[3]}

Replacement rule: transforming each subpart of an expression      /.

In[1]:=

a /.{a→ b}

Out[1]=

b

In[2]:=

f[a] /.{a→ b}

Out[2]=

f[b]

In[14]:=

a == b/.Equal→ Less

Out[14]=

a<b

Pattern: for representing forms of Mathematica expressions       x_

In[15]:=

{a^b, c^d, e^f}/.{x_^y_→ y^x}

Out[15]=

{b^a, d^c, f^e}

In[16]:=

{a * b == c, c * d == e, e * f == g}/.{x_ * y_ → y}

Out[16]=

{b == c, d == e, f == g}

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

Mathematica Preliminaries

Map: apply a function to each element of an expression        Map[ Function, List]

In[12]:=

Map[ Sin, {a, b, c}]

Out[12]=

{Sin[a], Sin[b], Sin[c]}

In[7]:=

Map[ 2^#&, {a, b, c}]

Out[7]=

{2^a, 2^b, 2^c}

Pick: take elements out of a list, according to some criterion         Pick[List, Criterion]

In[23]:=

NumberList = {2, 3, 4} ;

Map[#>2&, NumberList]

Pick[NumberList, % ]

Out[24]=

{False, True, True}

Out[25]=

{3, 4}

Solve: solve system of algebraic equations, for selected variables       Solve[ LHS==RHS, SolvedVariables ]

In[9]:=

Solve[a x^2 + b x + c == 0, x]

Out[9]=

{{x→ (-b - (b^2 - 4 a c)^(1/2))/(2 a)}, {x→ (-b + (b^2 - 4 a c)^(1/2))/(2 a)}}

 ◂ ▸  2 of 4

Mass-Action Kinetics

● What is a flux?
    flux of a chemical species is the rate of change in the concentration per unit time
        Units of concerntration: mole/liter ≈ 6.022*10^23 entities/liter, ...
        Units of time: second, minute
             ∴ Units of flux: mole/(liter * second), ...
             Mathematical expression for expressing the flux of species A: d [A]/(d t)       

● Rate laws: how does the flux depend on the concentration of the reactants?

 ◂ ▸  3 of 4

Mass-Action Kinetics

● Mass-action rate law: reaction flux is proportional to the probability of finding reacting molecules in a small volume.

In[906]:=

NumMolecule1 = 40 ;

Molecule1Location = Table[{Random[], Random[]}, {i, 1, NumMolecule1}] ;

PlotCircleRadius = 0.03 ;

PlotOffMargin = 0.05 ;

Plot1 = PlotMolecule[Hue[0.9],   Molecule1Location] ;

Show[Plot1, DisplayFunction→ $DisplayFunction, ImageSize→ {250, 250}, Frame→ True, FrameTicks→ None] ;

[Graphics:HTMLFiles/Lecture2_58.gif]

 ◂ ▸  3 of 4

Mass-Action Kinetics

● Mass-action rate law: reaction flux is proportional to the probability of finding reacting molecules in a small volume.

In[835]:=

NumMolecule2 = 5 ;

Molecule2Location = Table[{Random[], Random[]}, {i, 1, NumMolecule2}] ; <br />

(* Plotting intersection lines *)

DistanceList = Tuples[{Molecule1Location, Molecule2Location}] ;

PickedElements = Map[ Norm[#[[1]] - #[[2]]] < PlotCircleRadius * 2&, DistanceList  ] ;

LineIntersection = Show[Graphics[Map[{ Thickness[0.006], Line[#]} &, Pick[DistanceList,   PickedElements]]], DisplayFunction→ Identity] ;

Plot2 = PlotMolecule[Hue[0.7],   Molecule2Location] ;

PlotCombinedSpecies—Density5 = Show[{Plot1, Plot2, LineIntersection}, DisplayFunction→ $DisplayFunction, ImageSize→ {250, 250}, Frame→ True, FrameTicks→ None] ;

[Graphics:HTMLFiles/Lecture2_72.gif]

 ◂ ▸  3 of 4

Mass-Action Kinetics

● Mass-action rate law: reaction flux is proportional to the probability of finding reacting molecules in a small volume.

In[863]:=

NumMolecule2 = 10 ;

Molecule2Location = Table[{Random[], Random[]}, {i, 1, NumMolecule2}] ;

(* Plotting intersection lines *)

DistanceList = Tuples[{Molecule1Location, Molecule2Location}] ;

PickedElements = Map[ Norm[#[[1]] - #[[2]]] < PlotCircleRadius * 2&, DistanceList  ] ;

LineIntersection = Show[Graphics[Map[{ Thickness[0.006], Line[#]} &, Pick[DistanceList,   PickedElements]]], DisplayFunction→ Identity] ;

Plot2 = PlotMolecule[Hue[0.7],   Molecule2Location] ;

PlotCombinedSpecies—Density10 = Show[{Plot1, Plot2, LineIntersection}, DisplayFunction→ $DisplayFunction, ImageSize→ {250, 250}, Frame→ True, FrameTicks→ None] ;

[Graphics:HTMLFiles/Lecture2_86.gif]

 ◂ ▸  3 of 4

Mass-Action Kinetics

● Mass-action rate law: reaction flux is proportional to the probability of finding reacting molecules in a small volume.

In[870]:=

NumMolecule2 = 40 ;

Molecule2Location = Table[{Random[], Random[]}, {i, 1, NumMolecule2}] ; <br />

(* Plotting intersection lines *)

DistanceList = Tuples[{Molecule1Location, Molecule2Location}] ;

PickedElements = Map[ Norm[#[[1]] - #[[2]]] < PlotCircleRadius * 2&, DistanceList  ] ;

LineIntersection = Show[Graphics[Map[{ Thickness[0.006], Line[#]} &, Pick[DistanceList,   PickedElements]]], DisplayFunction→ Identity] ; <br />

Plot2 = PlotMolecule[Hue[0.7],   Molecule2Location] ;

PlotCombinedSpecies—Density40 = Show[{Plot1, Plot2, LineIntersection}, DisplayFunction→Identity, ImageSize→ {300, 300}, Frame→ True, FrameTicks→ None] ;

[Graphics:HTMLFiles/Lecture2_101.gif]

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

Mass-Action Kinetics

Rate law:

"The rate of the single chemical reaction (A+B→ C) is directly proportional to the product of the concentrations of the participating species, A and B"
↕
d [A]/(d t)=d [B]/(d t)= - k_AB * [A] * [B]
d [C]/(d t)=   k_AB * [A] * [B]

"The rate of the reversible chemical reaction (A+B⇔C+D) is directly proportional to the product of the concentrations of the participating species, A and B"
↕
d [A]/(d t)=d [B]/(d t)= -   k_AB* [A] * [B] + k_CD* [C] * [D]
d [C]/(d t) = d [D]/(d t)=   k_AB* [A] * [B] - k_CD * [C] * [D]

Chemical Equilibrium:

● Forward rate = backward rate, i.e.,    - k_AB* [A] * [B] + k_CD* [C] * [D] = 0

● Characterizing equilirium solutions: equilibrium constant ≡ [C]   [D]/[A]   [B] = k_AB/k_CD

In[581]:=

EquilCondition = (-KAB * Conc[A] * Conc[B] + KCD * Conc[C] * Conc[D] == 0) ;

ConcSoln = Solve[EquilCondition, Conc[C]]//Flatten

Conc[C] * Conc[D]/(Conc[A] * Conc[B])     /.ConcSoln

Out[582]=

{Conc[C] → (KAB Conc[A] Conc[B])/(KCD Conc[D])}

Out[583]=

KAB/KCD

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

Michaelis-Menten Kinetics

Conversion of substrate (S) to product (P), via enzyme-substrate (EnzS) formation:

        Underoverscript[Enz + S⇔EnzS, k_1^b, arg3] ; <br />Overscript[EnzS→P + Enz, k_2^f] ;

Get the equations using the Cellerator^TM package (automatic equation generation for biological modelling; www.cellerator.org)

In[1]:=

<<cellerator.m ;

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

In[2]:=

reaction1 = {{Enz + S⇔EnzS, k1f, k1b}} ;

(reaction1ODE = interpret[reaction1]//First)//ColumnForm

Out[3]=

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

In[4]:=

reaction2 = {{EnzS→P + Enz, k2f}} ;

(reaction2ODE = interpret[reaction2]//First)//ColumnForm

Out[5]=

Enz^′[t] == k2f EnzS[t]
EnzS^′[t] == -k2f EnzS[t]
P^′[t] == k2f EnzS[t]

 ◂ ▸  1 of 1

Michaelis-Menten Kinetics

Let's check conservation relations:

In[6]:=

EnzTotal[t] = Enz[t] + EnzS[t] ;

TimeDerivEnzTotal = D[EnzTotal[t], t] ;

reaction2Rule = reaction2ODE/.{Equal→ Rule}

Out[8]=

{Enz^′[t] →k2f EnzS[t], EnzS^′[t] → -k2f EnzS[t], P^′[t] →k2f EnzS[t]}

In[9]:=

TimeDerivEnzTotal/. reaction2Rule

Out[9]=

0

In[10]:=

STotal[t] = S[t] + EnzS[t] + P[t] ;

TimeDerivSTotal = D[STotal[t], t] ;

reactionsCombined = ({reaction1//First, reaction2//First}//interpret//First) /.{Equal→ Rule} ;

(TimeDerivSTotal/.reactionsCombined) == 0

Out[13]=

True

4 species, 2 conservation relations, therefore obtain 2 independent ODES:

In[14]:=

Clear[EnzTotal] ;

EnzTotal[t_] := EnzTotal[0] ;

EnzReplaceRule = Solve[EnzTotal[t] == Enz[t] + EnzS[t], Enz[t]]//Flatten ;

SelectedVar = {EnzS^′[t], S^′[t]} ;

TwoDimEnzymeODE = Pick[reactionsCombined,   Map[MemberQ[SelectedVar, #] & , Map[#[[1]] &, reactionsCombined]]]//.EnzReplaceRule/.{Rule→ Equal} ;

%//ColumnForm

Out[19]=

EnzS^′[t] == -k1b EnzS[t] - k2f EnzS[t] + k1f (-EnzS[t] + EnzTotal[0]) S[t]
S^′[t] == k1b EnzS[t] - k1f (-EnzS[t] + EnzTotal[0]) S[t]

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

Michaelis-Menten Kinetics

In[20]:=

STotal[t_] := STotal[0] ;

TransformEqn = {Overscript[S,^][t] == S[t]/STotal[0], Overscript[EnzS,^][t] == EnzS[t]/EnzTotal[0]} ;

TransformVariables = Solve[%, SelectedVar/.{Derivative[1][Species_][t_] → Species[t]}]//Flatten ;

ODETransformRule = Append[%, Map[ D[#[[1]], t] → D[#[[2]], t] &, %]]//Flatten

Out[23]=

In[24]:=

ModTwoDimEnzymeODE = Simplify[TwoDimEnzymeODE/.%, Assumptions→ {STotal[0] > 0, EnzTotal[0] >0}] ;

%//ColumnForm

Out[25]=

Overscript[EnzS,^][t] (k1b + k2f + k1f STotal[0] Overscript[S,^][t]) + Overscript[EnzS,^]^′[t] == k1f STotal[0] Overscript[S,^][t]
EnzTotal[0] Overscript[EnzS,^][t] (k1b + k1f STotal[0] Overscript[S,^][t]) == STotal[0] (k1f EnzTotal[0] Overscript[S,^][t] + Overscript[S,^]^′[t])

In[26]:=

TimeChainRule = {t→Overscript[t,^] , Derivative[1][Species_][t] → Derivative[1][Species][Overscript[t,^]] * (EnzTotal[0] * k1f)} ;

ModTwoDimEnzymeODE = Simplify[ModTwoDimEnzymeODE/.TimeChainRule, Assumptions→ {STotal[0] > 0, EnzTotal[0] >0}]//. {x_[temp_] → x[Overscript[t,^]]}

Out[27]=

In[28]:=

ParamEnzymeODE = Simplify[ModTwoDimEnzymeODE/.{EnzTotal[0] → ε  STotal[0]}, Assumptions→ {ε> 0, k1f> 0}] ;

%//ColumnForm

Out[29]=

Overscript[EnzS,^][Overscript[t,^]] (k1b + k1f STotal[0] Overscript[S,^][Overscript[t,^]]) == k1f STotal[0] (Overscript[S,^][Overscript[t,^]] + Overscript[S,^]^′[Overscript[t,^]])

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

Michaelis-Menten Kinetics: Asymptotic Analysis

In[30]:=

ZeroOrderApprox = ParamEnzymeODE/.{ε→0}

Out[30]=

In[31]:=

EnzSEquilibrium = Solve[%[[1]], Overscript[EnzS,^][Overscript[t,^]] ]//Flatten

Out[31]=

{Overscript[EnzS,^][Overscript[t,^]] → (k1f STotal[0] Overscript[S,^][Overscript[t,^]])/(k1b + k2f + k1f STotal[0] Overscript[S,^][Overscript[t,^]])}

In[32]:=

ReducedODE = ZeroOrderApprox[[2]]/.EnzSEquilibrium

Out[32]=

In[33]:=

MMKinetic = ReducedODE//Simplify[#, Assumptions→ {Overscript[S,^][Overscript[t,^]] > 0}] & //Solve[#, Overscript[S,^]^′[Overscript[t,^]]] &//Flatten

Out[33]=

{Overscript[S,^]^′[Overscript[t,^]] → -(k2f Overscript[S,^][Overscript[t,^]])/(k1b + k2f + k1f STotal[0] Overscript[S,^][Overscript[t,^]])}

In[34]:=

(%[[1, 2]]//Numerator)/((%[[1, 2]]//Denominator) / (k1f  STotal[0])//Expand)

Out[34]=

-(k2f Overscript[S,^][Overscript[t,^]])/(k1b/(k1f STotal[0]) + k2f/(k1f STotal[0]) + Overscript[S,^][Overscript[t,^]])

In[35]:=

Overscript[S,^]^′[Overscript[t,^]]/.MMKinetic/.{k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1} ;

[Graphics:HTMLFiles/Lecture2_211.gif]

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

Michaelis-Menten Kinetics: Asymptotic Analysis

Out[441]=

-(k2f Overscript[S,^][Overscript[t,^]])/(k1b/(k1f STotal[0]) + k2f/(k1f STotal[0]) + Overscript[S,^][Overscript[t,^]])

Changing the Michaelis constant,  k1b/(k1f STotal[0])+k2f/(k1f STotal[0])

In[519]:=

Overscript[S,^]^′[Overscript[t,^]]/.MMKinetic/.{k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1} ;

Overscript[S,^]^′[Overscript[t,^]]/.MMKinetic/.{k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 2} ;

Show[GraphicsArray[{{Plot1, Plot2}}], ImageSize-> {600, 200}, DisplayFunction-> $DisplayFunction] ;

[Graphics:HTMLFiles/Lecture2_225.gif]

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

Michaelis-Menten Kinetics: Asymptotic Analysis

What happens to the solution in the limit  ε → 0?

In[37]:=

KineticParamReplaceList = {k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1, ε→10^(-1)} ;

Sol = NDSolve[Append[ParamEnzymeODE/.%, {Overscript[S,^][0] == 1, Overscript[EnzS,^][0] == 0}], {Overscript[S,^], Overscript[EnzS,^]}, {Overscript[t,^], 0, 10}] ;

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

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

Show[GraphicsArray[{{EnzymeSubstratePlotLarge, EnzymeSubstratePlotExpand}, {SubstratePlotLarge, SubstratePlotExpand}}], ImageSize→ {600, 300}] ;

[Graphics:HTMLFiles/Lecture2_240.gif]

 ◂ ▸  1 of 1

Michaelis-Menten Kinetics: Asymptotic Analysis

What happens to the solution in the limit ε→0?

In[88]:=

KineticParamReplaceList = {k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1, ε→10^(-2)} ;

Sol = NDSolve[Append[ParamEnzymeODE/.%, {Overscript[S,^][0] == 1, Overscript[EnzS,^][0] == 0}], {Overscript[S,^], Overscript[EnzS,^]}, {Overscript[t,^], 0, 10}] ;

Show[GraphicsArray[{{EnzymeSubstratePlotLarge, EnzymeSubstratePlotExpand}, {SubstratePlotLarge, SubstratePlotExpand}}], ImageSize→ {600, 300}] ;

[Graphics:HTMLFiles/Lecture2_253.gif]

 ◂ ▸  1 of 1

Michaelis-Menten Kinetics: Asymptotic Analysis

What happens to the solution in the limit ε→0? Need to use singular perturbation analysis, to match outer and boundary layer solution

In[95]:=

KineticParamReplaceList = {k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1, ε→10^(-3)} ;

Sol = NDSolve[Append[ParamEnzymeODE/.%, {Overscript[S,^][0] == 1, Overscript[EnzS,^][0] == 0}], {Overscript[S,^], Overscript[EnzS,^]}, {Overscript[t,^], 0, 10}] ;

Show[GraphicsArray[{{EnzymeSubstratePlotLarge, EnzymeSubstratePlotExpand}, {SubstratePlotLarge, SubstratePlotExpand}}], ImageSize→ {600, 300}] ;

[Graphics:HTMLFiles/Lecture2_266.gif]

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

Regular Perturbation Analysis

Singular perturbation involves more complex ideas; start off with regular perturbation analysis.
Suppose we want to solve the ODE:

In[1]:=

ODESys = {x^″[t] + ε x^′[t]^2 + x[t] == 0, x[0] == 1, x^′[0] == 0} ;

where we know that ε is very small (0 < ε ≪ 1). Can we use the smallness of ε to obtain approximate solutions?
Since ε enters the equation as a small parameter, consider:

In[2]:=

MaxOrder = 2 ;    xtest[t_] :=    Underoverscript[∑, i = 0, arg3] ε^i x[i][t] ;

In[10]:=

ExpansionODE = ODESys/.{x[t_] → xtest[t], Derivative[n_][x][t_] →  Derivative[n][xtest][t]} ; %//ColumnForm

Out[10]=

x[0][0] + ε x[1][0] + ε^2 x[2][0] == 1
x[0]^′[0] + ε x[1]^′[0] + ε^2 x[2]^′[0] == 0

Now lets look at equation for each order of ε:

In[125]:=

εOrder = 0 ;

Order0ODE = Thread[ Equal[ Coefficient[ Map[#//First &,   ExpansionODE], ε, εOrder] , Map[#//Last&,   ExpansionODE] ]   ]

Out[126]=

{x[0][t] + x[0]^′′[t] == 0, x[0][0] == 1, x[0]^′[0] == 0}

In[127]:=

εOrder = 1 ;

Order1ODE = Thread[ Equal[ Coefficient[ Map[#//First &,   ExpansionODE], ε, εOrder] , Map[#//Last&,   ExpansionODE] ]   ]

Out[128]=

{x[1][t] + x[0]^′[t]^2 + x[1]^′′[t] == 0, x[1][0] == 1, x[1]^′[0] == 0}

In[129]:=

εOrder = 2 ;

Order2ODE = Thread[ Equal[ Coefficient[ Map[#//First &,   ExpansionODE], ε, εOrder] , Map[#//Last&,   ExpansionODE] ]   ]

Out[130]=

{x[2][t] + 2 x[0]^′[t] x[1]^′[t] + x[2]^′′[t] == 0, x[2][0] == 1, x[2]^′[0] == 0}

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

Regular Perturbation Analysis

Now we can solve the expansion equations at each order:

In[131]:=

Order0Soln = NDSolve[Order0ODE, x[0], {t, 0, 1}] ;

Plot[x[0][t]/.Order0Soln, {t, 0, 1}, PlotStyle→ {Hue[0.8], Thickness[0.01]}, Frame→ True, PlotRange→All, Axes→ None] ;

[Graphics:HTMLFiles/Lecture2_294.gif]

In[133]:=

Order1Soln = NDSolve[Order1ODE/.Order0Soln, x[1], {t, 0, 1}] ;

Plot[x[1][t]/.Order1Soln, {t, 0, 1}, PlotStyle→ {Hue[0.8], Thickness[0.01]}, Frame→ True, PlotRange→ All, Axes→ None] ;

[Graphics:HTMLFiles/Lecture2_297.gif]

In[135]:=

Order2Soln = NDSolve[Order2ODE/.Order0Soln /.Order1Soln, x[2], {t, 0, 1}] ;

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

Regular Perturbation Analysis

Lets compare the expansion solution with the solution for the original ODE system:

In[136]:=

ODESoln = NDSolve[ODESys/.{ε-> 0.1}, x, {t, 0, 1}] ;

ODESolnPlot = Plot[x[t]/.ODESoln, {t, 0, 1}, PlotStyle→ {Hue[0.3], Thickness[0.01]}, Frame→ True, PlotRange→All, DisplayFunction→ Identity, Axes→ None] ;

Show[{ODESolnPlot, ExpansionODESolnPlot}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/Lecture2_308.gif]

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

Regular Perturbation Analysis

Lets compare the expansion solution with the solution for the original ODE system:

In[140]:=

ODESoln = NDSolve[ODESys/.{ε-> 0.01}, x, {t, 0, 1}] ;

ODESolnPlot = Plot[x[t]/.ODESoln, {t, 0, 1}, PlotStyle→ {Hue[0.3], Thickness[0.01]}, Frame→ True, PlotRange→All, DisplayFunction→ Identity, Axes→ None] ;

Show[{ODESolnPlot, ExpansionODESolnPlot}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/Lecture2_318.gif]

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

Regular Perturbation Analysis

Lets compare the expansion solution with the solution for the original ODE system:

In[144]:=

ODESoln = NDSolve[ODESys/.{ε-> 0.001}, x, {t, 0, 1}] ;

ODESolnPlot = Plot[x[t]/.ODESoln, {t, 0, 1}, PlotStyle→ {Hue[0.3], Thickness[0.01]}, Frame→ True, PlotRange→All, DisplayFunction→ Identity, Axes→ None] ;

Show[{ODESolnPlot, ExpansionODESolnPlot}, DisplayFunction→$DisplayFunction] ;

[Graphics:HTMLFiles/Lecture2_328.gif]

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

Michaelis-Menten Kinetics: Asymptotic Analysis

In[1073]:=

ODESys = Append[ParamEnzymeODE/.{k1b→ 1, k1f→ 1, k2f→ 1, STotal[0] → 1}, {Overscript[S,^][0] == 1, Overscript[EnzS,^][0] == 0}]//Flatten

Out[1073]=

Read in the PerturbationODE package, that does the perturbation analysis automatically

In[1075]:=

<<PerturbationODE.m

In[1099]:=

OrderList = PolyOrderList[ODESys, {Overscript[S,^][Overscript[t,^]], Overscript[EnzS,^][Overscript[t,^]]}, ε, 1] ;

(* ε : zeroth order equation *)

OrderList[[1]]//ColumnForm

Out[1100]=

2 Overscript[EnzS,^][0][Overscript[t,^]] + Overscript[EnzS,^][0][Overscript[t,^]] Overscript[S,^][0][Overscript[t,^]] == Overscript[S,^][0][Overscript[t,^]]
Overscript[S,^][0][0] == 1
Overscript[EnzS,^][0][0] == 0

In[1101]:=

(* ε : first order equation *)OrderList[[2]]//ColumnForm

Out[1101]=

Overscript[S,^][1][0] == 0
Overscript[EnzS,^][1][0] == 0

Conclusion: the expansion equations are not well-formed, ODE systems!
In fact, need to do asymptotic matching

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

Singular Perturbation

-Construct inner, or singular solution: boundary layer
-Construct outer, or quasi-steady state solution

-Matching condition: solution is continuous going from inner to outer solutions
-Obtain the composite solution, which is uniformly valid approximation for all time

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