Page 274
% make up an m-file, enzymeRate.m, with
%   function Yprime=enzymeRate(t,Y);
%   k1=0.1; km1=0.1; k2=0.1; e0=0.4;
%   Yprime=[k1*Y(2)*(e0-Y(1))-(km1+k2)*Y(1);
%     -k1*Y(2)*(e0-Y(1))+km1*Y(1)];
[t,Y]=ode23('enzymeRate',[0 100],[0; 0.8]);
plot(t,Y)

               Page 275
vmax=10; kM=15;
s=0:.1:150;
v=vmax.*s./(kM+s);
plot(s,v)
hold on
asy=vmax*ones(size(s));
plot(s,asy)
x=[0 13]; y=[0 13*vmax/kM];
plot(x,y)
x=[0 15.3]; y=[vmax/2 vmax/2];
plot(x,y)
x=[15.3 15.3]; y=[0 vmax/2];
plot(x,y)

             Exercises/Experiments
1a.
% Matlab cannot symbolically solve the system but we can proceed
% this way: add the two equations and notice that the first terms
% cancel and the second terms nearly cancel leaving -k2*c = 0.
% This shows that c=0.  With c=0 in either equation it is easy
% to see that s=0 too. Now we find the Jacobian numerically.
% Make an m-file, enzyme96.m, with
%   function csPrime=enzyme96(c,s);
%   % set values for k1, k2, km1, and e0
%   k1=1; k2=2; km1=1.5; e0=5;
%   csPrime=[k1*s.*(e0-c)-(km1+k2)*c; -k1*s.*(e0-c)+km1*c];
%
% The Jacobian = the matrix whose first column is the derivative
% of the component functions with respect to c and the second
% column is with respect to s.  Take derivatives at c=s=0.
J1=(enzyme96(0+eps,0) - enzyme96(0,0))/eps;
J2=(enzyme96(0,0+eps) - enzyme96(0,0))/eps;
J=[J1 J2]; % Jacobian at (0,0)
eig(J) % both values neg. real, so (0,0) stable

2.
% For problem 2 and 3
% Make an m-file, exer962.m, with
%   function Yprime=exer962(t,Y);
%   % Y(1)=c, Y(2)=s, Y(3)=e, Y(4)=p
%   k1=0.1; k2=0.1; km1=0.1;
%   % k1=1; k2=0.1; km1=0.025;
%   Yprime=[ k1*Y(3).*Y(2)-(km1+k2)*Y(1);
%     -k1*Y(3).*Y(2)+km1*Y(1);
%     -k1*Y(3).*Y(2)+(km1+k2)*Y(1);
%      k2*Y(1)];
%
s0=0.8; e0=0.4;
[t,Y]=ode23('exer962',[0 100],[0;s0; e0; 0]);
plot(t,Y) 

3.
% continued from previous problem and rerun
% for the second part change k1 and km1 above to k1=1; km1=0.025;
s0=10; e0=0.4;
[t,Y]=ode23('exer962',[0 100],[0;s0; e0; 0]);
plot(t,Y(:,1))  % graph of c

5b.
% contents of m-file exer964.m:
%   function Yprime=exer964(t,Y);
%   km1=1; k2=1;
%   k1=1;
%   Yprime=[-k1*Y(1)+k1*Y(2); k1*Y(1)-(km1+k2)*Y(2); k2*Y(2)];
%
a0=1; b0=0; c0=0;
[t1,Y1]=ode23('exer964',[0 7],[a0; b0; c0]);
plot(t1,Y1);
% now change k1 above to k1=10; note the y-axis scale
[t10,Y10]=ode23('exer964',[0 7],[a0; b0; c0]);
plot(t10,Y10)
% compare product c directly, be sure to check the y-axis scale
plot(t1,Y1(:,3),'b')
hold on
plot(t10,Y10(:,3),'r')