Location: FCU_adenylylCyclase @ 57f217602228 / OLDFCU_adenylylCyclase.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2021-11-10 14:07:37+13:00
Desc:
f
Permanent Source URI:
https://models.cellml.org/workspace/705/rawfile/57f217602228ae3574f3cb94da1455343ef79829/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