- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-11-10 14:40:25+13:00
- Desc:
- Fixing exposure files
- Permanent Source URI:
- https://models.cellml.org/workspace/705/rawfile/2fe62170795cead474fd15120998d64bb1d6e493/OLDFCU_adenylylCyclase.m
tic;
[voi, states, algebraic, constants, rates] = solveModel();
toc;
function [algebraicVariableCount] = getalgebraicVariableCount()
% Used later when setting a global variable with the number of algebraic variables.
% Note: This is not the "main method".
algebraicVariableCount =110;
end
% There are a total of 46 enp.stimDuries in each of the rate and state variable arrays.
% There are a total of 116 enp.stimDuries in the constant variable array.
%
function [voi, states, algebraic, constants, rates] = solveModel()
[legend_states, legend_algebraic, legend_voi, legend_constants] = createlegends();
% Create algebraic of correct size
global algebraicVariableCount; algebraicVariableCount = getalgebraicVariableCount();
% change amoiunt of total givben chemical species available (labels) by a
% multiplication factor (Gmult)
labels = {'q_AC_init in comp env', ...
'q_cAMP_init in comp env', ...
'q_R_B1_init in comp env', ...
'q_Gs_init in comp env', ...
'q_R_M2_init in comp env', ...
'q_Gi_init in comp env'};
% labels = {'q_R_B1_init in comp env', ...
% 'q_R_M2_init in comp env'};
[igs,~,~] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
Gmult = 1; %1e3;
[INIT_states, constants] = initConsts(Gmult, igs);
% ------------------------------------------
% Set timespan to solve over
voi = [0,1e-5]; % fast timescales due to FKC
% voi = [0 1e-3];
% ------------------------------------------
% Set numerical accuracy options for ODE solver
options = odeset('RelTol', 1e-07/3, 'AbsTol', 1e-07/3, 'MaxStep', 1e-9); % step: 1e-4 % no need for step=1e-9. same as 1e-5
% extra options for stimuli [ L_B1, L_M2 ] settings
p = struct;
% p.stimTimes = [[1,20];[3,40]];
p.stimThreshold = [-0.95, -0.95]; %[1e-3,1e-3];
p.stimHolding = [1e-9,1e-9];
p.stimMag = [1.888e1,1.888e1]; %[1e6,1e6];
p.stimStart = [1.2e-6,3e0];
p.tRamp = [1e-6,1]; % seconds for ramp start
p.stimDur = [1,1];
if false
p.stimMag = [0,0];
end
p.freq = [1,1]; %[0.25,0.25]; % per second
% Solve model with ODE solver
[voi, states] = ode15s(@(voi, states)computeRates(voi, states, constants,p), voi, INIT_states, options);
% Compute algebraic variables
[rates, algebraic] = computeRates(voi, states, constants,p);
algebraic = computeAlgebraic(algebraic, constants, states, voi,p);
labels = {'q_L_B1_init in comp env', ...
'q_L_M2_init in comp env', ...
'q_LR_B1Gs in comp env', ...
'q_LR_M2Gi in comp env', ...
'q_aGs_GTP in comp env', ...
'q_aGi_GTP in comp env', ...
'q_aGs_GTP_AC in comp cAMP', ...
'q_AC in comp env', ...
'q_cAMP in comp env', ...
'q_ATP in comp env'};
[i_con, i_st, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
if true
[~, i_st, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('LR_B1\nGmult=%g', Gmult),ceil(sqrt(length(i_st))))
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('LR_B1\nGmult=%g', Gmult),ceil(sqrt(length(i_alg))))
end
% plot chemostat concentrations: they should be constant
labels = {'L_B1_T in comp env',...
'L_M2_T in comp env',...
'R_B1_T in comp env',...
'Gs_T in comp env',...
'R_M2_T in comp env',...
'Gi_T in comp env',...
'adenosine_T in comp env'
};
labels_init = {'q_L_B1_init in comp env',...
'q_L_M2_init in comp env',...
'q_R_B1_init in comp env',...
'q_Gs_init in comp env',...
'q_R_M2_init in comp env',...
'q_Gi_init in comp env',...
'q_ATP_init in comp env'
};
[~, ~, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
[i_con, ~, i_algStim] = find_indices(labels_init, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
plot_2per(i_algStim,i_con, i_alg,voi,constants,algebraic,legend_algebraic,legend_constants,sprintf('CHEMOSTATS\nGmult=%g', Gmult),3)
% debug R_B1_T
if false
labels = {'q_R_B1 in comp LRGbinding_B1AR','q_R_B1 in comp GsProtein'};
[~, i_st, ~] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('R_B1\nGmult=%g', Gmult),ceil(sqrt(length(i_st))))
end
% debug adenosine_T
if false
labels = {'q_cAMP in comp env','q_PDE_cAMP in comp env','q_five_AMP in comp env','adenosine_T in comp env'};
[~, ~, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('AMP\nGmult=%g', Gmult),ceil(sqrt(length(i_alg))))
end
% debug aGs_GTP and aGiGTP
if false
labels = {'q_aGs_GTP in comp env','q_aGi_GTP in comp env',...
'q_aGs_GTP in comp GsProtein','q_aGi_GTP in comp GiProtein'...
'q_Gs in comp env','q_Gi in comp env'...
};
[~, i_st, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('aG_GTP\nGmult=%g', Gmult),ceil(sqrt(length(i_alg))))
plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('aG_GTP\nGmult=%g', Gmult),ceil(sqrt(length(i_st))))
end
% plot all q variables to confirm if we are SS
if false
figure();
n = ceil(sqrt(size(states,2)));
for i=1:size(states,2)
subplot(n,n,i)
plot(voi, states(:,i));
xlabel(legend_voi);
title(legend_states(i,:));
end
end
end
function [] = plot_selected(ids,x,y,legend_x,legend_y,titlestr,ns)
istart = 30;
figure();
% plot stimuli
for i = 1:length(ids)
subplot(ns,ns,i)
plot(x(istart:end), y(istart:end,ids(i)));
xlabel('time (s)');
l = legend(legend_y(ids(i),:));
set(l,'Interpreter','none');
end
suptitle(titlestr)
end
function [] = plot_2per(i_alg0,id0,ids,x,y0,y,legend_y,legend_y_con,titlestr,ns)
istart = 30;
figure();
% plot 2 per subplot
for i = 1:2
subplot(ns,ns,i)
plot(x(istart:end), y(istart:end,i_alg0(i)),'x-'); hold on;
plot(x(istart:end), y(istart:end,ids(i)));
xlabel('time (s)');
diff = mean(abs(y(istart:end,i_alg0(i)) - y(istart:end,ids(i))));
str = split(legend_y(i_alg0(i),:),' in');
title([str(1), strcat('error = ', num2str(diff))]);
% legend('init','sumSpecies');
legend(legend_y(i_alg0(i),:),legend_y(ids(i),:));
end
for i = 3:length(ids)
subplot(ns,ns,i)
plot(x(istart:end), ones(length(x)-istart+1,1)*y0(id0(i-2)),'x-'); hold on;
plot(x(istart:end), y(istart:end,ids(i)));
diff = abs(mean(y(istart:end,ids(i))) - y0(id0(i-2)));
xlabel('time (s)');
str = split(legend_y(ids(i),:),' in');
title([str(1), strcat('error = ', num2str(diff))]);
legend(legend_y_con(id0(i-2),:),legend_y(ids(i),:));
end
suptitle(strcat(titlestr,' eps = ',num2str(eps)))
end
function [i_con, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
% return the indices for the selected labels
all_legends = [legend_constants; legend_states; legend_algebraic];
i_con = [];
for i = 1:length(labels)
i_con = [i_con; find(strcmp(labels{i},legend_constants))];
end
i_st = [];
for i = 1:length(labels)
i_st = [i_st; find(strcmp(labels{i},legend_states))];
end
i_alg = [];
for i = 1:length(labels)
i_alg = [i_alg; find(strcmp(labels{i},legend_algebraic))];
end
end
function [legend_states, legend_algebraic, legend_voi, legend_constants] = createlegends()
legend_states = ''; legend_algebraic = ''; legend_voi = ''; legend_constants = '';
legend_constants(:,1) = strpad('kappa_1a in comp BG_parameters (fmol_per_sec)');
legend_constants(:,2) = strpad('kappa_1b in comp BG_parameters (fmol_per_sec)');
legend_constants(:,3) = strpad('kappa_2a in comp BG_parameters (fmol_per_sec)');
legend_constants(:,4) = strpad('kappa_2b in comp BG_parameters (fmol_per_sec)');
legend_constants(:,5) = strpad('kappa_3a in comp BG_parameters (fmol_per_sec)');
legend_constants(:,6) = strpad('kappa_3b in comp BG_parameters (fmol_per_sec)');
legend_constants(:,7) = strpad('kappa_4a in comp BG_parameters (fmol_per_sec)');
legend_constants(:,8) = strpad('kappa_4b in comp BG_parameters (fmol_per_sec)');
legend_constants(:,9) = strpad('kappa_5 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,10) = strpad('kappa_6 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,11) = strpad('kappa_7 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,12) = strpad('kappa_GiAC in comp BG_parameters (fmol_per_sec)');
legend_constants(:,13) = strpad('kappa_sig1_B1 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,14) = strpad('kappa_sig2_B1 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,15) = strpad('kappa_sig3_B1 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,16) = strpad('kappa_sig4_B1 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,17) = strpad('kappa_sig1_M2 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,18) = strpad('kappa_sig2_M2 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,19) = strpad('kappa_sig3_M2 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,20) = strpad('kappa_sig4_M2 in comp BG_parameters (fmol_per_sec)');
legend_constants(:,21) = strpad('kappa_act1_Gs in comp BG_parameters (fmol_per_sec)');
legend_constants(:,22) = strpad('kappa_act2_Gs in comp BG_parameters (fmol_per_sec)');
legend_constants(:,23) = strpad('kappa_hyd_Gs in comp BG_parameters (fmol_per_sec)');
legend_constants(:,24) = strpad('kappa_reassoc_Gs in comp BG_parameters (fmol_per_sec)');
legend_constants(:,25) = strpad('kappa_act1_Gi in comp BG_parameters (fmol_per_sec)');
legend_constants(:,26) = strpad('kappa_act2_Gi in comp BG_parameters (fmol_per_sec)');
legend_constants(:,27) = strpad('kappa_hyd_Gi in comp BG_parameters (fmol_per_sec)');
legend_constants(:,28) = strpad('kappa_reassoc_Gi in comp BG_parameters (fmol_per_sec)');
legend_constants(:,29) = strpad('K_ATP in comp BG_parameters (per_fmol)');
legend_constants(:,30) = strpad('K_cAMP in comp BG_parameters (per_fmol)');
legend_constants(:,31) = strpad('K_AC in comp BG_parameters (per_fmol)');
legend_constants(:,32) = strpad('K_AC_ATP in comp BG_parameters (per_fmol)');
legend_constants(:,33) = strpad('K_aGs_GTP_AC in comp BG_parameters (per_fmol)');
legend_constants(:,34) = strpad('K_aGs_GTP_AC_ATP in comp BG_parameters (per_fmol)');
legend_constants(:,35) = strpad('K_FSK_AC in comp BG_parameters (per_fmol)');
legend_constants(:,36) = strpad('K_FSK_AC_ATP in comp BG_parameters (per_fmol)');
legend_constants(:,37) = strpad('K_PDE in comp BG_parameters (per_fmol)');
legend_constants(:,38) = strpad('K_PDE_cAMP in comp BG_parameters (per_fmol)');
legend_constants(:,39) = strpad('K_five_AMP in comp BG_parameters (per_fmol)');
legend_constants(:,40) = strpad('K_IBMX in comp BG_parameters (per_fmol)');
legend_constants(:,41) = strpad('K_PDEinh in comp BG_parameters (per_fmol)');
legend_constants(:,42) = strpad('K_aGs_GTP in comp BG_parameters (per_fmol)');
legend_constants(:,43) = strpad('K_FSK in comp BG_parameters (per_fmol)');
legend_constants(:,44) = strpad('K_aGi_GTP in comp BG_parameters (per_fmol)');
legend_constants(:,45) = strpad('K_ACinh in comp BG_parameters (per_fmol)');
legend_constants(:,46) = strpad('K_PPi in comp BG_parameters (per_fmol)');
legend_constants(:,47) = strpad('K_L_B1 in comp BG_parameters (per_fmol)');
legend_constants(:,48) = strpad('K_R_B1 in comp BG_parameters (per_fmol)');
legend_constants(:,49) = strpad('K_Gs in comp BG_parameters (per_fmol)');
legend_constants(:,50) = strpad('K_LR_B1 in comp BG_parameters (per_fmol)');
legend_constants(:,51) = strpad('K_R_B1Gs in comp BG_parameters (per_fmol)');
legend_constants(:,52) = strpad('K_LR_B1Gs in comp BG_parameters (per_fmol)');
legend_constants(:,53) = strpad('K_L_M2 in comp BG_parameters (per_fmol)');
legend_constants(:,54) = strpad('K_R_M2 in comp BG_parameters (per_fmol)');
legend_constants(:,55) = strpad('K_Gi in comp BG_parameters (per_fmol)');
legend_constants(:,56) = strpad('K_LR_M2 in comp BG_parameters (per_fmol)');
legend_constants(:,57) = strpad('K_R_M2Gi in comp BG_parameters (per_fmol)');
legend_constants(:,58) = strpad('K_LR_M2Gi in comp BG_parameters (per_fmol)');
legend_constants(:,59) = strpad('K_beta_gamma_Gs in comp BG_parameters (per_fmol)');
legend_constants(:,60) = strpad('K_aGs_GDP in comp BG_parameters (per_fmol)');
legend_constants(:,61) = strpad('K_Pi in comp BG_parameters (per_fmol)');
legend_constants(:,62) = strpad('K_beta_gamma_Gi in comp BG_parameters (per_fmol)');
legend_constants(:,63) = strpad('K_aGi_GDP in comp BG_parameters (per_fmol)');
legend_voi = strpad('time in comp env (second)');
legend_constants(:,64) = strpad('vol_myo in comp env (pL)');
legend_constants(:,65) = strpad('freq in comp env (dimensionless)');
legend_constants(:,66) = strpad('stimSt in comp env (second)');
legend_constants(:,67) = strpad('stimSt2 in comp env (second)');
legend_constants(:,68) = strpad('stimDur in comp env (second)');
legend_constants(:,69) = strpad('tR in comp env (second)');
legend_constants(:,70) = strpad('stimMag in comp env');
legend_constants(:,71) = strpad('stimHolding in comp env');
legend_constants(:,113) = strpad('m in comp env (fmol_per_sec)');
legend_constants(:,72) = strpad('q_ATP_init in comp env');
legend_constants(:,73) = strpad('q_AC_init in comp env');
legend_constants(:,74) = strpad('q_cAMP_init in comp env');
legend_constants(:,75) = strpad('q_AC_ATP_init in comp env');
legend_constants(:,76) = strpad('q_FSK_init in comp env');
legend_constants(:,77) = strpad('q_FSK_AC_init in comp env');
legend_constants(:,78) = strpad('q_FSK_AC_ATP_init in comp env');
legend_constants(:,79) = strpad('q_aGs_GTP_init in comp env');
legend_constants(:,80) = strpad('q_aGs_GTP_AC_init in comp env');
legend_constants(:,81) = strpad('q_aGs_GTP_AC_ATP_init in comp env');
legend_constants(:,82) = strpad('q_PDE_init in comp env');
legend_constants(:,83) = strpad('q_PDEinh_init in comp env');
legend_constants(:,84) = strpad('q_PDE_cAMP_init in comp env');
legend_constants(:,85) = strpad('q_IBMX_init in comp env');
legend_constants(:,86) = strpad('q_five_AMP_init in comp env');
legend_constants(:,87) = strpad('q_aGi_GTP_init in comp env');
legend_constants(:,88) = strpad('q_ACinh_init in comp env');
legend_constants(:,89) = strpad('q_PPi_init in comp env');
legend_constants(:,90) = strpad('q_R_B1_init in comp env');
legend_constants(:,91) = strpad('q_Gs_init in comp env');
legend_algebraic(:,1) = strpad('q_L_B1_init in comp env');
legend_constants(:,92) = strpad('q_LR_B1_init in comp env');
legend_constants(:,93) = strpad('q_R_B1Gs_init in comp env');
legend_constants(:,94) = strpad('q_LR_B1Gs_init in comp env');
legend_algebraic(:,2) = strpad('q_L_M2_init in comp env');
legend_constants(:,95) = strpad('q_R_M2_init in comp env');
legend_constants(:,96) = strpad('q_Gi_init in comp env');
legend_constants(:,97) = strpad('q_LR_M2_init in comp env');
legend_constants(:,98) = strpad('q_R_M2Gi_init in comp env');
legend_constants(:,99) = strpad('q_LR_M2Gi_init in comp env');
legend_constants(:,100) = strpad('q_beta_gamma_Gs_init in comp env');
legend_constants(:,101) = strpad('q_aGs_GDP_init in comp env');
legend_constants(:,102) = strpad('q_Pi_init in comp env');
legend_constants(:,103) = strpad('q_beta_gamma_Gi_init in comp env');
legend_constants(:,104) = strpad('q_aGi_GDP_init in comp env');
legend_algebraic(:,3) = strpad('q_ATP in comp env');
legend_algebraic(:,4) = strpad('q_cAMP in comp env');
legend_algebraic(:,5) = strpad('q_AC in comp env');
legend_algebraic(:,6) = strpad('q_AC_ATP in comp env');
legend_algebraic(:,7) = strpad('q_aGs_GTP_AC in comp env');
legend_algebraic(:,8) = strpad('q_aGs_GTP_AC_ATP in comp env');
legend_algebraic(:,9) = strpad('q_FSK_AC in comp env');
legend_algebraic(:,10) = strpad('q_FSK_AC_ATP in comp env');
legend_algebraic(:,11) = strpad('q_PDE in comp env');
legend_algebraic(:,14) = strpad('q_PDE_cAMP in comp env');
legend_algebraic(:,15) = strpad('q_five_AMP in comp env');
legend_algebraic(:,17) = strpad('q_IBMX in comp env');
legend_algebraic(:,19) = strpad('q_PDEinh in comp env');
legend_algebraic(:,20) = strpad('q_aGs_GTP in comp env');
legend_algebraic(:,22) = strpad('q_FSK in comp env');
legend_algebraic(:,25) = strpad('q_aGi_GTP in comp env');
legend_algebraic(:,27) = strpad('q_ACinh in comp env');
legend_algebraic(:,28) = strpad('q_PPi in comp env');
legend_algebraic(:,21) = strpad('q_L_B1 in comp env');
legend_algebraic(:,24) = strpad('q_R_B1 in comp env');
legend_algebraic(:,26) = strpad('q_Gs in comp env');
legend_algebraic(:,30) = strpad('q_LR_B1 in comp env');
legend_algebraic(:,31) = strpad('q_R_B1Gs in comp env');
legend_algebraic(:,34) = strpad('q_LR_B1Gs in comp env');
legend_algebraic(:,29) = strpad('q_L_M2 in comp env');
legend_algebraic(:,32) = strpad('q_R_M2 in comp env');
legend_algebraic(:,35) = strpad('q_Gi in comp env');
legend_algebraic(:,39) = strpad('q_LR_M2 in comp env');
legend_algebraic(:,49) = strpad('q_R_M2Gi in comp env');
legend_algebraic(:,53) = strpad('q_LR_M2Gi in comp env');
legend_algebraic(:,38) = strpad('q_beta_gamma_Gs in comp env');
legend_algebraic(:,52) = strpad('q_aGs_GDP in comp env');
legend_constants(:,114) = strpad('q_Pi in comp env');
legend_algebraic(:,60) = strpad('q_beta_gamma_Gi in comp env');
legend_algebraic(:,76) = strpad('q_aGi_GDP in comp env');
legend_states(:,1) = strpad('q_ATP in comp cAMP');
legend_states(:,2) = strpad('q_cAMP in comp cAMP');
legend_states(:,3) = strpad('q_AC in comp cAMP');
legend_states(:,4) = strpad('q_AC_ATP in comp cAMP');
legend_states(:,5) = strpad('q_aGs_GTP_AC in comp cAMP');
legend_states(:,6) = strpad('q_aGs_GTP_AC_ATP in comp cAMP');
legend_states(:,7) = strpad('q_FSK_AC in comp cAMP');
legend_states(:,8) = strpad('q_FSK_AC_ATP in comp cAMP');
legend_states(:,9) = strpad('q_PDE in comp cAMP');
legend_states(:,10) = strpad('q_PDE_cAMP in comp cAMP');
legend_states(:,11) = strpad('q_five_AMP in comp cAMP');
legend_states(:,12) = strpad('q_IBMX in comp cAMP');
legend_states(:,13) = strpad('q_PDEinh in comp cAMP');
legend_states(:,14) = strpad('q_aGs_GTP in comp cAMP');
legend_states(:,15) = strpad('q_FSK in comp cAMP');
legend_states(:,16) = strpad('q_aGi_GTP in comp cAMP');
legend_states(:,17) = strpad('q_ACinh in comp cAMP');
legend_states(:,18) = strpad('q_PPi in comp cAMP');
legend_states(:,19) = strpad('q_L_B1 in comp LRGbinding_B1AR');
legend_states(:,20) = strpad('q_R_B1 in comp LRGbinding_B1AR');
legend_states(:,21) = strpad('q_Gs in comp LRGbinding_B1AR');
legend_states(:,22) = strpad('q_LR_B1 in comp LRGbinding_B1AR');
legend_states(:,23) = strpad('q_R_B1Gs in comp LRGbinding_B1AR');
legend_states(:,24) = strpad('q_LR_B1Gs in comp LRGbinding_B1AR');
legend_states(:,25) = strpad('q_L_M2 in comp LRGbinding_M2');
legend_states(:,26) = strpad('q_R_M2 in comp LRGbinding_M2');
legend_states(:,27) = strpad('q_Gi in comp LRGbinding_M2');
legend_states(:,28) = strpad('q_LR_M2 in comp LRGbinding_M2');
legend_states(:,29) = strpad('q_R_M2Gi in comp LRGbinding_M2');
legend_states(:,30) = strpad('q_LR_M2Gi in comp LRGbinding_M2');
legend_states(:,31) = strpad('q_R_B1 in comp GsProtein');
legend_states(:,32) = strpad('q_Gs in comp GsProtein');
legend_states(:,33) = strpad('q_R_B1Gs in comp GsProtein');
legend_states(:,34) = strpad('q_LR_B1 in comp GsProtein');
legend_states(:,35) = strpad('q_LR_B1Gs in comp GsProtein');
legend_states(:,36) = strpad('q_aGs_GTP in comp GsProtein');
legend_states(:,37) = strpad('q_beta_gamma_Gs in comp GsProtein');
legend_states(:,38) = strpad('q_aGs_GDP in comp GsProtein');
legend_constants(:,105) = strpad('q_Pi in comp GsProtein');
legend_states(:,39) = strpad('q_R_M2 in comp GiProtein');
legend_states(:,40) = strpad('q_Gi in comp GiProtein');
legend_states(:,41) = strpad('q_R_M2Gi in comp GiProtein');
legend_states(:,42) = strpad('q_LR_M2 in comp GiProtein');
legend_states(:,43) = strpad('q_LR_M2Gi in comp GiProtein');
legend_states(:,44) = strpad('q_aGi_GTP in comp GiProtein');
legend_states(:,45) = strpad('q_beta_gamma_Gi in comp GiProtein');
legend_states(:,46) = strpad('q_aGi_GDP in comp GiProtein');
legend_constants(:,106) = strpad('q_Pi in comp GiProtein');
legend_algebraic(:,40) = strpad('L_B1_T in comp env');
legend_algebraic(:,41) = strpad('R_B1_T in comp env');
legend_algebraic(:,54) = strpad('Gs_T in comp env');
legend_algebraic(:,61) = strpad('L_M2_T in comp env');
legend_algebraic(:,62) = strpad('R_M2_T in comp env');
legend_algebraic(:,77) = strpad('Gi_T in comp env');
legend_algebraic(:,18) = strpad('adenosine_T in comp env');
legend_constants(:,107) = strpad('R in comp constants (J_per_K_per_mol)');
legend_constants(:,108) = strpad('T in comp constants (kelvin)');
legend_constants(:,109) = strpad('F in comp constants (C_per_mol)');
legend_algebraic(:,120) = strpad('v1a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,121) = strpad('v1b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,122) = strpad('v2a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,123) = strpad('v2b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,124) = strpad('v3a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,125) = strpad('v3b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,126) = strpad('v4a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,128) = strpad('v4b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,130) = strpad('v5 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,127) = strpad('v6 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,129) = strpad('v7 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,131) = strpad('vGiAC in comp cAMP (fmol_per_sec)');
legend_algebraic(:,33) = strpad('mu_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,42) = strpad('mu_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,36) = strpad('mu_cAMP in comp cAMP (J_per_mol)');
legend_algebraic(:,50) = strpad('mu_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,112) = strpad('mu_FSK in comp cAMP (J_per_mol)');
legend_algebraic(:,72) = strpad('mu_FSK_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,78) = strpad('mu_FSK_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,108) = strpad('mu_aGs_GTP in comp cAMP (J_per_mol)');
legend_algebraic(:,55) = strpad('mu_aGs_GTP_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,63) = strpad('mu_aGs_GTP_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,84) = strpad('mu_PDE in comp cAMP (J_per_mol)');
legend_algebraic(:,104) = strpad('mu_PDEinh in comp cAMP (J_per_mol)');
legend_algebraic(:,89) = strpad('mu_PDE_cAMP in comp cAMP (J_per_mol)');
legend_algebraic(:,99) = strpad('mu_IBMX in comp cAMP (J_per_mol)');
legend_algebraic(:,94) = strpad('mu_five_AMP in comp cAMP (J_per_mol)');
legend_algebraic(:,115) = strpad('mu_aGi_GTP in comp cAMP (J_per_mol)');
legend_algebraic(:,117) = strpad('mu_ACinh in comp cAMP (J_per_mol)');
legend_algebraic(:,119) = strpad('mu_PPi in comp cAMP (J_per_mol)');
legend_constants(:,110) = strpad('vol in comp cAMP (pL)');
legend_algebraic(:,12) = strpad('ATP_T in comp cAMP');
legend_algebraic(:,13) = strpad('AC_T in comp cAMP');
legend_algebraic(:,23) = strpad('Gs_T in comp cAMP');
legend_algebraic(:,16) = strpad('cadenosine_T in comp cAMP');
legend_algebraic(:,43) = strpad('mu_L_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,51) = strpad('mu_R_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,56) = strpad('mu_Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,64) = strpad('mu_LR_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,73) = strpad('mu_R_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,79) = strpad('mu_LR_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,85) = strpad('vsig1_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,90) = strpad('vsig2_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,95) = strpad('vsig3_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,100) = strpad('vsig4_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_constants(:,111) = strpad('vol in comp LRGbinding_B1AR (pL)');
legend_algebraic(:,44) = strpad('L_T in comp LRGbinding_B1AR');
legend_algebraic(:,45) = strpad('R_T in comp LRGbinding_B1AR');
legend_algebraic(:,46) = strpad('Gs_T in comp LRGbinding_B1AR');
legend_algebraic(:,65) = strpad('mu_L_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,74) = strpad('mu_R_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,80) = strpad('mu_Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,86) = strpad('mu_LR_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,91) = strpad('mu_R_M2Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,96) = strpad('mu_LR_M2Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,101) = strpad('vsig1_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,105) = strpad('vsig2_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,109) = strpad('vsig3_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,113) = strpad('vsig4_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_constants(:,112) = strpad('vol in comp LRGbinding_M2 (pL)');
legend_algebraic(:,66) = strpad('L_T in comp LRGbinding_M2');
legend_algebraic(:,67) = strpad('R_T in comp LRGbinding_M2');
legend_algebraic(:,68) = strpad('Gi_T in comp LRGbinding_M2');
legend_algebraic(:,97) = strpad('vact1_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,102) = strpad('vact2_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,106) = strpad('vhyd_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,110) = strpad('vreassoc_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,57) = strpad('mu_R_B1 in comp GsProtein (J_per_mol)');
legend_algebraic(:,69) = strpad('mu_Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,37) = strpad('mu_R_B1Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,75) = strpad('mu_LR_B1 in comp GsProtein (J_per_mol)');
legend_algebraic(:,47) = strpad('mu_LR_B1Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,81) = strpad('mu_aGs_GTP in comp GsProtein (J_per_mol)');
legend_algebraic(:,87) = strpad('mu_beta_gamma_Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,92) = strpad('mu_aGs_GDP in comp GsProtein (J_per_mol)');
legend_constants(:,115) = strpad('mu_Pi in comp GsProtein (J_per_mol)');
legend_algebraic(:,48) = strpad('R_T in comp GsProtein');
legend_algebraic(:,58) = strpad('Gs_T in comp GsProtein');
legend_algebraic(:,111) = strpad('vact1_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,114) = strpad('vact2_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,116) = strpad('vhyd_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,118) = strpad('vreassoc_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,82) = strpad('mu_R_M2 in comp GiProtein (J_per_mol)');
legend_algebraic(:,88) = strpad('mu_Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,59) = strpad('mu_R_M2Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,93) = strpad('mu_LR_M2 in comp GiProtein (J_per_mol)');
legend_algebraic(:,70) = strpad('mu_LR_M2Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,98) = strpad('mu_aGi_GTP in comp GiProtein (J_per_mol)');
legend_algebraic(:,103) = strpad('mu_beta_gamma_Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,107) = strpad('mu_aGi_GDP in comp GiProtein (J_per_mol)');
legend_constants(:,116) = strpad('mu_Pi in comp GiProtein (J_per_mol)');
legend_algebraic(:,71) = strpad('R_T in comp GiProtein');
legend_algebraic(:,83) = strpad('Gi_T in comp GiProtein');
legend_RATES(:,1) = strpad('d/dt q_ATP in comp cAMP');
legend_RATES(:,3) = strpad('d/dt q_AC in comp cAMP');
legend_RATES(:,4) = strpad('d/dt q_AC_ATP in comp cAMP');
legend_RATES(:,2) = strpad('d/dt q_cAMP in comp cAMP');
legend_RATES(:,15) = strpad('d/dt q_FSK in comp cAMP');
legend_RATES(:,7) = strpad('d/dt q_FSK_AC in comp cAMP');
legend_RATES(:,8) = strpad('d/dt q_FSK_AC_ATP in comp cAMP');
legend_RATES(:,14) = strpad('d/dt q_aGs_GTP in comp cAMP');
legend_RATES(:,5) = strpad('d/dt q_aGs_GTP_AC in comp cAMP');
legend_RATES(:,6) = strpad('d/dt q_aGs_GTP_AC_ATP in comp cAMP');
legend_RATES(:,10) = strpad('d/dt q_PDE_cAMP in comp cAMP');
legend_RATES(:,9) = strpad('d/dt q_PDE in comp cAMP');
legend_RATES(:,12) = strpad('d/dt q_IBMX in comp cAMP');
legend_RATES(:,13) = strpad('d/dt q_PDEinh in comp cAMP');
legend_RATES(:,11) = strpad('d/dt q_five_AMP in comp cAMP');
legend_RATES(:,16) = strpad('d/dt q_aGi_GTP in comp cAMP');
legend_RATES(:,17) = strpad('d/dt q_ACinh in comp cAMP');
legend_RATES(:,18) = strpad('d/dt q_PPi in comp cAMP');
legend_RATES(:,19) = strpad('d/dt q_L_B1 in comp LRGbinding_B1AR');
legend_RATES(:,20) = strpad('d/dt q_R_B1 in comp LRGbinding_B1AR');
legend_RATES(:,21) = strpad('d/dt q_Gs in comp LRGbinding_B1AR');
legend_RATES(:,22) = strpad('d/dt q_LR_B1 in comp LRGbinding_B1AR');
legend_RATES(:,23) = strpad('d/dt q_R_B1Gs in comp LRGbinding_B1AR');
legend_RATES(:,24) = strpad('d/dt q_LR_B1Gs in comp LRGbinding_B1AR');
legend_RATES(:,25) = strpad('d/dt q_L_M2 in comp LRGbinding_M2');
legend_RATES(:,26) = strpad('d/dt q_R_M2 in comp LRGbinding_M2');
legend_RATES(:,27) = strpad('d/dt q_Gi in comp LRGbinding_M2');
legend_RATES(:,28) = strpad('d/dt q_LR_M2 in comp LRGbinding_M2');
legend_RATES(:,29) = strpad('d/dt q_R_M2Gi in comp LRGbinding_M2');
legend_RATES(:,30) = strpad('d/dt q_LR_M2Gi in comp LRGbinding_M2');
legend_RATES(:,31) = strpad('d/dt q_R_B1 in comp GsProtein');
legend_RATES(:,33) = strpad('d/dt q_R_B1Gs in comp GsProtein');
legend_RATES(:,32) = strpad('d/dt q_Gs in comp GsProtein');
legend_RATES(:,34) = strpad('d/dt q_LR_B1 in comp GsProtein');
legend_RATES(:,35) = strpad('d/dt q_LR_B1Gs in comp GsProtein');
legend_RATES(:,36) = strpad('d/dt q_aGs_GTP in comp GsProtein');
legend_RATES(:,37) = strpad('d/dt q_beta_gamma_Gs in comp GsProtein');
legend_RATES(:,38) = strpad('d/dt q_aGs_GDP in comp GsProtein');
legend_RATES(:,39) = strpad('d/dt q_R_M2 in comp GiProtein');
legend_RATES(:,41) = strpad('d/dt q_R_M2Gi in comp GiProtein');
legend_RATES(:,40) = strpad('d/dt q_Gi in comp GiProtein');
legend_RATES(:,42) = strpad('d/dt q_LR_M2 in comp GiProtein');
legend_RATES(:,43) = strpad('d/dt q_LR_M2Gi in comp GiProtein');
legend_RATES(:,44) = strpad('d/dt q_aGi_GTP in comp GiProtein');
legend_RATES(:,45) = strpad('d/dt q_beta_gamma_Gi in comp GiProtein');
legend_RATES(:,46) = strpad('d/dt q_aGi_GDP in comp GiProtein');
legend_states = legend_states';
legend_algebraic = legend_algebraic';
legend_RATES = legend_RATES';
legend_constants = legend_constants';
end
function [states, constants] = initConsts(Gmult,igs)
voi = 0; constants = []; states = []; algebraic = [];
% ---------------------- begin kappa and K
constants(:,1) = 3.92757e+06;
constants(:,2) = 0.000820038;
constants(:,3) = 310777;
constants(:,4) = 0.0838605;
constants(:,5) = 2.41805e+08;
constants(:,6) = 2.81169e-17;
constants(:,7) = 30855;
constants(:,8) = 0.118673;
constants(:,9) = 825.588;
constants(:,10) = 5856.59;
constants(:,11) = 455682;
constants(:,12) = 1648.62;
constants(:,13) = 1083.42;
constants(:,14) = 16.7671;
constants(:,15) = 3.64758;
constants(:,16) = 3.34378e+06;
constants(:,17) = 2855.54;
constants(:,18) = 68.8777;
constants(:,19) = 14.9839;
constants(:,20) = 2.83927e+06;
constants(:,21) = 2.19297;
constants(:,22) = 0.00738313e6;
constants(:,23) = 0.235482;
constants(:,24) = 2.56289e+06;
constants(:,25) = 0.903114;
constants(:,26) = 9.47787e-05;
constants(:,27) = 1.65719;
constants(:,28) = 1.93732e+06;
constants(:,29) = 5.4782e-05;
constants(:,30) = 0.0159048;
constants(:,31) = 3.65262;
constants(:,32) = 7.08986;
constants(:,33) = 4.9636;
constants(:,34) = 2.94648;
constants(:,35) = 0.063794;
constants(:,36) = 0.103389;
constants(:,37) = 1.72199;
constants(:,38) = 1.22478;
constants(:,39) = 0.0159048;
constants(:,40) = 0.0198138;
constants(:,41) = 35.211;
constants(:,42) = 0.0987584;
constants(:,43) = 1.15389e-05;
constants(:,44) = 0.0140332;
constants(:,45) = 17.6328;
constants(:,46) = 0.000100127;
constants(:,47) = 0.0569202;
constants(:,48) = 0.00443997;
constants(:,49) = 175.673;
constants(:,50) = 1.31878;
constants(:,51) = 885.438;
constants(:,52) = 494.115;
constants(:,53) = 0.0365205;
constants(:,54) = 0.00814968;
constants(:,55) = 36.3124;
constants(:,56) = 1.55311;
constants(:,57) = 335.945;
constants(:,58) = 120.284;
constants(:,59) = 2.55469e-05;
constants(:,60) = 0.015617;
constants(:,61) = 0.000229788;
constants(:,62) = 0.000237839;
constants(:,63) = 0.00221913;
% ---------------------- end kappa and K
constants(:,64) = 34.4;
constants(:,65) = 500;
constants(:,66) = 1.5e-4;
constants(:,67) = 3e-4;
constants(:,68) = 0.25e-4;
constants(:,69) = 1.8e-4;
constants(:,70) = 0e-1;
constants(:,71) = 1e-5;
% **************************q_init concentrations begin
constants(:,72) = 190; %190 % q_ATP_init
constants(:,73) = 1.889E-03;
constants(:,74) = 1e-18; % cAMP_init 3.212E-02;
constants(:,75) = 0;
constants(:,76) = 3.800E-05;
constants(:,77) = 0;
constants(:,78) = 0;
constants(:,79) = 0;
constants(:,80) = 0;
constants(:,81) = 0;
constants(:,82) = 1.482E-03;
constants(:,83) = 0;
constants(:,84) = 0;
constants(:,85) = 3.80E-02;
constants(:,86) = 0;
constants(:,87) = 0;
constants(:,88) = 0;
constants(:,89) = 0;
constants(:,90) = 0.0004579000; % q_R_B1_init
constants(:,91) = 0.1455400000e0; % q_Gs_init
constants(:,92) = 0;
constants(:,93) = 0;
constants(:,94) = 0;
constants(:,95) = 0.00072;
constants(:,96) = 0.00836e0; % q_Gi_init
constants(:,97) = 0;
constants(:,98) = 0;
constants(:,99) = 0;
constants(:,100) = 0;
constants(:,101) = 0;
constants(:,102) = 570;
constants(:,103) = 0;
constants(:,104) = 0;
% ************************q_init concentrations end
states(:,1) = 1e-16;
states(:,2) = 1e-16;
states(:,3) = 1e-16;
states(:,4) = 1e-16;
states(:,5) = 1e-16;
states(:,6) = 1e-16;
states(:,7) = 1e-16;
states(:,8) = 1e-16;
states(:,9) = 1e-16;
states(:,10) = 1e-16;
states(:,11) = 1e-16;
states(:,12) = 1e-16;
states(:,13) = 1e-16;
states(:,14) = 1e-16;
states(:,15) = 1e-16;
states(:,16) = 1e-16;
states(:,17) = 1e-16;
states(:,18) = 1e-16;
states(:,19) = 0;
states(:,20) = 0;
states(:,21) = 0;
states(:,22) = 0;
states(:,23) = 0;
states(:,24) = 0;
states(:,25) = 1e-16;
states(:,26) = 1e-16;
states(:,27) = 1e-16;
states(:,28) = 1e-16;
states(:,29) = 1e-16;
states(:,30) = 1e-16;
states(:,31) = 1e-16;
states(:,32) = 1e-16;
states(:,33) = 1e-16;
states(:,34) = 1e-16;
states(:,35) = 1e-16;
states(:,36) = 1e-16;
states(:,37) = 1e-16;
states(:,38) = 1e-16;
constants(:,105) = 1e-16;
states(:,39) = 1e-16;
states(:,40) = 1e-16;
states(:,41) = 1e-16;
states(:,42) = 1e-16;
states(:,43) = 1e-16;
states(:,44) = 1e-16;
states(:,45) = 1e-16;
states(:,46) = 1e-16;
constants(:,106) = 1e-16;
constants(:,107) = 8.31;
constants(:,108) = 310;
constants(:,109) = 96485;
constants(:,110) = 38.0;
constants(:,111) = 34.4;
constants(:,112) = 34.4;
constants(:,113) = constants(:,70)./constants(:,69);
constants(:,114) = constants(:,102)+constants(:,105)+constants(:,106);
% constants(:,116) = 0.000000;
% constants(:,117) = 0.000000;
% constants(:,118) = 0.000000;
% constants(:,119) = 0.000000;
% constants(:,120) = 0.000000;
% constants(:,121) = 0.000000;
% constants(:,122) = 0.000000;
% constants(:,123) = 0.000000;
constants(:,115) = constants(:,107).*constants(:,108).*log( constants(:,61).*constants(:,114));
constants(:,116) = constants(:,107).*constants(:,108).*log( constants(:,61).*constants(:,114));
if (isempty(states)), warning('Initial values for states not set'); end
% multiply constant values by Gmult, with indices given in igs
for j=1:length(igs)
constants(:,igs(j)) = constants(:,igs(j))*Gmult;
end
end
function [rates, algebraic] = computeRates(voi, states, constants,p)
disp(voi)
global algebraicVariableCount;
statesSize = size(states);
statesColumnCount = statesSize(2);
if ( statesColumnCount == 1)
states = states';
algebraic = zeros(1, algebraicVariableCount);
utilOnes = 1;
else
statesRowCount = statesSize(1);
algebraic = zeros(statesRowCount, algebraicVariableCount);
rates = zeros(statesRowCount, statesColumnCount);
utilOnes = ones(statesRowCount, 1);
end
rates(:,31) = 0; % avoid confusion. %constants(:,117);
rates(:,33) = 0; % avoid confusion. %constants(:,117);
rates(:,34) = 0; % avoid confusion. %constants(:,118);
rates(:,35) = 0; % avoid confusion. %constants(:,119);
rates(:,39) = 0; % avoid confusion. %constants(:,120);
rates(:,41) = 0; % avoid confusion. %constants(:,121);
rates(:,42) = 0; % avoid confusion. %constants(:,122);
rates(:,43) = 0; % avoid confusion. %constants(:,123);
algebraic(:,24) = constants(:,90)+states(:,20)+states(:,31);
algebraic(:,51) = constants(:,107).*constants(:,108).*log( constants(:,48).*algebraic(:,24));
algebraic(:,26) = constants(:,91)+states(:,21)+states(:,32);
algebraic(:,56) = constants(:,107).*constants(:,108).*log( constants(:,49).*algebraic(:,26));
algebraic(:,31) = constants(:,93)+states(:,23)+states(:,33);
algebraic(:,73) = constants(:,107).*constants(:,108).*log( constants(:,51).*algebraic(:,31));
algebraic(:,85) = constants(:,13).*exp((algebraic(:,51)+algebraic(:,56))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,73)./( constants(:,107).*constants(:,108)));
%algebraic(:,1) = piecewise({voi<constants(:,66)&voi>constants(:,66) - constants(:,69), constants(:,71)+ constants(:,113).*((voi - constants(:,66))+constants(:,69)) , voi>=constants(:,66)&voi<constants(:,66)+constants(:,68), constants(:,70)+constants(:,71) , voi<constants(:,66)+constants(:,69)+constants(:,68)&voi>=constants(:,66)+constants(:,68), constants(:,71)+ - constants(:,113).*(((voi - constants(:,66)) - constants(:,69)) - constants(:,68)) }, constants(:,71));
if voi > p.stimStart(1)
j = 10;
end
m = p.stimMag(1)/p.tRamp(1);
algebraic(:,1) = piecewise({voi<p.stimStart(1)&voi>p.stimStart(1) - p.tRamp(1), p.stimHolding(1)+ m.*((voi - p.stimStart(1))+p.tRamp(1)) , voi>=p.stimStart(1)&voi<p.stimStart(1)+p.stimDur(1), p.stimMag(1)+p.stimHolding(1) , voi<p.stimStart(1)+p.tRamp(1)+p.stimDur(1)&voi>=p.stimStart(1)+p.stimDur(1), p.stimHolding(1)+ - m.*(((voi - p.stimStart(1)) - p.tRamp(1)) - p.stimDur(1)) }, p.stimHolding(1));
algebraic(:,21) = algebraic(:,1)+states(:,19);
algebraic(:,43) = constants(:,107).*constants(:,108).*log( constants(:,47).*algebraic(:,21));
algebraic(:,34) = constants(:,94)+states(:,24)+states(:,35);
algebraic(:,79) = constants(:,107).*constants(:,108).*log( constants(:,52).*algebraic(:,34));
algebraic(:,90) = constants(:,14).*exp((algebraic(:,73)+algebraic(:,43))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,79)./( constants(:,107).*constants(:,108)));
rates(:,23) = algebraic(:,85) - algebraic(:,90);
algebraic(:,30) = constants(:,92)+states(:,22)+states(:,34);
algebraic(:,64) = constants(:,107).*constants(:,108).*log( constants(:,50).*algebraic(:,30));
algebraic(:,95) = constants(:,15).*exp((algebraic(:,64)+algebraic(:,56))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,79)./( constants(:,107).*constants(:,108)));
rates(:,21) = - algebraic(:,85) - algebraic(:,95);
rates(:,24) = algebraic(:,90)+algebraic(:,95);
algebraic(:,100) = constants(:,16).*exp((algebraic(:,51)+algebraic(:,43))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,64)./( constants(:,107).*constants(:,108)));
rates(:,19) = - algebraic(:,90) - algebraic(:,100);
rates(:,20) = - algebraic(:,85) - algebraic(:,100);
rates(:,22) = - algebraic(:,95)+algebraic(:,100);
algebraic(:,32) = constants(:,95)+states(:,26)+states(:,39);
algebraic(:,74) = constants(:,107).*constants(:,108).*log( constants(:,54).*algebraic(:,32));
algebraic(:,35) = constants(:,96)+states(:,27)+states(:,40);
algebraic(:,80) = constants(:,107).*constants(:,108).*log( constants(:,55).*algebraic(:,35));
algebraic(:,49) = constants(:,98)+states(:,29)+states(:,41);
algebraic(:,91) = constants(:,107).*constants(:,108).*log( constants(:,57).*algebraic(:,49));
algebraic(:,101) = constants(:,17).*exp((algebraic(:,74)+algebraic(:,80))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,91)./( constants(:,107).*constants(:,108)));
%algebraic(:,2) = piecewise({voi<constants(:,67)&voi>constants(:,67) - constants(:,69), constants(:,71)+ constants(:,113).*((voi - constants(:,67))+constants(:,69)) , voi>=constants(:,67)&voi<constants(:,67)+constants(:,68), constants(:,70)+constants(:,71) , voi<constants(:,67)+constants(:,69)+constants(:,68)&voi>=constants(:,67)+constants(:,68), constants(:,71)+ - constants(:,113).*(((voi - constants(:,67)) - constants(:,69)) - constants(:,68)) }, constants(:,71));
m = p.stimMag(2)/p.tRamp(2);
algebraic(:,2) = piecewise({voi<p.stimStart(2)&voi>p.stimStart(2) - p.tRamp(2), p.stimHolding(2)+ m.*((voi - p.stimStart(2))+p.tRamp(2)) , voi>=p.stimStart(2)&voi<p.stimStart(2)+p.stimDur(2), p.stimMag(2)+p.stimHolding(2) , voi<p.stimStart(2)+p.tRamp(2)+p.stimDur(2)&voi>=p.stimStart(2)+p.stimDur(2), p.stimHolding(2)+ - m.*(((voi - p.stimStart(2)) - p.tRamp(2)) - p.stimDur(2)) }, p.stimHolding(2));
algebraic(:,29) = algebraic(:,2)+states(:,25);
algebraic(:,65) = constants(:,107).*constants(:,108).*log( constants(:,53).*algebraic(:,29));
algebraic(:,53) = constants(:,99)+states(:,30)+states(:,43);
algebraic(:,96) = constants(:,107).*constants(:,108).*log( constants(:,58).*algebraic(:,53));
algebraic(:,105) = constants(:,18).*exp((algebraic(:,91)+algebraic(:,65))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,96)./( constants(:,107).*constants(:,108)));
rates(:,29) = algebraic(:,101) - algebraic(:,105);
algebraic(:,57) = constants(:,107).*constants(:,108).*log( constants(:,48).*algebraic(:,24));
algebraic(:,69) = constants(:,107).*constants(:,108).*log( constants(:,49).*algebraic(:,26));
algebraic(:,20) = constants(:,79)+states(:,14)+states(:,36);
algebraic(:,81) = constants(:,107).*constants(:,108).*log( constants(:,42).*algebraic(:,20));
algebraic(:,38) = constants(:,100)+states(:,37);
algebraic(:,87) = constants(:,107).*constants(:,108).*log( constants(:,59).*algebraic(:,38));
algebraic(:,97) = constants(:,21).*(exp((algebraic(:,57)+algebraic(:,69)+constants(:,115))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,81)+algebraic(:,87)+algebraic(:,57))./( constants(:,107).*constants(:,108))));
algebraic(:,75) = constants(:,107).*constants(:,108).*log( constants(:,50).*algebraic(:,30));
algebraic(:,102) = constants(:,22).*(exp((algebraic(:,75)+algebraic(:,69)+constants(:,115))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,81)+algebraic(:,87)+algebraic(:,75))./( constants(:,107).*constants(:,108))));
algebraic(:,52) = constants(:,101)+states(:,38);
algebraic(:,92) = constants(:,107).*constants(:,108).*log( constants(:,60).*algebraic(:,52));
algebraic(:,106) = constants(:,23).*(exp(algebraic(:,81)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,92)+constants(:,115))./( constants(:,107).*constants(:,108))));
rates(:,36) = (algebraic(:,97)+algebraic(:,102)) - algebraic(:,106);
algebraic(:,39) = constants(:,97)+states(:,28)+states(:,42);
algebraic(:,86) = constants(:,107).*constants(:,108).*log( constants(:,56).*algebraic(:,39));
algebraic(:,109) = constants(:,19).*exp((algebraic(:,86)+algebraic(:,80))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,96)./( constants(:,107).*constants(:,108)));
rates(:,27) = - algebraic(:,101) - algebraic(:,109);
rates(:,30) = algebraic(:,105)+algebraic(:,109);
algebraic(:,110) = constants(:,24).*(exp((algebraic(:,92)+algebraic(:,87))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,69)./( constants(:,107).*constants(:,108))));
rates(:,32) = ( - algebraic(:,97) - algebraic(:,102))+algebraic(:,110);
rates(:,37) = (algebraic(:,97)+algebraic(:,102)) - algebraic(:,110);
rates(:,38) = algebraic(:,106) - algebraic(:,110);
algebraic(:,113) = constants(:,20).*exp((algebraic(:,74)+algebraic(:,65))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,86)./( constants(:,107).*constants(:,108)));
rates(:,25) = - algebraic(:,105) - algebraic(:,113);
rates(:,26) = - algebraic(:,101) - algebraic(:,113);
rates(:,28) = - algebraic(:,109)+algebraic(:,113);
algebraic(:,82) = constants(:,107).*constants(:,108).*log( constants(:,54).*algebraic(:,32));
algebraic(:,88) = constants(:,107).*constants(:,108).*log( constants(:,55).*algebraic(:,35));
algebraic(:,25) = constants(:,87)+states(:,16)+states(:,44);
algebraic(:,98) = constants(:,107).*constants(:,108).*log( constants(:,44).*algebraic(:,25));
algebraic(:,60) = constants(:,103)+states(:,45);
algebraic(:,103) = constants(:,107).*constants(:,108).*log( constants(:,62).*algebraic(:,60));
algebraic(:,111) = constants(:,25).*(exp((algebraic(:,82)+algebraic(:,88)+constants(:,116))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,98)+algebraic(:,103)+algebraic(:,82))./( constants(:,107).*constants(:,108))));
algebraic(:,93) = constants(:,107).*constants(:,108).*log( constants(:,56).*algebraic(:,39));
algebraic(:,114) = constants(:,26).*(exp((algebraic(:,93)+algebraic(:,88)+constants(:,116))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,98)+algebraic(:,103)+algebraic(:,93))./( constants(:,107).*constants(:,108))));
algebraic(:,76) = constants(:,104)+states(:,46);
algebraic(:,107) = constants(:,107).*constants(:,108).*log( constants(:,63).*algebraic(:,76));
algebraic(:,116) = constants(:,27).*(exp(algebraic(:,98)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,107)+constants(:,116))./( constants(:,107).*constants(:,108))));
rates(:,44) = (algebraic(:,111)+algebraic(:,114)) - algebraic(:,116);
algebraic(:,118) = constants(:,28).*(exp((algebraic(:,107)+algebraic(:,103))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,88)./( constants(:,107).*constants(:,108))));
rates(:,40) = ( - algebraic(:,111) - algebraic(:,114))+algebraic(:,118);
rates(:,45) = (algebraic(:,111)+algebraic(:,114)) - algebraic(:,118);
rates(:,46) = algebraic(:,116) - algebraic(:,118);
algebraic(:,3) = constants(:,72)+states(:,1);
algebraic(:,33) = constants(:,107).*constants(:,108).*log( constants(:,29).*algebraic(:,3));
algebraic(:,5) = constants(:,73)+states(:,3);
algebraic(:,42) = constants(:,107).*constants(:,108).*log( constants(:,31).*algebraic(:,5));
algebraic(:,6) = constants(:,75)+states(:,4);
algebraic(:,50) = constants(:,107).*constants(:,108).*log( constants(:,32).*algebraic(:,6));
algebraic(:,120) = constants(:,1).*(exp((algebraic(:,42)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,50)./( constants(:,107).*constants(:,108))));
algebraic(:,4) = constants(:,74)+states(:,2);
algebraic(:,36) = constants(:,107).*constants(:,108).*log( constants(:,30).*algebraic(:,4));
algebraic(:,28) = constants(:,89)+states(:,18);
algebraic(:,119) = constants(:,107).*constants(:,108).*log( constants(:,46).*algebraic(:,28));
algebraic(:,121) = constants(:,2).*(exp(algebraic(:,50)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,42)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
rates(:,4) = algebraic(:,120) - algebraic(:,121);
algebraic(:,7) = constants(:,80)+states(:,5);
algebraic(:,55) = constants(:,107).*constants(:,108).*log( constants(:,33).*algebraic(:,7));
algebraic(:,8) = constants(:,81)+states(:,6);
algebraic(:,63) = constants(:,107).*constants(:,108).*log( constants(:,34).*algebraic(:,8));
algebraic(:,122) = constants(:,3).*(exp((algebraic(:,55)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,63)./( constants(:,107).*constants(:,108))));
algebraic(:,123) = constants(:,4).*(exp(algebraic(:,63)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,55)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
rates(:,6) = algebraic(:,122) - algebraic(:,123);
algebraic(:,9) = constants(:,77)+states(:,7);
algebraic(:,72) = constants(:,107).*constants(:,108).*log( constants(:,35).*algebraic(:,9));
algebraic(:,10) = constants(:,78)+states(:,8);
algebraic(:,78) = constants(:,107).*constants(:,108).*log( constants(:,36).*algebraic(:,10));
algebraic(:,124) = constants(:,5).*(exp((algebraic(:,72)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,78)./( constants(:,107).*constants(:,108))));
rates(:,1) = ( - algebraic(:,120) - algebraic(:,124)) - algebraic(:,122);
algebraic(:,125) = constants(:,6).*(exp(algebraic(:,78)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,72)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
rates(:,8) = algebraic(:,124) - algebraic(:,125);
algebraic(:,11) = constants(:,82)+states(:,9);
algebraic(:,84) = constants(:,107).*constants(:,108).*log( constants(:,37).*algebraic(:,11));
algebraic(:,14) = constants(:,84)+states(:,10);
algebraic(:,89) = constants(:,107).*constants(:,108).*log( constants(:,38).*algebraic(:,14));
algebraic(:,126) = constants(:,7).*(exp((algebraic(:,84)+algebraic(:,36))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,89)./( constants(:,107).*constants(:,108))));
rates(:,2) = (algebraic(:,121)+algebraic(:,125)+algebraic(:,123)) - algebraic(:,126);
algebraic(:,108) = constants(:,107).*constants(:,108).*log( constants(:,42).*algebraic(:,20));
algebraic(:,127) = constants(:,10).*(exp((algebraic(:,42)+algebraic(:,108))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,55)./( constants(:,107).*constants(:,108))));
rates(:,14) = - algebraic(:,127);
rates(:,5) = (algebraic(:,127) - algebraic(:,122))+algebraic(:,123);
algebraic(:,22) = constants(:,76)+states(:,15);
algebraic(:,112) = constants(:,107).*constants(:,108).*log( constants(:,43).*algebraic(:,22));
algebraic(:,129) = constants(:,11).*(exp((algebraic(:,112)+algebraic(:,42))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,72)./( constants(:,107).*constants(:,108))));
rates(:,15) = - algebraic(:,129);
rates(:,7) = (algebraic(:,129)+algebraic(:,125)) - algebraic(:,124);
algebraic(:,15) = constants(:,86)+states(:,11);
algebraic(:,94) = constants(:,107).*constants(:,108).*log( constants(:,39).*algebraic(:,15));
algebraic(:,128) = constants(:,8).*(exp(algebraic(:,89)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,84)+algebraic(:,94))./( constants(:,107).*constants(:,108))));
rates(:,10) = algebraic(:,126) - algebraic(:,128);
rates(:,11) = algebraic(:,128);
algebraic(:,115) = constants(:,107).*constants(:,108).*log( constants(:,44).*algebraic(:,25));
algebraic(:,27) = constants(:,88)+states(:,17);
algebraic(:,117) = constants(:,107).*constants(:,108).*log( constants(:,45).*algebraic(:,27));
algebraic(:,131) = constants(:,12).*(exp((algebraic(:,42)+algebraic(:,115))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,117)./( constants(:,107).*constants(:,108))));
rates(:,3) = (((algebraic(:,121) - algebraic(:,120)) - algebraic(:,127)) - algebraic(:,129)) - algebraic(:,131);
algebraic(:,19) = constants(:,83)+states(:,13);
algebraic(:,104) = constants(:,107).*constants(:,108).*log( constants(:,41).*algebraic(:,19));
algebraic(:,17) = constants(:,85)+states(:,12);
algebraic(:,99) = constants(:,107).*constants(:,108).*log( constants(:,40).*algebraic(:,17));
algebraic(:,130) = constants(:,9).*(exp((algebraic(:,84)+algebraic(:,99))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,104)./( constants(:,107).*constants(:,108))));
rates(:,9) = (algebraic(:,128) - algebraic(:,126)) - algebraic(:,130);
rates(:,12) = - algebraic(:,130);
rates(:,13) = algebraic(:,130);
rates(:,16) = - algebraic(:,131);
rates(:,17) = algebraic(:,131);
rates(:,18) = algebraic(:,131);
rates = rates';
end
% Calculate algebraic variables
function algebraic = computeAlgebraic(algebraic, constants, states, voi, p)
statesSize = size(states);
statesColumnCount = statesSize(2);
if ( statesColumnCount == 1)
states = states';
utilOnes = 1;
else
statesRowCount = statesSize(1);
utilOnes = ones(statesRowCount, 1);
end
algebraic(:,24) = constants(:,90)+states(:,20)+states(:,31);
algebraic(:,51) = constants(:,107).*constants(:,108).*log( constants(:,48).*algebraic(:,24));
algebraic(:,26) = constants(:,91)+states(:,21)+states(:,32);
algebraic(:,56) = constants(:,107).*constants(:,108).*log( constants(:,49).*algebraic(:,26));
algebraic(:,31) = constants(:,93)+states(:,23)+states(:,33);
algebraic(:,73) = constants(:,107).*constants(:,108).*log( constants(:,51).*algebraic(:,31));
algebraic(:,85) = constants(:,13).*exp((algebraic(:,51)+algebraic(:,56))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,73)./( constants(:,107).*constants(:,108)));
%algebraic(:,1) = piecewise({voi<constants(:,66)&voi>constants(:,66) - constants(:,69), constants(:,71)+ constants(:,113).*((voi - constants(:,66))+constants(:,69)) , voi>=constants(:,66)&voi<constants(:,66)+constants(:,68), constants(:,70)+constants(:,71) , voi<constants(:,66)+constants(:,69)+constants(:,68)&voi>=constants(:,66)+constants(:,68), constants(:,71)+ - constants(:,113).*(((voi - constants(:,66)) - constants(:,69)) - constants(:,68)) }, constants(:,71));
m = p.stimMag(1)/p.tRamp(1);
algebraic(:,1) = piecewise({voi<p.stimStart(1)&voi>p.stimStart(1) - p.tRamp(1), p.stimHolding(1)+ m.*((voi - p.stimStart(1))+p.tRamp(1)) , voi>=p.stimStart(1)&voi<p.stimStart(1)+p.stimDur(1), p.stimMag(1)+p.stimHolding(1) , voi<p.stimStart(1)+p.tRamp(1)+p.stimDur(1)&voi>=p.stimStart(1)+p.stimDur(1), p.stimHolding(1)+ - m.*(((voi - p.stimStart(1)) - p.tRamp(1)) - p.stimDur(1)) }, p.stimHolding(1));
algebraic(:,21) = algebraic(:,1)+states(:,19);
algebraic(:,43) = constants(:,107).*constants(:,108).*log( constants(:,47).*algebraic(:,21));
algebraic(:,34) = constants(:,94)+states(:,24)+states(:,35);
algebraic(:,79) = constants(:,107).*constants(:,108).*log( constants(:,52).*algebraic(:,34));
algebraic(:,90) = constants(:,14).*exp((algebraic(:,73)+algebraic(:,43))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,79)./( constants(:,107).*constants(:,108)));
algebraic(:,30) = constants(:,92)+states(:,22)+states(:,34);
algebraic(:,64) = constants(:,107).*constants(:,108).*log( constants(:,50).*algebraic(:,30));
algebraic(:,95) = constants(:,15).*exp((algebraic(:,64)+algebraic(:,56))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,79)./( constants(:,107).*constants(:,108)));
algebraic(:,100) = constants(:,16).*exp((algebraic(:,51)+algebraic(:,43))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,64)./( constants(:,107).*constants(:,108)));
algebraic(:,32) = constants(:,95)+states(:,26)+states(:,39);
algebraic(:,74) = constants(:,107).*constants(:,108).*log( constants(:,54).*algebraic(:,32));
algebraic(:,35) = constants(:,96)+states(:,27)+states(:,40);
algebraic(:,80) = constants(:,107).*constants(:,108).*log( constants(:,55).*algebraic(:,35));
algebraic(:,49) = constants(:,98)+states(:,29)+states(:,41);
algebraic(:,91) = constants(:,107).*constants(:,108).*log( constants(:,57).*algebraic(:,49));
algebraic(:,101) = constants(:,17).*exp((algebraic(:,74)+algebraic(:,80))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,91)./( constants(:,107).*constants(:,108)));
% algebraic(:,2) = piecewise({voi<constants(:,67)&voi>constants(:,67) - constants(:,69), constants(:,71)+ constants(:,113).*((voi - constants(:,67))+constants(:,69)) , voi>=constants(:,67)&voi<constants(:,67)+constants(:,68), constants(:,70)+constants(:,71) , voi<constants(:,67)+constants(:,69)+constants(:,68)&voi>=constants(:,67)+constants(:,68), constants(:,71)+ - constants(:,113).*(((voi - constants(:,67)) - constants(:,69)) - constants(:,68)) }, constants(:,71));
m = p.stimMag(2)/p.tRamp(2);
algebraic(:,2) = piecewise({voi<p.stimStart(2)&voi>p.stimStart(2) - p.tRamp(2), p.stimHolding(2)+ m.*((voi - p.stimStart(2))+p.tRamp(2)) , voi>=p.stimStart(2)&voi<p.stimStart(2)+p.stimDur(2), p.stimMag(2)+p.stimHolding(2) , voi<p.stimStart(2)+p.tRamp(2)+p.stimDur(2)&voi>=p.stimStart(2)+p.stimDur(2), p.stimHolding(2)+ - m.*(((voi - p.stimStart(2)) - p.tRamp(2)) - p.stimDur(2)) }, p.stimHolding(2));
algebraic(:,29) = algebraic(:,2)+states(:,25);
algebraic(:,65) = constants(:,107).*constants(:,108).*log( constants(:,53).*algebraic(:,29));
algebraic(:,53) = constants(:,99)+states(:,30)+states(:,43);
algebraic(:,96) = constants(:,107).*constants(:,108).*log( constants(:,58).*algebraic(:,53));
algebraic(:,105) = constants(:,18).*exp((algebraic(:,91)+algebraic(:,65))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,96)./( constants(:,107).*constants(:,108)));
algebraic(:,57) = constants(:,107).*constants(:,108).*log( constants(:,48).*algebraic(:,24));
algebraic(:,69) = constants(:,107).*constants(:,108).*log( constants(:,49).*algebraic(:,26));
algebraic(:,20) = constants(:,79)+states(:,14)+states(:,36);
algebraic(:,81) = constants(:,107).*constants(:,108).*log( constants(:,42).*algebraic(:,20));
algebraic(:,38) = constants(:,100)+states(:,37);
algebraic(:,87) = constants(:,107).*constants(:,108).*log( constants(:,59).*algebraic(:,38));
algebraic(:,97) = constants(:,21).*(exp((algebraic(:,57)+algebraic(:,69)+constants(:,115))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,81)+algebraic(:,87)+algebraic(:,57))./( constants(:,107).*constants(:,108))));
algebraic(:,75) = constants(:,107).*constants(:,108).*log( constants(:,50).*algebraic(:,30));
algebraic(:,102) = constants(:,22).*(exp((algebraic(:,75)+algebraic(:,69)+constants(:,115))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,81)+algebraic(:,87)+algebraic(:,75))./( constants(:,107).*constants(:,108))));
algebraic(:,52) = constants(:,101)+states(:,38);
algebraic(:,92) = constants(:,107).*constants(:,108).*log( constants(:,60).*algebraic(:,52));
algebraic(:,106) = constants(:,23).*(exp(algebraic(:,81)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,92)+constants(:,115))./( constants(:,107).*constants(:,108))));
algebraic(:,39) = constants(:,97)+states(:,28)+states(:,42);
algebraic(:,86) = constants(:,107).*constants(:,108).*log( constants(:,56).*algebraic(:,39));
algebraic(:,109) = constants(:,19).*exp((algebraic(:,86)+algebraic(:,80))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,96)./( constants(:,107).*constants(:,108)));
algebraic(:,110) = constants(:,24).*(exp((algebraic(:,92)+algebraic(:,87))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,69)./( constants(:,107).*constants(:,108))));
algebraic(:,113) = constants(:,20).*exp((algebraic(:,74)+algebraic(:,65))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,86)./( constants(:,107).*constants(:,108)));
algebraic(:,82) = constants(:,107).*constants(:,108).*log( constants(:,54).*algebraic(:,32));
algebraic(:,88) = constants(:,107).*constants(:,108).*log( constants(:,55).*algebraic(:,35));
algebraic(:,25) = constants(:,87)+states(:,16)+states(:,44);
algebraic(:,98) = constants(:,107).*constants(:,108).*log( constants(:,44).*algebraic(:,25));
algebraic(:,60) = constants(:,103)+states(:,45);
algebraic(:,103) = constants(:,107).*constants(:,108).*log( constants(:,62).*algebraic(:,60));
algebraic(:,111) = constants(:,25).*(exp((algebraic(:,82)+algebraic(:,88)+constants(:,116))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,98)+algebraic(:,103)+algebraic(:,82))./( constants(:,107).*constants(:,108))));
algebraic(:,93) = constants(:,107).*constants(:,108).*log( constants(:,56).*algebraic(:,39));
algebraic(:,114) = constants(:,26).*(exp((algebraic(:,93)+algebraic(:,88)+constants(:,116))./( constants(:,107).*constants(:,108))) - exp((algebraic(:,98)+algebraic(:,103)+algebraic(:,93))./( constants(:,107).*constants(:,108))));
algebraic(:,76) = constants(:,104)+states(:,46);
algebraic(:,107) = constants(:,107).*constants(:,108).*log( constants(:,63).*algebraic(:,76));
algebraic(:,116) = constants(:,27).*(exp(algebraic(:,98)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,107)+constants(:,116))./( constants(:,107).*constants(:,108))));
algebraic(:,118) = constants(:,28).*(exp((algebraic(:,107)+algebraic(:,103))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,88)./( constants(:,107).*constants(:,108))));
algebraic(:,3) = constants(:,72)+states(:,1);
algebraic(:,33) = constants(:,107).*constants(:,108).*log( constants(:,29).*algebraic(:,3));
algebraic(:,5) = constants(:,73)+states(:,3);
algebraic(:,42) = constants(:,107).*constants(:,108).*log( constants(:,31).*algebraic(:,5));
algebraic(:,6) = constants(:,75)+states(:,4);
algebraic(:,50) = constants(:,107).*constants(:,108).*log( constants(:,32).*algebraic(:,6));
algebraic(:,120) = constants(:,1).*(exp((algebraic(:,42)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,50)./( constants(:,107).*constants(:,108))));
algebraic(:,4) = constants(:,74)+states(:,2);
algebraic(:,36) = constants(:,107).*constants(:,108).*log( constants(:,30).*algebraic(:,4));
algebraic(:,28) = constants(:,89)+states(:,18);
algebraic(:,119) = constants(:,107).*constants(:,108).*log( constants(:,46).*algebraic(:,28));
algebraic(:,121) = constants(:,2).*(exp(algebraic(:,50)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,42)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
algebraic(:,7) = constants(:,80)+states(:,5);
algebraic(:,55) = constants(:,107).*constants(:,108).*log( constants(:,33).*algebraic(:,7));
algebraic(:,8) = constants(:,81)+states(:,6);
algebraic(:,63) = constants(:,107).*constants(:,108).*log( constants(:,34).*algebraic(:,8));
algebraic(:,122) = constants(:,3).*(exp((algebraic(:,55)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,63)./( constants(:,107).*constants(:,108))));
algebraic(:,123) = constants(:,4).*(exp(algebraic(:,63)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,55)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
algebraic(:,9) = constants(:,77)+states(:,7);
algebraic(:,72) = constants(:,107).*constants(:,108).*log( constants(:,35).*algebraic(:,9));
algebraic(:,10) = constants(:,78)+states(:,8);
algebraic(:,78) = constants(:,107).*constants(:,108).*log( constants(:,36).*algebraic(:,10));
algebraic(:,124) = constants(:,5).*(exp((algebraic(:,72)+algebraic(:,33))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,78)./( constants(:,107).*constants(:,108))));
algebraic(:,125) = constants(:,6).*(exp(algebraic(:,78)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,72)+algebraic(:,36)+algebraic(:,119))./( constants(:,107).*constants(:,108))));
algebraic(:,11) = constants(:,82)+states(:,9);
algebraic(:,84) = constants(:,107).*constants(:,108).*log( constants(:,37).*algebraic(:,11));
algebraic(:,14) = constants(:,84)+states(:,10);
algebraic(:,89) = constants(:,107).*constants(:,108).*log( constants(:,38).*algebraic(:,14));
algebraic(:,126) = constants(:,7).*(exp((algebraic(:,84)+algebraic(:,36))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,89)./( constants(:,107).*constants(:,108))));
algebraic(:,108) = constants(:,107).*constants(:,108).*log( constants(:,42).*algebraic(:,20));
algebraic(:,127) = constants(:,10).*(exp((algebraic(:,42)+algebraic(:,108))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,55)./( constants(:,107).*constants(:,108))));
algebraic(:,22) = constants(:,76)+states(:,15);
algebraic(:,112) = constants(:,107).*constants(:,108).*log( constants(:,43).*algebraic(:,22));
algebraic(:,129) = constants(:,11).*(exp((algebraic(:,112)+algebraic(:,42))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,72)./( constants(:,107).*constants(:,108))));
algebraic(:,15) = constants(:,86)+states(:,11);
algebraic(:,94) = constants(:,107).*constants(:,108).*log( constants(:,39).*algebraic(:,15));
algebraic(:,128) = constants(:,8).*(exp(algebraic(:,89)./( constants(:,107).*constants(:,108))) - exp((algebraic(:,84)+algebraic(:,94))./( constants(:,107).*constants(:,108))));
algebraic(:,115) = constants(:,107).*constants(:,108).*log( constants(:,44).*algebraic(:,25));
algebraic(:,27) = constants(:,88)+states(:,17);
algebraic(:,117) = constants(:,107).*constants(:,108).*log( constants(:,45).*algebraic(:,27));
algebraic(:,131) = constants(:,12).*(exp((algebraic(:,42)+algebraic(:,115))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,117)./( constants(:,107).*constants(:,108))));
algebraic(:,19) = constants(:,83)+states(:,13);
algebraic(:,104) = constants(:,107).*constants(:,108).*log( constants(:,41).*algebraic(:,19));
algebraic(:,17) = constants(:,85)+states(:,12);
algebraic(:,99) = constants(:,107).*constants(:,108).*log( constants(:,40).*algebraic(:,17));
algebraic(:,130) = constants(:,9).*(exp((algebraic(:,84)+algebraic(:,99))./( constants(:,107).*constants(:,108))) - exp(algebraic(:,104)./( constants(:,107).*constants(:,108))));
algebraic(:,12) = algebraic(:,3)+algebraic(:,6)+algebraic(:,10)+algebraic(:,8);
algebraic(:,13) = algebraic(:,5)+algebraic(:,6)+algebraic(:,9)+algebraic(:,10)+algebraic(:,7)+algebraic(:,8);
algebraic(:,16) = algebraic(:,4)+algebraic(:,14)+states(:,11);
algebraic(:,18) = algebraic(:,4)+algebraic(:,14)+algebraic(:,15) + algebraic(:,3)+algebraic(:,6)+algebraic(:,8) + algebraic(:,10);
algebraic(:,23) = algebraic(:,20)+algebraic(:,7)+algebraic(:,8);
algebraic(:,37) = constants(:,107).*constants(:,108).*log( constants(:,51).*algebraic(:,31));
algebraic(:,40) = algebraic(:,21)+algebraic(:,30)+algebraic(:,34);
algebraic(:,41) = algebraic(:,24)+algebraic(:,31)+algebraic(:,30)+algebraic(:,34);
algebraic(:,44) = algebraic(:,21)+algebraic(:,30)+algebraic(:,34);
algebraic(:,45) = algebraic(:,24)+algebraic(:,30)+algebraic(:,31)+algebraic(:,34);
algebraic(:,46) = algebraic(:,26)+algebraic(:,31)+algebraic(:,34);
algebraic(:,47) = constants(:,107).*constants(:,108).*log( constants(:,52).*algebraic(:,34));
algebraic(:,48) = algebraic(:,24)+algebraic(:,31)+algebraic(:,30)+algebraic(:,34);
algebraic(:,54) = algebraic(:,26)+algebraic(:,31)+algebraic(:,34)+algebraic(:,20)+algebraic(:,52)+algebraic(:,7)+algebraic(:,8);
algebraic(:,58) = algebraic(:,26)+algebraic(:,31)+algebraic(:,34)+algebraic(:,20)+algebraic(:,52);
algebraic(:,59) = constants(:,107).*constants(:,108).*log( constants(:,57).*algebraic(:,49));
algebraic(:,61) = algebraic(:,29)+algebraic(:,39)+algebraic(:,53);
algebraic(:,62) = algebraic(:,32)+algebraic(:,49)+algebraic(:,39)+algebraic(:,53);
algebraic(:,66) = algebraic(:,29)+algebraic(:,39)+algebraic(:,53);
algebraic(:,67) = algebraic(:,32)+algebraic(:,39)+algebraic(:,49)+algebraic(:,53);
algebraic(:,68) = algebraic(:,35)+algebraic(:,49)+algebraic(:,53);
algebraic(:,70) = constants(:,107).*constants(:,108).*log( constants(:,58).*algebraic(:,53));
algebraic(:,71) = algebraic(:,32)+algebraic(:,49)+algebraic(:,39)+algebraic(:,53);
algebraic(:,77) = algebraic(:,35)+algebraic(:,49)+algebraic(:,53)+algebraic(:,25)+algebraic(:,76)+algebraic(:,27);
algebraic(:,83) = algebraic(:,35)+algebraic(:,49)+algebraic(:,53)+algebraic(:,25)+algebraic(:,76);
end
% Compute result of a piecewise function
function x = piecewise(cases, default)
set = [0];
for i = 1:2:length(cases)
if (length(cases{i+1}) == 1)
x(cases{i} & ~set,:) = cases{i+1};
else
x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set);
end
set = set | cases{i};
if(set), break, end
end
if (length(default) == 1)
x(~set,:) = default;
else
x(~set,:) = default(~set);
end
end
% Pad out or shorten strings to a set length
function strout = strpad(strin)
req_length = 60;
insize = size(strin,2);
if insize > req_length
strout = strin(1:req_length);
else
strout = [strin, blanks(req_length - insize)];
end
end
