I am modeling a biological model that includes several (27) stiff ordinary differential equations using C ++. My program works fine under the MS C ++ 2010 expression compiler, but it does not run under the g ++ compiler (NetBeans 6.8, Ubuntu 10.04 LTS). The problem is that some of the variables become NaN . The following are the values of the Vm variable after each step of the program under the g ++ compiler:
.......... Nanon nanometer nanon nanoman nanon nanom nanom nanom nanom nanom nanom nanom nanocomponent
And here is the output of the same code without any changes in the MS C ++ compiler:
-59.4 -59.3993 -59.3986 -59.3979 -59.3972 -59.3966 -59.3959 -59.3952 -59.3946 -59.3939 -59.3933 -59.3926 -59.392 -59.3914 -59.3907 -59.3901 -59.3895 -59.3889 -59.3883 -59.3877 -59.3871 -59.3865 -59.3859 -59.3853 -59.3847 -59.3841 -59.3836 -59.383 -59.3824 -59.3819 -59.3813 -59.3808 -59.3802 -59.3797 -59.3791 -59.3786 -59.3781 -59.3775 -59.377 -
I used the libraries "cmath" and "fstream".
Where is the problem? The code in both scenarios is exactly the same.
Change 1:
OK guys, here is the whole code:
#include <iostream> #include <fstream> #include <cmath> using namespace std; void FCN(void); const int TMAX = 10000; //[ms] simulation time const int TSTEP = 1; const int MXSTEP = TMAX / TSTEP; int ISTEPPRINT = MXSTEP / 5000; //time step for storing on disc int ISTEP = 0; const double R = 8341.0; const double temp = 293.0; const double F = 96487.0; const double RT_F = R*temp / F; const double z_K = 1; const double z_Na = 1; const double z_Ca = 2; const double z_Cl = -1; const double N_Av = 6.022e23; double Ca_o = 2.0; double Na_o = 140.0; double Cl_o = 129.0; double K_o = 5; double NE = 0; double NO = 0; double cGMP = 0; //[mM] [cGMP]i double cGMPprime = 0; //var double IP3 = 0; //[mM] [IP3]i double IP3d = 0; //var double IP3prime = 0; //var double DAG = 0; //[mM] double DAGprime = 0; //var double Ca_u = 0.66; double Ca_r = 0.57; double Ca_i = 68e-6; double Na_i = 8.4; double K_i = 140; double Cl_i = 59.4; double V_m = -59.4; double V_mprime; double Na_iprime; double K_iprime; double Cl_iprime; double Ca_iprime; double vol_i = 1; double I_Natotm; double I_Ktotm; double I_Cltotm; double I_Catotm; //[pA] total membrane Ca current //Reversal potentials double E_Ca; //[mV] double E_Na; //[mV] double E_K; //[mV] double E_Cl; //Membrane capacitance and area double C_m = 25.0; double A_m = C_m / 1e6; //Voltage dependent calcium current I_CaL double I_CaL; double P_CaL = 1.88e-5; double d_L; double d_Lprime; double d_Lbar; double tau_d_L; double f_L; double f_Lbar; double tau_f_L; double f_Lprime; //Delayed rectifier current I_K double I_K; double g_K = 1.35; double p_K; double p_Kbar; double V_1_2 = -11.0; double k = 15.0; double tau_p_K; double p_Kprime; double q_1 = 1; double q_2 = 1; double q_bar; double q_1prime; double q_2prime; double Pmin_NSC = 0.4344; double Po_NSC; double PNa_NSC = (5.11e-7); double PCa_NSC = (5.11e-7)*4.54; double PK_NSC = (5.11e-7)*1.06; // double d_NSCmin = 0.0244; double K_NSC = 3.0e-3; double INa_NSC; double ICa_NSC; double IK_NSC; double I_NSC; //KATP current I_KATP double I_KATP; //[pA] background K current double g_KATP = 0.067; //[nS] max. background K current conductance //Inward rectifier current I_K_i double I_K_i; //[pA] double g_maxK_i; //[nS] max. slope conductance const double G_K_i = 0; // inward rectifier constant const double n_K_i = 0.5; // inward rectifier constant //Calcium-activated potassium current I_KCa double I_KCa; double i1_KCa; double P_BKCa = 3.9e-13; double N_BKCa = 6.6e6; double P_KCa; double p_obar; double V_1_2_KCa; double p_f; double p_s; double p_fprime; double p_sprime; double tau_pf = 0.84; double tau_ps = 35.9; double dV_1_2_KCa_NO = 46.3; double R_NO; double dV_1_2_KCa_cGMP = 76; double R_cGMP; double k_leak = 1; double R_00; double R_01 = 0.9955; double R_10 = 0.0033; double R_11 = 4.0e-6; double R_01prime; double R_10prime; double R_11prime; const double Kr1 = 2500.0; const double Kr2 = 1.05; const double K_r1 = 0.0076; const double K_r2 = 0.084; double I_up; const double I_upbar = 3.34 * (k_leak + 1); const double K_mup = 0.001; double I_tr; const double vol_u = 0.07; double tau_tr = 1000.0; double I_rel; const double vol_r = 0.007; const double tau_rel = 0.0333; //[ms] const double R_leak = 1.07e-5 * (k_leak); ////equal to R_10^2 during concentration clamp // time constant of SR release double Ca_uprime; // dCa_u/dt double Ca_rprime; // dCa_r/dt double S_CM; const double K_d = 2.6e-4; const double S_CMbar = 0.1; double CaCM; const double K_dB = 5.298e-4; const double B_Fbar = 0.1; const double vol_Ca = 0.7; const double CSQNbar = 15; const double K_CSQN = 0.8; double I_PMCA; double I_PMCAbar = 5.37; double K_mPMCA = 170e-6; double I_NaK; ////[pA] Na/K pump double I_NaKbar = 2.3083; const double K_mK = 1.6; const double K_mNa = 22; const double Q_10_NaK = 1.87; double I_NCX; const double gamma2 = 0.45; // double g_NCX = 0.000487; //[nS] double d_NCX = 0.0003; // double Fi_F; // double Fi_R; // double I_NaKCl_Cl; //[pA] double L_cotr = 1.79e-8; double I_M = 0; //[pA] double I_MCa = 0; double I_MNa = 0; //[pA] Na component double I_MK = 0; double I_SOC; //[pA] double I_SOCCa; double I_SOCNa; //[pA] Na component const double g_SOCCa = 0.0083; //[nS] const double g_SOCNa = 0.0575; //[nS] const double H_SOC = 1; const double K_SOC = 0.0001; const double tau_SOC = 100; double P_SOCbar; double P_SOC = 0; double P_SOCprime; //Chloride currents double I_Cl; double g_Cl = 0.23; double alpha_Cl; double P_Cl; //Stimulation current double I_stim = 0; //[pA] //IP3 receptor double I_IP3; double I_IP3bar = 2880e-6; //[1/ms] double K_IP3 = 0.12e-3; double K_actIP3 = 0.17e-3; double K_inhIP3 = 0.1e-3; //[mM] double h_IP3; double k_onIP3 = 1.4; double h_IP3prime; double R_TG = 2e4; double K_1G = 0.01; double K_2G = 0.2; double k_rG = 1.75e-7; double k_pG = 0.1e-3; double k_eG = 6e-6; double ksi_G = 0.85; double G_TG = 1e5; double k_degG = 1.25e-3; double k_aG = 0.17e-3; double k_dG = 1.5e-3; double PIP2_T = 5e7; double r_rG = 0.015e-3; double K_cG = 0.4e-3; double alpha_G = 2.781e-8; double vol_IP3 = vol_i; double gamma_G = N_Av*vol_IP3 * 1e-15; double RS_G = R_TG*ksi_G; double RS_PG = 0; double PIP2; // double r_hG; double G; double delta_G; // double RS_Gprime; double RS_PGprime; double Gprime; double PIP2prime; double rho_rG; //cGMP formation double k1sGC = 2e3; //[1/mM/ms] double k_1sGC = 15e-3; //[1/ms] double k2sGC = 0.64e-5; //[1/ms] double k_2sGC = 0.1e-6; //[1/ms] double k3sGC = 4.2; //[1/mM/ms] double kDsGC = 0.4e-3; double kDact_deactsGC = 0.1e-3; //[1/ms] double V_cGMP = 0; // double V_cGMPprime; double V_cGMPmax = 0.1 * 1.26e-6; //[mM/ms] double V_cGMPbar; double B5sGC = k2sGC / k3sGC; double A0sGC = ((k_1sGC + k2sGC) * kDsGC + k_1sGC*k_2sGC) / (k1sGC*k3sGC); double A1sGC = ((k1sGC + k3sGC) * kDsGC + (k2sGC + k_2sGC) * k1sGC) / (k1sGC*k3sGC); double kpde_cGMP = 0.0695e-3; //[1/ms] double tausGC; const int N = 27; double Y[N], Y1[N], YPRIME[N]; int N1 = 33; double T = 0; int main(void) { ofstream fileY, fileY1, fileT; // initial conditions SMC //ICaL d_Lbar = 1.0 / (1 + exp(-(V_m) / 8.3)); d_L = d_Lbar; f_Lbar = 1.0 / (1 + exp((V_m + 42.0) / 9.1)); f_L = f_Lbar; //IKCa R_NO = NO / (NO + 200e-6); R_cGMP = pow(cGMP, 2) / (pow(cGMP, 2) + pow(0.55e-3, 2)); V_1_2_KCa = -41.7 * log10(Ca_i) - 128.2 - dV_1_2_KCa_NO * R_NO - dV_1_2_KCa_cGMP*R_cGMP; p_obar = 1 / (1 + exp(-(V_m - V_1_2_KCa) / 18.25)); p_f = p_obar; p_s = p_obar; //I_K p_Kbar = 1 / (1 + exp(-(V_m - V_1_2) / k)); p_K = p_Kbar; q_bar = 1.0 / (1 + exp((V_m + 40) / 14)); q_1 = q_bar; q_2 = q_bar; //IP3 receptor h_IP3 = K_inhIP3 / (Ca_i + K_inhIP3); //norepinephrine receptor PIP2 = PIP2_T - (1 + k_degG / r_rG) * gamma_G*IP3; r_hG = k_degG * gamma_G * IP3 / PIP2; G = (K_cG + Ca_i) / (alpha_G * Ca_i) * r_hG; delta_G = k_dG * G / (k_aG * (G_TG - G)); Y[0] = V_m; Y[1] = d_L; Y[2] = f_L; Y[3] = p_K; Y[4] = q_1; Y[5] = p_f; Y[6] = p_s; Y[7] = R_01; Y[8] = R_10; Y[9] = R_11; Y[10] = Ca_u; Y[11] = Ca_r; Y[12] = Ca_i; Y[13] = Na_i; Y[14] = K_i; Y[15] = q_2; Y[16] = P_SOC; Y[17] = Cl_i; Y[18] = h_IP3; Y[19] = RS_G; Y[20] = RS_PG; Y[21] = G; Y[22] = IP3; Y[23] = PIP2; Y[24] = DAG; Y[25] = V_cGMP; Y[26] = cGMP; ISTEP = -1; T = 0.0 - TSTEP; fileY.open("Y.txt"); fileY1.open("Y1.txt"); fileT.open("T.txt"); for (;;) { ISTEP = ISTEP + 1; T = T + TSTEP; //Norepinephrine if (T > 10000) NE = 1e-3; //NE [mM] beginning of stimulation if (T > 70000) NE = 0; //end of stimulation //Nitric oxide //IF (T>30000) NO = 1D-3 //NO [mM] //IF (T>70000) NO = 0 //Extracellular potassium //IF (T>10000) K_o = 30 //IF (T>70000) K_o = 5 //Current //IF (T>10000) I_stim = 5 //I_stim [pA] current injection //IF (T>40000) I_stim = -5 //IF (T>70000) I_stim = 0 // For the time being, I just interested in Y[0] values (which is V_m actually) fileY << Y[0]; fileY << "\t"; if ((ISTEP % ISTEPPRINT) == 0) { // for (int i=0; i< N; i++) { // fileY << Y[i]; // fileY << "\t"; // } // fileY << endl; // for (int i=0; i< N1; i++) { // fileY1 << Y1[i]; // fileY1 << "\t"; // } // fileY1 << endl; // // // // fileT << T; // fileT << "\t"; } FCN(); for (int i = 0; i < N; i++) { Y[i] = Y[i] + TSTEP * YPRIME[i]; } // disp(Yconcat(1)) if (ISTEP == MXSTEP) break; } cout << "It is done!" << endl; cout << Y[0] << endl; fileY.close(); fileY1.close(); fileT.close(); return 0; } void FCN(void) { V_m = Y[0]; d_L = Y[1]; f_L = Y[2]; p_K = Y[3]; q_1 = Y[4]; p_f = Y[5]; p_s = Y[6]; R_01 = Y[7]; R_10 = Y[8]; R_11 = Y[9]; Ca_u = Y[10]; Ca_r = Y[11]; Ca_i = Y[12]; Na_i = Y[13]; K_i = Y[14]; q_2 = Y[15]; P_SOC = Y[16]; Cl_i = Y[17]; h_IP3 = Y[18]; RS_G = Y[19]; RS_PG = Y[20]; G = Y[21]; IP3 = Y[22]; PIP2 = Y[23]; DAG = Y[24]; V_cGMP = Y[25]; cGMP = Y[26]; //-------------------------------------- Model equations --------------------------------------------- //cGMP formation V_cGMPbar = V_cGMPmax * (B5sGC * NO + pow(NO, 2)) / (A0sGC + A1sGC * NO + pow(NO, 2)); if ((V_cGMPbar - V_cGMP) >= 0) { tausGC = 1 / (k3sGC * NO + kDact_deactsGC); //kDact_deactsGC different from original kDsGC to uncouple Km from time constants } else { tausGC = 1 / (kDact_deactsGC + k_2sGC); } V_cGMPprime = (V_cGMPbar - V_cGMP) / tausGC; cGMPprime = V_cGMP - kpde_cGMP * cGMP * cGMP / (1e-3 + cGMP); //Norepinephrine receptor RS_Gprime = (k_rG * ksi_G * R_TG - (k_rG + k_pG * NE / (K_1G + NE)) * RS_G - k_rG * RS_PG); RS_PGprime = NE * (k_pG * RS_G / (K_1G + NE) - k_eG * RS_PG / (K_2G + NE)); rho_rG = NE * RS_G / (ksi_G * R_TG * (K_1G + NE)); Gprime = k_aG * (delta_G + rho_rG)*(G_TG - G) - k_dG*G; r_hG = alpha_G * Ca_i / (K_cG + Ca_i) * G; IP3prime = r_hG / gamma_G * PIP2 - k_degG*IP3; PIP2prime = -(r_hG + r_rG) * PIP2 - r_rG * gamma_G * IP3 + r_rG*PIP2_T; DAGprime = r_hG / gamma_G * PIP2 - k_degG*DAG; //Reversal potentials E_Ca = RT_F / z_Ca * log(Ca_o / Ca_i); E_Na = RT_F * log(Na_o / Na_i); E_K = RT_F * log(K_o / K_i); E_Cl = RT_F / z_Cl * log(Cl_o / Cl_i); //Voltage dependent calcium current I_CaL tau_d_L = 2.5 * exp(-pow((V_m + 40) / 30, 2)) + 1.15; d_Lbar = 1.0 / (1 + exp(-(V_m) / 8.3)); d_Lprime = (d_Lbar - d_L) / tau_d_L; f_Lbar = 1.0 / (1 + exp((V_m + 42.0) / 9.1)); tau_f_L = 65 * exp(-pow((V_m + 35) / 25, 2)) + 45; f_Lprime = (f_Lbar - f_L) / tau_f_L; if (V_m == 0) { I_CaL = d_L * f_L * P_CaL * A_m * 1e6 * z_Ca * F * (Ca_i - Ca_o); //[pA] } else { I_CaL = d_L * f_L * P_CaL * A_m * 1e6 * V_m * pow(z_Ca * F, 2) / (R * temp)*(Ca_o - Ca_i * exp(V_m * z_Ca / (RT_F))) / (1 - exp(V_m * z_Ca / (RT_F))); //[pA] } //Delayed rectifier current I_K p_Kbar = 1 / (1 + exp(-(V_m - V_1_2) / k)); tau_p_K = 61.49 * exp(-0.0268 * V_m); p_Kprime = (p_Kbar - p_K) / tau_p_K; q_bar = 1.0 / (1 + exp((V_m + 40) / 14)); q_1prime = (q_bar - q_1) / 371; q_2prime = (q_bar - q_2) / 2884; I_K = 1 * g_K * p_K * (0.45 * q_1 + 0.55 * q_2) * (V_m - E_K); //Alpha-adrenoceptor-activated nonselective cation channel NSC Po_NSC = Pmin_NSC + (1 - Pmin_NSC) / (1 + exp(-(V_m - 47.12) / 24.24)); if (V_m == 0) { INa_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PNa_NSC * A_m * 1e6 * F * (Na_i - Na_o); ICa_NSC = 1 * (0 * DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PCa_NSC * A_m * 1e6 * z_Ca * F * (Ca_i - Ca_o); IK_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PK_NSC * A_m * 1e6 * F * (K_i - K_o); } else { INa_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PNa_NSC * A_m * 1e6 * V_m * pow(F, 2) / (R * temp)*(Na_o - Na_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F)); ICa_NSC = 1 * (0 * DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PCa_NSC * A_m * 1e6 * V_m * pow(z_Ca * F, 2) / (R * temp)*(Ca_o - Ca_i * exp(V_m * z_Ca / RT_F)) / (1 - exp(V_m * z_Ca / RT_F)); IK_NSC = 1 * (DAG / (DAG + K_NSC) + d_NSCmin) * Po_NSC * PK_NSC * A_m * 1e6 * V_m * pow(F, 2) / (R * temp)*(K_o - K_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F)); } I_NSC = ICa_NSC + INa_NSC + IK_NSC; //KATP current I_KATP I_KATP = g_KATP * (V_m - E_K); //Inward rectifier current I_K_i g_maxK_i = G_K_i * pow(K_o, n_K_i); I_K_i = g_maxK_i * (V_m - E_K) / (1 + exp((V_m - E_K) / 28.89)); //Calcium-activated potassium current I_KCa if (V_m == 0) { i1_KCa = 1e6 * P_BKCa * F * (K_i - K_o); //[pA] } else { i1_KCa = 1e6 * P_BKCa * V_m * F / RT_F * (K_o - K_i * exp(V_m / RT_F)) / (1 - exp(V_m / RT_F)); //[pA] } //Mistry and Garland 1998 R_NO = NO / (NO + 200e-6); R_cGMP = pow(cGMP, 2) / (pow(cGMP, 2) + pow(1.5e-3, 2)); V_1_2_KCa = -41.7 * log10(Ca_i) - 128.2 - dV_1_2_KCa_NO * R_NO - dV_1_2_KCa_cGMP*R_cGMP; p_obar = 1 / (1 + exp(-(V_m - V_1_2_KCa) / 18.25)); p_fprime = (p_obar - p_f) / tau_pf; p_sprime = (p_obar - p_s) / tau_ps; P_KCa = 0.17 * p_f + 0.83 * p_s; I_KCa = A_m * N_BKCa * i1_KCa * P_KCa; //Store operated non-selective cation channel P_SOCbar = 1 / (1 + pow(Ca_u / K_SOC, H_SOC)); P_SOCprime = (P_SOCbar - P_SOC) / tau_SOC; I_SOCCa = 1 * P_SOC * g_SOCCa * (V_m - E_Ca); I_SOCNa = 1 * P_SOC * g_SOCNa * (V_m - E_Na); I_SOC = I_SOCCa + I_SOCNa; //Chloride currents alpha_Cl = pow(cGMP, 3.3) / (pow(cGMP, 3.3) + pow(6.4e-3, 3.3)); P_Cl = pow(Ca_i, 2) / (pow(Ca_i, 2) + pow(365e-6, 2)) * 0.0132 + pow(Ca_i, 2) / (pow(Ca_i, 2) + pow(400e-6 * (1 - alpha_Cl * 0.9), 2)) * alpha_Cl; I_Cl = P_Cl * g_Cl * C_m * (V_m - E_Cl); //IP3 receptor current h_IP3prime = k_onIP3 * (K_inhIP3 - (Ca_i + K_inhIP3) * h_IP3); I_IP3 = I_IP3bar * pow(IP3 / (IP3 + K_IP3) * Ca_i / (Ca_i + K_actIP3) * h_IP3, 3)*(Ca_u - Ca_i) * z_Ca * F*vol_Ca; //Calcium-induced calcium release R_00 = 1 - R_01 - R_10 - R_11; R_10prime = Kr1 * pow(Ca_i, 2) * R_00 - (K_r1 + Kr2 * Ca_i) * R_10 + K_r2*R_11; R_11prime = Kr2 * Ca_i * R_10 - (K_r1 + K_r2) * R_11 + Kr1 * pow(Ca_i, 2) * R_01; R_01prime = Kr2 * Ca_i * R_00 + K_r1 * R_11 - (K_r2 + Kr1 * pow(Ca_i, 2)) * R_01; I_up = I_upbar * Ca_i / (Ca_i + K_mup); I_tr = (Ca_u - Ca_r) * (2 * F * vol_u) / tau_tr; I_rel = (pow(R_10, 2) + R_leak) * (Ca_r - Ca_i) * (2 * F * vol_r) / tau_rel; Ca_uprime = (I_up - I_tr - I_IP3) / (2 * F * vol_u); Ca_rprime = (I_tr - I_rel) / (2 * F * vol_r) / (1 + CSQNbar * K_CSQN / pow((K_CSQN + Ca_r), 2)); //Ca buffering and cytosolic material balance S_CM = S_CMbar * K_d / (K_d + Ca_i); CaCM = S_CMbar - S_CM; I_PMCA = I_PMCAbar * Ca_i / (Ca_i + K_mPMCA); //NaK pump I_NaK = pow(Q_10_NaK, ((temp - 309.15) / 10)) * C_m * I_NaKbar * ((pow(K_o, 1.1)) / (pow(K_o, 1.1) + (pow(K_mK, 1.1))) *(pow(Na_i, 1.7)) / ((pow(Na_i, 1.7))+(pow(K_mNa, 1.7)))) *(V_m + 150) / (V_m + 200); Fi_F = exp(gamma2 * V_m * F / (R * temp)); Fi_R = exp((gamma2 - 1) * V_m * F / (R * temp)); I_NCX = 1 * (1 + 0.55 * cGMP / (cGMP + (45e-3))) * g_NCX * (pow(Na_i, 3) * Ca_o * Fi_F - pow(Na_o, 3) * Ca_i * Fi_R) / (1 + d_NCX * (pow(Na_o, 3) * Ca_i + pow(Na_i, 3) * Ca_o)); I_NaKCl_Cl = (1 + 7 / 2 * cGMP / (cGMP + 6.4e-3))*(-A_m * L_cotr * R * temp * z_Cl * F * log(Na_o / Na_i * K_o / K_i * pow(Cl_o / Cl_i, 2))); I_Catotm = I_SOCCa + I_CaL - 2 * I_NCX + I_PMCA + ICa_NSC + I_MCa; Ca_iprime = -(I_Catotm + I_up - I_rel - I_IP3) / (2 * F * vol_Ca) / (1 + S_CMbar * K_d / (pow(K_d + Ca_i, 2)) + B_Fbar * K_dB / (pow(K_dB + Ca_i, 2))); I_Natotm = -0.5 * I_NaKCl_Cl + I_SOCNa + 3 * I_NaK + 3 * I_NCX + INa_NSC + I_MNa; Na_iprime = -(I_Natotm) / (F * vol_i); I_Ktotm = -0.5 * I_NaKCl_Cl + I_K + I_KCa + I_K_i + IK_NSC + I_KATP - 2 * I_NaK + I_MK; K_iprime = -(I_Ktotm) / (F * vol_i); I_Cltotm = I_NaKCl_Cl + I_Cl; Cl_iprime = -(I_Cltotm) / (z_Cl * F * vol_i); //Transmembrane potential V_mprime = -1 / C_m * (-I_stim + I_Cl + I_SOC + I_CaL + I_K + I_KCa + I_K_i + I_M + I_NCX + I_NaK + I_PMCA + I_NSC + I_KATP); //YPRIME = zeros(1, N); YPRIME[0] = V_mprime; YPRIME[1] = d_Lprime; YPRIME[2] = f_Lprime; YPRIME[3] = p_Kprime; YPRIME[4] = q_1prime; YPRIME[5] = p_fprime; YPRIME[6] = p_sprime; YPRIME[7] = R_01prime; YPRIME[8] = R_10prime; YPRIME[9] = R_11prime; YPRIME[10] = Ca_uprime; YPRIME[11] = Ca_rprime; YPRIME[12] = Ca_iprime; YPRIME[13] = Na_iprime; YPRIME[14] = K_iprime; YPRIME[15] = q_2prime; YPRIME[16] = P_SOCprime; YPRIME[17] = Cl_iprime; YPRIME[18] = h_IP3prime; YPRIME[19] = RS_Gprime; YPRIME[20] = RS_PGprime; YPRIME[21] = Gprime; YPRIME[22] = IP3prime; YPRIME[23] = PIP2prime; YPRIME[24] = DAGprime; YPRIME[25] = V_cGMPprime; YPRIME[26] = cGMPprime; //Non state variables Y1[0] = I_CaL; Y1[1] = I_K; Y1[2] = I_K_i; Y1[3] = I_NSC; Y1[4] = I_KCa; Y1[5] = I_up; Y1[6] = I_rel; Y1[7] = I_PMCA; Y1[8] = I_NCX; Y1[9] = I_NaK; Y1[10] = INa_NSC; Y1[11] = ICa_NSC; Y1[12] = IK_NSC; Y1[13] = I_SOCCa; Y1[14] = I_SOCNa; Y1[15] = I_Cl; Y1[16] = I_NaKCl_Cl; Y1[17] = I_IP3; Y1[18] = I_tr; Y1[19] = NE; Y1[20] = I_KATP; Y1[21] = I_stim; Y1[22] = V_cGMPbar; Y1[23] = NO; Y1[29] = I_Catotm; Y1[30] = I_Natotm; Y1[31] = I_Cltotm; Y1[32] = I_Ktotm; }