| sol = NDSolve[{ |
| x'[t]==-x[t]+x[t]*y[t], |
| y'[t]==y[t]-x[t]*y[t], |
| x[0]==0.65, y[0]==0.65}, {x,y},{t,0,10}] |
| g1 = ParametricPlot[Evaluate[{x[t], y[t]}/.sol],{t,0,10}] |
| g2 = VectorPlot[{-x + x*y, y - x*y}, {x,0,2}, {y,0,2}] |
| Show[{g1, g2}, PlotRange ->{{0,2}, {0,2}}] |
| B.8 Predator-prey models. |
| First, generate the solution. |
| a = 1; b = 1; c = 1; d = 1; xinit = 0.65; yinit = 0.65; tmax = 10; |
| sol = NDSolve[{ |
| x'[t]==-a*x[t] + b*x[t]*y[t], |
| y'[t]==c*y[t] - d*x[t]*y[t], | |
| x[0]==xinit, y[0]==yinit}, {x,y},{t,0,tmax}] |
|
| Now plot the x(t) and y(t) curves, |
| Plot[{Evaluate[x[t]/.sol],Evaluate[y[t]/.sol]},{t,0,tmax}] |
|
| and plot the trajectory (x(t), y(t))). |
| ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tmax}] |
|
| B.9 The Fitzhugh-Nagumo equations. |
| The square pulse function is defined by |
| z[t_]:=If[t < 10,0,If[t < 20,1,0]] |
|
| To solve the Fitzhugh-Nagumo equations, use this code. Modify the parameter values as you wish. |
| a = 0.8; b = 0.35; c = 1; tmax = 30; |
| sol = NDSolve[{ |
| x'[t]==c*(y[t] + x[t] - x[t]^3/3 + z[t]), |
| y'[t]==(-x[t] + a - b*y[t])/c, |
| x[0]==0.6, y[0]==0.6}, {x,y}, {t,0,tmax}, MaxSteps->1000] |
|
| To plot the population curves in a single graph use this code. |
| Plot[{Evaluate[x[t]/.sol],Evaluate[y[t]/.sol]},{t,0,tmax},PlotRange->All] |
|
| and to plot the trajectory use |
| ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tmax},PlotRange->All] |
|
| B.10 Eigenvalues and eigenvectors. |
| To find the eigenvalues and eigenvectors of the matrix m, entered in this example as a list of lists,
use these commands. Enter your own matrix, of coure. |
| m = {{1, 1},{4,1}}; |
| Eigenvalues[m] |
| Eigenvectors[m] |
|
Here we'll present Mathematica code from Appendix B of Mathematical Models in the Biosciences, 2 by
Michael Frame, Yale University Press, 2021. Explanations for the code, and what the code is written to do,
can be found in that volume.
| B.11 Eigenvalues in higher dimensions. |
| We can find eigenvalues of any square matrix. But while for 2 × 2 matrices the
eigenvalues are the solutions of a quadratic equation, for larger matrices the exact expresion for the
eigenvalues may be far too complicated, or impossible to find. If the output of the Eigenvalues[m] command is
too complex, the N[Eigenvalues[m]] command gives numerical approximations. |
| m = {{ 1, 0, 2 },{ 0, -1, 1 },{0, -1, 2}}; |
| Eigenvalues[m] |
| N[Eigenvalues[m]] |
|
| B.12 SIR calculations. |
| First, generate the solution with this code. |
| b=2; e=3; d=1.0; g=0.5; n=0.75; tmax=45; |
| sol = NDSolve[{ |
| x'[t]==b*x[t] - e*x[t]*z[t] - d*x[t] + g*z[t], |
| y'[t]==e*x[t]*z[t] - n*y[t] - d*y[t], |
| z'[t]==n*y[t] - d*z[t] - g*z[t], |
| x[0]==0.5,y[0]==0.25,z[0]==0.25},{x,y,z},{t,0,tmax}] |
|
| To plot the three population curves on the same graph, use |
| Plot[{Evaluate[x[t] /. sol], Evaluate[y[t] /. sol], Evaluate[z[t] /. sol]}, {t, 0, tmax}] |
|
| and to plot the trajectory use |
| ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tmax}] |
|
| B.13 Michaelis-Menten calculations. |
| We generate the solution with this code |
| a = 0.2; b = 0.03; c = 0.05; tmax = 100; |
| sol=NDSolve[{ |
| x'[t]==-a*x[t]*y[t] + b*z[t], |
| y'[t]==-a*x[t]*y[t] + (b+c)*z[t], |
| z'[t]==a*x[t]*y[t] - (b+c)*z[t],w'[t]==c*z[t], |
| x[0]==1,y[0]==.5,z[0]==0.0,w[0]==0},{x,y,z,w},{t,0,tmax}] |
|
| To plot the four population curves individually, use |
| Plot[Evaluate[x[t] /. sol], {t, 0, tmax}] |
| Plot[Evaluate[y[t] /. sol], {t, 0, tmax}] |
| Plot[Evaluate[z[t] /. sol], {t, 0, tmax}] |
| Plot[Evaluate[w[t] /. sol], {t, 0, tmax}] |
|
| To plot the four population curves on the same graph, use |
| Plot[{Evaluate[w[t] /. sol], Evaluate[x[t] /. sol], Evaluate[y[t] /. sol], Evaluate[z[t] /. sol]}, {t, 0, tmax}] |
|
| B.14 Virus dynamics. |
| We generate the solution with this code. |
| l=100000; d=0.1; a=0.5; b=.0000002; k=100; u=5; tmax=20; |
| sol = NDSolve[{ |
| x'[t] == l - d*x[t] - b*x[t]*z[t], |
| y'[t] == b*x[t]*z[t] - a*y[t], |
| z'[t] == k*y[t] - u*z[t], |
| w'[t] == -u*w[t], |
| x[0] == 1000, y[0] == 1000, z[0] == 0, w[0] == 5000000}, {x, y, z, w}, {t, 0, tmax}] |
|
| To plot the w and z curves in the same graph, use this code. |
| Plot[{Evaluate[z[t] /. sol], Evaluate[w[t] /. sol]}, {t, 0, tmax}] |
|
| B.15 Immune system dynamics. |
| We generate the solution with this code. |
| b=100000; dx=0.1; g=0.0000002; dy=0.5; e=0.002; k=80; dz=5; h=0.002; dw=50; tmax=10; |
| sol = NDSolve[ |
| {x'[t] == b - dx*x[t]- g*x[t]*z[t], |
| y'[t] == g*x[t]*z[t] - dy*y[t] - e*y[t]*w[t], |
| z'[t] == k*y[t] - dz*z[t], |
| w'[t] == h*y[t]*w[t] - dw*w[t], |
| x[0] == 1000000, y[0] == 0, z[0] == 2000, w[0] == 0.005}, |
| {x, y, z, w}, {t, 0, tmax},MaxSteps->2000] |
|
| To plot the y and w curves in the same graph, use this code. |
| Plot[{Evaluate[y[t] /. sol], Evaluate[w[t] /. sol],},{t, 0, tmax}] |
|
| B.16 The Hodgkin-Huxley equations. |
| We'll plot trajectories, find the nullclines, and locate the fixed points. First, generate
two solutions with slightly different initial points. |
| tmin=12; tmax=40; |
| c=1.9; vna=-115; vk=12; vl=-10.6; cgna=120; cgk=36; cgl=0.3; |
| i[t_]:=Sin[Pi*t] + Sin[Pi*t/Sqrt[2]]; |
| gna[m_,h_]:=cgna*(m^3)*h; |
| gk[n_]:=cgk*n^4; |
| an[v_]:=0.01*(v + 10)/(Exp[(v + 10)/10] - 1); |
| bn[v_]:=0.125*Exp[v/80]; |
| am[v_]:=0.1*(v + 25)/(Exp[(v + 25)/10] - 1); |
| bm[v_]:=4*Exp[v/18]; |
| ah[v_]:=0.07*Exp[v/20]; |
| bh[v_]:=1/(Exp[(v + 30)/10] + 1); |
| sol1 = NDSolve[{ |
| n'[t]==an[v[t]]*(1 - n[t]) - bn[v[t]]*n[t], |
| m'[t]==am[v[t]]*(1 - m[t]) - bm[v[t]]*m[t], |
| h'[t]==ah[v[t]]*(1 - h[t]) - bh[v[t]]*h[t], |
| v'[t]==(-1/c)*(gna[m[t],h[t]]*(v[t] - vna) + gk[n[t]]*(v[t] - vk) + cgl*(v[t] - vl) + i[t]), |
| n[0]==0.577, m[0]==0.046, h[0]==0.874, v[0]==0.78}, |
| {n, m, h, v}, {t, 0, tmax},MaxSteps->6000] |
| sol2 = NDSolve[{ |
| n'[t]==an[v[t]]*(1 - n[t]) - bn[v[t]]*n[t], |
| m'[t]==am[v[t]]*(1 - m[t]) - bm[v[t]]*m[t], |
| h'[t]==ah[v[t]]*(1 - h[t]) - bh[v[t]]*h[t], |
| v'[t]==(-1/c)*(gna[m[t],h[t]]*(v[t] - vna) + gk[n[t]]*(v[t] - vk) + cgl*(v[t] - vl) + i[t]), |
| n[0]==0.577, m[0]==0.045, h[0]==0.874, v[0]==0.78}, |
| {n, m, h, v}, {t, 0, tmax},MaxSteps->6000] |
|
| To plot the (m(t), V(t)) trajectories for both solutions, use |
| ParametricPlot[{Evaluate[{m[t], v[t]} /. sol1],Evaluate[{m[t], v[t]} /. sol2]}, |
| {t, tmin, tmax}, PlotRange -> All, PlotPoints -> 100, AspectRatio -> 1] |
|
| To plot the m and V nullclines, use |
| gk0 = 36; ninf = 0.318; vk = 12; gl = 0.3; vl = -10.6; gna = 120; hinf = 0.596; vna = -115; |
| am[v_]:=0.1*(v + 25)/(Exp[(v + 25)/10] - 1); |
| bm[v_]:=4*Exp[v/18]; |
| mncl[v_]:=am[v]/(am[v] + bm[v]); |
| vncl[v_]:=((-gk0*(ninf^4)*(v - vk) - gl*(v - vl))/(gna*hinf*(v - vna)))^(1/3); |
| Plot[{mncl[v],vncl[v]}, {v,-122,2}] |
|
| Code to generate periodic rectangular pulses |
| i[t_]:=If[t<1,1,If[t<2,0,i[t-2]]] |
|
| To find approximate coordinates of the fixed points, use |
| FindRoot[mncl[v]==vncl[v],{v,-100},MaxIterations->30] |
|
| B.17 Chaos in predator-prey models. |
| First, generate a solution with this code. |
| r1=1; r2=1; r3=-1; a11=0.001; a12=0.001; a13=0.01; a21=0.0016; |
| a22=0.001; a23=0.001; a31=-0.005; a32=-0.0005; a33=0; tmax=2000; |
| sol = NDSolve[{ |
| x'[t] == r1*x[t] - x[t]*(a11*x[t] + a12*y[t] + a13*z[t]), |
| y'[t] == r2*y[t] - y[t]*(a21*x[t] + a22*y[t] + a23*z[t]), |
| z'[t] == r3*z[t] - z[t]*(a31*x[t] + a32*y[t] + a33*z[t]), |
| x[0]==0.1, y[0]==0.1, z[0]==0.1},{x, y, z}, {t, 0, tmax},MaxSteps->20000] |
|
| To plot the trajectory in 3 dimensions, use this code. |
| ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tmax}, |
| PlotPoints->12000, BoxRatios->{1,1,1}, PlotRange->All] |
|
| B.18 The Duffing equation. |
| First we'll generate solutions for two slightly different initial points. |
| a=0.005; b=-0.5; c=0.5; d=0.1; w=1; tmax=100; |
| sol1 = NDSolve[{ |
| x'[t] == y[t], |
| y'[t] == -a*y[t] - b*x[t] - c*x[t]^3 + d*Cos[w*z[t]], |
| z'[t] == 1, |
| x[0]==1.0, y[0]==0, z[0]==0}, {x, y, z}, {t, 0, tmax}, MaxSteps -> 20000] |
| sol2 = NDSolve[{ |
| x'[t] == y[t], |
| y'[t] == -a*y[t] - b*x[t] - c*x[t]^3 + d*Cos[w*z[t]], |
| z'[t] == 1, |
| x[0]==1.1, y[0]==0, z[0]==0}, {x, y, z}, {t, 0, tmax}, MaxSteps -> 20000] |
|
| To plot the x(t) curves for both systems, use |
| Plot[{Evaluate[x[t] /. sol1], Evaluate[x[t] /. sol2]},{t, 0, tmax}] |
|
| To plot the differences of the x(t) curves, use |
| Plot[Evaluate[x[t] /. sol1] - Evaluate[x[t] /. sol2], {t, tmin, tmax}] |
|
| B.19 Control of chaos. |
| Here are codes to generate two figures from this section. |
| l[s_, x_] := s*x*(1 - x) |
| SeedRandom[1234]; s = 4; num = 200; del = 0.02; rd = 0.000001; |
| x = Random[Real, {0, 1}]; |
| xlst = {x}; |
| Do[{If[Abs[x - 0.25] < del, s = Solve[0.25 == t*x*(1 - x), t][[1, 1, 2]] + Random[Real, {-rd, rd}]], |
| x = l[s, x], AppendTo[xlst, x], s = 4, x = l[s, x], AppendTo[xlst, x]},{i, 1, num}] |
| ptlst = {}; |
| Do[AppendTo[ptlst, Point[{i, xlst[[i]]}]], {i, 1, Length[xlst]}]; |
| Show[Graphics[{Line[{{0, 1}, {0, 0}, {Length[xlst], 0}}], ptlst}], |
| PlotRange -> {0, 1}, AspectRatio -> 1/2] |
|
| l[s_, x_] := s*x*(1 - x) |
| SeedRandom[3142]; s = 4; num = 200; del = 0.02; eps = 0.05; rd = 0.000001; |
| x = Random[Real,{0,1}]; |
| xlst = {x}; |
| Do[{If[Abs[x - 0.25] < del, {s = Solve[0.25 == t*x*(1 - x), t][[1, 1, 2]] |
| + Random[Real, {-rd, rd}], x = l[s, x], AppendTo[xlst, x]}], |
| s = 4, x = l[s, x], AppendTo[xlst, x], |
| If[Abs[x - 0.75] > eps, {s = Solve[0.75 == t*x*(1 - x), t][[1, 1, 2]] |
| + Random[Real, {-rd, rd}], x = l[s, x], AppendTo[xlst, x]}], |
| s = 4, x = l[s, x], AppendTo[xlst, x]}, {i, 1, num}] |
| ptlst = {}; |
| Do[AppendTo[ptlst, Point[{i, xlst[[i]]}]], {i, 1, Length[xlst]}]; |
| Show[Graphics[{Line[{{0, 1}, {0, 0}, {Length[xlst], 0}}], ptlst}], |
| PlotRange -> {0, 1}, AspectRatio -> 1/2] |
|
| B.20 Sensitivity analysis. |
| Probably this one you don't need to copy and paste this. |
| MatrixForm[Transpose[{{1,2},{3,4}}]] |
|
| B.21 The Clancy-Rudy models. |
| Run these bits of code in the order presented. |
| p = {{0.8, 0.5, 0, 0, 0}, {0.2, 0.2, 0.2, 0.2, 0}, {0, 0.2, 0.7, 0.1, 0}, {0, 0.1, 0.1, 0.6, 0.1}, {0, 0, 0, 0.1, 0.9}}; |
|
| Now we generate the cumulative probability matrix. |
| cp = {{}, {}, {}, {}, {}}; |
| cp[[1]] = p[[1]]; |
| cp[[2]] = cp[[1]] + p[[2]]; |
| cp[[3]] = cp[[2]] + p[[3]]; |
| cp[[4]] = cp[[3]] + p[[4]]; |
| cp[[5]] = cp[[4]] + p[[5]]; |
|
| Start with a random state. Initialize the list of states. Initialize the state counters ot 0. |
| s = Random[Integer, {1, 5}]; |
| slst = {s}; ct[1] = 0; ct[2] = 0; ct[3] = 0; ct[4] = 0; ct[5] = 0; |
| num = 100000; |
|
| Next we generate a sequence of random numbers, and based on the current state and the
cumulative probability matrix, determine the next state and augment the appropriate state counter.
This may take a few moments to run. |
| Do[{x = Random[Real, {0, 1}], |
| If[x < cp[[1, s]], t = 1, |
| If[x < cp[[2, s]], t = 2, |
| If[x < cp[[3, s]], t = 3, |
| If[x < cp[[4, s]], t = 4, t = 5]]]], |
| s = t, ct[s] = ct[s] + 1, AppendTo[slst, s]}, {i, 1, num}] |
|
| Finally, to compare the counts with the fractions calculated in
Practice Problem 16.4.1 (b) use this command. |
| N[{ct[1]/num, ct[2]/num, ct[3]/num, ct[4]/num, ct[5]/num}] |
|
| B.22 A genetic toggle switch |
| Here is code to plot the curves (t,x(t)) and (t,y(t)), and to
plot the trajectory (x(t), y(t)). |
| a1 = 2.5; a2 = 2.5; b1 = 2; b2 = 2; c = 0.08; tmax = 300; |
| sol = NDSolve[{ |
| x'[t] == -x[t] + a1/(1 + y[t]^b1) + c*Sin[t/10], |
| y'[t] == -y[t] + a2/(1 + x[t]^b2)+ c*Sin[t/15], |
| x[0] == 0.5, y[0] == 1.5},{x, y}, {t, 0, tmax},MaxSteps->8000] |
| Plot[{Evaluate[x[t] /. sol],Evaluate[y[t]/.sol]}, {t, 0, tmax},PlotRange->All] |
| ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tmax}, PlotRange->All,PlotPoints->100] |
|
| B.23 Transcription networks. |
| ay=1; by=1; az=1; bz=2; kxy=0.5; kxz=1.0; kyz=0.75; tmax=10; |
| x[t_]:=If[t<5,1.25,0.55]; |
| sol = NDSolve[{ |
| y'[t]==If[x[t]>kxy, by - ay*y[t],-ay*y[t]], |
| z'[t]==If[x[t]>kxz && y[t]>kyz,bz - az*z[t],-az*z[t]], |
| y[0]==0,z[0]==0}, {y,z}, {t, 0, tmax},MaxSteps->8000, |
| Method->{"DiscontinuityProcessing"->False}] |
| Plot[{Evaluate[y[t] /. sol], Evaluate[z[t] /. sol], x[t]}, {t, 0, tmax}, PlotRange -> All] |
|
| B.24 Topological methods |
| First, symbol-driven IFS. Paste a sequence of c, t, g, and a between the brackets of dlst = {}; |
| dlst={}; |
| tr[c,{x_,y_}]:={0.5*x,0.5*y}; |
| tr[a,{x_,y_}]:={0.5*x + 0.5,0.5*y}; |
| tr[t,{x_,y_}]:={0.5*x,0.5*y + 0.5}; |
| tr[g,{x_,y_}]:={0.5*x + 0.5,0.5*y + 0.5}; |
| ptlst={}; |
| {x,y}={0.5,0.5}; |
| Do[{{nx,ny}=tr[dlst[[i]],{x,y}],{x,y}={nx,ny}, |
| AppendTo[ptlst,Point[{x,y}]]},{i,1,Length[dlst]}] |
| Show[Graphics[{Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],PointSize[0.01],ptlst}], |
| PlotRange->{{-.01,1.01},{-.01,1.01}},AspectRatio->Automatic] |
|
| Here is code to plot the driven IFS for numerical data, entered in as comma-separated numbers in dlst,
and to calculate the address occupancies. |
| For equal-size bins the bin boundaries are |
| rng = Max[dlst] - Min[dlst]; |
| B1 = Min[dlst] + 0.25*rng; |
| B2 = Min[dlst] + 0.50*rng; |
| B3 = Min[dlst] + 0.75*rng; |
|
| For equal-size bins the bin boundaries are |
| qlst = Quartiles[dlst]; |
| B1 = qlst[[1]]; |
| B2 = qlst[[2]]; |
| B3 = qlst[[3]]; |
|
| For median-centered bins with parameter d the bin boundaries are |
| rng = Max[dlst] - Min[dlst]; |
| d=0.3; |
| B2 = Median[dlst]; |
| B1 = B2 - d*rng; |
| B3 = B2 + d*rng; |
|
| Paste the choice of bin boundaries immediately after the dlst = {}; line. |
| To generate the driven IFS image use this code. |
| dlst = {}; |
| (* Remember to paste the choice of bin boundaries here *) |
| tr[1,{x_,y_}]:={0.5*x, 0.5*y}; |
| tr[2,{x_,y_}]:={0.5*x + 0.5, 0.5*y}; |
| tr[3,{x_,y_}]:={0.5*x, 0.5*y + 0.5}; |
| tr[4,{x_,y_}]:={0.5*x + 0.5, 0.5*y + 0.5}; |
| ind[z_]:=If[z < B1,1,If[z < B2,2,If[z < B3,3,4]]]; |
| Do[ad1[i]=0.,{i,1,4}]; |
| Do[Do[ad2[i,j]=0,{i,1,4}],{j,1,4}]; |
| Do[Do[Do[ad3[i,j,k]=0,{i,1,4}],{j,1,4}],{k,1,4}]; |
| {x,y}={0.5,0.5}; |
| drlst={Point[{x,y}]}; |
| w=ind[dlst[[1]]]; {x,y}=tr[w,{x,y}]; |
| AppendTo[drlst,Point[{x,y}]];ad1[w]=ad1[w]+1; |
| v=ind[dlst[[2]]]; {x,y}=tr[v,{x,y}]; |
| AppendTo[drlst,Point[{x,y}]]; ad1[v]=ad1[v]+1; adr[v,w]=ad2[v,w]+1; |
| Do[{u=ind[dlst[[i]]], {x,y}=tr[u,{x,y}], |
| AppendTo[drlst,Point[{x,y}]], ad1[u]=ad1[u]+1, |
| ad2[u,v]=ad2[u,v]+1, ad3[u,v,w]=ad3[u,v,w]+1, w=v,v=u},{i,3,Length[dlst]}]; |
| lnlst1={}; Do[{AppendTo[lnlst1,Line[{{x,0},{x,1}}]],AppendTo[lnlst1,Line[{{0,x},{1,x}}]]},{x,0,1,1/4}]; |
| Show[Graphics[{{GrayLevel[.6],lnlst1},PointSize[.015],drlst}], |
| PlotRange->{{-.01,1.01},{-.01,1.01}},AspectRatio->1] |
|
| To find the number of points with address 124, for example, after the program has run use the command |
|
|