我的C++代码在MS C++编译器上运行得很好,但是我给了Non on g++编译器。为什么?

我用C++模拟了一个涉及27个刚性常微分方程的生物模型。我的程序在MS C++ 2010表达式编译器下运行得很好,但是在G++编译器(NETBeaS 6.8,Ubuntu 10.04 LTS)下失败了。问题是有些变量 NaN . 以下是变量的值 Vm

这里是MSC C++编译器下没有任何更改的相同代码的输出:

-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 -

我只使用了“cmath”和“fstream”库。

编辑1:

好了,伙计们,下面是全部代码:

#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;
}