Generated Code
The following is matlab code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction() % This is the "main function". In Matlab, things work best if you rename this function to match the filename. [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel(); end function [algebraicVariableCount] = getAlgebraicVariableCount() % Used later when setting a global variable with the number of algebraic variables. % Note: This is not the "main method". algebraicVariableCount =119; end % There are a total of 31 entries in each of the rate and state variable arrays. % There are a total of 130 entries in the constant variable array. % function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel() % Create ALGEBRAIC of correct size global algebraicVariableCount; algebraicVariableCount = getAlgebraicVariableCount(); % Initialise constants and state variables [INIT_STATES, CONSTANTS] = initConsts; % Set timespan to solve over tspan = [0, 10]; % Set numerical accuracy options for ODE solver options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1); % Solve model with ODE solver [VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options); % Compute algebraic variables [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS); ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI); % Plot state variables against variable of integration [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends(); figure(); plot(VOI, STATES); xlabel(LEGEND_VOI); l = legend(LEGEND_STATES); set(l,'Interpreter','none'); end function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends() LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = ''; LEGEND_VOI = strpad('time in component environment (second)'); LEGEND_CONSTANTS(:,1) = strpad('L_totmax in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,2) = strpad('sum_b1_AR in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,3) = strpad('Gs_tot in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,4) = strpad('Kl in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,5) = strpad('Kr in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,6) = strpad('Kc in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,7) = strpad('k_bar_kp in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,8) = strpad('k_bar_km in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,9) = strpad('k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)'); LEGEND_CONSTANTS(:,10) = strpad('k_p_kam in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,11) = strpad('k_g_act in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,12) = strpad('k_hyd in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,13) = strpad('k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)'); LEGEND_CONSTANTS(:,14) = strpad('AC_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,15) = strpad('ATP in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,16) = strpad('PDE_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,17) = strpad('IBMX_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,18) = strpad('Fsk_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,19) = strpad('k_ac_basal in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,20) = strpad('k_ac_gsa in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,21) = strpad('k_ac_fsk in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,22) = strpad('k_pde in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,23) = strpad('Km_basal in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,24) = strpad('Km_gsa in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,25) = strpad('Km_fsk in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,26) = strpad('Km_pde in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,27) = strpad('K_gsa in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,28) = strpad('K_fsk in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,29) = strpad('Ki_ibmx in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,30) = strpad('PKAI_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,31) = strpad('PKAII_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,32) = strpad('PKI_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,33) = strpad('K_a in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,34) = strpad('K_b in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,35) = strpad('K_d in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,36) = strpad('Ki_pki in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,37) = strpad('PLB_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,38) = strpad('PP1_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,39) = strpad('Inhib1_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,40) = strpad('k_pka_plb in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,41) = strpad('Km_pka_plb in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,42) = strpad('k_pp1_plb in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,43) = strpad('Km_pp1_plb in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,44) = strpad('k_pka_i1 in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,45) = strpad('Km_pka_i1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,46) = strpad('Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)'); LEGEND_CONSTANTS(:,47) = strpad('Km_pp2a_i1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,48) = strpad('Ki_inhib1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,49) = strpad('epsilon in component LCC_parameters (dimensionless)'); LEGEND_CONSTANTS(:,50) = strpad('LCC_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,51) = strpad('PP1_lcc_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,52) = strpad('PP2A_lcc_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,53) = strpad('k_pka_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,54) = strpad('Km_pka_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,55) = strpad('k_pp1_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,56) = strpad('Km_pp1_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,57) = strpad('k_pp2a_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,58) = strpad('Km_pp2a_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,59) = strpad('V_myo in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,60) = strpad('V_nsr in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,61) = strpad('V_jsr in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,62) = strpad('A_Cap in component EC_Coupling_Parameters (cm2)'); LEGEND_CONSTANTS(:,63) = strpad('Temp in component EC_Coupling_Parameters (kelvin)'); LEGEND_CONSTANTS(:,64) = strpad('Na_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,65) = strpad('K_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,66) = strpad('Ca_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,67) = strpad('G_Na in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,68) = strpad('G_to in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,69) = strpad('G_ss in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,70) = strpad('G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,71) = strpad('G_Kp in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,72) = strpad('f in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,73) = strpad('g in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,74) = strpad('gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)'); LEGEND_CONSTANTS(:,75) = strpad('omega in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,76) = strpad('p_Ca in component EC_Coupling_Parameters (cm_per_sec)'); LEGEND_CONSTANTS(:,77) = strpad('p_K in component EC_Coupling_Parameters (cm_per_sec)'); LEGEND_CONSTANTS(:,78) = strpad('N_lcc in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,79) = strpad('I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,80) = strpad('k_NaCa in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,81) = strpad('Km_Na in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,82) = strpad('Km_Ca in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,83) = strpad('k_sat in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,84) = strpad('eta in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,85) = strpad('I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,86) = strpad('Km_Nai in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,87) = strpad('Km_Ko in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,88) = strpad('I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,89) = strpad('Km_PCa in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,90) = strpad('G_CaB in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,91) = strpad('G_NaB in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,92) = strpad('Pns in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,93) = strpad('Km_NS in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,94) = strpad('I_up_bar in component EC_Coupling_Parameters (mM_per_sec)'); LEGEND_CONSTANTS(:,95) = strpad('Km_up0 in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,96) = strpad('NSR_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,97) = strpad('tau_on in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,98) = strpad('tau_off in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,99) = strpad('G_max_rel in component EC_Coupling_Parameters (mM_per_sec)'); LEGEND_CONSTANTS(:,100) = strpad('d_Cai_th in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,101) = strpad('Km_rel in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,102) = strpad('CSQN_th in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,103) = strpad('CSQN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,104) = strpad('Km_CSQN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,105) = strpad('tau_tr in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,106) = strpad('TRPN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,107) = strpad('CMDN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,108) = strpad('INDO_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,109) = strpad('Km_TRPN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,110) = strpad('Km_CMDN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,111) = strpad('Km_INDO in component EC_Coupling_Parameters (mM)'); LEGEND_ALGEBRAIC(:,63) = strpad('LR in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,64) = strpad('LRG in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,65) = strpad('RG in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,81) = strpad('BARK_DESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,12) = strpad('BARK_RESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,89) = strpad('PKA_DESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,24) = strpad('PKA_RESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,82) = strpad('G_ACT in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,26) = strpad('HYD in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,28) = strpad('REASSOC in component b1_AR_module (per_sec)'); LEGEND_ALGEBRAIC(:,66) = strpad('L in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,67) = strpad('R in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,68) = strpad('Gs in component b1_AR_module (uM)'); LEGEND_STATES(:,1) = strpad('b1_AR_tot in component b1_AR_module (uM)'); LEGEND_STATES(:,2) = strpad('b1_AR_d in component b1_AR_module (uM)'); LEGEND_STATES(:,3) = strpad('b1_AR_p in component b1_AR_module (uM)'); LEGEND_STATES(:,4) = strpad('Gs_agtp_tot in component b1_AR_module (uM)'); LEGEND_STATES(:,5) = strpad('Gs_agdp in component b1_AR_module (uM)'); LEGEND_STATES(:,6) = strpad('Gs_bg in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,69) = strpad('PKAC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,70) = strpad('cAMP in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,43) = strpad('Gsa_GTP in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,51) = strpad('Fsk in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,44) = strpad('AC in component cAMP_module (uM)'); LEGEND_CONSTANTS(:,128) = strpad('PDE in component cAMP_module (uM)'); LEGEND_CONSTANTS(:,129) = strpad('IBMX in component cAMP_module (uM)'); LEGEND_STATES(:,7) = strpad('cAMP_tot in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,45) = strpad('Gsa_GTP_AC in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,52) = strpad('Fsk_AC in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,49) = strpad('AC_ACT_GSA in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,47) = strpad('AC_ACT_BASAL in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,54) = strpad('AC_ACT_FSK in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,83) = strpad('PDE_ACT in component cAMP_module (uM)'); LEGEND_CONSTANTS(:,130) = strpad('PDE_IBMX in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,71) = strpad('PKI in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,72) = strpad('A2RC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,73) = strpad('A2R_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,74) = strpad('A2RC_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,75) = strpad('A2R_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,76) = strpad('ARC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,77) = strpad('ARC_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,78) = strpad('PKA_temp in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,79) = strpad('PKAC_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,30) = strpad('PLB in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,84) = strpad('PLB_PHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,61) = strpad('PLB_DEPHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,31) = strpad('Inhib1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,56) = strpad('Inhib1p_PP1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,85) = strpad('Inhib1_PHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,33) = strpad('Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)'); LEGEND_STATES(:,8) = strpad('PLB_p in component PLB_module (uM)'); LEGEND_STATES(:,9) = strpad('Inhib1p_tot in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,57) = strpad('Inhib1p in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,58) = strpad('PP1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,1) = strpad('frac_PLB_p in component PLB_module (dimensionless)'); LEGEND_ALGEBRAIC(:,32) = strpad('frac_PLB in component PLB_module (dimensionless)'); LEGEND_CONSTANTS(:,120) = strpad('frac_PLB_o in component PLB_module (dimensionless)'); LEGEND_ALGEBRAIC(:,35) = strpad('LCCa in component LCC_module (uM)'); LEGEND_ALGEBRAIC(:,86) = strpad('LCCa_PHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,37) = strpad('LCCa_DEPHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,39) = strpad('LCCb in component LCC_module (uM)'); LEGEND_ALGEBRAIC(:,87) = strpad('LCCb_PHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,41) = strpad('LCCb_DEPHOSPH in component LCC_module (uM_per_sec)'); LEGEND_STATES(:,10) = strpad('LCCa_p in component LCC_module (uM)'); LEGEND_STATES(:,11) = strpad('LCCb_p in component LCC_module (uM)'); LEGEND_ALGEBRAIC(:,2) = strpad('frac_LCCa_p in component LCC_module (dimensionless)'); LEGEND_CONSTANTS(:,121) = strpad('frac_LCCa_po in component LCC_module (dimensionless)'); LEGEND_ALGEBRAIC(:,34) = strpad('frac_LCCb_p in component LCC_module (dimensionless)'); LEGEND_CONSTANTS(:,122) = strpad('frac_LCCb_po in component LCC_module (dimensionless)'); LEGEND_ALGEBRAIC(:,36) = strpad('E_Na in component Nernst_Potentials (mV)'); LEGEND_ALGEBRAIC(:,38) = strpad('E_K in component Nernst_Potentials (mV)'); LEGEND_ALGEBRAIC(:,40) = strpad('E_Ca in component Nernst_Potentials (mV)'); LEGEND_CONSTANTS(:,123) = strpad('E_Cl in component Nernst_Potentials (mV)'); LEGEND_CONSTANTS(:,112) = strpad('R in component Nernst_Potentials (joules_per_mole_kelvin)'); LEGEND_CONSTANTS(:,113) = strpad('Frdy in component Nernst_Potentials (coulombs_per_mole)'); LEGEND_CONSTANTS(:,124) = strpad('FoRT in component Nernst_Potentials (per_mV)'); LEGEND_CONSTANTS(:,114) = strpad('z_Na in component Nernst_Potentials (dimensionless)'); LEGEND_CONSTANTS(:,115) = strpad('z_K in component Nernst_Potentials (dimensionless)'); LEGEND_CONSTANTS(:,116) = strpad('z_Ca in component Nernst_Potentials (dimensionless)'); LEGEND_STATES(:,12) = strpad('Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_STATES(:,13) = strpad('K_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_STATES(:,14) = strpad('Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_ALGEBRAIC(:,3) = strpad('am in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,14) = strpad('bm in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,4) = strpad('ah in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,5) = strpad('aj in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,15) = strpad('bh in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,16) = strpad('bj in component Fast_Na_Current (per_sec)'); LEGEND_STATES(:,15) = strpad('m in component Fast_Na_Current (dimensionless)'); LEGEND_STATES(:,16) = strpad('h in component Fast_Na_Current (dimensionless)'); LEGEND_STATES(:,17) = strpad('j in component Fast_Na_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,42) = strpad('I_Na in component Fast_Na_Current (uA_per_uF)'); LEGEND_STATES(:,18) = strpad('V_m in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_ALGEBRAIC(:,6) = strpad('a_lcc in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,17) = strpad('b_lcc in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,19) = strpad('f_lcc in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,7) = strpad('y_lcc_inf in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,18) = strpad('tau_y_lcc in component L_Type_Calcium_Current (second)'); LEGEND_ALGEBRAIC(:,25) = strpad('gamma in component L_Type_Calcium_Current (per_mM_per_sec)'); LEGEND_ALGEBRAIC(:,27) = strpad('v_gamma in component L_Type_Calcium_Current (per_mM_per_sec)'); LEGEND_ALGEBRAIC(:,29) = strpad('v_omega in component L_Type_Calcium_Current (dimensionless)'); LEGEND_STATES(:,19) = strpad('v in component L_Type_Calcium_Current (dimensionless)'); LEGEND_STATES(:,20) = strpad('w in component L_Type_Calcium_Current (dimensionless)'); LEGEND_STATES(:,21) = strpad('x in component L_Type_Calcium_Current (dimensionless)'); LEGEND_STATES(:,22) = strpad('y in component L_Type_Calcium_Current (dimensionless)'); LEGEND_STATES(:,23) = strpad('z in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,46) = strpad('i_bar_Ca in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,48) = strpad('i_bar_K in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,50) = strpad('f_avail in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,53) = strpad('I_Ca in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,55) = strpad('I_CaK in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,59) = strpad('I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,8) = strpad('r_toss in component Transient_Outward_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,9) = strpad('s_toss in component Transient_Outward_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,20) = strpad('tau_r_to in component Transient_Outward_K_Current (second)'); LEGEND_ALGEBRAIC(:,21) = strpad('tau_s_to in component Transient_Outward_K_Current (second)'); LEGEND_ALGEBRAIC(:,22) = strpad('tau_ss_to in component Transient_Outward_K_Current (second)'); LEGEND_STATES(:,24) = strpad('r_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_STATES(:,25) = strpad('s_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_STATES(:,26) = strpad('ss_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,60) = strpad('I_to in component Transient_Outward_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,10) = strpad('r_ss_inf in component Steady_State_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,23) = strpad('tau_r_ss in component Steady_State_K_Current (second)'); LEGEND_ALGEBRAIC(:,11) = strpad('s_ss_inf in component Steady_State_K_Current (dimensionless)'); LEGEND_CONSTANTS(:,125) = strpad('tau_s_ss in component Steady_State_K_Current (second)'); LEGEND_STATES(:,27) = strpad('r_ss in component Steady_State_K_Current (dimensionless)'); LEGEND_STATES(:,28) = strpad('s_ss in component Steady_State_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,62) = strpad('I_ss in component Steady_State_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,80) = strpad('a_Ki in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,88) = strpad('b_Ki in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,90) = strpad('Ki_ss in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,91) = strpad('I_Ki in component Time_Independent_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,92) = strpad('Kp in component Plateau_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,93) = strpad('I_Kp in component Plateau_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,94) = strpad('s4 in component Na_Ca_Exchanger_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,95) = strpad('s5 in component Na_Ca_Exchanger_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,96) = strpad('I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)'); LEGEND_CONSTANTS(:,126) = strpad('sigma in component Na_K_Pump_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,97) = strpad('f_NaK in component Na_K_Pump_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,98) = strpad('I_NaK in component Na_K_Pump_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,99) = strpad('I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,100) = strpad('I_CaB in component Ca_Background_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,101) = strpad('I_NaB in component Na_Background_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,102) = strpad('I_Na_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,103) = strpad('I_K_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,104) = strpad('I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,105) = strpad('t_rel in component Calcium_Induced_Calcium_Release (second)'); LEGEND_ALGEBRAIC(:,106) = strpad('ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)'); LEGEND_ALGEBRAIC(:,107) = strpad('ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)'); LEGEND_ALGEBRAIC(:,108) = strpad('g_rel in component Calcium_Induced_Calcium_Release (per_sec)'); LEGEND_ALGEBRAIC(:,109) = strpad('I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)'); LEGEND_STATES(:,29) = strpad('Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_STATES(:,30) = strpad('trel in component Ion_Concentrations_and_Membrane_Potential (second)'); LEGEND_ALGEBRAIC(:,110) = strpad('Km_up in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_ALGEBRAIC(:,111) = strpad('I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,112) = strpad('I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,113) = strpad('I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,115) = strpad('B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)'); LEGEND_STATES(:,31) = strpad('Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_ALGEBRAIC(:,117) = strpad('SR_content in component Other_SR_Fluxes_and_Concentrations (uM)'); LEGEND_ALGEBRAIC(:,114) = strpad('b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,116) = strpad('b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,118) = strpad('b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,119) = strpad('B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_CONSTANTS(:,127) = strpad('I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)'); LEGEND_CONSTANTS(:,117) = strpad('V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_CONSTANTS(:,118) = strpad('V_test in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_ALGEBRAIC(:,13) = strpad('V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_CONSTANTS(:,119) = strpad('R_clamp in component Ion_Concentrations_and_Membrane_Potential (dimensionless)'); LEGEND_RATES(:,1) = strpad('d/dt b1_AR_tot in component b1_AR_module (uM)'); LEGEND_RATES(:,2) = strpad('d/dt b1_AR_d in component b1_AR_module (uM)'); LEGEND_RATES(:,3) = strpad('d/dt b1_AR_p in component b1_AR_module (uM)'); LEGEND_RATES(:,4) = strpad('d/dt Gs_agtp_tot in component b1_AR_module (uM)'); LEGEND_RATES(:,5) = strpad('d/dt Gs_agdp in component b1_AR_module (uM)'); LEGEND_RATES(:,6) = strpad('d/dt Gs_bg in component b1_AR_module (uM)'); LEGEND_RATES(:,7) = strpad('d/dt cAMP_tot in component cAMP_module (uM)'); LEGEND_RATES(:,8) = strpad('d/dt PLB_p in component PLB_module (uM)'); LEGEND_RATES(:,9) = strpad('d/dt Inhib1p_tot in component PLB_module (uM)'); LEGEND_RATES(:,10) = strpad('d/dt LCCa_p in component LCC_module (uM)'); LEGEND_RATES(:,11) = strpad('d/dt LCCb_p in component LCC_module (uM)'); LEGEND_RATES(:,15) = strpad('d/dt m in component Fast_Na_Current (dimensionless)'); LEGEND_RATES(:,16) = strpad('d/dt h in component Fast_Na_Current (dimensionless)'); LEGEND_RATES(:,17) = strpad('d/dt j in component Fast_Na_Current (dimensionless)'); LEGEND_RATES(:,19) = strpad('d/dt v in component L_Type_Calcium_Current (dimensionless)'); LEGEND_RATES(:,20) = strpad('d/dt w in component L_Type_Calcium_Current (dimensionless)'); LEGEND_RATES(:,21) = strpad('d/dt x in component L_Type_Calcium_Current (dimensionless)'); LEGEND_RATES(:,22) = strpad('d/dt y in component L_Type_Calcium_Current (dimensionless)'); LEGEND_RATES(:,23) = strpad('d/dt z in component L_Type_Calcium_Current (dimensionless)'); LEGEND_RATES(:,24) = strpad('d/dt r_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_RATES(:,25) = strpad('d/dt s_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_RATES(:,26) = strpad('d/dt ss_to in component Transient_Outward_K_Current (dimensionless)'); LEGEND_RATES(:,27) = strpad('d/dt r_ss in component Steady_State_K_Current (dimensionless)'); LEGEND_RATES(:,28) = strpad('d/dt s_ss in component Steady_State_K_Current (dimensionless)'); LEGEND_RATES(:,31) = strpad('d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_RATES(:,29) = strpad('d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_RATES(:,12) = strpad('d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_RATES(:,13) = strpad('d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_RATES(:,14) = strpad('d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)'); LEGEND_RATES(:,18) = strpad('d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_RATES(:,30) = strpad('d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)'); LEGEND_STATES = LEGEND_STATES'; LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC'; LEGEND_RATES = LEGEND_RATES'; LEGEND_CONSTANTS = LEGEND_CONSTANTS'; end function [STATES, CONSTANTS] = initConsts() VOI = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = []; CONSTANTS(:,1) = 1; CONSTANTS(:,2) = 0.0132; CONSTANTS(:,3) = 3.83; CONSTANTS(:,4) = 0.285; CONSTANTS(:,5) = 0.062; CONSTANTS(:,6) = 33; CONSTANTS(:,7) = 1.1e-3; CONSTANTS(:,8) = 2.2e-3; CONSTANTS(:,9) = 3.6e-3; CONSTANTS(:,10) = 2.2e-3; CONSTANTS(:,11) = 16; CONSTANTS(:,12) = 0.8; CONSTANTS(:,13) = 1.21e3; CONSTANTS(:,14) = 49.7e-3; CONSTANTS(:,15) = 5e3; CONSTANTS(:,16) = 38.9e-3; CONSTANTS(:,17) = 0; CONSTANTS(:,18) = 0; CONSTANTS(:,19) = 0.2; CONSTANTS(:,20) = 8.5; CONSTANTS(:,21) = 7.3; CONSTANTS(:,22) = 5; CONSTANTS(:,23) = 1.03e3; CONSTANTS(:,24) = 315; CONSTANTS(:,25) = 860; CONSTANTS(:,26) = 1.3; CONSTANTS(:,27) = 0.4; CONSTANTS(:,28) = 44; CONSTANTS(:,29) = 30; CONSTANTS(:,30) = 0.59; CONSTANTS(:,31) = 0.025; CONSTANTS(:,32) = 0.18; CONSTANTS(:,33) = 9.14; CONSTANTS(:,34) = 1.64; CONSTANTS(:,35) = 4.375; CONSTANTS(:,36) = 0.2e-3; CONSTANTS(:,37) = 106; CONSTANTS(:,38) = 0.89; CONSTANTS(:,39) = 0.3; CONSTANTS(:,40) = 54; CONSTANTS(:,41) = 21; CONSTANTS(:,42) = 8.5; CONSTANTS(:,43) = 7; CONSTANTS(:,44) = 60; CONSTANTS(:,45) = 1; CONSTANTS(:,46) = 14; CONSTANTS(:,47) = 1; CONSTANTS(:,48) = 1e-3; CONSTANTS(:,49) = 10; CONSTANTS(:,50) = 25e-3; CONSTANTS(:,51) = 25e-3; CONSTANTS(:,52) = 25e-3; CONSTANTS(:,53) = 54; CONSTANTS(:,54) = 21; CONSTANTS(:,55) = 8.52; CONSTANTS(:,56) = 3; CONSTANTS(:,57) = 10.1; CONSTANTS(:,58) = 3; CONSTANTS(:,59) = 20.8e-6; CONSTANTS(:,60) = 9.88e-7; CONSTANTS(:,61) = 9.3e-8; CONSTANTS(:,62) = 1.534e-4; CONSTANTS(:,63) = 310; CONSTANTS(:,64) = 140; CONSTANTS(:,65) = 5.4; CONSTANTS(:,66) = 1.8; CONSTANTS(:,67) = 8; CONSTANTS(:,68) = 0.35; CONSTANTS(:,69) = 0.07; CONSTANTS(:,70) = 0.24; CONSTANTS(:,71) = 0.008; CONSTANTS(:,72) = 300; CONSTANTS(:,73) = 2e3; CONSTANTS(:,74) = 5187.5; CONSTANTS(:,75) = 10; CONSTANTS(:,76) = 1.7469e-8; CONSTANTS(:,77) = 3.234e-11; CONSTANTS(:,78) = 3e5; CONSTANTS(:,79) = -0.458; CONSTANTS(:,80) = 1483; CONSTANTS(:,81) = 87.5; CONSTANTS(:,82) = 1.38; CONSTANTS(:,83) = 0.1; CONSTANTS(:,84) = 0.35; CONSTANTS(:,85) = 1.1; CONSTANTS(:,86) = 10; CONSTANTS(:,87) = 1.5; CONSTANTS(:,88) = 1.15; CONSTANTS(:,89) = 0.5e-3; CONSTANTS(:,90) = 2.8e-3; CONSTANTS(:,91) = 1.18e-3; CONSTANTS(:,92) = 0; CONSTANTS(:,93) = 1.2e-3; CONSTANTS(:,94) = 4.7; CONSTANTS(:,95) = 3e-4; CONSTANTS(:,96) = 15; CONSTANTS(:,97) = 2e-3; CONSTANTS(:,98) = 2e-3; CONSTANTS(:,99) = 60e3; CONSTANTS(:,100) = 0.18e-3; CONSTANTS(:,101) = 0.8e-3; CONSTANTS(:,102) = 8.75; CONSTANTS(:,103) = 15; CONSTANTS(:,104) = 0.8; CONSTANTS(:,105) = 5.7e-4; CONSTANTS(:,106) = 0.07; CONSTANTS(:,107) = 0.05; CONSTANTS(:,108) = 0.07; CONSTANTS(:,109) = 0.5128e-3; CONSTANTS(:,110) = 2.38e-3; CONSTANTS(:,111) = 8.44e-4; STATES(:,1) = 0.01205; STATES(:,2) = 0; STATES(:,3) = 1.154e-3; STATES(:,4) = 0.02505; STATES(:,5) = 6.446e-4; STATES(:,6) = 0.02569; STATES(:,7) = 0.8453; STATES(:,8) = 4.105; STATES(:,9) = 0.0526; STATES(:,10) = 5.103e-3; STATES(:,11) = 5.841e-3; CONSTANTS(:,112) = 8314; CONSTANTS(:,113) = 96485; CONSTANTS(:,114) = 1; CONSTANTS(:,115) = 1; CONSTANTS(:,116) = 2; STATES(:,12) = 16; STATES(:,13) = 145; STATES(:,14) = 1.58e-4; STATES(:,15) = 1.4e-3; STATES(:,16) = 0.99; STATES(:,17) = 0.99; STATES(:,18) = -85.66; STATES(:,19) = 0; STATES(:,20) = 0; STATES(:,21) = 0.13; STATES(:,22) = 0.96; STATES(:,23) = 0.92; STATES(:,24) = 1.4e-3; STATES(:,25) = 1; STATES(:,26) = 0.613; STATES(:,27) = 198e-3; STATES(:,28) = 0.43; STATES(:,29) = 1.92; STATES(:,30) = 0.9; STATES(:,31) = 1.92; CONSTANTS(:,117) = -40; CONSTANTS(:,118) = -10; CONSTANTS(:,119) = 0.02; CONSTANTS(:,120) = 0.961300; CONSTANTS(:,121) = 0.204100; CONSTANTS(:,122) = 0.233600; CONSTANTS(:,123) = - 40.0000; CONSTANTS(:,124) = (CONSTANTS(:,113)./CONSTANTS(:,112))./CONSTANTS(:,63); CONSTANTS(:,125) = 2.10000; CONSTANTS(:,126) = (exp(CONSTANTS(:,64)./67.3000) - 1.00000)./7.00000; CONSTANTS(:,127) = 0.00000; [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS, STATES, ALGEBRAIC); if (isempty(STATES)), warning('Initial values for states not set');, end end function [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS) 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 ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000)); RATES(:,28) = (ALGEBRAIC(:,11) - STATES(:,28))./CONSTANTS(:,125); ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp( - 0.100000.*(STATES(:,18)+47.1300))); ALGEBRAIC(:,14) = 0.0800000.*exp( - STATES(:,18)./11.0000); RATES(:,15) = 1000.00.*( ALGEBRAIC(:,3).*(1.00000 - STATES(:,15)) - ALGEBRAIC(:,14).*STATES(:,15)); ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, 0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000)); ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) }, 3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18))); RATES(:,16) = 1000.00.*( ALGEBRAIC(:,4).*(1.00000 - STATES(:,16)) - ALGEBRAIC(:,15).*STATES(:,16)); ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( ( - 127140..*exp( 0.244400.*STATES(:,18)) - 3.47400e-05.*exp( - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300)))); ALGEBRAIC(:,16) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp( - 2.57500e-07.*STATES(:,18)))./(1.00000+exp( - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp( - 0.0105200.*STATES(:,18)))./(1.00000+exp( - 0.137800.*(STATES(:,18)+40.1400)))); RATES(:,17) = 1000.00.*( ALGEBRAIC(:,5).*(1.00000 - STATES(:,17)) - ALGEBRAIC(:,16).*STATES(:,17)); ALGEBRAIC(:,6) = 400.000.*exp((STATES(:,18)+2.00000)./10.0000); ALGEBRAIC(:,17) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000); RATES(:,19) = ALGEBRAIC(:,6).*(1.00000 - STATES(:,19)) - ALGEBRAIC(:,17).*STATES(:,19); RATES(:,20) = 2.00000.*ALGEBRAIC(:,6).*(1.00000 - STATES(:,20)) - (ALGEBRAIC(:,17)./2.00000).*STATES(:,20); ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,50); ALGEBRAIC(:,19) = CONSTANTS(:,72).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,121)+0.625000); RATES(:,21) = ALGEBRAIC(:,19).*(1.00000 - STATES(:,21)) - CONSTANTS(:,73).*STATES(:,21); ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000)); ALGEBRAIC(:,18) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000)); RATES(:,22) = (ALGEBRAIC(:,7) - STATES(:,22))./ALGEBRAIC(:,18); ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200)); ALGEBRAIC(:,20) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000))); RATES(:,24) = (ALGEBRAIC(:,8) - STATES(:,24))./ALGEBRAIC(:,20); ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410)); ALGEBRAIC(:,21) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000; RATES(:,25) = (ALGEBRAIC(:,9) - STATES(:,25))./ALGEBRAIC(:,21); ALGEBRAIC(:,22) = 3.70000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000; RATES(:,26) = (ALGEBRAIC(:,9) - STATES(:,26))./ALGEBRAIC(:,22); ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200)); ALGEBRAIC(:,23) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000))); RATES(:,27) = (ALGEBRAIC(:,10) - STATES(:,27))./ALGEBRAIC(:,23); ALGEBRAIC(:,26) = CONSTANTS(:,12).*STATES(:,4); ALGEBRAIC(:,28) = CONSTANTS(:,13).*STATES(:,5).*STATES(:,6); RATES(:,5) = ALGEBRAIC(:,26) - ALGEBRAIC(:,28); ALGEBRAIC(:,25) = CONSTANTS(:,74).*STATES(:,14); ALGEBRAIC(:,27) = ALGEBRAIC(:,25).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,19)./CONSTANTS(:,73))); ALGEBRAIC(:,29) = CONSTANTS(:,75).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000)); RATES(:,23) = ALGEBRAIC(:,29).*(1.00000 - STATES(:,23)) - ALGEBRAIC(:,27).*STATES(:,23); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,81) = CONSTANTS(:,7).*(ALGEBRAIC(:,63)+ALGEBRAIC(:,64)); ALGEBRAIC(:,12) = CONSTANTS(:,8).*STATES(:,2); RATES(:,2) = ALGEBRAIC(:,81) - ALGEBRAIC(:,12); ALGEBRAIC(:,82) = CONSTANTS(:,11).*(ALGEBRAIC(:,65)+ALGEBRAIC(:,64)); RATES(:,4) = ALGEBRAIC(:,82) - ALGEBRAIC(:,26); RATES(:,6) = ALGEBRAIC(:,82) - ALGEBRAIC(:,28); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,49) = ( CONSTANTS(:,20).*ALGEBRAIC(:,45).*CONSTANTS(:,15))./(CONSTANTS(:,24)+CONSTANTS(:,15)); ALGEBRAIC(:,47) = ( CONSTANTS(:,19).*ALGEBRAIC(:,44).*CONSTANTS(:,15))./(CONSTANTS(:,23)+CONSTANTS(:,15)); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,54) = ( CONSTANTS(:,21).*ALGEBRAIC(:,52).*CONSTANTS(:,15))./(CONSTANTS(:,25)+CONSTANTS(:,15)); ALGEBRAIC(:,83) = ( CONSTANTS(:,22).*CONSTANTS(:,128).*ALGEBRAIC(:,70))./(CONSTANTS(:,26)+ALGEBRAIC(:,70)); RATES(:,7) = (ALGEBRAIC(:,47)+ALGEBRAIC(:,49)+ALGEBRAIC(:,54)) - ALGEBRAIC(:,83); ALGEBRAIC(:,30) = CONSTANTS(:,37) - STATES(:,8); ALGEBRAIC(:,84) = ( CONSTANTS(:,40).*ALGEBRAIC(:,69).*ALGEBRAIC(:,30))./(CONSTANTS(:,41)+ALGEBRAIC(:,30)); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,61) = ( CONSTANTS(:,42).*ALGEBRAIC(:,58).*STATES(:,8))./(CONSTANTS(:,43)+STATES(:,8)); RATES(:,8) = ALGEBRAIC(:,84) - ALGEBRAIC(:,61); ALGEBRAIC(:,31) = CONSTANTS(:,39) - STATES(:,9); ALGEBRAIC(:,85) = ( CONSTANTS(:,44).*ALGEBRAIC(:,69).*ALGEBRAIC(:,31))./(CONSTANTS(:,45)+ALGEBRAIC(:,31)); ALGEBRAIC(:,33) = ( CONSTANTS(:,46).*STATES(:,9))./(CONSTANTS(:,47)+STATES(:,9)); RATES(:,9) = ALGEBRAIC(:,85) - ALGEBRAIC(:,33); ALGEBRAIC(:,35) = CONSTANTS(:,50) - STATES(:,10); ALGEBRAIC(:,86) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,35))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,35)); ALGEBRAIC(:,37) = ( CONSTANTS(:,49).*CONSTANTS(:,57).*CONSTANTS(:,52).*STATES(:,10))./(CONSTANTS(:,58)+ CONSTANTS(:,49).*STATES(:,10)); RATES(:,10) = ALGEBRAIC(:,86) - ALGEBRAIC(:,37); ALGEBRAIC(:,39) = CONSTANTS(:,50) - STATES(:,11); ALGEBRAIC(:,87) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,39))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,39)); ALGEBRAIC(:,41) = ( CONSTANTS(:,49).*CONSTANTS(:,55).*CONSTANTS(:,51).*STATES(:,11))./(CONSTANTS(:,56)+ CONSTANTS(:,49).*STATES(:,11)); RATES(:,11) = ALGEBRAIC(:,87) - ALGEBRAIC(:,41); ALGEBRAIC(:,89) = CONSTANTS(:,9).*ALGEBRAIC(:,69).*STATES(:,1); ALGEBRAIC(:,24) = CONSTANTS(:,10).*STATES(:,3); RATES(:,1) = (ALGEBRAIC(:,12) - ALGEBRAIC(:,81))+(ALGEBRAIC(:,24) - ALGEBRAIC(:,89)); RATES(:,3) = ALGEBRAIC(:,89) - ALGEBRAIC(:,24); ALGEBRAIC(:,36) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,114))).*log(CONSTANTS(:,64)./STATES(:,12)); ALGEBRAIC(:,42) = CONSTANTS(:,67).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,36)); ALGEBRAIC(:,94) = exp( CONSTANTS(:,84).*STATES(:,18).*CONSTANTS(:,124)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,66); ALGEBRAIC(:,95) = exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124)).*power(CONSTANTS(:,64), 3.00000).*STATES(:,14); ALGEBRAIC(:,96) = (CONSTANTS(:,80)./( (power(CONSTANTS(:,81), 3.00000)+power(CONSTANTS(:,64), 3.00000)).*(CONSTANTS(:,82)+CONSTANTS(:,66)).*(1.00000+ CONSTANTS(:,83).*exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124))))).*(ALGEBRAIC(:,94) - ALGEBRAIC(:,95)); ALGEBRAIC(:,97) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,124))+ 0.0365000.*CONSTANTS(:,126).*exp( - STATES(:,18).*CONSTANTS(:,124))); ALGEBRAIC(:,98) = ( (( CONSTANTS(:,85).*ALGEBRAIC(:,97))./(1.00000+power(CONSTANTS(:,86)./STATES(:,12), 1.50000))).*CONSTANTS(:,65))./(CONSTANTS(:,65)+CONSTANTS(:,87)); ALGEBRAIC(:,101) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,36)); ALGEBRAIC(:,102) = ALGEBRAIC(:,42)+ALGEBRAIC(:,101)+ 3.00000.*ALGEBRAIC(:,96)+ 3.00000.*ALGEBRAIC(:,98); RATES(:,12) = ( - 1000.00.*ALGEBRAIC(:,102).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,114).*CONSTANTS(:,113)); ALGEBRAIC(:,48) = ( CONSTANTS(:,77).*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,124)) - CONSTANTS(:,65)))./(exp( STATES(:,18).*CONSTANTS(:,124)) - 1.00000); ALGEBRAIC(:,34) = STATES(:,11)./CONSTANTS(:,50); ALGEBRAIC(:,50) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,34))./CONSTANTS(:,122)+0.600000); ALGEBRAIC(:,46) = ( CONSTANTS(:,76).*4.00000.*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 0.341000.*CONSTANTS(:,66)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 1.00000); ALGEBRAIC(:,53) = ALGEBRAIC(:,46).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,55) = (ALGEBRAIC(:,48)./(1.00000+ALGEBRAIC(:,53)./CONSTANTS(:,79))).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,38) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,13)); ALGEBRAIC(:,60) = CONSTANTS(:,68).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,62) = CONSTANTS(:,69).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,80) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 59.2150))); ALGEBRAIC(:,88) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,38)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,38))+4.75300))); ALGEBRAIC(:,90) = ALGEBRAIC(:,80)./(ALGEBRAIC(:,80)+ALGEBRAIC(:,88)); ALGEBRAIC(:,91) = CONSTANTS(:,70).*power((CONSTANTS(:,65)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,90).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,92) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000)); ALGEBRAIC(:,93) = CONSTANTS(:,71).*ALGEBRAIC(:,92).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,103) = ((ALGEBRAIC(:,60)+ALGEBRAIC(:,62)+ALGEBRAIC(:,91)+ALGEBRAIC(:,93)) - 2.00000.*ALGEBRAIC(:,98))+ALGEBRAIC(:,55); RATES(:,13) = ( - 1000.00.*ALGEBRAIC(:,103).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,115).*CONSTANTS(:,113)); ALGEBRAIC(:,99) = ( CONSTANTS(:,88).*STATES(:,14))./(CONSTANTS(:,89)+STATES(:,14)); ALGEBRAIC(:,40) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,14)); ALGEBRAIC(:,100) = CONSTANTS(:,90).*(STATES(:,18) - ALGEBRAIC(:,40)); ALGEBRAIC(:,104) = (ALGEBRAIC(:,53)+ALGEBRAIC(:,100)+ALGEBRAIC(:,99)) - 2.00000.*ALGEBRAIC(:,96); RATES(:,18) = - 1000.00.*((ALGEBRAIC(:,104)+ALGEBRAIC(:,103)+ALGEBRAIC(:,102)) - CONSTANTS(:,127)); RATES(:,30) = piecewise({ - 1000.00.*((ALGEBRAIC(:,104)+ALGEBRAIC(:,103)+ALGEBRAIC(:,102)) - CONSTANTS(:,127))>30000.0, 1.00000 - 10000.0.*STATES(:,30) }, 1.00000); ALGEBRAIC(:,32) = ALGEBRAIC(:,30)./CONSTANTS(:,37); ALGEBRAIC(:,110) = ( CONSTANTS(:,95).*(1.00000+ 2.00000.*ALGEBRAIC(:,32)))./(1.00000+ 2.00000.*CONSTANTS(:,120)); ALGEBRAIC(:,111) = ( CONSTANTS(:,94).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,110), 2.00000)+power(STATES(:,14), 2.00000)); ALGEBRAIC(:,112) = ( CONSTANTS(:,94).*STATES(:,31))./CONSTANTS(:,96); ALGEBRAIC(:,113) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,105); RATES(:,31) = (ALGEBRAIC(:,111) - ALGEBRAIC(:,112)) - ( ALGEBRAIC(:,113).*CONSTANTS(:,61))./CONSTANTS(:,60); ALGEBRAIC(:,105) = STATES(:,30)+0.00200000; ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,105)./CONSTANTS(:,97)); ALGEBRAIC(:,107) = exp( - ALGEBRAIC(:,105)./CONSTANTS(:,98)); ALGEBRAIC(:,108) = CONSTANTS(:,99)./(1.00000+exp((ALGEBRAIC(:,104)+5.00000)./0.900000)); ALGEBRAIC(:,109) = ALGEBRAIC(:,108).*ALGEBRAIC(:,106).*ALGEBRAIC(:,107).*(STATES(:,29) - STATES(:,14)); ALGEBRAIC(:,115) = 1.00000./(1.00000+( CONSTANTS(:,103).*CONSTANTS(:,104))./power(CONSTANTS(:,104)+STATES(:,29), 2.00000)); RATES(:,29) = ALGEBRAIC(:,115).*(ALGEBRAIC(:,113) - ALGEBRAIC(:,109)); ALGEBRAIC(:,114) = ( CONSTANTS(:,106).*CONSTANTS(:,109))./power(CONSTANTS(:,109)+STATES(:,14), 2.00000); ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000); ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000); ALGEBRAIC(:,119) = 1.00000./(1.00000+ALGEBRAIC(:,116)+ALGEBRAIC(:,114)+ALGEBRAIC(:,114)+ALGEBRAIC(:,118)); RATES(:,14) = - ALGEBRAIC(:,119).*((( 1000.00.*ALGEBRAIC(:,104).*CONSTANTS(:,62))./( CONSTANTS(:,59).*CONSTANTS(:,116).*CONSTANTS(:,113))+( (ALGEBRAIC(:,111) - ALGEBRAIC(:,112)).*CONSTANTS(:,60))./CONSTANTS(:,59)) - ( ALGEBRAIC(:,109).*CONSTANTS(:,61))./CONSTANTS(:,59)); RATES = RATES'; end % Calculate algebraic variables function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI) statesSize = size(STATES); statesColumnCount = statesSize(2); if ( statesColumnCount == 1) STATES = STATES'; utilOnes = 1; else statesRowCount = statesSize(1); utilOnes = ones(statesRowCount, 1); end ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000)); ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp( - 0.100000.*(STATES(:,18)+47.1300))); ALGEBRAIC(:,14) = 0.0800000.*exp( - STATES(:,18)./11.0000); ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, 0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000)); ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) }, 3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18))); ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( ( - 127140..*exp( 0.244400.*STATES(:,18)) - 3.47400e-05.*exp( - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300)))); ALGEBRAIC(:,16) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp( - 2.57500e-07.*STATES(:,18)))./(1.00000+exp( - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp( - 0.0105200.*STATES(:,18)))./(1.00000+exp( - 0.137800.*(STATES(:,18)+40.1400)))); ALGEBRAIC(:,6) = 400.000.*exp((STATES(:,18)+2.00000)./10.0000); ALGEBRAIC(:,17) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000); ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,50); ALGEBRAIC(:,19) = CONSTANTS(:,72).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,121)+0.625000); ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000)); ALGEBRAIC(:,18) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000)); ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200)); ALGEBRAIC(:,20) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000))); ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410)); ALGEBRAIC(:,21) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000; ALGEBRAIC(:,22) = 3.70000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000; ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200)); ALGEBRAIC(:,23) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000))); ALGEBRAIC(:,26) = CONSTANTS(:,12).*STATES(:,4); ALGEBRAIC(:,28) = CONSTANTS(:,13).*STATES(:,5).*STATES(:,6); ALGEBRAIC(:,25) = CONSTANTS(:,74).*STATES(:,14); ALGEBRAIC(:,27) = ALGEBRAIC(:,25).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,19)./CONSTANTS(:,73))); ALGEBRAIC(:,29) = CONSTANTS(:,75).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000)); ALGEBRAIC(:,81) = CONSTANTS(:,7).*(ALGEBRAIC(:,63)+ALGEBRAIC(:,64)); ALGEBRAIC(:,12) = CONSTANTS(:,8).*STATES(:,2); ALGEBRAIC(:,82) = CONSTANTS(:,11).*(ALGEBRAIC(:,65)+ALGEBRAIC(:,64)); ALGEBRAIC(:,49) = ( CONSTANTS(:,20).*ALGEBRAIC(:,45).*CONSTANTS(:,15))./(CONSTANTS(:,24)+CONSTANTS(:,15)); ALGEBRAIC(:,47) = ( CONSTANTS(:,19).*ALGEBRAIC(:,44).*CONSTANTS(:,15))./(CONSTANTS(:,23)+CONSTANTS(:,15)); ALGEBRAIC(:,54) = ( CONSTANTS(:,21).*ALGEBRAIC(:,52).*CONSTANTS(:,15))./(CONSTANTS(:,25)+CONSTANTS(:,15)); ALGEBRAIC(:,83) = ( CONSTANTS(:,22).*CONSTANTS(:,128).*ALGEBRAIC(:,70))./(CONSTANTS(:,26)+ALGEBRAIC(:,70)); ALGEBRAIC(:,30) = CONSTANTS(:,37) - STATES(:,8); ALGEBRAIC(:,84) = ( CONSTANTS(:,40).*ALGEBRAIC(:,69).*ALGEBRAIC(:,30))./(CONSTANTS(:,41)+ALGEBRAIC(:,30)); ALGEBRAIC(:,61) = ( CONSTANTS(:,42).*ALGEBRAIC(:,58).*STATES(:,8))./(CONSTANTS(:,43)+STATES(:,8)); ALGEBRAIC(:,31) = CONSTANTS(:,39) - STATES(:,9); ALGEBRAIC(:,85) = ( CONSTANTS(:,44).*ALGEBRAIC(:,69).*ALGEBRAIC(:,31))./(CONSTANTS(:,45)+ALGEBRAIC(:,31)); ALGEBRAIC(:,33) = ( CONSTANTS(:,46).*STATES(:,9))./(CONSTANTS(:,47)+STATES(:,9)); ALGEBRAIC(:,35) = CONSTANTS(:,50) - STATES(:,10); ALGEBRAIC(:,86) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,35))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,35)); ALGEBRAIC(:,37) = ( CONSTANTS(:,49).*CONSTANTS(:,57).*CONSTANTS(:,52).*STATES(:,10))./(CONSTANTS(:,58)+ CONSTANTS(:,49).*STATES(:,10)); ALGEBRAIC(:,39) = CONSTANTS(:,50) - STATES(:,11); ALGEBRAIC(:,87) = ( CONSTANTS(:,49).*CONSTANTS(:,53).*ALGEBRAIC(:,79).*ALGEBRAIC(:,39))./(CONSTANTS(:,54)+ CONSTANTS(:,49).*ALGEBRAIC(:,39)); ALGEBRAIC(:,41) = ( CONSTANTS(:,49).*CONSTANTS(:,55).*CONSTANTS(:,51).*STATES(:,11))./(CONSTANTS(:,56)+ CONSTANTS(:,49).*STATES(:,11)); ALGEBRAIC(:,89) = CONSTANTS(:,9).*ALGEBRAIC(:,69).*STATES(:,1); ALGEBRAIC(:,24) = CONSTANTS(:,10).*STATES(:,3); ALGEBRAIC(:,36) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,114))).*log(CONSTANTS(:,64)./STATES(:,12)); ALGEBRAIC(:,42) = CONSTANTS(:,67).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,36)); ALGEBRAIC(:,94) = exp( CONSTANTS(:,84).*STATES(:,18).*CONSTANTS(:,124)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,66); ALGEBRAIC(:,95) = exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124)).*power(CONSTANTS(:,64), 3.00000).*STATES(:,14); ALGEBRAIC(:,96) = (CONSTANTS(:,80)./( (power(CONSTANTS(:,81), 3.00000)+power(CONSTANTS(:,64), 3.00000)).*(CONSTANTS(:,82)+CONSTANTS(:,66)).*(1.00000+ CONSTANTS(:,83).*exp( (CONSTANTS(:,84) - 1.00000).*STATES(:,18).*CONSTANTS(:,124))))).*(ALGEBRAIC(:,94) - ALGEBRAIC(:,95)); ALGEBRAIC(:,97) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,124))+ 0.0365000.*CONSTANTS(:,126).*exp( - STATES(:,18).*CONSTANTS(:,124))); ALGEBRAIC(:,98) = ( (( CONSTANTS(:,85).*ALGEBRAIC(:,97))./(1.00000+power(CONSTANTS(:,86)./STATES(:,12), 1.50000))).*CONSTANTS(:,65))./(CONSTANTS(:,65)+CONSTANTS(:,87)); ALGEBRAIC(:,101) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,36)); ALGEBRAIC(:,102) = ALGEBRAIC(:,42)+ALGEBRAIC(:,101)+ 3.00000.*ALGEBRAIC(:,96)+ 3.00000.*ALGEBRAIC(:,98); ALGEBRAIC(:,48) = ( CONSTANTS(:,77).*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,124)) - CONSTANTS(:,65)))./(exp( STATES(:,18).*CONSTANTS(:,124)) - 1.00000); ALGEBRAIC(:,34) = STATES(:,11)./CONSTANTS(:,50); ALGEBRAIC(:,50) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,34))./CONSTANTS(:,122)+0.600000); ALGEBRAIC(:,46) = ( CONSTANTS(:,76).*4.00000.*STATES(:,18).*CONSTANTS(:,113).*CONSTANTS(:,124).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 0.341000.*CONSTANTS(:,66)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,124)) - 1.00000); ALGEBRAIC(:,53) = ALGEBRAIC(:,46).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,55) = (ALGEBRAIC(:,48)./(1.00000+ALGEBRAIC(:,53)./CONSTANTS(:,79))).*CONSTANTS(:,78).*ALGEBRAIC(:,50).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,38) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,13)); ALGEBRAIC(:,60) = CONSTANTS(:,68).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,62) = CONSTANTS(:,69).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,80) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 59.2150))); ALGEBRAIC(:,88) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,38)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,38)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,38))+4.75300))); ALGEBRAIC(:,90) = ALGEBRAIC(:,80)./(ALGEBRAIC(:,80)+ALGEBRAIC(:,88)); ALGEBRAIC(:,91) = CONSTANTS(:,70).*power((CONSTANTS(:,65)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,90).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,92) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000)); ALGEBRAIC(:,93) = CONSTANTS(:,71).*ALGEBRAIC(:,92).*(STATES(:,18) - ALGEBRAIC(:,38)); ALGEBRAIC(:,103) = ((ALGEBRAIC(:,60)+ALGEBRAIC(:,62)+ALGEBRAIC(:,91)+ALGEBRAIC(:,93)) - 2.00000.*ALGEBRAIC(:,98))+ALGEBRAIC(:,55); ALGEBRAIC(:,99) = ( CONSTANTS(:,88).*STATES(:,14))./(CONSTANTS(:,89)+STATES(:,14)); ALGEBRAIC(:,40) = (1.00000./( CONSTANTS(:,124).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,14)); ALGEBRAIC(:,100) = CONSTANTS(:,90).*(STATES(:,18) - ALGEBRAIC(:,40)); ALGEBRAIC(:,104) = (ALGEBRAIC(:,53)+ALGEBRAIC(:,100)+ALGEBRAIC(:,99)) - 2.00000.*ALGEBRAIC(:,96); ALGEBRAIC(:,32) = ALGEBRAIC(:,30)./CONSTANTS(:,37); ALGEBRAIC(:,110) = ( CONSTANTS(:,95).*(1.00000+ 2.00000.*ALGEBRAIC(:,32)))./(1.00000+ 2.00000.*CONSTANTS(:,120)); ALGEBRAIC(:,111) = ( CONSTANTS(:,94).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,110), 2.00000)+power(STATES(:,14), 2.00000)); ALGEBRAIC(:,112) = ( CONSTANTS(:,94).*STATES(:,31))./CONSTANTS(:,96); ALGEBRAIC(:,113) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,105); ALGEBRAIC(:,105) = STATES(:,30)+0.00200000; ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,105)./CONSTANTS(:,97)); ALGEBRAIC(:,107) = exp( - ALGEBRAIC(:,105)./CONSTANTS(:,98)); ALGEBRAIC(:,108) = CONSTANTS(:,99)./(1.00000+exp((ALGEBRAIC(:,104)+5.00000)./0.900000)); ALGEBRAIC(:,109) = ALGEBRAIC(:,108).*ALGEBRAIC(:,106).*ALGEBRAIC(:,107).*(STATES(:,29) - STATES(:,14)); ALGEBRAIC(:,115) = 1.00000./(1.00000+( CONSTANTS(:,103).*CONSTANTS(:,104))./power(CONSTANTS(:,104)+STATES(:,29), 2.00000)); ALGEBRAIC(:,114) = ( CONSTANTS(:,106).*CONSTANTS(:,109))./power(CONSTANTS(:,109)+STATES(:,14), 2.00000); ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000); ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000); ALGEBRAIC(:,119) = 1.00000./(1.00000+ALGEBRAIC(:,116)+ALGEBRAIC(:,114)+ALGEBRAIC(:,114)+ALGEBRAIC(:,118)); ALGEBRAIC(:,1) = STATES(:,8)./CONSTANTS(:,37); ALGEBRAIC(:,13) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,118) }, CONSTANTS(:,117)); ALGEBRAIC(:,59) = ALGEBRAIC(:,53)+ALGEBRAIC(:,55); ALGEBRAIC(:,117) = 1000.00.*(( (STATES(:,29)+STATES(:,29)./ALGEBRAIC(:,115)).*CONSTANTS(:,61))./CONSTANTS(:,59)+( STATES(:,31).*CONSTANTS(:,60))./CONSTANTS(:,59)); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_0; if (length(initialGuess_0) ~= 3), initialGuess_0 = [0.0389,0,0.1];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_0, options); initialGuess_0 = soln; CONSTANTS(:,128) = soln(1); CONSTANTS(:,129) = soln(2); CONSTANTS(:,130) = soln(3); else SET_CONSTANTS(:,128) = logical(1); SET_CONSTANTS(:,129) = logical(1); SET_CONSTANTS(:,130) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_0, options); initialGuess_0 = soln; TEMP_CONSTANTS(:,128) = soln(1); TEMP_CONSTANTS(:,129) = soln(2); TEMP_CONSTANTS(:,130) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) CONSTANTS(:,128) = algebraicCandidate(1); CONSTANTS(:,129) = algebraicCandidate(2); CONSTANTS(:,130) = algebraicCandidate(3); resid(1) = CONSTANTS(:,130) - ( CONSTANTS(:,128).*CONSTANTS(:,129))./CONSTANTS(:,29); resid(2) = CONSTANTS(:,128) - (CONSTANTS(:,16) - CONSTANTS(:,130)); resid(3) = CONSTANTS(:,129) - (CONSTANTS(:,17) - CONSTANTS(:,130)); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_1; if (length(initialGuess_1) ~= 17), initialGuess_1 = [0.1,0.1,0.1,0.988,5.5258e-05,3.8182,0.0588,0.227,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.0083];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_1, options); initialGuess_1 = soln; ALGEBRAIC(:,63) = soln(1); ALGEBRAIC(:,64) = soln(2); ALGEBRAIC(:,65) = soln(3); ALGEBRAIC(:,66) = soln(4); ALGEBRAIC(:,67) = soln(5); ALGEBRAIC(:,68) = soln(6); ALGEBRAIC(:,69) = soln(7); ALGEBRAIC(:,70) = soln(8); ALGEBRAIC(:,71) = soln(9); ALGEBRAIC(:,72) = soln(10); ALGEBRAIC(:,73) = soln(11); ALGEBRAIC(:,74) = soln(12); ALGEBRAIC(:,75) = soln(13); ALGEBRAIC(:,76) = soln(14); ALGEBRAIC(:,77) = soln(15); ALGEBRAIC(:,78) = soln(16); ALGEBRAIC(:,79) = soln(17); else SET_ALGEBRAIC(:,63) = logical(1); SET_ALGEBRAIC(:,64) = logical(1); SET_ALGEBRAIC(:,65) = logical(1); SET_ALGEBRAIC(:,66) = logical(1); SET_ALGEBRAIC(:,67) = logical(1); SET_ALGEBRAIC(:,68) = logical(1); SET_ALGEBRAIC(:,69) = logical(1); SET_ALGEBRAIC(:,70) = logical(1); SET_ALGEBRAIC(:,71) = logical(1); SET_ALGEBRAIC(:,72) = logical(1); SET_ALGEBRAIC(:,73) = logical(1); SET_ALGEBRAIC(:,74) = logical(1); SET_ALGEBRAIC(:,75) = logical(1); SET_ALGEBRAIC(:,76) = logical(1); SET_ALGEBRAIC(:,77) = logical(1); SET_ALGEBRAIC(:,78) = logical(1); SET_ALGEBRAIC(:,79) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_1, options); initialGuess_1 = soln; TEMP_ALGEBRAIC(:,63) = soln(1); TEMP_ALGEBRAIC(:,64) = soln(2); TEMP_ALGEBRAIC(:,65) = soln(3); TEMP_ALGEBRAIC(:,66) = soln(4); TEMP_ALGEBRAIC(:,67) = soln(5); TEMP_ALGEBRAIC(:,68) = soln(6); TEMP_ALGEBRAIC(:,69) = soln(7); TEMP_ALGEBRAIC(:,70) = soln(8); TEMP_ALGEBRAIC(:,71) = soln(9); TEMP_ALGEBRAIC(:,72) = soln(10); TEMP_ALGEBRAIC(:,73) = soln(11); TEMP_ALGEBRAIC(:,74) = soln(12); TEMP_ALGEBRAIC(:,75) = soln(13); TEMP_ALGEBRAIC(:,76) = soln(14); TEMP_ALGEBRAIC(:,77) = soln(15); TEMP_ALGEBRAIC(:,78) = soln(16); TEMP_ALGEBRAIC(:,79) = soln(17); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,63) = algebraicCandidate(1); ALGEBRAIC(:,64) = algebraicCandidate(2); ALGEBRAIC(:,65) = algebraicCandidate(3); ALGEBRAIC(:,66) = algebraicCandidate(4); ALGEBRAIC(:,67) = algebraicCandidate(5); ALGEBRAIC(:,68) = algebraicCandidate(6); ALGEBRAIC(:,69) = algebraicCandidate(7); ALGEBRAIC(:,70) = algebraicCandidate(8); ALGEBRAIC(:,71) = algebraicCandidate(9); ALGEBRAIC(:,72) = algebraicCandidate(10); ALGEBRAIC(:,73) = algebraicCandidate(11); ALGEBRAIC(:,74) = algebraicCandidate(12); ALGEBRAIC(:,75) = algebraicCandidate(13); ALGEBRAIC(:,76) = algebraicCandidate(14); ALGEBRAIC(:,77) = algebraicCandidate(15); ALGEBRAIC(:,78) = algebraicCandidate(16); ALGEBRAIC(:,79) = algebraicCandidate(17); resid(1) = ALGEBRAIC(:,63) - ( ALGEBRAIC(:,66).*ALGEBRAIC(:,67))./CONSTANTS(:,4); resid(2) = ALGEBRAIC(:,64) - ( ALGEBRAIC(:,63).*ALGEBRAIC(:,68))./CONSTANTS(:,5); resid(3) = ALGEBRAIC(:,65) - ( ALGEBRAIC(:,67).*ALGEBRAIC(:,68))./CONSTANTS(:,6); resid(4) = ALGEBRAIC(:,66) - ((CONSTANTS(:,1) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64)); resid(5) = ALGEBRAIC(:,67) - (((STATES(:,1) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64)) - ALGEBRAIC(:,65)); resid(6) = ALGEBRAIC(:,68) - ((CONSTANTS(:,3) - ALGEBRAIC(:,64)) - ALGEBRAIC(:,65)); resid(7) = ALGEBRAIC(:,71) - ( CONSTANTS(:,32).*CONSTANTS(:,36))./(CONSTANTS(:,36)+ALGEBRAIC(:,69)+ALGEBRAIC(:,79)); resid(8) = ALGEBRAIC(:,72) - (ALGEBRAIC(:,69)./CONSTANTS(:,35)).*ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)); resid(9) = ALGEBRAIC(:,73) - ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)); resid(10) = ALGEBRAIC(:,74) - (ALGEBRAIC(:,79)./CONSTANTS(:,35)).*ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)); resid(11) = ALGEBRAIC(:,75) - ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)); resid(12) = ALGEBRAIC(:,76) - (CONSTANTS(:,33)./ALGEBRAIC(:,70)).*ALGEBRAIC(:,72); resid(13) = ALGEBRAIC(:,77) - (CONSTANTS(:,33)./ALGEBRAIC(:,70)).*ALGEBRAIC(:,74); resid(14) = ALGEBRAIC(:,78) - (( CONSTANTS(:,33).*CONSTANTS(:,34))./CONSTANTS(:,35)+( CONSTANTS(:,33).*ALGEBRAIC(:,70))./CONSTANTS(:,35)+power(ALGEBRAIC(:,70), 2.00000)./CONSTANTS(:,35)); resid(15) = ALGEBRAIC(:,70) - ((STATES(:,7) - (ALGEBRAIC(:,76)+ 2.00000.*ALGEBRAIC(:,72)+ 2.00000.*ALGEBRAIC(:,73))) - (ALGEBRAIC(:,77)+ 2.00000.*ALGEBRAIC(:,74)+ 2.00000.*ALGEBRAIC(:,75))); resid(16) = 0.00000 - ( 2.00000.*CONSTANTS(:,30).*power(ALGEBRAIC(:,70), 2.00000) - ALGEBRAIC(:,69).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)).*( ALGEBRAIC(:,78).*ALGEBRAIC(:,69)+power(ALGEBRAIC(:,70), 2.00000))); resid(17) = 0.00000 - ( 2.00000.*CONSTANTS(:,31).*power(ALGEBRAIC(:,70), 2.00000) - ALGEBRAIC(:,79).*(1.00000+ALGEBRAIC(:,71)./CONSTANTS(:,36)).*( ALGEBRAIC(:,78).*ALGEBRAIC(:,79)+power(ALGEBRAIC(:,70), 2.00000))); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_2; if (length(initialGuess_2) ~= 3), initialGuess_2 = [0.0224,0.0471,0.1];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_2, options); initialGuess_2 = soln; ALGEBRAIC(:,43) = soln(1); ALGEBRAIC(:,44) = soln(2); ALGEBRAIC(:,45) = soln(3); else SET_ALGEBRAIC(:,43) = logical(1); SET_ALGEBRAIC(:,44) = logical(1); SET_ALGEBRAIC(:,45) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_2, options); initialGuess_2 = soln; TEMP_ALGEBRAIC(:,43) = soln(1); TEMP_ALGEBRAIC(:,44) = soln(2); TEMP_ALGEBRAIC(:,45) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,43) = algebraicCandidate(1); ALGEBRAIC(:,44) = algebraicCandidate(2); ALGEBRAIC(:,45) = algebraicCandidate(3); resid(1) = ALGEBRAIC(:,45) - ( ALGEBRAIC(:,43).*ALGEBRAIC(:,44))./CONSTANTS(:,27); resid(2) = ALGEBRAIC(:,43) - (STATES(:,4) - ALGEBRAIC(:,45)); resid(3) = ALGEBRAIC(:,44) - (CONSTANTS(:,14) - ALGEBRAIC(:,45)); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_3; if (length(initialGuess_3) ~= 2), initialGuess_3 = [0,0.1];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_3, options); initialGuess_3 = soln; ALGEBRAIC(:,51) = soln(1); ALGEBRAIC(:,52) = soln(2); else SET_ALGEBRAIC(:,51) = logical(1); SET_ALGEBRAIC(:,52) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_3, options); initialGuess_3 = soln; TEMP_ALGEBRAIC(:,51) = soln(1); TEMP_ALGEBRAIC(:,52) = soln(2); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,51) = algebraicCandidate(1); ALGEBRAIC(:,52) = algebraicCandidate(2); resid(1) = ALGEBRAIC(:,52) - ( ALGEBRAIC(:,51).*ALGEBRAIC(:,44))./CONSTANTS(:,28); resid(2) = ALGEBRAIC(:,51) - (CONSTANTS(:,18) - ALGEBRAIC(:,52)); end % Functions required for solving differential algebraic equation function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN) ALGEBRAIC = ALGEBRAIC_IN; CONSTANTS = CONSTANTS_IN; STATES = STATES_IN; global initialGuess_4; if (length(initialGuess_4) ~= 3), initialGuess_4 = [0.1,6.2734e-05,0.8375];, end options = optimset('Display', 'off', 'TolX', 1E-6); if length(VOI) == 1 residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES); soln = fsolve(residualfn, initialGuess_4, options); initialGuess_4 = soln; ALGEBRAIC(:,56) = soln(1); ALGEBRAIC(:,57) = soln(2); ALGEBRAIC(:,58) = soln(3); else SET_ALGEBRAIC(:,56) = logical(1); SET_ALGEBRAIC(:,57) = logical(1); SET_ALGEBRAIC(:,58) = logical(1); for i=1:length(VOI) residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:)); soln = fsolve(residualfn, initialGuess_4, options); initialGuess_4 = soln; TEMP_ALGEBRAIC(:,56) = soln(1); TEMP_ALGEBRAIC(:,57) = soln(2); TEMP_ALGEBRAIC(:,58) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,56) = algebraicCandidate(1); ALGEBRAIC(:,57) = algebraicCandidate(2); ALGEBRAIC(:,58) = algebraicCandidate(3); resid(1) = ALGEBRAIC(:,56) - ( ALGEBRAIC(:,57).*ALGEBRAIC(:,58))./CONSTANTS(:,48); resid(2) = ALGEBRAIC(:,57) - (STATES(:,9) - ALGEBRAIC(:,56)); resid(3) = ALGEBRAIC(:,58) - (CONSTANTS(:,38) - ALGEBRAIC(:,56)); 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