Figure 5.9 (page 214):

Comparison of equilibrium assumption to full model for k_1=1, k_{-1}=0.5, k_2=k_{-2}=10.

Code for Figure 5.9

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

global k1 k_1 k2 k_2 nscheme
k1      = 1;
k_1     = 1/2;
k2      = 10;
k_2     = 10;
c_A0    = 1;
c_B0    = 0.4;
c_C0    = 0;
%
% nscheme: 0 full model, 1 1st is rls, 2 2nd is rls; 3 qssa on B 
%          4 incorrect rls on 1, 5 correct rls on 1 but wrong IC
nscheme = 1;


if (nscheme == 0 || nscheme == 4 || nscheme ==5)
  %
  % IC for full model
  %
  x0=[c_A0; c_B0; c_C0];
elseif (nscheme == 1)
  %
  % IC for rx 1 is rls
  %
  K2    = k2/k_2;
  x0=[c_A0; 1/(1+K2)*(c_B0+c_C0); K2/(1+K2)*(c_B0+c_C0)];
elseif (nscheme == 2)
  %
  % IC for rx 2 is rls
  %
  K1    = k1/k_1;
  x0=[1/(1+K1)*(c_A0+c_B0); K1/(1+K1)*(c_A0+c_B0); c_C0];
elseif (nscheme == 3)
  %
  % IC for QSSA model
  %
  x0=[c_A0; k1/(k_1+k2)*c_A0 + k_2/(k_1+k2)*c_C0; c_C0];
else
  error ('nscheme out of range');
end

tfinal = 5;
ntimes = 200;
tout   = linspace(0, tfinal, ntimes)';
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s (@rhs, tout, x0, opts);

table = [tout x ];

%%
%% compare to full model
%%

nscheme=0;
if (nscheme == 0 || nscheme == 4 || nscheme ==5)
  %
  % IC for full model
  %
  x0=[c_A0; c_B0; c_C0];
elseif (nscheme == 1)
  %
  % IC for rx 1 is rls
  %
  K2    = k2/k_2;
  x0=[c_A0; 1/(1+K2)*(c_B0+c_C0); K2/(1+K2)*(c_B0+c_C0)];
elseif (nscheme == 2)
  %
  % IC for rx 2 is rls
  %
  K1    = k1/k_1;
  x0=[1/(1+K1)*(c_A0+c_B0); K1/(1+K1)*(c_A0+c_B0); c_C0];
elseif (nscheme == 3)
  %
  % IC for QSSA model
  %
  x0=[c_A0; k1/(k_1+k2)*c_A0 + k_2/(k_1+k2)*c_C0; c_C0];
else
  error ('nscheme out of range');
end

tfinal = 5;
ntimes = 200;
tout   = linspace(0, tfinal, ntimes)';
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s (@rhs, tout, x0, opts);

table = [table x ];


plot (table(:,1),table(:,2:7));
axis ([0,3,0,1]);
title ('Figure 5.9')

rhs.m

function retval = rhs(t,x)
global k1 k_1 k2 k_2 nscheme
c_A   = x(1);
c_B   = x(2);
c_C   = x(3);
r1  = k1*c_A - k_1*c_B;
r2  = k2*c_B - k_2*c_C;
retval = zeros (3, 1);
if (nscheme == 0)
   %
   % full model
   %
   retval(1) = -r1; 
   retval(2) = r1 - r2;
   retval(3) = r2;
elseif (nscheme == 1 || nscheme == 5)
   %
   % reaction 1 is rls, reaction 2 at equilibrium
   %
   K2    = k2/k_2;
   retval(1) = -r1; 
   retval(2) =  1/(1+K2)*r1;
   retval(3) =  K2/(1+K2)*r1;
elseif (nscheme == 2)
  %
  % reaction 2 is rls, reaction 1 at equilibrium
  %
  K1    = k1/k_1;
  retval(1) = -1/(1+K1) * r2; 
  retval(2) = -K1/(1+K1) * r2;
  retval(3) = r2;
elseif (nscheme == 3)
  %
  % QSSA on B, redefine the c_B concentration under QSSA
  %
  c_B= k1/(k_1+k2)*c_A + k_2/(k_1+k2)*c_C;
  r1  = k1*c_A - k_1*c_B;
  r2  = k2*c_B - k_2*c_C;
  retval(1) = -r1; 
  %retval(2) = r1 - r2; don't use this one, it is always zero and
  %therefore incorrect
  retval(2) = k1*k2*(k_2-k1)/(k_1+k2)^2 * c_A +  ...
              k_1*k_2*(k1-k_2)/(k_1+k2)^2*c_C;
  retval(3) = r2;
elseif (nscheme == 4)
  %
  % reaction 1 is rls, but incorrectly setting r2=0 everywhere
  %
  K2 = k2/k_2;
  retval(1) = -r1; 
  retval(2) =  r1;
  retval(3) =  0;
else
  error ('nscheme out of range');
end