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 =121; end % There are a total of 31 entries in each of the rate and state variable arrays. % There are a total of 134 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('protocol in component environment (dimensionless)'); LEGEND_CONSTANTS(:,2) = strpad('L_totmax in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,3) = strpad('sum_b1_AR in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,4) = strpad('Gs_tot in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,5) = strpad('Kl in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,6) = strpad('Kr in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,7) = strpad('Kc in component b1_AR_Gs_parameters (uM)'); LEGEND_CONSTANTS(:,8) = strpad('k_bar_kp in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,9) = strpad('k_bar_km in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,10) = strpad('k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)'); LEGEND_CONSTANTS(:,11) = strpad('k_p_kam in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,12) = strpad('k_g_act in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,13) = strpad('k_hyd in component b1_AR_Gs_parameters (per_sec)'); LEGEND_CONSTANTS(:,14) = strpad('k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)'); LEGEND_CONSTANTS(:,15) = strpad('AC_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,16) = strpad('ATP in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,17) = strpad('PDE_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,18) = strpad('IBMX_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,19) = strpad('Fsk_tot in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,20) = strpad('k_ac_basal in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,21) = strpad('k_ac_gsa in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,22) = strpad('k_ac_fsk in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,23) = strpad('k_pde in component cAMP_parameters (per_sec)'); LEGEND_CONSTANTS(:,24) = strpad('Km_basal in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,25) = strpad('Km_gsa in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,26) = strpad('Km_fsk in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,27) = strpad('Km_pde in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,28) = strpad('K_gsa in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,29) = strpad('K_fsk in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,30) = strpad('Ki_ibmx in component cAMP_parameters (uM)'); LEGEND_CONSTANTS(:,31) = strpad('PKAI_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,32) = strpad('PKAII_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,33) = strpad('PKI_tot in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,34) = strpad('K_a in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,35) = strpad('K_b in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,36) = strpad('K_d in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,37) = strpad('Ki_pki in component PKA_parameters (uM)'); LEGEND_CONSTANTS(:,38) = strpad('PLB_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,39) = strpad('PP1_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,40) = strpad('Inhib1_tot in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,41) = strpad('k_pka_plb in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,42) = strpad('Km_pka_plb in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,43) = strpad('k_pp1_plb in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,44) = strpad('Km_pp1_plb in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,45) = strpad('k_pka_i1 in component PLB_parameters (per_sec)'); LEGEND_CONSTANTS(:,46) = strpad('Km_pka_i1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,47) = strpad('Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)'); LEGEND_CONSTANTS(:,48) = strpad('Km_pp2a_i1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,49) = strpad('Ki_inhib1 in component PLB_parameters (uM)'); LEGEND_CONSTANTS(:,50) = strpad('epsilon in component LCC_parameters (dimensionless)'); LEGEND_CONSTANTS(:,51) = strpad('LCC_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,52) = strpad('PP1_lcc_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,53) = strpad('PP2A_lcc_tot in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,54) = strpad('k_pka_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,55) = strpad('Km_pka_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,56) = strpad('k_pp1_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,57) = strpad('Km_pp1_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,58) = strpad('k_pp2a_lcc in component LCC_parameters (per_sec)'); LEGEND_CONSTANTS(:,59) = strpad('Km_pp2a_lcc in component LCC_parameters (uM)'); LEGEND_CONSTANTS(:,60) = strpad('V_myo in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,61) = strpad('V_nsr in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,62) = strpad('V_jsr in component EC_Coupling_Parameters (uL)'); LEGEND_CONSTANTS(:,63) = strpad('A_Cap in component EC_Coupling_Parameters (cm2)'); LEGEND_CONSTANTS(:,64) = strpad('Temp in component EC_Coupling_Parameters (kelvin)'); LEGEND_CONSTANTS(:,65) = strpad('Na_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,66) = strpad('K_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,67) = strpad('Ca_ext in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,68) = strpad('G_Na in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,69) = strpad('G_to in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,70) = strpad('G_ss in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,71) = strpad('G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,72) = strpad('G_Kp in component EC_Coupling_Parameters (mS_per_uF)'); LEGEND_CONSTANTS(:,73) = strpad('f in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,74) = strpad('g in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,75) = strpad('gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)'); LEGEND_CONSTANTS(:,76) = strpad('omega in component EC_Coupling_Parameters (per_sec)'); LEGEND_CONSTANTS(:,77) = strpad('p_Ca in component EC_Coupling_Parameters (cm_per_sec)'); LEGEND_CONSTANTS(:,78) = strpad('p_K in component EC_Coupling_Parameters (cm_per_sec)'); LEGEND_CONSTANTS(:,79) = strpad('N_lcc in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,80) = strpad('I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,81) = strpad('k_NaCa in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,82) = strpad('Km_Na in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,83) = strpad('Km_Ca in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,84) = strpad('k_sat in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,85) = strpad('eta in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,86) = strpad('I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,87) = strpad('Km_Nai in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,88) = strpad('Km_Ko in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,89) = strpad('I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,90) = strpad('Km_PCa in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,91) = strpad('G_CaB in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,92) = strpad('G_NaB in component EC_Coupling_Parameters (uA_per_uF)'); LEGEND_CONSTANTS(:,93) = strpad('Pns in component EC_Coupling_Parameters (dimensionless)'); LEGEND_CONSTANTS(:,94) = strpad('Km_NS in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,95) = strpad('I_up_bar in component EC_Coupling_Parameters (mM_per_sec)'); LEGEND_CONSTANTS(:,96) = strpad('Km_up0 in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,97) = strpad('NSR_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,98) = strpad('tau_on in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,99) = strpad('tau_off in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,100) = strpad('G_max_rel in component EC_Coupling_Parameters (mM_per_sec)'); LEGEND_CONSTANTS(:,101) = strpad('d_Cai_th in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,102) = strpad('Km_rel in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,103) = strpad('CSQN_th in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,104) = strpad('CSQN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,105) = strpad('Km_CSQN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,106) = strpad('tau_tr in component EC_Coupling_Parameters (second)'); LEGEND_CONSTANTS(:,107) = strpad('TRPN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,108) = strpad('CMDN_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,109) = strpad('INDO_bar in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,110) = strpad('Km_TRPN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,111) = strpad('Km_CMDN in component EC_Coupling_Parameters (mM)'); LEGEND_CONSTANTS(:,112) = strpad('Km_INDO in component EC_Coupling_Parameters (mM)'); LEGEND_ALGEBRAIC(:,62) = strpad('LR in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,63) = strpad('LRG in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,64) = strpad('RG in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,80) = 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(:,88) = strpad('PKA_DESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,23) = strpad('PKA_RESENS in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,81) = strpad('G_ACT in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,25) = strpad('HYD in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,27) = strpad('REASSOC in component b1_AR_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,65) = strpad('L in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,66) = strpad('R in component b1_AR_module (uM)'); LEGEND_ALGEBRAIC(:,67) = 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(:,68) = strpad('PKAC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,69) = strpad('cAMP in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,42) = strpad('Gsa_GTP in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,50) = strpad('Fsk in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,43) = strpad('AC in component cAMP_module (uM)'); LEGEND_CONSTANTS(:,132) = strpad('PDE in component cAMP_module (uM)'); LEGEND_CONSTANTS(:,133) = strpad('IBMX in component cAMP_module (uM)'); LEGEND_STATES(:,7) = strpad('cAMP_tot in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,44) = strpad('Gsa_GTP_AC in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,51) = strpad('Fsk_AC in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,48) = strpad('AC_ACT_GSA in component cAMP_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,46) = strpad('AC_ACT_BASAL in component cAMP_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,53) = strpad('AC_ACT_FSK in component cAMP_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,82) = strpad('PDE_ACT in component cAMP_module (uM_per_sec)'); LEGEND_CONSTANTS(:,134) = strpad('PDE_IBMX in component cAMP_module (uM)'); LEGEND_ALGEBRAIC(:,70) = strpad('PKI in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,71) = strpad('A2RC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,72) = strpad('A2R_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,73) = strpad('A2RC_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,74) = strpad('A2R_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,75) = strpad('ARC_I in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,76) = strpad('ARC_II in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,77) = strpad('PKA_temp in component PKA_module (uM)'); LEGEND_ALGEBRAIC(:,78) = strpad('PKAC_II in component PKA_module (uM)'); LEGEND_CONSTANTS(:,113) = strpad('zed in component PKA_module (dimensionless)'); LEGEND_CONSTANTS(:,114) = strpad('zed1 in component PKA_module (dimensionless)'); LEGEND_ALGEBRAIC(:,29) = strpad('PLB in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,83) = strpad('PLB_PHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,60) = strpad('PLB_DEPHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,30) = strpad('Inhib1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,55) = strpad('Inhib1p_PP1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,84) = strpad('Inhib1_PHOSPH in component PLB_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,32) = 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(:,56) = strpad('Inhib1p in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,57) = strpad('PP1 in component PLB_module (uM)'); LEGEND_ALGEBRAIC(:,1) = strpad('frac_PLB_p in component PLB_module (dimensionless)'); LEGEND_ALGEBRAIC(:,31) = strpad('frac_PLB in component PLB_module (dimensionless)'); LEGEND_CONSTANTS(:,125) = strpad('frac_PLB_o in component PLB_module (dimensionless)'); LEGEND_ALGEBRAIC(:,34) = strpad('LCCa in component LCC_module (uM)'); LEGEND_ALGEBRAIC(:,85) = strpad('LCCa_PHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,36) = strpad('LCCa_DEPHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,38) = strpad('LCCb in component LCC_module (uM)'); LEGEND_ALGEBRAIC(:,86) = strpad('LCCb_PHOSPH in component LCC_module (uM_per_sec)'); LEGEND_ALGEBRAIC(:,40) = 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(:,126) = strpad('frac_LCCa_po in component LCC_module (dimensionless)'); LEGEND_ALGEBRAIC(:,33) = strpad('frac_LCCb_p in component LCC_module (dimensionless)'); LEGEND_CONSTANTS(:,127) = strpad('frac_LCCb_po in component LCC_module (dimensionless)'); LEGEND_ALGEBRAIC(:,35) = strpad('E_Na in component Nernst_Potentials (mV)'); LEGEND_ALGEBRAIC(:,37) = strpad('E_K in component Nernst_Potentials (mV)'); LEGEND_ALGEBRAIC(:,39) = strpad('E_Ca in component Nernst_Potentials (mV)'); LEGEND_CONSTANTS(:,128) = strpad('E_Cl in component Nernst_Potentials (mV)'); LEGEND_CONSTANTS(:,115) = strpad('R in component Nernst_Potentials (joules_per_kmole_kelvin)'); LEGEND_CONSTANTS(:,116) = strpad('Frdy in component Nernst_Potentials (coulombs_per_mole)'); LEGEND_CONSTANTS(:,129) = strpad('FoRT in component Nernst_Potentials (per_mV)'); LEGEND_CONSTANTS(:,117) = strpad('z_Na in component Nernst_Potentials (dimensionless)'); LEGEND_CONSTANTS(:,118) = strpad('z_K in component Nernst_Potentials (dimensionless)'); LEGEND_CONSTANTS(:,119) = 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(:,13) = 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(:,14) = strpad('bh in component Fast_Na_Current (per_sec)'); LEGEND_ALGEBRAIC(:,15) = 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(:,41) = 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(:,16) = strpad('b_lcc in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,18) = strpad('f_lcc in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,7) = strpad('y_lcc_inf in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,17) = strpad('tau_y_lcc in component L_Type_Calcium_Current (second)'); LEGEND_ALGEBRAIC(:,24) = strpad('gamma in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,26) = strpad('v_gamma in component L_Type_Calcium_Current (per_sec)'); LEGEND_ALGEBRAIC(:,28) = strpad('v_omega in component L_Type_Calcium_Current (per_sec)'); 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(:,45) = strpad('i_bar_Ca in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,47) = strpad('i_bar_K in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,49) = strpad('f_avail in component L_Type_Calcium_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,52) = strpad('I_Ca in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,54) = strpad('I_CaK in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,58) = strpad('I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)'); LEGEND_CONSTANTS(:,120) = strpad('alpha in component L_Type_Calcium_Current (ibar_to_I)'); 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(:,19) = strpad('tau_r_to in component Transient_Outward_K_Current (second)'); LEGEND_ALGEBRAIC(:,20) = strpad('tau_s_to in component Transient_Outward_K_Current (second)'); LEGEND_ALGEBRAIC(:,21) = 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(:,59) = 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(:,22) = 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(:,130) = 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(:,61) = strpad('I_ss in component Steady_State_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,79) = strpad('a_Ki in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,87) = strpad('b_Ki in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,89) = strpad('Ki_ss in component Time_Independent_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,90) = strpad('I_Ki in component Time_Independent_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,91) = strpad('Kp in component Plateau_K_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,92) = strpad('I_Kp in component Plateau_K_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,93) = strpad('s4 in component Na_Ca_Exchanger_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,94) = strpad('s5 in component Na_Ca_Exchanger_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,95) = strpad('I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)'); LEGEND_CONSTANTS(:,131) = strpad('sigma in component Na_K_Pump_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,96) = strpad('f_NaK in component Na_K_Pump_Current (dimensionless)'); LEGEND_ALGEBRAIC(:,97) = strpad('I_NaK in component Na_K_Pump_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,98) = strpad('I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,99) = strpad('I_CaB in component Ca_Background_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,100) = strpad('I_NaB in component Na_Background_Current (uA_per_uF)'); LEGEND_ALGEBRAIC(:,101) = strpad('I_Na_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,102) = strpad('I_K_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,103) = strpad('I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)'); LEGEND_ALGEBRAIC(:,104) = strpad('t_rel_sub in component Calcium_Induced_Calcium_Release (second)'); LEGEND_ALGEBRAIC(:,106) = strpad('ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)'); LEGEND_ALGEBRAIC(:,108) = strpad('ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)'); LEGEND_ALGEBRAIC(:,110) = strpad('g_rel in component Calcium_Induced_Calcium_Release (per_sec)'); LEGEND_ALGEBRAIC(:,111) = 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(:,112) = strpad('Km_up in component Other_SR_Fluxes_and_Concentrations (mM)'); LEGEND_ALGEBRAIC(:,113) = strpad('I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,114) = strpad('I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,115) = strpad('I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)'); LEGEND_ALGEBRAIC(:,117) = 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(:,119) = strpad('SR_content in component Other_SR_Fluxes_and_Concentrations (uM)'); LEGEND_ALGEBRAIC(:,116) = strpad('b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,118) = strpad('b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,120) = strpad('b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,121) = strpad('B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)'); LEGEND_ALGEBRAIC(:,109) = strpad('I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)'); LEGEND_CONSTANTS(:,121) = strpad('V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_CONSTANTS(:,122) = strpad('V_test in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_ALGEBRAIC(:,105) = strpad('V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)'); LEGEND_CONSTANTS(:,123) = strpad('R_clamp in component Ion_Concentrations_and_Membrane_Potential (ms)'); LEGEND_CONSTANTS(:,124) = strpad('lambda in component Ion_Concentrations_and_Membrane_Potential (nFC_per_sAcm2)'); LEGEND_ALGEBRAIC(:,107) = strpad('I_pace in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)'); 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.01; CONSTANTS(:,3) = 0.0132; CONSTANTS(:,4) = 3.83; CONSTANTS(:,5) = 0.285; CONSTANTS(:,6) = 0.062; CONSTANTS(:,7) = 33; CONSTANTS(:,8) = 1.1e-3; CONSTANTS(:,9) = 2.2e-3; CONSTANTS(:,10) = 3.6e-3; CONSTANTS(:,11) = 2.2e-3; CONSTANTS(:,12) = 16; CONSTANTS(:,13) = 0.8; CONSTANTS(:,14) = 1.21e3; CONSTANTS(:,15) = 49.7e-3; CONSTANTS(:,16) = 5e3; CONSTANTS(:,17) = 38.9e-3; CONSTANTS(:,18) = 0; CONSTANTS(:,19) = 0; CONSTANTS(:,20) = 0.2; CONSTANTS(:,21) = 8.5; CONSTANTS(:,22) = 7.3; CONSTANTS(:,23) = 5; CONSTANTS(:,24) = 1.03e3; CONSTANTS(:,25) = 315; CONSTANTS(:,26) = 860; CONSTANTS(:,27) = 1.3; CONSTANTS(:,28) = 0.4; CONSTANTS(:,29) = 44; CONSTANTS(:,30) = 30; CONSTANTS(:,31) = 0.59; CONSTANTS(:,32) = 0.025; CONSTANTS(:,33) = 0.18; CONSTANTS(:,34) = 9.14; CONSTANTS(:,35) = 1.64; CONSTANTS(:,36) = 4.375; CONSTANTS(:,37) = 0.2e-3; CONSTANTS(:,38) = 106; CONSTANTS(:,39) = 0.89; CONSTANTS(:,40) = 0.3; CONSTANTS(:,41) = 54; CONSTANTS(:,42) = 21; CONSTANTS(:,43) = 8.5; CONSTANTS(:,44) = 7; CONSTANTS(:,45) = 60; CONSTANTS(:,46) = 1; CONSTANTS(:,47) = 14; CONSTANTS(:,48) = 1; CONSTANTS(:,49) = 1e-3; CONSTANTS(:,50) = 10; CONSTANTS(:,51) = 25e-3; CONSTANTS(:,52) = 25e-3; CONSTANTS(:,53) = 25e-3; CONSTANTS(:,54) = 54; CONSTANTS(:,55) = 21; CONSTANTS(:,56) = 8.52; CONSTANTS(:,57) = 3; CONSTANTS(:,58) = 10.1; CONSTANTS(:,59) = 3; CONSTANTS(:,60) = 20.8e-6; CONSTANTS(:,61) = 9.88e-7; CONSTANTS(:,62) = 9.3e-8; CONSTANTS(:,63) = 1.534e-4; CONSTANTS(:,64) = 310; CONSTANTS(:,65) = 140; CONSTANTS(:,66) = 5.4; CONSTANTS(:,67) = 1.8; CONSTANTS(:,68) = 8; CONSTANTS(:,69) = 0.35; CONSTANTS(:,70) = 0.07; CONSTANTS(:,71) = 0.24; CONSTANTS(:,72) = 0.008; CONSTANTS(:,73) = 300; CONSTANTS(:,74) = 2e3; CONSTANTS(:,75) = 5187.5; CONSTANTS(:,76) = 10; CONSTANTS(:,77) = 1.7469e-8; CONSTANTS(:,78) = 3.234e-11; CONSTANTS(:,79) = 3e5; CONSTANTS(:,80) = -0.458; CONSTANTS(:,81) = 1483; CONSTANTS(:,82) = 87.5; CONSTANTS(:,83) = 1.38; CONSTANTS(:,84) = 0.1; CONSTANTS(:,85) = 0.35; CONSTANTS(:,86) = 1.1; CONSTANTS(:,87) = 10; CONSTANTS(:,88) = 1.5; CONSTANTS(:,89) = 1.15; CONSTANTS(:,90) = 0.5e-3; CONSTANTS(:,91) = 2.8e-3; CONSTANTS(:,92) = 1.18e-3; CONSTANTS(:,93) = 0; CONSTANTS(:,94) = 1.2e-3; CONSTANTS(:,95) = 4.7; CONSTANTS(:,96) = 3e-4; CONSTANTS(:,97) = 15; CONSTANTS(:,98) = 2e-3; CONSTANTS(:,99) = 2e-3; CONSTANTS(:,100) = 60e3; CONSTANTS(:,101) = 0.18e-3; CONSTANTS(:,102) = 0.8e-3; CONSTANTS(:,103) = 8.75; CONSTANTS(:,104) = 15; CONSTANTS(:,105) = 0.8; CONSTANTS(:,106) = 5.7e-4; CONSTANTS(:,107) = 0.07; CONSTANTS(:,108) = 0.05; CONSTANTS(:,109) = 0.07; CONSTANTS(:,110) = 0.5128e-3; CONSTANTS(:,111) = 2.38e-3; CONSTANTS(:,112) = 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; CONSTANTS(:,113) = 0; CONSTANTS(:,114) = 0; STATES(:,8) = 4.105; STATES(:,9) = 0.0526; STATES(:,10) = 5.103e-3; STATES(:,11) = 5.841e-3; CONSTANTS(:,115) = 8314; CONSTANTS(:,116) = 96485; CONSTANTS(:,117) = 1; CONSTANTS(:,118) = 1; CONSTANTS(:,119) = 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; CONSTANTS(:,120) = 10e5; 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(:,121) = -40; CONSTANTS(:,122) = -10; CONSTANTS(:,123) = 0.02; CONSTANTS(:,124) = -1e3; CONSTANTS(:,125) = 0.961300; CONSTANTS(:,126) = 0.204100; CONSTANTS(:,127) = 0.233600; CONSTANTS(:,128) = - 40.0000; CONSTANTS(:,129) = (CONSTANTS(:,116)./CONSTANTS(:,115))./CONSTANTS(:,64); CONSTANTS(:,130) = 2.10000; CONSTANTS(:,131) = (exp(CONSTANTS(:,65)./67.3000) - 1.00000)./7.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(:,130); ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp( - 0.100000.*(STATES(:,18)+47.1300))); ALGEBRAIC(:,13) = 0.0800000.*exp( - STATES(:,18)./11.0000); RATES(:,15) = 1000.00.*( ALGEBRAIC(:,3).*(1.00000 - STATES(:,15)) - ALGEBRAIC(:,13).*STATES(:,15)); ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, 0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000)); ALGEBRAIC(:,14) = 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(:,14).*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(:,15) = 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(:,15).*STATES(:,17)); ALGEBRAIC(:,6) = 400.000.*exp((STATES(:,18)+2.00000)./10.0000); ALGEBRAIC(:,16) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000); RATES(:,19) = ALGEBRAIC(:,6).*(1.00000 - STATES(:,19)) - ALGEBRAIC(:,16).*STATES(:,19); RATES(:,20) = 2.00000.*ALGEBRAIC(:,6).*(1.00000 - STATES(:,20)) - (ALGEBRAIC(:,16)./2.00000).*STATES(:,20); ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,51); ALGEBRAIC(:,18) = CONSTANTS(:,73).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,126)+0.625000); RATES(:,21) = ALGEBRAIC(:,18).*(1.00000 - STATES(:,21)) - CONSTANTS(:,74).*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(:,17) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000)); RATES(:,22) = (ALGEBRAIC(:,7) - STATES(:,22))./ALGEBRAIC(:,17); ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200)); ALGEBRAIC(:,19) = 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(:,19); ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410)); ALGEBRAIC(:,20) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000; RATES(:,25) = (ALGEBRAIC(:,9) - STATES(:,25))./ALGEBRAIC(:,20); ALGEBRAIC(:,21) = 3.70000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000; RATES(:,26) = (ALGEBRAIC(:,9) - STATES(:,26))./ALGEBRAIC(:,21); ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200)); ALGEBRAIC(:,22) = 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(:,22); ALGEBRAIC(:,25) = CONSTANTS(:,13).*STATES(:,4); ALGEBRAIC(:,27) = CONSTANTS(:,14).*STATES(:,5).*STATES(:,6); RATES(:,5) = ALGEBRAIC(:,25) - ALGEBRAIC(:,27); ALGEBRAIC(:,24) = CONSTANTS(:,75).*STATES(:,14); ALGEBRAIC(:,26) = ALGEBRAIC(:,24).*(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(:,18)./CONSTANTS(:,74))); ALGEBRAIC(:,28) = CONSTANTS(:,76).*(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(:,28).*(1.00000 - STATES(:,23)) - ALGEBRAIC(:,26).*STATES(:,23); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,80) = CONSTANTS(:,8).*(ALGEBRAIC(:,62)+ALGEBRAIC(:,63)); ALGEBRAIC(:,12) = CONSTANTS(:,9).*STATES(:,2); RATES(:,2) = ALGEBRAIC(:,80) - ALGEBRAIC(:,12); ALGEBRAIC(:,81) = CONSTANTS(:,12).*(ALGEBRAIC(:,64)+ALGEBRAIC(:,63)); RATES(:,4) = ALGEBRAIC(:,81) - ALGEBRAIC(:,25); RATES(:,6) = ALGEBRAIC(:,81) - ALGEBRAIC(:,27); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,48) = ( CONSTANTS(:,21).*ALGEBRAIC(:,44).*CONSTANTS(:,16))./(CONSTANTS(:,25)+CONSTANTS(:,16)); ALGEBRAIC(:,46) = ( CONSTANTS(:,20).*ALGEBRAIC(:,43).*CONSTANTS(:,16))./(CONSTANTS(:,24)+CONSTANTS(:,16)); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,53) = ( CONSTANTS(:,22).*ALGEBRAIC(:,51).*CONSTANTS(:,16))./(CONSTANTS(:,26)+CONSTANTS(:,16)); ALGEBRAIC(:,82) = ( CONSTANTS(:,23).*CONSTANTS(:,132).*ALGEBRAIC(:,69))./(CONSTANTS(:,27)+ALGEBRAIC(:,69)); RATES(:,7) = (ALGEBRAIC(:,46)+ALGEBRAIC(:,48)+ALGEBRAIC(:,53)) - ALGEBRAIC(:,82); ALGEBRAIC(:,29) = CONSTANTS(:,38) - STATES(:,8); ALGEBRAIC(:,83) = ( CONSTANTS(:,41).*ALGEBRAIC(:,68).*ALGEBRAIC(:,29))./(CONSTANTS(:,42)+ALGEBRAIC(:,29)); [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS, STATES, ALGEBRAIC); ALGEBRAIC(:,60) = ( CONSTANTS(:,43).*ALGEBRAIC(:,57).*STATES(:,8))./(CONSTANTS(:,44)+STATES(:,8)); RATES(:,8) = ALGEBRAIC(:,83) - ALGEBRAIC(:,60); ALGEBRAIC(:,30) = CONSTANTS(:,40) - STATES(:,9); ALGEBRAIC(:,84) = ( CONSTANTS(:,45).*ALGEBRAIC(:,68).*ALGEBRAIC(:,30))./(CONSTANTS(:,46)+ALGEBRAIC(:,30)); ALGEBRAIC(:,32) = ( CONSTANTS(:,47).*STATES(:,9))./(CONSTANTS(:,48)+STATES(:,9)); RATES(:,9) = ALGEBRAIC(:,84) - ALGEBRAIC(:,32); ALGEBRAIC(:,34) = CONSTANTS(:,51) - STATES(:,10); ALGEBRAIC(:,85) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,34))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,34)); ALGEBRAIC(:,36) = ( CONSTANTS(:,50).*CONSTANTS(:,58).*CONSTANTS(:,53).*STATES(:,10))./(CONSTANTS(:,59)+ CONSTANTS(:,50).*STATES(:,10)); RATES(:,10) = ALGEBRAIC(:,85) - ALGEBRAIC(:,36); ALGEBRAIC(:,38) = CONSTANTS(:,51) - STATES(:,11); ALGEBRAIC(:,86) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,38))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,38)); ALGEBRAIC(:,40) = ( CONSTANTS(:,50).*CONSTANTS(:,56).*CONSTANTS(:,52).*STATES(:,11))./(CONSTANTS(:,57)+ CONSTANTS(:,50).*STATES(:,11)); RATES(:,11) = ALGEBRAIC(:,86) - ALGEBRAIC(:,40); ALGEBRAIC(:,88) = CONSTANTS(:,10).*ALGEBRAIC(:,68).*STATES(:,1); ALGEBRAIC(:,23) = CONSTANTS(:,11).*STATES(:,3); RATES(:,1) = ((ALGEBRAIC(:,12) - ALGEBRAIC(:,80))+ALGEBRAIC(:,23)) - ALGEBRAIC(:,88); RATES(:,3) = ALGEBRAIC(:,88) - ALGEBRAIC(:,23); ALGEBRAIC(:,35) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,117))).*log(CONSTANTS(:,65)./STATES(:,12)); ALGEBRAIC(:,41) = CONSTANTS(:,68).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,35)); ALGEBRAIC(:,93) = 1.00000.*exp( CONSTANTS(:,85).*STATES(:,18).*CONSTANTS(:,129)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,67); ALGEBRAIC(:,94) = 1.00000.*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,129)).*power(CONSTANTS(:,65), 3.00000).*STATES(:,14); ALGEBRAIC(:,95) = (CONSTANTS(:,81)./( (power(CONSTANTS(:,82), 3.00000)+power(CONSTANTS(:,65), 3.00000)).*(CONSTANTS(:,83)+CONSTANTS(:,67)).*(1.00000+ CONSTANTS(:,84).*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,129))))).*(ALGEBRAIC(:,93) - ALGEBRAIC(:,94)); ALGEBRAIC(:,96) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,129))+ 0.0365000.*CONSTANTS(:,131).*exp( - STATES(:,18).*CONSTANTS(:,129))); ALGEBRAIC(:,97) = ( (( CONSTANTS(:,86).*ALGEBRAIC(:,96))./(1.00000+power(CONSTANTS(:,87)./STATES(:,12), 1.50000))).*CONSTANTS(:,66))./(CONSTANTS(:,66)+CONSTANTS(:,88)); ALGEBRAIC(:,100) = CONSTANTS(:,92).*(STATES(:,18) - ALGEBRAIC(:,35)); ALGEBRAIC(:,101) = ALGEBRAIC(:,41)+ALGEBRAIC(:,100)+ 3.00000.*ALGEBRAIC(:,95)+ 3.00000.*ALGEBRAIC(:,97); RATES(:,12) = ( CONSTANTS(:,124).*ALGEBRAIC(:,101).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,117).*CONSTANTS(:,116)); ALGEBRAIC(:,47) = ( CONSTANTS(:,120).*CONSTANTS(:,78).*STATES(:,18).*CONSTANTS(:,116).*CONSTANTS(:,129).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,129)) - CONSTANTS(:,66)))./(exp( STATES(:,18).*CONSTANTS(:,129)) - 1.00000); ALGEBRAIC(:,33) = STATES(:,11)./CONSTANTS(:,51); ALGEBRAIC(:,49) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,33))./CONSTANTS(:,127)+0.600000); ALGEBRAIC(:,45) = ( CONSTANTS(:,120).*CONSTANTS(:,77).*4.00000.*STATES(:,18).*CONSTANTS(:,116).*CONSTANTS(:,129).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,129)) - 0.341000.*CONSTANTS(:,67)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,129)) - 1.00000); ALGEBRAIC(:,52) = ALGEBRAIC(:,45).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,54) = (ALGEBRAIC(:,47)./(1.00000+ALGEBRAIC(:,52)./CONSTANTS(:,80))).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,37) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,118))).*log(CONSTANTS(:,66)./STATES(:,13)); ALGEBRAIC(:,59) = CONSTANTS(:,69).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,61) = CONSTANTS(:,70).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,79) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 59.2150))); ALGEBRAIC(:,87) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,37)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,37))+4.75300))); ALGEBRAIC(:,89) = ALGEBRAIC(:,79)./(ALGEBRAIC(:,79)+ALGEBRAIC(:,87)); ALGEBRAIC(:,90) = CONSTANTS(:,71).*power((CONSTANTS(:,66)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,89).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,91) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000)); ALGEBRAIC(:,92) = CONSTANTS(:,72).*ALGEBRAIC(:,91).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,102) = ((ALGEBRAIC(:,59)+ALGEBRAIC(:,61)+ALGEBRAIC(:,90)+ALGEBRAIC(:,92)) - 2.00000.*ALGEBRAIC(:,97))+ALGEBRAIC(:,54); RATES(:,13) = ( CONSTANTS(:,124).*ALGEBRAIC(:,102).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,118).*CONSTANTS(:,116)); ALGEBRAIC(:,98) = ( CONSTANTS(:,89).*STATES(:,14))./(CONSTANTS(:,90)+STATES(:,14)); ALGEBRAIC(:,39) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,119))).*log(CONSTANTS(:,67)./STATES(:,14)); ALGEBRAIC(:,99) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,39)); ALGEBRAIC(:,103) = (ALGEBRAIC(:,52)+ALGEBRAIC(:,99)+ALGEBRAIC(:,98)) - 2.00000.*ALGEBRAIC(:,95); ALGEBRAIC(:,105) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,122) }, CONSTANTS(:,121)); ALGEBRAIC(:,107) = piecewise({ sin( 2.00000.* pi.*VOI)>0.999900, 10.0000 }, 0.00000); ALGEBRAIC(:,109) = piecewise({CONSTANTS(:,1)==0.00000, 0.00000 , CONSTANTS(:,1)==1.00000, ALGEBRAIC(:,107) }, (ALGEBRAIC(:,105) - STATES(:,18))./CONSTANTS(:,123)); RATES(:,18) = - 1000.00.*((ALGEBRAIC(:,103)+ALGEBRAIC(:,102)+ALGEBRAIC(:,101)) - ALGEBRAIC(:,109)); RATES(:,30) = piecewise({ - 1000.00.*((ALGEBRAIC(:,103)+ALGEBRAIC(:,102)+ALGEBRAIC(:,101)) - ALGEBRAIC(:,109))>30000.0, 1.00000 - 10000.0.*STATES(:,30) }, 1.00000); ALGEBRAIC(:,31) = ALGEBRAIC(:,29)./CONSTANTS(:,38); ALGEBRAIC(:,112) = ( CONSTANTS(:,96).*(1.00000+ 2.00000.*ALGEBRAIC(:,31)))./(1.00000+ 2.00000.*CONSTANTS(:,125)); ALGEBRAIC(:,113) = ( CONSTANTS(:,95).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,112), 2.00000)+power(STATES(:,14), 2.00000)); ALGEBRAIC(:,114) = ( CONSTANTS(:,95).*STATES(:,31))./CONSTANTS(:,97); ALGEBRAIC(:,115) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,106); RATES(:,31) = (ALGEBRAIC(:,113) - ALGEBRAIC(:,114)) - ( ALGEBRAIC(:,115).*CONSTANTS(:,62))./CONSTANTS(:,61); ALGEBRAIC(:,104) = STATES(:,30)+0.00200000; ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,104)./CONSTANTS(:,98)); ALGEBRAIC(:,108) = exp( - ALGEBRAIC(:,104)./CONSTANTS(:,99)); ALGEBRAIC(:,110) = CONSTANTS(:,100)./( 1.00000.*(1.00000+exp((ALGEBRAIC(:,103)+5.00000)./0.900000))); ALGEBRAIC(:,111) = ALGEBRAIC(:,110).*ALGEBRAIC(:,106).*ALGEBRAIC(:,108).*(STATES(:,29) - STATES(:,14)); ALGEBRAIC(:,117) = 1.00000./(1.00000+( CONSTANTS(:,104).*CONSTANTS(:,105))./power(CONSTANTS(:,105)+STATES(:,29), 2.00000)); RATES(:,29) = ALGEBRAIC(:,117).*(ALGEBRAIC(:,115) - ALGEBRAIC(:,111)); 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(:,120) = ( CONSTANTS(:,109).*CONSTANTS(:,112))./power(CONSTANTS(:,112)+STATES(:,14), 2.00000); ALGEBRAIC(:,121) = 1.00000./(1.00000+ALGEBRAIC(:,118)+ 2.00000.*ALGEBRAIC(:,116)+ALGEBRAIC(:,120)); RATES(:,14) = ALGEBRAIC(:,121).*((( CONSTANTS(:,124).*ALGEBRAIC(:,103).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,119).*CONSTANTS(:,116)) - ( (ALGEBRAIC(:,113) - ALGEBRAIC(:,114)).*CONSTANTS(:,61))./CONSTANTS(:,60))+( ALGEBRAIC(:,111).*CONSTANTS(:,62))./CONSTANTS(:,60)); 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(:,13) = 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(:,14) = 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(:,15) = 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(:,16) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000); ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,51); ALGEBRAIC(:,18) = CONSTANTS(:,73).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,126)+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(:,17) = 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(:,19) = 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(:,20) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000; ALGEBRAIC(:,21) = 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(:,22) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000))); ALGEBRAIC(:,25) = CONSTANTS(:,13).*STATES(:,4); ALGEBRAIC(:,27) = CONSTANTS(:,14).*STATES(:,5).*STATES(:,6); ALGEBRAIC(:,24) = CONSTANTS(:,75).*STATES(:,14); ALGEBRAIC(:,26) = ALGEBRAIC(:,24).*(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(:,18)./CONSTANTS(:,74))); ALGEBRAIC(:,28) = CONSTANTS(:,76).*(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(:,80) = CONSTANTS(:,8).*(ALGEBRAIC(:,62)+ALGEBRAIC(:,63)); ALGEBRAIC(:,12) = CONSTANTS(:,9).*STATES(:,2); ALGEBRAIC(:,81) = CONSTANTS(:,12).*(ALGEBRAIC(:,64)+ALGEBRAIC(:,63)); ALGEBRAIC(:,48) = ( CONSTANTS(:,21).*ALGEBRAIC(:,44).*CONSTANTS(:,16))./(CONSTANTS(:,25)+CONSTANTS(:,16)); ALGEBRAIC(:,46) = ( CONSTANTS(:,20).*ALGEBRAIC(:,43).*CONSTANTS(:,16))./(CONSTANTS(:,24)+CONSTANTS(:,16)); ALGEBRAIC(:,53) = ( CONSTANTS(:,22).*ALGEBRAIC(:,51).*CONSTANTS(:,16))./(CONSTANTS(:,26)+CONSTANTS(:,16)); ALGEBRAIC(:,82) = ( CONSTANTS(:,23).*CONSTANTS(:,132).*ALGEBRAIC(:,69))./(CONSTANTS(:,27)+ALGEBRAIC(:,69)); ALGEBRAIC(:,29) = CONSTANTS(:,38) - STATES(:,8); ALGEBRAIC(:,83) = ( CONSTANTS(:,41).*ALGEBRAIC(:,68).*ALGEBRAIC(:,29))./(CONSTANTS(:,42)+ALGEBRAIC(:,29)); ALGEBRAIC(:,60) = ( CONSTANTS(:,43).*ALGEBRAIC(:,57).*STATES(:,8))./(CONSTANTS(:,44)+STATES(:,8)); ALGEBRAIC(:,30) = CONSTANTS(:,40) - STATES(:,9); ALGEBRAIC(:,84) = ( CONSTANTS(:,45).*ALGEBRAIC(:,68).*ALGEBRAIC(:,30))./(CONSTANTS(:,46)+ALGEBRAIC(:,30)); ALGEBRAIC(:,32) = ( CONSTANTS(:,47).*STATES(:,9))./(CONSTANTS(:,48)+STATES(:,9)); ALGEBRAIC(:,34) = CONSTANTS(:,51) - STATES(:,10); ALGEBRAIC(:,85) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,34))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,34)); ALGEBRAIC(:,36) = ( CONSTANTS(:,50).*CONSTANTS(:,58).*CONSTANTS(:,53).*STATES(:,10))./(CONSTANTS(:,59)+ CONSTANTS(:,50).*STATES(:,10)); ALGEBRAIC(:,38) = CONSTANTS(:,51) - STATES(:,11); ALGEBRAIC(:,86) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,38))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,38)); ALGEBRAIC(:,40) = ( CONSTANTS(:,50).*CONSTANTS(:,56).*CONSTANTS(:,52).*STATES(:,11))./(CONSTANTS(:,57)+ CONSTANTS(:,50).*STATES(:,11)); ALGEBRAIC(:,88) = CONSTANTS(:,10).*ALGEBRAIC(:,68).*STATES(:,1); ALGEBRAIC(:,23) = CONSTANTS(:,11).*STATES(:,3); ALGEBRAIC(:,35) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,117))).*log(CONSTANTS(:,65)./STATES(:,12)); ALGEBRAIC(:,41) = CONSTANTS(:,68).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,35)); ALGEBRAIC(:,93) = 1.00000.*exp( CONSTANTS(:,85).*STATES(:,18).*CONSTANTS(:,129)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,67); ALGEBRAIC(:,94) = 1.00000.*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,129)).*power(CONSTANTS(:,65), 3.00000).*STATES(:,14); ALGEBRAIC(:,95) = (CONSTANTS(:,81)./( (power(CONSTANTS(:,82), 3.00000)+power(CONSTANTS(:,65), 3.00000)).*(CONSTANTS(:,83)+CONSTANTS(:,67)).*(1.00000+ CONSTANTS(:,84).*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,129))))).*(ALGEBRAIC(:,93) - ALGEBRAIC(:,94)); ALGEBRAIC(:,96) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,129))+ 0.0365000.*CONSTANTS(:,131).*exp( - STATES(:,18).*CONSTANTS(:,129))); ALGEBRAIC(:,97) = ( (( CONSTANTS(:,86).*ALGEBRAIC(:,96))./(1.00000+power(CONSTANTS(:,87)./STATES(:,12), 1.50000))).*CONSTANTS(:,66))./(CONSTANTS(:,66)+CONSTANTS(:,88)); ALGEBRAIC(:,100) = CONSTANTS(:,92).*(STATES(:,18) - ALGEBRAIC(:,35)); ALGEBRAIC(:,101) = ALGEBRAIC(:,41)+ALGEBRAIC(:,100)+ 3.00000.*ALGEBRAIC(:,95)+ 3.00000.*ALGEBRAIC(:,97); ALGEBRAIC(:,47) = ( CONSTANTS(:,120).*CONSTANTS(:,78).*STATES(:,18).*CONSTANTS(:,116).*CONSTANTS(:,129).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,129)) - CONSTANTS(:,66)))./(exp( STATES(:,18).*CONSTANTS(:,129)) - 1.00000); ALGEBRAIC(:,33) = STATES(:,11)./CONSTANTS(:,51); ALGEBRAIC(:,49) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,33))./CONSTANTS(:,127)+0.600000); ALGEBRAIC(:,45) = ( CONSTANTS(:,120).*CONSTANTS(:,77).*4.00000.*STATES(:,18).*CONSTANTS(:,116).*CONSTANTS(:,129).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,129)) - 0.341000.*CONSTANTS(:,67)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,129)) - 1.00000); ALGEBRAIC(:,52) = ALGEBRAIC(:,45).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,54) = (ALGEBRAIC(:,47)./(1.00000+ALGEBRAIC(:,52)./CONSTANTS(:,80))).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23); ALGEBRAIC(:,37) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,118))).*log(CONSTANTS(:,66)./STATES(:,13)); ALGEBRAIC(:,59) = CONSTANTS(:,69).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,61) = CONSTANTS(:,70).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,79) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 59.2150))); ALGEBRAIC(:,87) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,37)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,37))+4.75300))); ALGEBRAIC(:,89) = ALGEBRAIC(:,79)./(ALGEBRAIC(:,79)+ALGEBRAIC(:,87)); ALGEBRAIC(:,90) = CONSTANTS(:,71).*power((CONSTANTS(:,66)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,89).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,91) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000)); ALGEBRAIC(:,92) = CONSTANTS(:,72).*ALGEBRAIC(:,91).*(STATES(:,18) - ALGEBRAIC(:,37)); ALGEBRAIC(:,102) = ((ALGEBRAIC(:,59)+ALGEBRAIC(:,61)+ALGEBRAIC(:,90)+ALGEBRAIC(:,92)) - 2.00000.*ALGEBRAIC(:,97))+ALGEBRAIC(:,54); ALGEBRAIC(:,98) = ( CONSTANTS(:,89).*STATES(:,14))./(CONSTANTS(:,90)+STATES(:,14)); ALGEBRAIC(:,39) = (1.00000./( CONSTANTS(:,129).*CONSTANTS(:,119))).*log(CONSTANTS(:,67)./STATES(:,14)); ALGEBRAIC(:,99) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,39)); ALGEBRAIC(:,103) = (ALGEBRAIC(:,52)+ALGEBRAIC(:,99)+ALGEBRAIC(:,98)) - 2.00000.*ALGEBRAIC(:,95); ALGEBRAIC(:,105) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,122) }, CONSTANTS(:,121)); ALGEBRAIC(:,107) = piecewise({ sin( 2.00000.* pi.*VOI)>0.999900, 10.0000 }, 0.00000); ALGEBRAIC(:,109) = piecewise({CONSTANTS(:,1)==0.00000, 0.00000 , CONSTANTS(:,1)==1.00000, ALGEBRAIC(:,107) }, (ALGEBRAIC(:,105) - STATES(:,18))./CONSTANTS(:,123)); ALGEBRAIC(:,31) = ALGEBRAIC(:,29)./CONSTANTS(:,38); ALGEBRAIC(:,112) = ( CONSTANTS(:,96).*(1.00000+ 2.00000.*ALGEBRAIC(:,31)))./(1.00000+ 2.00000.*CONSTANTS(:,125)); ALGEBRAIC(:,113) = ( CONSTANTS(:,95).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,112), 2.00000)+power(STATES(:,14), 2.00000)); ALGEBRAIC(:,114) = ( CONSTANTS(:,95).*STATES(:,31))./CONSTANTS(:,97); ALGEBRAIC(:,115) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,106); ALGEBRAIC(:,104) = STATES(:,30)+0.00200000; ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,104)./CONSTANTS(:,98)); ALGEBRAIC(:,108) = exp( - ALGEBRAIC(:,104)./CONSTANTS(:,99)); ALGEBRAIC(:,110) = CONSTANTS(:,100)./( 1.00000.*(1.00000+exp((ALGEBRAIC(:,103)+5.00000)./0.900000))); ALGEBRAIC(:,111) = ALGEBRAIC(:,110).*ALGEBRAIC(:,106).*ALGEBRAIC(:,108).*(STATES(:,29) - STATES(:,14)); ALGEBRAIC(:,117) = 1.00000./(1.00000+( CONSTANTS(:,104).*CONSTANTS(:,105))./power(CONSTANTS(:,105)+STATES(:,29), 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(:,120) = ( CONSTANTS(:,109).*CONSTANTS(:,112))./power(CONSTANTS(:,112)+STATES(:,14), 2.00000); ALGEBRAIC(:,121) = 1.00000./(1.00000+ALGEBRAIC(:,118)+ 2.00000.*ALGEBRAIC(:,116)+ALGEBRAIC(:,120)); ALGEBRAIC(:,1) = STATES(:,8)./CONSTANTS(:,38); ALGEBRAIC(:,58) = ALGEBRAIC(:,52)+ALGEBRAIC(:,54); ALGEBRAIC(:,119) = 1000.00.*(( (STATES(:,29)+STATES(:,29)./ALGEBRAIC(:,117)).*CONSTANTS(:,62))./CONSTANTS(:,60)+( STATES(:,31).*CONSTANTS(:,61))./CONSTANTS(:,60)); 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(:,132) = soln(1); CONSTANTS(:,133) = soln(2); CONSTANTS(:,134) = soln(3); else SET_CONSTANTS(:,132) = logical(1); SET_CONSTANTS(:,133) = logical(1); SET_CONSTANTS(:,134) = 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(:,132) = soln(1); TEMP_CONSTANTS(:,133) = soln(2); TEMP_CONSTANTS(:,134) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) CONSTANTS(:,132) = algebraicCandidate(1); CONSTANTS(:,133) = algebraicCandidate(2); CONSTANTS(:,134) = algebraicCandidate(3); resid(1) = CONSTANTS(:,134) - ( CONSTANTS(:,132).*CONSTANTS(:,133))./CONSTANTS(:,30); resid(2) = CONSTANTS(:,132) - (CONSTANTS(:,17) - CONSTANTS(:,134)); resid(3) = CONSTANTS(:,133) - (CONSTANTS(:,18) - CONSTANTS(:,134)); 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.000191562,0.0117971,6.39347e-06,0.988,5.5258e-05,3.8182,0.0588,0.227,0.000534718,0.00290213,0.21597,5.83396e-05,0.0306208,0.116856,0.00234908,3.91219,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(:,62) = soln(1); ALGEBRAIC(:,63) = soln(2); ALGEBRAIC(:,64) = soln(3); ALGEBRAIC(:,65) = soln(4); ALGEBRAIC(:,66) = soln(5); ALGEBRAIC(:,67) = soln(6); ALGEBRAIC(:,68) = soln(7); ALGEBRAIC(:,69) = soln(8); ALGEBRAIC(:,70) = soln(9); ALGEBRAIC(:,71) = soln(10); ALGEBRAIC(:,72) = soln(11); ALGEBRAIC(:,73) = soln(12); ALGEBRAIC(:,74) = soln(13); ALGEBRAIC(:,75) = soln(14); ALGEBRAIC(:,76) = soln(15); ALGEBRAIC(:,77) = soln(16); ALGEBRAIC(:,78) = soln(17); else SET_ALGEBRAIC(:,62) = logical(1); 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); 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(:,62) = soln(1); TEMP_ALGEBRAIC(:,63) = soln(2); TEMP_ALGEBRAIC(:,64) = soln(3); TEMP_ALGEBRAIC(:,65) = soln(4); TEMP_ALGEBRAIC(:,66) = soln(5); TEMP_ALGEBRAIC(:,67) = soln(6); TEMP_ALGEBRAIC(:,68) = soln(7); TEMP_ALGEBRAIC(:,69) = soln(8); TEMP_ALGEBRAIC(:,70) = soln(9); TEMP_ALGEBRAIC(:,71) = soln(10); TEMP_ALGEBRAIC(:,72) = soln(11); TEMP_ALGEBRAIC(:,73) = soln(12); TEMP_ALGEBRAIC(:,74) = soln(13); TEMP_ALGEBRAIC(:,75) = soln(14); TEMP_ALGEBRAIC(:,76) = soln(15); TEMP_ALGEBRAIC(:,77) = soln(16); TEMP_ALGEBRAIC(:,78) = soln(17); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,62) = algebraicCandidate(1); ALGEBRAIC(:,63) = algebraicCandidate(2); ALGEBRAIC(:,64) = algebraicCandidate(3); ALGEBRAIC(:,65) = algebraicCandidate(4); ALGEBRAIC(:,66) = algebraicCandidate(5); ALGEBRAIC(:,67) = algebraicCandidate(6); ALGEBRAIC(:,68) = algebraicCandidate(7); ALGEBRAIC(:,69) = algebraicCandidate(8); ALGEBRAIC(:,70) = algebraicCandidate(9); ALGEBRAIC(:,71) = algebraicCandidate(10); ALGEBRAIC(:,72) = algebraicCandidate(11); ALGEBRAIC(:,73) = algebraicCandidate(12); ALGEBRAIC(:,74) = algebraicCandidate(13); ALGEBRAIC(:,75) = algebraicCandidate(14); ALGEBRAIC(:,76) = algebraicCandidate(15); ALGEBRAIC(:,77) = algebraicCandidate(16); ALGEBRAIC(:,78) = algebraicCandidate(17); resid(1) = ALGEBRAIC(:,62) - ( ALGEBRAIC(:,65).*ALGEBRAIC(:,66))./CONSTANTS(:,5); resid(2) = ALGEBRAIC(:,63) - ( ALGEBRAIC(:,62).*ALGEBRAIC(:,67))./CONSTANTS(:,6); resid(3) = ALGEBRAIC(:,64) - ( ALGEBRAIC(:,66).*ALGEBRAIC(:,67))./CONSTANTS(:,7); resid(4) = ALGEBRAIC(:,65) - ((CONSTANTS(:,2) - ALGEBRAIC(:,62)) - ALGEBRAIC(:,63)); resid(5) = ALGEBRAIC(:,66) - (((STATES(:,1) - ALGEBRAIC(:,62)) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64)); resid(6) = ALGEBRAIC(:,67) - ((CONSTANTS(:,4) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64)); resid(7) = ALGEBRAIC(:,70) - ( CONSTANTS(:,33).*CONSTANTS(:,37))./(CONSTANTS(:,37)+ALGEBRAIC(:,68)+ALGEBRAIC(:,78)); resid(8) = ALGEBRAIC(:,71) - (ALGEBRAIC(:,68)./CONSTANTS(:,36)).*ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)); resid(9) = ALGEBRAIC(:,72) - ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)); resid(10) = ALGEBRAIC(:,73) - (ALGEBRAIC(:,78)./CONSTANTS(:,36)).*ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)); resid(11) = ALGEBRAIC(:,74) - ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)); resid(12) = ALGEBRAIC(:,75) - (CONSTANTS(:,34)./ALGEBRAIC(:,69)).*ALGEBRAIC(:,71); resid(13) = ALGEBRAIC(:,76) - (CONSTANTS(:,34)./ALGEBRAIC(:,69)).*ALGEBRAIC(:,73); resid(14) = ALGEBRAIC(:,77) - (( CONSTANTS(:,34).*CONSTANTS(:,35))./CONSTANTS(:,36)+( CONSTANTS(:,34).*ALGEBRAIC(:,69))./CONSTANTS(:,36)+power(ALGEBRAIC(:,69), 2.00000)./CONSTANTS(:,36)); resid(15) = ALGEBRAIC(:,69) - ((STATES(:,7) - (ALGEBRAIC(:,75)+ 2.00000.*ALGEBRAIC(:,71)+ 2.00000.*ALGEBRAIC(:,72))) - (ALGEBRAIC(:,76)+ 2.00000.*ALGEBRAIC(:,73)+ 2.00000.*ALGEBRAIC(:,74))); resid(16) = CONSTANTS(:,113) - ( 2.00000.*CONSTANTS(:,31).*power(ALGEBRAIC(:,69), 2.00000) - ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)).*( ALGEBRAIC(:,77).*ALGEBRAIC(:,68)+power(ALGEBRAIC(:,69), 2.00000))); resid(17) = CONSTANTS(:,114) - ( 2.00000.*CONSTANTS(:,32).*power(ALGEBRAIC(:,69), 2.00000) - ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)).*( ALGEBRAIC(:,77).*ALGEBRAIC(:,78)+power(ALGEBRAIC(:,69), 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.00263705];, 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(:,42) = soln(1); ALGEBRAIC(:,43) = soln(2); ALGEBRAIC(:,44) = soln(3); else SET_ALGEBRAIC(:,42) = logical(1); SET_ALGEBRAIC(:,43) = logical(1); SET_ALGEBRAIC(:,44) = 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(:,42) = soln(1); TEMP_ALGEBRAIC(:,43) = soln(2); TEMP_ALGEBRAIC(:,44) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,42) = algebraicCandidate(1); ALGEBRAIC(:,43) = algebraicCandidate(2); ALGEBRAIC(:,44) = algebraicCandidate(3); resid(1) = ALGEBRAIC(:,44) - ( ALGEBRAIC(:,42).*ALGEBRAIC(:,43))./CONSTANTS(:,28); resid(2) = ALGEBRAIC(:,42) - (STATES(:,4) - ALGEBRAIC(:,44)); resid(3) = ALGEBRAIC(:,43) - (CONSTANTS(:,15) - ALGEBRAIC(:,44)); 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];, 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(:,50) = soln(1); ALGEBRAIC(:,51) = soln(2); else SET_ALGEBRAIC(:,50) = logical(1); SET_ALGEBRAIC(:,51) = 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(:,50) = soln(1); TEMP_ALGEBRAIC(:,51) = soln(2); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,50) = algebraicCandidate(1); ALGEBRAIC(:,51) = algebraicCandidate(2); resid(1) = ALGEBRAIC(:,51) - ( ALGEBRAIC(:,50).*ALGEBRAIC(:,43))./CONSTANTS(:,29); resid(2) = ALGEBRAIC(:,50) - (CONSTANTS(:,19) - ALGEBRAIC(:,51)); 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.0525373,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(:,55) = soln(1); ALGEBRAIC(:,56) = soln(2); ALGEBRAIC(:,57) = soln(3); else SET_ALGEBRAIC(:,55) = logical(1); SET_ALGEBRAIC(:,56) = logical(1); SET_ALGEBRAIC(:,57) = 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(:,55) = soln(1); TEMP_ALGEBRAIC(:,56) = soln(2); TEMP_ALGEBRAIC(:,57) = soln(3); ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC); end end end function resid = residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES) ALGEBRAIC(:,55) = algebraicCandidate(1); ALGEBRAIC(:,56) = algebraicCandidate(2); ALGEBRAIC(:,57) = algebraicCandidate(3); resid(1) = ALGEBRAIC(:,55) - ( ALGEBRAIC(:,56).*ALGEBRAIC(:,57))./CONSTANTS(:,49); resid(2) = ALGEBRAIC(:,56) - (STATES(:,9) - ALGEBRAIC(:,55)); resid(3) = ALGEBRAIC(:,57) - (CONSTANTS(:,39) - ALGEBRAIC(:,55)); 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