(************** Content-type: application/mathematica ************** CreatedBy='Mathematica 5.2' Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 22965, 691]*) (*NotebookOutlinePosition[ 23692, 717]*) (* CellTagsIndexPosition[ 23648, 713]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell[TextData[{ StyleBox["Math. Model & Sci. Comput. in Biosciences", FontSize->24, FontSlant->"Italic"], StyleBox["\nExercise 1: ", FontSlant->"Italic"], "Singular Perturbation Analysis" }], "Title"], Cell[TextData[{ StyleBox["Instruction: complete and hand-in electronically by ", FontSize->18], StyleBox["Friday, May 18", FontSize->18, FontWeight->"Bold"], StyleBox[".\n\nComplete the ", FontSize->18], StyleBox["areas high-lighted as such", FontSize->18, Background->RGBColor[0, 1, 1]], StyleBox[". In particular, \"XXX\" denote fields that need to be filled in. \ Discussions are allowed, but each person must complete on their own.", FontSize->18] }], "Subtitle"], Cell[CellGroupData[{ Cell["Numerical Exploration", "Section"], Cell["\<\ Here we will extend the singular perturbation analysis that was \ carried out in Lecture3.nb. In the lecture, singular perturbation analysis \ for the Michaelis-Menten kinetics was analyzed to order 0; we will obtain \ equations up to order 1. Remember, in Lecture 2 we created the ODE system using the Cellerator package \ and simplified to the following 2-dimensional ODE system: \ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ StyleBox["SimpEnzymeODE", FontColor->RGBColor[0, 0, 1]], StyleBox["=", FontSize->10], RowBox[{ StyleBox["{", FontSize->12], RowBox[{ RowBox[{ StyleBox[ RowBox[{ StyleBox["\[Epsilon]", FontSize->12, FontColor->RGBColor[1, 0, 0]], StyleBox[" ", FontSize->12], StyleBox[ RowBox[{ SuperscriptBox[\(EnzS\&^\), "\[Prime]", MultilineFunction->None], "[", \(t\&^\), "]"}], FontSize->12]}], FontSize->14], StyleBox[" ", FontSize->12], StyleBox["==", FontSize->12], StyleBox[\(S\&^[ t\&^] - \((\((k1b + k2f)\)/\((k1f\ STotal[0])\) + S\&^[t\&^])\) EnzS\&^[t\&^]\), FontSize->12]}], StyleBox[",", FontSize->12], StyleBox[ " \t\t\t \ ", FontSize->12], StyleBox[ RowBox[{ RowBox[{ SuperscriptBox[\(S\&^\), "\[Prime]", MultilineFunction->None], "[", \(t\&^\), "]"}], "==", " ", \(\(-S\&^[t\&^]\) + \ EnzS\&^[ t\&^]\ \((k1b/\((k1f\ STotal[0])\) + \ S\&^[t\&^])\)\)}], FontSize->12]}], StyleBox["}", FontSize->12]}]}], ";"}]], "Input", FontSize->14], Cell[TextData[{ "where the \"^\" denotes dimensionless, normalized values, ", StyleBox["STotal[0]", "InlineInput"], " the initial total substrate concentration, ", StyleBox["\[Epsilon]", "InlineInput", FontColor->RGBColor[1, 0, 0]], StyleBox[" = EnzTotal[0]/STotal[0]", "InlineInput"], " a small dimensionless parameter describing the small, enzyme-to-substrate \ concentration ratio." }], "Text"], Cell[TextData[{ "For convenience, let's remove the \"^\" notation in the equations. ", StyleBox["Hint", FontSlant->"Italic"], ": apply the replacement rule, ", StyleBox["OverHat[x_]\[Rule] x", "InlineInput"], ", to ", StyleBox["SimpEnzymeODE", "InlineInput", FontColor->RGBColor[0, 0, 1]], StyleBox[".", FontColor->RGBColor[0, 0, 1]] }], "Text"], Cell[BoxData[ RowBox[{ StyleBox["SimpEnzymeODE", FontColor->RGBColor[0, 0, 1]], "=", RowBox[{ StyleBox["SimpEnzymeODE", FontColor->RGBColor[0, 0, 1]], "/.", \({\ \ XXX\ \ \ }\)}]}]], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "\nTypically, prior to doing mathematical analysis, it is a good idea to \ numerically explore the solution behavior.\n\nLet's look at the list of \ variables and functions that appear in the ODE system. ", StyleBox["Hint", FontSlant->"Italic"], ": use ", StyleBox["Cases", "InlineInput"], ", with selection criterion ", StyleBox["_Symbol|_Symbol[x_]", "InlineInput"], ", and apply ", StyleBox["Union ", "InlineInput"], "as demonstrated in Tutorial1.nb. \n" }], "Text"], Cell[BoxData[ \(Cases[SimpEnzymeODE, \ \ \ \ \ \ , \ \ Infinity\ ]\)], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "\nTo numerically solve the ODE, we need to provide initial condition as \ well as choose parameter values. \nLet's append the ODE list with the initial \ condition, ", StyleBox["EnzS[0]==0, S[0]==1", "InlineInput"], ".\n" }], "Text"], Cell[BoxData[{\(IC\ = \ {\ \ \ };\), "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"ODEIC", "=", RowBox[{"Join", "[", StyleBox[\(\(SimpEnzymeODE\)\(,\)\), FontColor->RGBColor[0, 0, 1]], StyleBox[" ", FontColor->RGBColor[0, 0, 1]], "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", \( (*\ just\ for\ display, \ in\ column\ form\ *) \)}], "\[IndentingNewLine]", \(% // ColumnForm\)}], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "\nFurthermore, to numerically solve ODES we need to provide numerical \ values of parameters. \nLets set ", StyleBox["k1b->1, k1f->1, k2f->1, STotal[0]->1, \[Epsilon]->0.1", "InlineInput"], ".\n" }], "Text"], Cell[BoxData[{ \(\(\(\(ParamReplaceRule\)\(=\)\)\ ;\)\[IndentingNewLine]\ \[IndentingNewLine] (*\ give\ the\ "\"\ ODE\ and\ initial\ condition\ system\ as\ \ NumODEIC\ *) \), "\[IndentingNewLine]", \(\(\(NumODEIC\)\(=\)\(\ \)\)\)}], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "When calling numerical ODE integrator ", StyleBox["NDSolve", "InlineInput"], ", one needs to specify the functions being solved for. We can extract this \ using again ", StyleBox["Cases", "InlineInput", Background->None], " command:" }], "Text"], Cell[BoxData[ \(SolvedFunctions = Cases[\(NumODEIC\)\(,\)\ \ \ \ \ \ \ \ ] // Union\)], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Now call ", StyleBox["NDSolve", "InlineInput"], " and ", StyleBox["Plot", "InlineInput"], " the result, as was done in Tutorial1." }], "Text"], Cell[BoxData[{ \(\(ODESoln = NDSolve[\ \ \ \ \(,\)\(\ \)\({t, 0, 1}\)];\)\[IndentingNewLine]\), "\[IndentingNewLine]", \(\(Plot[SolvedFunctions[\([1]\)] /. \ \ XXX\ \ , \ {t, 0, 1}, \ AxesLabel \[Rule] \ \((ToString /@ SolvedFunctions[\([1]\)])\)];\)\), "\[IndentingNewLine]", \(\(Plot[SolvedFunctions[\([2]\)] /. \ \ XXX, \ \ \ {t, 0, 1}, \ AxesLabel \[Rule] \ \((ToString /@ SolvedFunctions[\([2]\)])\)];\)\)}], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Now it seems ", StyleBox["EnzS[t]", "InlineInput"], " reaches a maximum at some t>0. Lets find out when it does. \n\nWe can do \ this by generating a uniformly distributed values of times from t=0 to t=1, \ evaluate ", StyleBox["EnzS[t]", "InlineInput"], " at this list of times, and find its maximum.\nFirst, generate the list of \ times using ", StyleBox["Table", "InlineInput"] }], "Text"], Cell[BoxData[{ \(\(TimeInterval = 10000;\)\), "\[IndentingNewLine]", \(\(timeList = Table[t \[Rule] \ XXX, \ {i, 0, TimeInterval}];\)\[IndentingNewLine]\[IndentingNewLine] (*\ write\ out\ the\ short\ form\ of\ the\ output\ with\ Short\ *) \), "\ \[IndentingNewLine]", \(% // Short\)}], "Input", Background->RGBColor[0, 1, 1]], Cell["\<\ Then, from ODESoln, obtain the EnzSSoln in the form of an \ interpolated function,\ \>", "Text"], Cell[BoxData[ \(EnzSSoln = XXX /. ODESoln\)], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Now evaluate by the numerical values of EnzSSoln at each of the time \ points in timeList. ", StyleBox["Hint", FontSlant->"Italic"], ": use commands ", StyleBox["Map", "InlineInput"], " and the fact that ", StyleBox["timeList", "InlineInput"], " is a list of replacement rules. Remember that one way to use ", StyleBox["Map", "InlineInput"], " is to create a ", StyleBox["pure function", FontSlant->"Italic"], ":" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(\( (*\ simple\ example\ on\ using\ Map, \ with\ a\ pure\ function\ defined\ by\ "\<#^5 &\>"\ \ *) \)\(\ \[IndentingNewLine]\)\( (*\ "\<#\>"\ shows\ where\ the\ input\ of\ the\ \ function\ goes, \ and\ "\<&\>"\ simply\ ends\ the\ function\ definition\ *) \)\(\ \[IndentingNewLine]\)\(\[IndentingNewLine]\)\(Map[#^5\ &, {a, b, c}]\)\)\)], "Input"], Cell[BoxData[ \({a\^5, b\^5, c\^5}\)], "Output"] }, Open ]], Cell[BoxData[{ \(\(\(\(NumList\)\(=\)\)\ \ \ \ \ \ \ ;\)\), "\[IndentingNewLine]", \(% // Short\)}], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "From this list of numerical values, find out the position of the maximum \ within the list using", StyleBox[" Ordering[XXX,-1]", "InlineInput"], ". \n\nThen get the time at the maximum, as well as ", StyleBox["EnzSSoln", "InlineInput"], " at maximum." }], "Text"], Cell[BoxData[{ \(\(\(MaxEnzSLocation = Ordering[NumList, \ \(-1\)];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(TimeOfMaxEnzS = timeList[\([MaxEnzSLocation]\)]\), "\[IndentingNewLine]", \(MaxEnzS = EnzSSoln /. %\)}], "Input", Background->RGBColor[0, 1, 1]], Cell[TextData[{ "Now, we can combine all the previous steps to write a routine (i.e., \ function) that takes as an input the \[Epsilon] value and does the above \ computation automatically.\n\nTo write this routine, we wrap everything \ within a ", StyleBox["Module", "InlineInput"], "; use the template below. After defining the function ", StyleBox["PlotPeakOfEnzS", "InlineInput"], ", we can then apply to the list of epsilon values ", StyleBox["EpsValList={0.1,0.01, 0.001, 0.0001}", "InlineInput"], ". We'd like to test: is Time_Of_Peak_Of_EnzS(\[Alpha]) = ", Cell[BoxData[ \(TraditionalForm\`\[Epsilon]\^\[Alpha]\)]], " for some value \[Alpha]?\n\nTo test this, we do a ", StyleBox["LogLogListPlot", "InlineInput"], " of the peak-times vs. \[Epsilon] values after computing the results." }], "Text"], Cell[BoxData[{ \(\(PlotPeakOfEnzS[EpsVal_] := Module[\({TimeInterval = 10^4}\)\(,\)\ \[IndentingNewLine]\[IndentingNewLine]];\)\ \[IndentingNewLine]\[IndentingNewLine]\), "\[IndentingNewLine]", \(\(EpsValList = {0.1, 0.01, \ 0.001, \ 0.0001};\)\), "\[IndentingNewLine]", \(\(\(PeakLocationList = \(N[t /. #\ , 10] &\)\ /@ \ \(PlotPeakOfEnzS\ /@ \ EpsValList\) // Flatten\ \ ;\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(<< Graphics`Graphics`\), "\[IndentingNewLine]", \(\(LogLogListPlot[\ \ \ \ \ \ XXX\ \ \ , \ PlotJoined \[Rule] \ True\ , \ AspectRatio \[Rule] \ Automatic, \ AxesLabel \[Rule] \ {"\<\[Epsilon]\>", \ "\