Figure 5.6:

Full model solution for k_1=1, k_{-1}=0.5, k_2=k_{-2}=1.

Code for Figure 5.6

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
% 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.

k1   = 1;
k_1  = 1/2;
k2   = 1;
k_2  = 1;
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
nscheme = 0;


if (nscheme == 0)
    % 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

p = struct(); % Create structure to pass parameters to ode15s function
p.k1 = k1;
p.k_1 = k_1;
p.k2 = k2;
p.k_2 = k_2;
p.nscheme = nscheme;

tfinal = 5;
ntimes = 200;
tout   = linspace(0, tfinal, ntimes)';

opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s (@(t,x) rhs(t,x,p), tout, x0, opts);

% compute RA RB RC production rates
RA = -k1*x(:,1) + k_1*x(:,2);
RC =  k2*x(:,2) - k_2*x(:,3);
RB = -RA - RC;

table = [tout x RA RB RC];

save -ascii rls_1.dat table;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    plot (table(:,1),table(:,2:4));
    % TITLE
end % PLOTTING

rhs.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function retval = rhs(t,x,p)
    c_A = x(1);
    c_B = x(2);
    c_C = x(3);
    r1  = p.k1*c_A - p.k_1*c_B;
    r2  = p.k2*c_B - p.k_2*c_C;
    retval = zeros (3, 1);

    if (p.nscheme == 0)
        % full model
        retval(1) = -r1;
        retval(2) = r1 - r2;
        retval(3) = r2;
    elseif (p.nscheme == 1)
        % reaction 1 is rls, reaction 2 at equilibrium
        K2 = p.k2/p.k_2;
        retval(1) = -r1;
        retval(2) = 1/(1+K2)*r1;
        retval(3) = K2/(1+K2)*r1;
    elseif (p.nscheme == 2)
        % reaction 2 is rls, reaction 1 at equilibrium
        K1 = p.k1/p.k_1;
        retval(1) = -1/(1+K1) * r2;
        retval(2) = -K1/(1+K1) * r2;
        retval(3) = r2;
    elseif (p.nscheme == 3)
        % QSSA on B, redefine the c_B concentration under QSSA
        c_B = p.k1/(p.k_1 + p.k2)*c_A + p.k_2/(p.k_1 + p.k2)*c_C;
        r1 = p.k1*c_A - p.k_1*c_B;
        r2 = p.k2*c_B - p.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) = p.k1*p.k2*(p.k_2 - p.k1)/...
        (p.k_1 + p.k2)^2*c_A + p.k_1*p.k_2*...
        (p.k1 - p.k_2)/(p.k_1 + p.k2)^2*c_C;
        retval(3) = r2;
    else
        error ('nscheme out of range');
    end