- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-11-10 10:33:08+13:00
- Desc:
- Fixing figure
- Permanent Source URI:
- https://models.cellml.org/workspace/705/rawfile/9229501f7d14360a0ca80dc890186901c48a0dbf/FCU_adenylylCyclase_Exported.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 = 1e3; %1e3;
[INIT_states, constants] = initConsts(Gmult, igs);
% Set timespan to solve over
voi = [0,6.5];%linspace(0, 0.5,2000);
% voi = [0 1e-3];
% Set numerical accuracy options for ODE solver
options = odeset('RelTol', 1e-07, 'AbsTol', 1e-07, 'MaxStep', 1e-4); % step: 1e-4
% 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,2];
% p.stimEnd = [3,4]; % not used
% p.stimOffset = 0.0; % seconds % not used
p.tRamp = [1,1];
p.stimDur = [1,1];
if false
p.stimMag = p.stimHolding(1);
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
plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('states\nGmult=%g', Gmult),1)
plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('algebraic\nGmult=%g', Gmult),ceil(sqrt(length(i_alg))))
end
% debugging for q_LR_B1
labels = {'q_LR_B1 in comp env',...
'q_LR_B1 in comp LRGbinding_B1AR',...
'q_LR_B1 in comp GsProtein'};
if false
[~, 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_st))))
end
ids = [80,85,44,39,55];
figure,
for i = 1:2
subplot(2,2,i)
plot(voi(30:end), algebraic(30:end,ids(i)));
xlabel('time (s)');
legend(legend_algebraic(ids(i),:));
end
subplot(2,2,3)
plot(voi(30:end), -algebraic(30:end,ids(1))+algebraic(30:end,ids(2)));
xlabel('time (s)');
title(legend_states(22,:));
subplot(2,2,4)
for j=1:3
plot(voi(30:end), algebraic(30:end,ids(j+2))); hold on;
end
legend(legend_algebraic(ids(3),:),legend_algebraic(ids(4),:),legend_algebraic(ids(5),:));
xlabel('time (s)');
% debugging aGsGTP spike
if true
istart = 1; %55250;
iend = length(voi); %57850;
figure,
subplot(2,2,1)
plot(voi(istart:iend), states(istart:iend,14),'.'); % hold on;
title('aGsGTP cAMP');
subplot(2,2,2)
istart2 = 1;
iend2 = length(voi);
plot(voi(istart2:iend2), states(istart2:iend2,36),'.');
% legend('comp cAMP','comp GsProtein');
title('aGsGTP GsProtein');
subplot(2,2,3)
plot(voi(istart:iend), rates(14,istart:iend),'.'); % hold on;
title('d/dtaGsGTP cAMP');
subplot(2,2,4)
plot(voi(istart:iend), rates(36,istart:iend),'.');
% legend('comp cAMP','comp GsProtein');
title('d/dt aGsGTP GsProtein');
end
end
function [] = plot_selected(ids,x,y,legend_x,legend_y,titlestr,ns)
figure();
% plot stimuli
for i = 1:length(ids)
subplot(ns,ns,i)
plot(x(30:end), y(30:end,ids(i)));
xlabel('time (s)');
l = legend(legend_y(ids(i),:));
set(l,'Interpreter','none');
end
suptitle(titlestr)
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('q_ATP_init in comp env');
legend_constants(:,66) = strpad('q_AC_init in comp env');
legend_constants(:,67) = strpad('q_cAMP_init in comp env');
legend_constants(:,68) = strpad('q_AC_ATP_init in comp env');
legend_constants(:,69) = strpad('q_FSK_init in comp env');
legend_constants(:,70) = strpad('q_FSK_AC_init in comp env');
legend_constants(:,71) = strpad('q_FSK_AC_ATP_init in comp env');
legend_constants(:,72) = strpad('q_aGs_GTP_init in comp env');
legend_constants(:,73) = strpad('q_aGs_GTP_AC_init in comp env');
legend_constants(:,74) = strpad('q_aGs_GTP_AC_ATP_init in comp env');
legend_constants(:,75) = strpad('q_PDE_init in comp env');
legend_constants(:,76) = strpad('q_PDEinh_init in comp env');
legend_constants(:,77) = strpad('q_PDE_cAMP_init in comp env');
legend_constants(:,78) = strpad('q_IBMX_init in comp env');
legend_constants(:,79) = strpad('q_five_AMP_init in comp env');
legend_constants(:,80) = strpad('q_aGi_GTP_init in comp env');
legend_constants(:,81) = strpad('q_ACinh_init in comp env');
legend_constants(:,82) = strpad('q_PPi_init in comp env');
legend_constants(:,83) = strpad('q_R_B1_init in comp env');
legend_constants(:,84) = strpad('q_Gs_init in comp env');
legend_algebraic(:,1) = strpad('q_L_B1_init in comp env');
legend_constants(:,85) = strpad('q_LR_B1_init in comp env');
legend_constants(:,86) = strpad('q_R_B1Gs_init in comp env');
legend_constants(:,87) = strpad('q_LR_B1Gs_init in comp env');
legend_algebraic(:,2) = strpad('q_L_M2_init in comp env');
legend_constants(:,88) = strpad('q_R_M2_init in comp env');
legend_constants(:,89) = strpad('q_Gi_init in comp env');
legend_constants(:,90) = strpad('q_LR_M2_init in comp env');
legend_constants(:,91) = strpad('q_R_M2Gi_init in comp env');
legend_constants(:,92) = strpad('q_LR_M2Gi_init in comp env');
legend_constants(:,93) = strpad('q_beta_gamma_Gs_init in comp env');
legend_constants(:,94) = strpad('q_aGs_GDP_init in comp env');
legend_constants(:,95) = strpad('q_Pi_init in comp env');
legend_constants(:,96) = strpad('q_beta_gamma_Gi_init in comp env');
legend_constants(:,97) = 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(:,12) = strpad('q_PDE_cAMP in comp env');
legend_algebraic(:,13) = strpad('q_five_AMP in comp env');
legend_algebraic(:,14) = strpad('q_IBMX in comp env');
legend_algebraic(:,15) = strpad('q_PDEinh in comp env');
legend_algebraic(:,16) = strpad('q_aGs_GTP in comp env');
legend_algebraic(:,17) = strpad('q_FSK in comp env');
legend_algebraic(:,19) = strpad('q_aGi_GTP in comp env');
legend_algebraic(:,21) = strpad('q_ACinh in comp env');
legend_algebraic(:,24) = strpad('q_PPi in comp env');
legend_algebraic(:,18) = strpad('q_L_B1 in comp env');
legend_algebraic(:,20) = strpad('q_R_B1 in comp env');
legend_algebraic(:,22) = strpad('q_Gs in comp env');
legend_algebraic(:,26) = strpad('q_LR_B1 in comp env');
legend_algebraic(:,27) = strpad('q_R_B1Gs in comp env');
legend_algebraic(:,31) = strpad('q_LR_B1Gs in comp env');
legend_algebraic(:,23) = strpad('q_L_M2 in comp env');
legend_algebraic(:,25) = strpad('q_R_M2 in comp env');
legend_algebraic(:,28) = strpad('q_Gi in comp env');
legend_algebraic(:,32) = strpad('q_LR_M2 in comp env');
legend_algebraic(:,36) = strpad('q_R_M2Gi in comp env');
legend_algebraic(:,42) = strpad('q_LR_M2Gi in comp env');
legend_algebraic(:,29) = strpad('q_beta_gamma_Gs in comp env');
legend_algebraic(:,34) = strpad('q_aGs_GDP in comp env');
legend_constants(:,106) = strpad('q_Pi in comp env');
legend_algebraic(:,37) = strpad('q_beta_gamma_Gi in comp env');
legend_algebraic(:,46) = 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(:,98) = 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(:,99) = strpad('q_Pi in comp GiProtein');
legend_constants(:,100) = strpad('R in comp constants (J_per_K_per_mol)');
legend_constants(:,101) = strpad('T in comp constants (kelvin)');
legend_constants(:,102) = strpad('F in comp constants (C_per_mol)');
legend_algebraic(:,99) = strpad('v1a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,100) = strpad('v1b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,101) = strpad('v2a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,102) = strpad('v2b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,103) = strpad('v3a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,104) = strpad('v3b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,105) = strpad('v4a in comp cAMP (fmol_per_sec)');
legend_algebraic(:,107) = strpad('v4b in comp cAMP (fmol_per_sec)');
legend_algebraic(:,109) = strpad('v5 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,106) = strpad('v6 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,108) = strpad('v7 in comp cAMP (fmol_per_sec)');
legend_algebraic(:,110) = strpad('vGiAC in comp cAMP (fmol_per_sec)');
legend_algebraic(:,30) = strpad('mu_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,38) = strpad('mu_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,33) = strpad('mu_cAMP in comp cAMP (J_per_mol)');
legend_algebraic(:,43) = strpad('mu_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,95) = strpad('mu_FSK in comp cAMP (J_per_mol)');
legend_algebraic(:,59) = strpad('mu_FSK_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,64) = strpad('mu_FSK_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,92) = strpad('mu_aGs_GTP in comp cAMP (J_per_mol)');
legend_algebraic(:,48) = strpad('mu_aGs_GTP_AC in comp cAMP (J_per_mol)');
legend_algebraic(:,54) = strpad('mu_aGs_GTP_AC_ATP in comp cAMP (J_per_mol)');
legend_algebraic(:,69) = strpad('mu_PDE in comp cAMP (J_per_mol)');
legend_algebraic(:,89) = strpad('mu_PDEinh in comp cAMP (J_per_mol)');
legend_algebraic(:,74) = strpad('mu_PDE_cAMP in comp cAMP (J_per_mol)');
legend_algebraic(:,84) = strpad('mu_IBMX in comp cAMP (J_per_mol)');
legend_algebraic(:,79) = strpad('mu_five_AMP in comp cAMP (J_per_mol)');
legend_algebraic(:,96) = strpad('mu_aGi_GTP in comp cAMP (J_per_mol)');
legend_algebraic(:,97) = strpad('mu_ACinh in comp cAMP (J_per_mol)');
legend_algebraic(:,98) = strpad('mu_PPi in comp cAMP (J_per_mol)');
legend_constants(:,103) = strpad('vol in comp cAMP (pL)');
legend_algebraic(:,39) = strpad('mu_L_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,44) = strpad('mu_R_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,49) = strpad('mu_Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,55) = strpad('mu_LR_B1 in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,60) = strpad('mu_R_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,65) = strpad('mu_LR_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
legend_algebraic(:,70) = strpad('vsig1_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,75) = strpad('vsig2_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,80) = strpad('vsig3_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_algebraic(:,85) = strpad('vsig4_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
legend_constants(:,104) = strpad('vol in comp LRGbinding_B1AR (pL)');
legend_algebraic(:,50) = strpad('mu_L_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,56) = strpad('mu_R_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,61) = strpad('mu_Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,66) = strpad('mu_LR_M2 in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,71) = strpad('mu_R_M2Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,76) = strpad('mu_LR_M2Gi in comp LRGbinding_M2 (J_per_mol)');
legend_algebraic(:,81) = strpad('vsig1_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,86) = strpad('vsig2_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,90) = strpad('vsig3_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_algebraic(:,93) = strpad('vsig4_M2 in comp LRGbinding_M2 (fmol_per_sec)');
legend_constants(:,105) = strpad('vol in comp LRGbinding_M2 (pL)');
legend_algebraic(:,72) = strpad('vact1_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,77) = strpad('vact2_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,82) = strpad('vhyd_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,87) = strpad('vreassoc_Gs in comp GsProtein (fmol_per_sec)');
legend_algebraic(:,40) = strpad('mu_R_B1 in comp GsProtein (J_per_mol)');
legend_algebraic(:,45) = strpad('mu_Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,35) = strpad('mu_R_B1Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,51) = strpad('mu_LR_B1 in comp GsProtein (J_per_mol)');
legend_algebraic(:,41) = strpad('mu_LR_B1Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,57) = strpad('mu_aGs_GTP in comp GsProtein (J_per_mol)');
legend_algebraic(:,62) = strpad('mu_beta_gamma_Gs in comp GsProtein (J_per_mol)');
legend_algebraic(:,67) = strpad('mu_aGs_GDP in comp GsProtein (J_per_mol)');
legend_constants(:,107) = strpad('mu_Pi in comp GsProtein (J_per_mol)');
legend_algebraic(:,83) = strpad('vact1_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,88) = strpad('vact2_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,91) = strpad('vhyd_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,94) = strpad('vreassoc_Gi in comp GiProtein (fmol_per_sec)');
legend_algebraic(:,52) = strpad('mu_R_M2 in comp GiProtein (J_per_mol)');
legend_algebraic(:,58) = strpad('mu_Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,47) = strpad('mu_R_M2Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,63) = strpad('mu_LR_M2 in comp GiProtein (J_per_mol)');
legend_algebraic(:,53) = strpad('mu_LR_M2Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,68) = strpad('mu_aGi_GTP in comp GiProtein (J_per_mol)');
legend_algebraic(:,73) = strpad('mu_beta_gamma_Gi in comp GiProtein (J_per_mol)');
legend_algebraic(:,78) = strpad('mu_aGi_GDP in comp GiProtein (J_per_mol)');
legend_constants(:,108) = strpad('mu_Pi in comp GiProtein (J_per_mol)');
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 = [];
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.00738313;
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;
constants(:,64) = 34.4;
constants(:,65) = 190;
constants(:,66) = 0.0006143333; %*Gmult; % q_AC;
constants(:,67) = 1.889E-03; %*Gmult; % q_camp_init
constants(:,68) = 1e-15;
constants(:,69) = 3.800E-05;
constants(:,70) = 1e-15;
constants(:,71) = 1e-15;
constants(:,72) = 0.0009519;
constants(:,73) = 1e-15;
constants(:,74) = 1e-15;
constants(:,75) = 1.434E-03;
constants(:,76) = 1e-15;
constants(:,77) = 1e-15;
constants(:,78) = 3.80E-02;
constants(:,79) = 1e-15;
constants(:,80) = 1e-15;
constants(:,81) = 1e-15;
constants(:,82) = 0.038;
constants(:,83) = 0.0004579000; %*Gmult; % q_R_B1_init
constants(:,84) = 0.1455400000; %*Gmult; % q_Gs_init
constants(:,85) = 1e-1; % q_LR_B1_init % init, to get a positive mu
constants(:,86) = 1e-15;
constants(:,87) = 1e-15;
constants(:,88) = 0.00072; %*Gmult; % q_R_M2_init
constants(:,89) = 0.00836; %*Gmult; % q_Gi_init
constants(:,90) = 1e-15; %q_LR_M2_init
constants(:,91) = 1e-15;
constants(:,92) = 1e-15;
constants(:,93) = 1e-15;
constants(:,94) = 1e-15;
constants(:,95) = 570;
constants(:,96) = 1e-15;
constants(:,97) = 1e-15;
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) = 1e-16;
states(:,20) = 1e-16;
states(:,21) = 1e-16;
states(:,22) = 1e-16;
states(:,23) = 1e-16;
states(:,24) = 1e-16;
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(:,98) = 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(:,99) = 1e-16;
constants(:,100) = 8.31;
constants(:,101) = 310;
constants(:,102) = 96485;
constants(:,103) = 38.0;
constants(:,104) = 34.4;
constants(:,105) = 34.4;
constants(:,106) = constants(:,95)+constants(:,98)+constants(:,99);
constants(:,108) = 0.000000;
constants(:,109) = 0.000000;
constants(:,110) = 0.000000;
constants(:,111) = 0.000000;
constants(:,112) = 0.000000;
constants(:,113) = 0.000000;
constants(:,114) = 0.000000;
constants(:,115) = 0.000000;
constants(:,107) = constants(:,100).*constants(:,101).*log( constants(:,61).*constants(:,106));
constants(:,108) = constants(:,100).*constants(:,101).*log( constants(:,61).*constants(:,106));
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) = constants(:,108);
rates(:,33) = constants(:,109);
rates(:,34) = constants(:,110);
rates(:,35) = constants(:,111);
rates(:,39) = constants(:,112);
rates(:,41) = constants(:,113);
rates(:,42) = constants(:,114);
rates(:,43) = constants(:,115);
algebraic(:,20) = constants(:,83)+states(:,20)+states(:,31);
algebraic(:,44) = constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
algebraic(:,22) = constants(:,84)+states(:,21)+states(:,32);
algebraic(:,49) = constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
algebraic(:,27) = constants(:,86)+states(:,23)+states(:,33);
algebraic(:,60) = constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
algebraic(:,70) = constants(:,13).*exp((algebraic(:,44)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,60)./( constants(:,100).*constants(:,101)));
if voi > 3
j = 10;
end
% oldStim = piecewise({voi>p.stimStart(1)-1 & voi<=p.stimStart(1), p.stimMag(1).*voi , voi>p.stimStart(1)&voi<p.stimEnd(1), p.stimMag(1)}, p.stimHolding(1)(1));
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(:,18) = algebraic(:,1)+states(:,19);
algebraic(:,39) = constants(:,100).*constants(:,101).*log( constants(:,47).*algebraic(:,18));
algebraic(:,31) = constants(:,87)+states(:,24)+states(:,35);
algebraic(:,65) = constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
algebraic(:,75) = constants(:,14).*exp((algebraic(:,60)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
rates(:,23) = algebraic(:,70) - algebraic(:,75);
algebraic(:,26) = constants(:,85)+states(:,22)+states(:,34);
algebraic(:,55) = constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
algebraic(:,80) = constants(:,15).*exp((algebraic(:,55)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
rates(:,21) = - algebraic(:,70) - algebraic(:,80);
rates(:,24) = algebraic(:,75)+algebraic(:,80);
algebraic(:,40) = constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
algebraic(:,45) = constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
algebraic(:,16) = constants(:,72)+states(:,14)+states(:,36);
algebraic(:,57) = constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
algebraic(:,29) = constants(:,93)+states(:,37);
algebraic(:,62) = constants(:,100).*constants(:,101).*log( constants(:,59).*algebraic(:,29));
algebraic(:,72) = constants(:,21).*(exp((algebraic(:,40)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,40))./( constants(:,100).*constants(:,101))));
algebraic(:,51) = constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
algebraic(:,77) = constants(:,22).*(exp((algebraic(:,51)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,51))./( constants(:,100).*constants(:,101))));
algebraic(:,34) = constants(:,94)+states(:,38);
algebraic(:,67) = constants(:,100).*constants(:,101).*log( constants(:,60).*algebraic(:,34));
algebraic(:,82) = constants(:,23).*(exp(algebraic(:,57)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,67)+constants(:,107))./( constants(:,100).*constants(:,101))));
rates(:,36) = (algebraic(:,72)+algebraic(:,77)) - algebraic(:,82);
algebraic(:,85) = constants(:,16).*exp((algebraic(:,44)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,55)./( constants(:,100).*constants(:,101)));
rates(:,19) = - algebraic(:,75) - algebraic(:,85);
rates(:,20) = - algebraic(:,70) - algebraic(:,85);
rates(:,22) = - algebraic(:,80)+algebraic(:,85);
algebraic(:,25) = constants(:,88)+states(:,26)+states(:,39);
algebraic(:,56) = constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
algebraic(:,28) = constants(:,89)+states(:,27)+states(:,40);
algebraic(:,61) = constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
algebraic(:,36) = constants(:,91)+states(:,29)+states(:,41);
algebraic(:,71) = constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
algebraic(:,81) = constants(:,17).*exp((algebraic(:,56)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,71)./( constants(:,100).*constants(:,101)));
% oldStim = piecewise({voi>p.stimStart(2)-1 & voi<=p.stimStart(2), p.stimMag(2).*voi , voi>p.stimStart(2)&voi<p.stimEnd(2), p.stimMag(2)}, p.stimHolding(1)(2));
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(:,23) = algebraic(:,2)+states(:,25);
algebraic(:,50) = constants(:,100).*constants(:,101).*log( constants(:,53).*algebraic(:,23));
algebraic(:,42) = constants(:,92)+states(:,30)+states(:,43);
algebraic(:,76) = constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
algebraic(:,86) = constants(:,18).*exp((algebraic(:,71)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
rates(:,29) = algebraic(:,81) - algebraic(:,86);
algebraic(:,87) = constants(:,24).*(exp((algebraic(:,67)+algebraic(:,62))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,45)./( constants(:,100).*constants(:,101))));
rates(:,32) = ( - algebraic(:,72) - algebraic(:,77))+algebraic(:,87);
rates(:,37) = (algebraic(:,72)+algebraic(:,77)) - algebraic(:,87);
rates(:,38) = algebraic(:,82) - algebraic(:,87);
algebraic(:,32) = constants(:,90)+states(:,28)+states(:,42);
algebraic(:,66) = constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
algebraic(:,90) = constants(:,19).*exp((algebraic(:,66)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
rates(:,27) = - algebraic(:,81) - algebraic(:,90);
rates(:,30) = algebraic(:,86)+algebraic(:,90);
algebraic(:,52) = constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
algebraic(:,58) = constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
algebraic(:,19) = constants(:,80)+states(:,16)+states(:,44);
algebraic(:,68) = constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
algebraic(:,37) = constants(:,96)+states(:,45);
algebraic(:,73) = constants(:,100).*constants(:,101).*log( constants(:,62).*algebraic(:,37));
algebraic(:,83) = constants(:,25).*(exp((algebraic(:,52)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,52))./( constants(:,100).*constants(:,101))));
algebraic(:,63) = constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
algebraic(:,88) = constants(:,26).*(exp((algebraic(:,63)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,63))./( constants(:,100).*constants(:,101))));
algebraic(:,46) = constants(:,97)+states(:,46);
algebraic(:,78) = constants(:,100).*constants(:,101).*log( constants(:,63).*algebraic(:,46));
algebraic(:,91) = constants(:,27).*(exp(algebraic(:,68)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,78)+constants(:,108))./( constants(:,100).*constants(:,101))));
rates(:,44) = (algebraic(:,83)+algebraic(:,88)) - algebraic(:,91);
algebraic(:,93) = constants(:,20).*exp((algebraic(:,56)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,66)./( constants(:,100).*constants(:,101)));
rates(:,25) = - algebraic(:,86) - algebraic(:,93);
rates(:,26) = - algebraic(:,81) - algebraic(:,93);
rates(:,28) = - algebraic(:,90)+algebraic(:,93);
algebraic(:,94) = constants(:,28).*(exp((algebraic(:,78)+algebraic(:,73))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,58)./( constants(:,100).*constants(:,101))));
rates(:,40) = ( - algebraic(:,83) - algebraic(:,88))+algebraic(:,94);
rates(:,45) = (algebraic(:,83)+algebraic(:,88)) - algebraic(:,94);
rates(:,46) = algebraic(:,91) - algebraic(:,94);
algebraic(:,3) = constants(:,65)+states(:,1);
algebraic(:,30) = constants(:,100).*constants(:,101).*log( constants(:,29).*algebraic(:,3));
algebraic(:,5) = constants(:,66)+states(:,3);
algebraic(:,38) = constants(:,100).*constants(:,101).*log( constants(:,31).*algebraic(:,5));
algebraic(:,6) = constants(:,68)+states(:,4);
algebraic(:,43) = constants(:,100).*constants(:,101).*log( constants(:,32).*algebraic(:,6));
algebraic(:,99) = constants(:,1).*(exp((algebraic(:,38)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,43)./( constants(:,100).*constants(:,101))));
algebraic(:,4) = constants(:,67)+states(:,2);
algebraic(:,33) = constants(:,100).*constants(:,101).*log( constants(:,30).*algebraic(:,4));
algebraic(:,24) = constants(:,82)+states(:,18);
algebraic(:,98) = constants(:,100).*constants(:,101).*log( constants(:,46).*algebraic(:,24));
algebraic(:,100) = constants(:,2).*(exp(algebraic(:,43)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,38)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
rates(:,4) = algebraic(:,99) - algebraic(:,100);
algebraic(:,7) = constants(:,73)+states(:,5);
algebraic(:,48) = constants(:,100).*constants(:,101).*log( constants(:,33).*algebraic(:,7));
algebraic(:,8) = constants(:,74)+states(:,6);
algebraic(:,54) = constants(:,100).*constants(:,101).*log( constants(:,34).*algebraic(:,8));
algebraic(:,101) = constants(:,3).*(exp((algebraic(:,48)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,54)./( constants(:,100).*constants(:,101))));
algebraic(:,102) = constants(:,4).*(exp(algebraic(:,54)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,48)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
rates(:,6) = algebraic(:,101) - algebraic(:,102);
algebraic(:,9) = constants(:,70)+states(:,7);
algebraic(:,59) = constants(:,100).*constants(:,101).*log( constants(:,35).*algebraic(:,9));
algebraic(:,10) = constants(:,71)+states(:,8);
algebraic(:,64) = constants(:,100).*constants(:,101).*log( constants(:,36).*algebraic(:,10));
algebraic(:,103) = constants(:,5).*(exp((algebraic(:,59)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,64)./( constants(:,100).*constants(:,101))));
rates(:,1) = ( - algebraic(:,99) - algebraic(:,103)) - algebraic(:,101);
algebraic(:,104) = constants(:,6).*(exp(algebraic(:,64)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,59)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
rates(:,8) = algebraic(:,103) - algebraic(:,104);
algebraic(:,11) = constants(:,75)+states(:,9);
algebraic(:,69) = constants(:,100).*constants(:,101).*log( constants(:,37).*algebraic(:,11));
algebraic(:,12) = constants(:,77)+states(:,10);
algebraic(:,74) = constants(:,100).*constants(:,101).*log( constants(:,38).*algebraic(:,12));
algebraic(:,105) = constants(:,7).*(exp((algebraic(:,69)+algebraic(:,33))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,74)./( constants(:,100).*constants(:,101))));
rates(:,2) = (algebraic(:,100)+algebraic(:,104)+algebraic(:,102)) - algebraic(:,105);
algebraic(:,92) = constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
algebraic(:,106) = constants(:,10).*(exp((algebraic(:,38)+algebraic(:,92))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,48)./( constants(:,100).*constants(:,101))));
rates(:,14) = - algebraic(:,106);
rates(:,5) = (algebraic(:,106) - algebraic(:,101))+algebraic(:,102);
algebraic(:,17) = constants(:,69)+states(:,15);
algebraic(:,95) = constants(:,100).*constants(:,101).*log( constants(:,43).*algebraic(:,17));
algebraic(:,108) = constants(:,11).*(exp((algebraic(:,95)+algebraic(:,38))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,59)./( constants(:,100).*constants(:,101))));
rates(:,15) = - algebraic(:,108);
rates(:,7) = (algebraic(:,108)+algebraic(:,104)) - algebraic(:,103);
algebraic(:,13) = constants(:,79)+states(:,11);
algebraic(:,79) = constants(:,100).*constants(:,101).*log( constants(:,39).*algebraic(:,13));
algebraic(:,107) = constants(:,8).*(exp(algebraic(:,74)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,69)+algebraic(:,79))./( constants(:,100).*constants(:,101))));
rates(:,10) = algebraic(:,105) - algebraic(:,107);
rates(:,11) = algebraic(:,107);
algebraic(:,96) = constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
algebraic(:,21) = constants(:,81)+states(:,17);
algebraic(:,97) = constants(:,100).*constants(:,101).*log( constants(:,45).*algebraic(:,21));
algebraic(:,110) = constants(:,12).*(exp((algebraic(:,38)+algebraic(:,96))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,97)./( constants(:,100).*constants(:,101))));
rates(:,3) = (((algebraic(:,100) - algebraic(:,99)) - algebraic(:,106)) - algebraic(:,108)) - algebraic(:,110);
algebraic(:,15) = constants(:,76)+states(:,13);
algebraic(:,89) = constants(:,100).*constants(:,101).*log( constants(:,41).*algebraic(:,15));
algebraic(:,14) = constants(:,78)+states(:,12);
algebraic(:,84) = constants(:,100).*constants(:,101).*log( constants(:,40).*algebraic(:,14));
algebraic(:,109) = constants(:,9).*(exp((algebraic(:,69)+algebraic(:,84))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,89)./( constants(:,100).*constants(:,101))));
rates(:,9) = (algebraic(:,107) - algebraic(:,105)) - algebraic(:,109);
rates(:,12) = - algebraic(:,109);
rates(:,13) = algebraic(:,109);
rates(:,16) = - algebraic(:,110);
rates(:,17) = algebraic(:,110);
rates(:,18) = algebraic(:,110);
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(:,20) = constants(:,83)+states(:,20)+states(:,31);
algebraic(:,44) = constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
algebraic(:,22) = constants(:,84)+states(:,21)+states(:,32);
algebraic(:,49) = constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
algebraic(:,27) = constants(:,86)+states(:,23)+states(:,33);
algebraic(:,60) = constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
algebraic(:,70) = constants(:,13).*exp((algebraic(:,44)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,60)./( constants(:,100).*constants(:,101)));
% oldStim = piecewise({voi>p.stimStart(1)-1 & voi<=p.stimStart(1), p.stimMag(1).*voi , voi>p.stimStart(1)&voi<p.stimEnd(1), p.stimMag(1)}, p.stimHolding(1)(1));
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(:,18) = algebraic(:,1)+states(:,19);
algebraic(:,39) = constants(:,100).*constants(:,101).*log( constants(:,47).*algebraic(:,18));
algebraic(:,31) = constants(:,87)+states(:,24)+states(:,35);
algebraic(:,65) = constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
algebraic(:,75) = constants(:,14).*exp((algebraic(:,60)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
algebraic(:,26) = constants(:,85)+states(:,22)+states(:,34);
algebraic(:,55) = constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
algebraic(:,80) = constants(:,15).*exp((algebraic(:,55)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
algebraic(:,40) = constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
algebraic(:,45) = constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
algebraic(:,16) = constants(:,72)+states(:,14)+states(:,36);
algebraic(:,57) = constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
algebraic(:,29) = constants(:,93)+states(:,37);
algebraic(:,62) = constants(:,100).*constants(:,101).*log( constants(:,59).*algebraic(:,29));
algebraic(:,72) = constants(:,21).*(exp((algebraic(:,40)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,40))./( constants(:,100).*constants(:,101))));
algebraic(:,51) = constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
algebraic(:,77) = constants(:,22).*(exp((algebraic(:,51)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,51))./( constants(:,100).*constants(:,101))));
algebraic(:,34) = constants(:,94)+states(:,38);
algebraic(:,67) = constants(:,100).*constants(:,101).*log( constants(:,60).*algebraic(:,34));
algebraic(:,82) = constants(:,23).*(exp(algebraic(:,57)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,67)+constants(:,107))./( constants(:,100).*constants(:,101))));
algebraic(:,85) = constants(:,16).*exp((algebraic(:,44)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,55)./( constants(:,100).*constants(:,101)));
algebraic(:,25) = constants(:,88)+states(:,26)+states(:,39);
algebraic(:,56) = constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
algebraic(:,28) = constants(:,89)+states(:,27)+states(:,40);
algebraic(:,61) = constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
algebraic(:,36) = constants(:,91)+states(:,29)+states(:,41);
algebraic(:,71) = constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
algebraic(:,81) = constants(:,17).*exp((algebraic(:,56)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,71)./( constants(:,100).*constants(:,101)));
% oldStim2 = piecewise({voi>p.stimStart(2)-1 & voi<=p.stimStart(2), p.stimMag(2).*voi , voi>p.stimStart(2)&voi<p.stimEnd(2), p.stimMag(2)}, p.stimHolding(1)(2));
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(:,23) = algebraic(:,2)+states(:,25);
algebraic(:,50) = constants(:,100).*constants(:,101).*log( constants(:,53).*algebraic(:,23));
algebraic(:,42) = constants(:,92)+states(:,30)+states(:,43);
algebraic(:,76) = constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
algebraic(:,86) = constants(:,18).*exp((algebraic(:,71)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
algebraic(:,87) = constants(:,24).*(exp((algebraic(:,67)+algebraic(:,62))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,45)./( constants(:,100).*constants(:,101))));
algebraic(:,32) = constants(:,90)+states(:,28)+states(:,42);
algebraic(:,66) = constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
algebraic(:,90) = constants(:,19).*exp((algebraic(:,66)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
algebraic(:,52) = constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
algebraic(:,58) = constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
algebraic(:,19) = constants(:,80)+states(:,16)+states(:,44);
algebraic(:,68) = constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
algebraic(:,37) = constants(:,96)+states(:,45);
algebraic(:,73) = constants(:,100).*constants(:,101).*log( constants(:,62).*algebraic(:,37));
algebraic(:,83) = constants(:,25).*(exp((algebraic(:,52)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,52))./( constants(:,100).*constants(:,101))));
algebraic(:,63) = constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
algebraic(:,88) = constants(:,26).*(exp((algebraic(:,63)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,63))./( constants(:,100).*constants(:,101))));
algebraic(:,46) = constants(:,97)+states(:,46);
algebraic(:,78) = constants(:,100).*constants(:,101).*log( constants(:,63).*algebraic(:,46));
algebraic(:,91) = constants(:,27).*(exp(algebraic(:,68)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,78)+constants(:,108))./( constants(:,100).*constants(:,101))));
algebraic(:,93) = constants(:,20).*exp((algebraic(:,56)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,66)./( constants(:,100).*constants(:,101)));
algebraic(:,94) = constants(:,28).*(exp((algebraic(:,78)+algebraic(:,73))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,58)./( constants(:,100).*constants(:,101))));
algebraic(:,3) = constants(:,65)+states(:,1);
algebraic(:,30) = constants(:,100).*constants(:,101).*log( constants(:,29).*algebraic(:,3));
algebraic(:,5) = constants(:,66)+states(:,3);
algebraic(:,38) = constants(:,100).*constants(:,101).*log( constants(:,31).*algebraic(:,5));
algebraic(:,6) = constants(:,68)+states(:,4);
algebraic(:,43) = constants(:,100).*constants(:,101).*log( constants(:,32).*algebraic(:,6));
algebraic(:,99) = constants(:,1).*(exp((algebraic(:,38)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,43)./( constants(:,100).*constants(:,101))));
algebraic(:,4) = constants(:,67)+states(:,2);
algebraic(:,33) = constants(:,100).*constants(:,101).*log( constants(:,30).*algebraic(:,4));
algebraic(:,24) = constants(:,82)+states(:,18);
algebraic(:,98) = constants(:,100).*constants(:,101).*log( constants(:,46).*algebraic(:,24));
algebraic(:,100) = constants(:,2).*(exp(algebraic(:,43)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,38)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
algebraic(:,7) = constants(:,73)+states(:,5);
algebraic(:,48) = constants(:,100).*constants(:,101).*log( constants(:,33).*algebraic(:,7));
algebraic(:,8) = constants(:,74)+states(:,6);
algebraic(:,54) = constants(:,100).*constants(:,101).*log( constants(:,34).*algebraic(:,8));
algebraic(:,101) = constants(:,3).*(exp((algebraic(:,48)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,54)./( constants(:,100).*constants(:,101))));
algebraic(:,102) = constants(:,4).*(exp(algebraic(:,54)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,48)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
algebraic(:,9) = constants(:,70)+states(:,7);
algebraic(:,59) = constants(:,100).*constants(:,101).*log( constants(:,35).*algebraic(:,9));
algebraic(:,10) = constants(:,71)+states(:,8);
algebraic(:,64) = constants(:,100).*constants(:,101).*log( constants(:,36).*algebraic(:,10));
algebraic(:,103) = constants(:,5).*(exp((algebraic(:,59)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,64)./( constants(:,100).*constants(:,101))));
algebraic(:,104) = constants(:,6).*(exp(algebraic(:,64)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,59)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
algebraic(:,11) = constants(:,75)+states(:,9);
algebraic(:,69) = constants(:,100).*constants(:,101).*log( constants(:,37).*algebraic(:,11));
algebraic(:,12) = constants(:,77)+states(:,10);
algebraic(:,74) = constants(:,100).*constants(:,101).*log( constants(:,38).*algebraic(:,12));
algebraic(:,105) = constants(:,7).*(exp((algebraic(:,69)+algebraic(:,33))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,74)./( constants(:,100).*constants(:,101))));
algebraic(:,92) = constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
algebraic(:,106) = constants(:,10).*(exp((algebraic(:,38)+algebraic(:,92))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,48)./( constants(:,100).*constants(:,101))));
algebraic(:,17) = constants(:,69)+states(:,15);
algebraic(:,95) = constants(:,100).*constants(:,101).*log( constants(:,43).*algebraic(:,17));
algebraic(:,108) = constants(:,11).*(exp((algebraic(:,95)+algebraic(:,38))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,59)./( constants(:,100).*constants(:,101))));
algebraic(:,13) = constants(:,79)+states(:,11);
algebraic(:,79) = constants(:,100).*constants(:,101).*log( constants(:,39).*algebraic(:,13));
algebraic(:,107) = constants(:,8).*(exp(algebraic(:,74)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,69)+algebraic(:,79))./( constants(:,100).*constants(:,101))));
algebraic(:,96) = constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
algebraic(:,21) = constants(:,81)+states(:,17);
algebraic(:,97) = constants(:,100).*constants(:,101).*log( constants(:,45).*algebraic(:,21));
algebraic(:,110) = constants(:,12).*(exp((algebraic(:,38)+algebraic(:,96))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,97)./( constants(:,100).*constants(:,101))));
algebraic(:,15) = constants(:,76)+states(:,13);
algebraic(:,89) = constants(:,100).*constants(:,101).*log( constants(:,41).*algebraic(:,15));
algebraic(:,14) = constants(:,78)+states(:,12);
algebraic(:,84) = constants(:,100).*constants(:,101).*log( constants(:,40).*algebraic(:,14));
algebraic(:,109) = constants(:,9).*(exp((algebraic(:,69)+algebraic(:,84))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,89)./( constants(:,100).*constants(:,101))));
algebraic(:,35) = constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
algebraic(:,41) = constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
algebraic(:,47) = constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
algebraic(:,53) = constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
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 = 160;
insize = size(strin,2);
if insize > req_length
strout = strin(1:req_length);
else
strout = [strin, blanks(req_length - insize)];
end
end
