%%%% solve ODE model to get output time profile % Modified version (by Peter Halling) of Hui's program % This version tries to omit calculation of sensitivities, and focus just % on progress of state variables. % Also calls ODE_ProchiralOrdered.m %clear all clear variables %following help guidance close all clc format long % Input values not normally to be changed %initrate = 0.005; %sets time scale, used to calculate EobySo % But now has to be set within function DoOneRun % Input values here may be changed showplot = false; D0 = 0.00333; % M units Q = 0.00167; % M units % These first are assumed to remain constant via cofactor recycling S0 = 0.2; % M units Keq = 5; Evalue = 39; % Next 2 values here are with concn units (based on M), then first order const % with s-1, then rest are dimensionless KMD = 0.001; KMS = 0.001; k4 = 1000; SCRrat = 1; k34rat = 1; k3rat = 1; k14rat = 1; k12Rrat = 1; k12Srat = 1; %dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, k3rat, ... % k14rat, k12Rrat, k12Srat, true); % For test if k4 < 0 % never for D0 = [0.005 0.01 0.02] % 0.001 gives EobySo too big dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end D0 = 0.01 for Q = [0.001 0.005 0.01 0.02] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end Q = 0.01 end if k4<0 for D0 = [0.005 0.00667 0.008 0.009] Q = 0.01 - D0; dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end for D0 = [0.0025 0.00333 0.004 0.0045] Q = 0.005 - D0; dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end for D0 = [0.00125 0.00167 0.002 0.00225] Q = 0.0025 - D0; dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end D0 = 0.00333 Q = 0.00167 for S0 = [0.05 0.1 0.2 0.4 1] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end S0 = 0.2 %end for Keq = [1 2.5 5 10 20] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end Keq = 5 end if k4 < 0 % never for Evalue = [9 19 39 100] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end Evalue = 19 end if k4 < 0 % never for KMD = [0.0005 0.001 0.002 0.01 0.1] for KMS = [0.001 0.01 0.1 1] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end KMD = 0.001 KMS = 0.1 %end %if k12Srat < 0 for KMD = [0.0005 0.001 0.002 0.01] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end KMD = 0.001 %end %if k34Srat < 0 % never for KMS = [0.0005 0.001 0.002] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end KMS = 0.001 %if k34Srat < 0 % i.e. never, just skip this bit for k4 = [500 1000 2000] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k4 = 1000 %end for SCRrat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end SCRrat = 1 for k34rat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k34rat = 1 for k3rat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k3rat = 1 %end %if k14rat < 0 for k14rat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k14rat = 1 %end %if k34Srat < 0 % i.e. never for k12Rrat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k12Rrat = 1 for k12Srat = [0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end k12Srat = 1 end if k4<0 for Keq = [1 2 5 10 20] for S0 = [0.05 0.1 0.2 0.5 1 2] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end KMD = 0.001 KMS = 0.001 end if k4 < 0 for bigrat = [0.0001 0.0002 0.0005 0.001 0.002 0.005 0.01] for KMS = [0.0005 0.001 0.002] for SCRrat = [1 2 5 10] for k34rat = [0.5 1 2] k3rat = bigrat*k34rat/(KMS*SCRrat); if (k3rat > 0.09) if (k3rat < 11) dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end end end end end end if k4 < 0 for SCRrat = [0.1 0.2 0.5 1 2 5 10] for KMS = [0.0005 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end SCRrat = 1 KMS = 0.001 for k3rat = [0.1 0.2 0.5 1 2 5 10] for k34rat = [0.1 0.2 0.5 1 2 5 10] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end SCRrat = 1 KMS = 0.001 %end %KMS = 0.1 for SCRrat = [0.1 0.2 0.5 1 2 5 10] for k3rat = [0.2 0.5 1 2 5] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end SCRrat = 1 k3rat = 1 end if k4 < 0 for S0 = [1] %[0.05 0.1 0.2 0.5 1] for Evalue = [19] % 39 99 199] dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); end end S0 = 0.2 Evalue = 39 end dummy = DoOneRun (D0, Q, 0.5, Keq, 19, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); dummy = DoOneRun (D0, Q, 1, Keq, Evalue, KMD, KMS, k4, SCRrat, k34rat, ... k3rat, k14rat, k12Rrat, k12Srat, false); function dummy = DoOneRun (D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, ... k34rat, k3rat, k14rat, k12Rrat, k12Srat, showplot) initrate = 0.005; %sets time scale, used to calculate EobySo % Now calculate rate constants, at first all on dimensioned basis k3R=k34rat*k4; k3S=k3R/(k3rat*Evalue); thisdenom = KMS*(2*Evalue*k3R*k3S + Evalue*k3R*k4 + Evalue*k3S*k4 + k3R*k3S + k3R*k4); SCRf=k4*k3R*k3S*Evalue/thisdenom; k1=SCRf*(Evalue+1)*KMS/(Evalue*KMD); SCSf=SCRf/Evalue; SCRb=SCRf/SCRrat; SCSb=SCRb/Evalue; km4=k1*k14rat; km1 = k1*k4*SCRf/(Keq*km4*SCRb); km2R=k12Rrat*km1; km2S = k12Srat*km1/Evalue; k2R = SCRf*(k3R + km2R)/k3R; km3R = SCRb*(k3R + km2R)/km2R; k2S = SCSf*(k3S + km2S)/k3S; km3S = SCSb*(k3S + km2S)/km2S; term1 = (1/k4 + 1/k3S)/Evalue + 1/k4 + 1/k3R + (Q*km4 + k4)*(1 + 1/Evalue)/(k1*D0*k4); term2 = (1 + km1*(Q*km4/k4 + 1)/(k1*D0))/(S0*SCRf); EobySo = initrate*(term1 + term2); if EobySo > 0.01 toobig = true; %error ('EobySo > 0.01') else toobig = false; end % Remove dimensions k1n = S0*k1; k2Rn = S0*k2R; k2Sn = S0*k2S; km3Rn = S0*km3R; km3Sn = S0*km3S; km4n = S0*km4; Sn0 = 1; % By definition of non-dimensional basis En0 = 1; % Same Dn0 = D0/S0; Qn0 = Q/S0; % All other state variables start at 0 % Set up the values of state variables and simulation parameters % initial condition of state variables Xinit = [En0 0 0 0 0 Dn0 Sn0 0 0 Qn0]'; % nominal parameter values k = [k1n k2Rn k2Sn k3R k3S k4 km1 km2R km2S km3Rn km3Sn km4n]'; % set up the simulation time and sampling interval %tspan = 0:0.001:0.1; % start with transient %tspan = 0:1:100; % initial rate period %tspan = 0:20:1000; % further progress tspan = 0:100:5000; % to get better conversion on slow reactions %tspan = 0:200:10000; % to get better conversion on slow reactions if ~toobig % With normalised non-dimensional values for all variables and parameters, % probably similar AbsTol for all, so try keeping default 1e-6. AbsTol_value = 1e-6; %AbsTol_value=[1e-11;1e-11;1e-13;1e-11;1e-12;(1e-6).*ones(5,1);(1e-13).*ones(11,1);... % (1e-14).*ones(11,1);(1e-15).*ones(11,1);(1e-13).*ones(11,1);(1e-14).*ones(11,1);... % (1e-9).*ones(55,1)]; %?????????? options = odeset('AbsTol',AbsTol_value); % Run the ODE simulation %[t, x_s] = feval(@ode15s, @ODE_Prochiral, tspan, xs_init, options, k, EobySo); % This doesn't work any more when have extra parameter EobySo, but managed % to avoid feval call completely after viewing ode15s help %[t,y] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0); Example from help [t, xout] = ode15s(@(t,x) ODE_ProchiralOrdered(t,x,k,EobySo), tspan, Xinit, options); % Help gives "[t,y] = ode15s(odefun,tspan,y0). Each row in the solution array %y corresponds to a value returned in column vector t." %disp(xout) %Yes, xout has a column for each state variable, and a row for each time, so x = transpose(xout); % x has one row for each state variable, one col for each time PRn = x(8,:); PSn = x(9,:); diff = PRn - PSn; conversion = PRn + PSn; ees = diff./conversion; else % is toobig showplot = false; t = '*'; conversion = 'EobySo too big'; ees = 'EobySo too big'; end % Prepare for output names = ["D0", "Q", "S0", "Keq", "Eval", "KMD", "KMS", "k4", "SCRrat", ... "k34rat", "k3rat", "k14rat", "k12Rrat", "k12Srat"]; values = [D0, Q, S0, Keq, Evalue, KMD, KMS, k4, SCRrat, ... k34rat, k3rat, k14rat, k12Rrat, k12Srat]; nonames = length(names); if length(values) ~= nonames error ('Different length names and values') end if showplot nolines = fix((nonames-1)/5) + 1; for lineno = 1:nolines firstctr = 5*(lineno-1) + 1; lastctr = 5*(lineno-1)+ 5; if lastctr > nonames lastctr = nonames; end lbl(lineno)=""; for ctr = firstctr:lastctr lbl(lineno) = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " "); %thislbl = strcat(lbl(lineno), names(ctr), "=", num2str(values(ctr)), " ") end end figure (1) %plot(t,x(1,:),'b--',t,x(2,:),'r-',t,x(3,:),'g-.',t,x(4,:),'m-',t,x(5,:),'c-.','LineWidth',2.5); %plot(t,x(6,:),'b--',t,x(7,:),'r-',t,x(8,:),'g-.',t,x(9,:),'m-',t,x(10,:),'c-.','LineWidth',2.5); plot(t,conversion,'b--',t,ees,'r-',t,x(1,:),'g-.','LineWidth',2.5); hold on; %plot(t(b),a,'g.','MarkerSize',24); %text(t(b),a+0.05,max_text);xlabel('\fontsize{16}time/s'); % From help %dim = [.2 .5 .3 .3]; %annotation('textbox',dim,'String',str,'FitBoxToText','on'); % By trial and error, x and y run from origin in bottom left %dim = [0.15 .6 .3 .3]; %annotation('textbox',dim,'String',lbl_text,'FitBoxToText','on'); %dim = [0.15 .53 .3 .3]; %annotation('textbox',dim,'String',lbl_text2,'FitBoxToText','on'); %dim = [0.15 .46 .3 .3]; %annotation('textbox',dim,'String',lbl_text3,'FitBoxToText','on'); for ctr = 1:nolines dim = [0.15 0.67-0.07*ctr .3 .3]; annotation('textbox',dim,'String',lbl(ctr), 'FitBoxToText', 'on'); end ylabel('\fontsize{16}non-dimensional concentrations'); %title('\fontsize{14}time profile of state variables'); %legend('En','EDn','EtRn','EtSn','EQn'); %legend('Dn','Qn','Sn','PRn','PSn'); legend('Conv','ee_P','En'); end % Output_progress filename = 'C:\Users\ph\My Documents\A Edited\ProgressCurves.xlsx'; %rowno = int8(5) startcol = 5; firstdatacol = startcol + nonames + 2; rowno = FindBlankRow(filename, 5, firstdatacol); for ctr = 1:nonames writematrix(names(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, startcol+ctr-1)); writematrix(values(ctr), filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+ 1, startcol+ctr-1)); end writematrix('time',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno, firstdatacol)); writematrix('conv',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+1, firstdatacol)); writematrix('ee_P',filename, 'Sheet', 1, 'Range', MakeCellRef(rowno+2, firstdatacol)); nextref = MakeCellRef(rowno,firstdatacol+1); writematrix(transpose(t),filename,'Sheet',1,'Range',nextref); nextref = MakeCellRef(rowno+1, firstdatacol+1); writematrix(conversion,filename,'Sheet',1,'Range',nextref); nextref = MakeCellRef(rowno+2, firstdatacol+1); writematrix(ees,filename,'Sheet',1,'Range',nextref); dummy = 1; %to prevent warning end function rowno = FindBlankRow(filename, startrow, colno) %rowno = startrow; %foundit = false; %while ~foundit % cellref = MakeCellRef(int8(rowno), int8(colno)) % entry = readmatrix(filename, 'Sheet', 1, 'Range', strcat(cellref, ':', cellref)) % Doesn't work, reads series of rows, then type mismatch on next % statement % if entry == '' % foundit = true % end % rowno = rowno + 1; cellref = MakeCellRef(int8(startrow), int8(colno)); entries = readmatrix(filename, 'Sheet', 1, 'Range', strcat(cellref, ':', cellref)); norows = length(entries); rowno = startrow + norows; end function cellref = MakeCellRef(rowno, colno) firstletno = idivide(colno - 1, int8(26)); secondletno = mod(colno - 1, int8(26)) + 1; if firstletno > 0 cellref = char(firstletno + 64); else cellref = ''; end cellref = strcat(cellref,char(secondletno+64)); cellref = strcat(cellref,num2str(rowno)); end