BackwardEulerMahajanModel2008.cpp

Go to the documentation of this file.

#include "BackwardEulerMahajanModel2008.hpp"
#include <cmath>
#include <cassert>
#include "CardiacNewtonSolver.hpp"
#include "Exception.hpp"
#include "OdeSystemInformation.hpp"

    BackwardEulerMahajanModel2008::BackwardEulerMahajanModel2008(boost::shared_ptr<AbstractIvpOdeSolver> /* unused; should be empty */, boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
        : AbstractBackwardEulerCardiacCell<15>(
                26,
                0,
                pIntracellularStimulus)
    {
        // Time units: ms
        //
        mpSystemInfo = OdeSystemInformation<BackwardEulerMahajanModel2008>::Instance();
        Init();

    }

    BackwardEulerMahajanModel2008::~BackwardEulerMahajanModel2008()
    {
    }

    void BackwardEulerMahajanModel2008::VerifyStateVariables()
    {}

    double BackwardEulerMahajanModel2008::GetIIonic()
    {
        std::vector<double>& rY = rGetStateVariables();
        double var_cell__V = rY[0];
        // Units: mV; Initial value: -87.169816169406
        double var_INa__xm = rY[1];
        // Units: dimensionless; Initial value: 0.001075453357
        double var_INa__xh = rY[2];
        // Units: dimensionless; Initial value: 0.990691306716
        double var_INa__xj = rY[3];
        // Units: dimensionless; Initial value: 0.993888937283
        double var_ICaL__c1 = rY[4];
        // Units: dimensionless; Initial value: 0.000018211252
        double var_ICaL__c2 = rY[5];
        // Units: dimensionless; Initial value: 0.979322592773
        double var_ICaL__xi1ca = rY[6];
        // Units: dimensionless; Initial value: 0.001208153482
        double var_ICaL__xi1ba = rY[7];
        // Units: dimensionless; Initial value: 0.000033616596
        double var_ICaL__xi2ca = rY[8];
        // Units: dimensionless; Initial value: 0.004173008466
        double var_ICaL__xi2ba = rY[9];
        // Units: dimensionless; Initial value: 0.015242594688
        double var_IKr__xr = rY[10];
        // Units: dimensionless; Initial value: 0.007074239331
        double var_IKs__xs1 = rY[11];
        // Units: dimensionless; Initial value: 0.048267587131
        double var_IKs__xs2 = rY[12];
        // Units: dimensionless; Initial value: 0.105468807033
        double var_Ito__xtos = rY[13];
        // Units: dimensionless; Initial value: 0.00364776906
        double var_Ito__ytos = rY[14];
        // Units: dimensionless; Initial value: 0.174403618112
        double var_Ito__xtof = rY[15];
        // Units: dimensionless; Initial value: 0.003643592594
        double var_Ito__ytof = rY[16];
        // Units: dimensionless; Initial value: 0.993331326442
        double var_Na__Na_i = rY[19];
        // Units: mM; Initial value: 11.441712311614
        double var_Ca__Ca_submem = rY[21];
        // Units: uM; Initial value: 0.226941113355
        double var_Ca__Ca_i = rY[22];
        // Units: uM; Initial value: 0.256752008084

        const double var_Environment__R = 8.314472;
        const double var_Environment__T = 308.0;
        const double var_Environment__F = 96.4853415;
        const double var_Environment__K_o = 5.4;
        const double var_Environment__Ca_o = 1.8;
        const double var_Environment__Na_o = 136.0;
        const double var_Environment__FonRT = var_Environment__F / (var_Environment__R * var_Environment__T);
        const double var_reversal_potentials__K_i = 140.0;
        const double var_reversal_potentials__FonRT = var_Environment__FonRT;
        const double var_reversal_potentials__K_o = var_Environment__K_o;
        const double var_reversal_potentials__ek = (1.0 / var_reversal_potentials__FonRT) * log(var_reversal_potentials__K_o / var_reversal_potentials__K_i);
        const double var_IK1__ek = var_reversal_potentials__ek;
        const double var_IK1__V = var_cell__V;
        const double var_IK1__aki = 1.02 / (1.0 + exp(0.2385 * ((var_IK1__V - var_IK1__ek) - 59.215)));
        const double var_IK1__bki = ((0.49124 * exp(0.08032 * ((var_IK1__V - var_IK1__ek) + 5.476))) + (1.0 * exp(0.06175 * ((var_IK1__V - var_IK1__ek) - 594.31)))) / (1.0 + exp((-0.5143) * ((var_IK1__V - var_IK1__ek) + 4.753)));
        const double var_IK1__xkin = var_IK1__aki / (var_IK1__aki + var_IK1__bki);
        const double var_IK1__K_o = var_Environment__K_o;
        const double var_IK1__gkix = 0.3;
        const double var_IK1__xik1 = var_IK1__gkix * sqrt(var_IK1__K_o / 5.4) * var_IK1__xkin * (var_IK1__V - var_IK1__ek);
        const double var_cell__xik1 = var_IK1__xik1;
        const double var_Ito__ek = var_reversal_potentials__ek;
        const double var_Ito__gtos = 0.04;
        const double var_Ito__V = var_cell__V;
        const double var_Ito__rt2 = (var_Ito__V + 33.5) / 10.0;
        const double var_Ito__rs_inf = 1.0 / (1.0 + exp(var_Ito__rt2));
        const double var_Ito__xitos = var_Ito__gtos * var_Ito__xtos * (var_Ito__ytos + (0.5 * var_Ito__rs_inf)) * (var_Ito__V - var_Ito__ek);
        const double var_Ito__gtof = 0.11;
        const double var_Ito__xitof = var_Ito__gtof * var_Ito__xtof * var_Ito__ytof * (var_Ito__V - var_Ito__ek);
        const double var_Ito__xito = var_Ito__xitos + var_Ito__xitof;
        const double var_cell__xito = var_Ito__xito;
        const double var_INaK__xkmko = 1.5;
        const double var_INaK__V = var_cell__V;
        const double var_INaK__FonRT = var_Environment__FonRT;
        const double var_INaK__Na_o = var_Environment__Na_o;
        const double var_INaK__sigma = (exp(var_INaK__Na_o / 67.3) - 1.0) / 7.0;
        const double var_INaK__fNaK = 1.0 / (1.0 + (0.1245 * exp((-0.1) * var_INaK__V * var_INaK__FonRT)) + (0.0365 * var_INaK__sigma * exp((-var_INaK__V) * var_INaK__FonRT)));
        const double var_INaK__Na_i = var_Na__Na_i;
        const double var_INaK__K_o = var_Environment__K_o;
        const double var_INaK__xkmnai = 12.0;
        const double var_INaK__gNaK = 1.5;
        const double var_INaK__xiNaK = (((var_INaK__gNaK * var_INaK__fNaK * var_INaK__Na_i) / (var_INaK__Na_i + var_INaK__xkmnai)) * var_INaK__K_o) / (var_INaK__K_o + var_INaK__xkmko);
        const double var_cell__xiNaK = var_INaK__xiNaK;
        const double var_cell__wca = 8.0;
        const double var_INaCa__gNaCa = 0.84;
        const double var_INaCa__V = var_cell__V;
        const double var_INaCa__FonRT = var_Environment__FonRT;
        const double var_INaCa__zw4 = 1.0 + (0.2 * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__xkdna = 0.3;
        const double var_INaCa__Ca_submem = var_Ca__Ca_submem;
        const double var_INaCa__aloss = 1.0 / (1.0 + pow(var_INaCa__xkdna / var_INaCa__Ca_submem, 3.0));
        const double var_INaCa__Na_o = var_Environment__Na_o;
        const double var_INaCa__Na_i = var_Na__Na_i;
        const double var_INaCa__Ca_o = var_Environment__Ca_o;
        const double var_Ca__csm = var_Ca__Ca_submem / 1000.0;
        const double var_INaCa__csm = var_Ca__csm;
        const double var_INaCa__yz4 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o) + (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcao = 1.3;
        const double var_INaCa__xmnao = 87.5;
        const double var_INaCa__yz1 = (var_INaCa__xmcao * pow(var_INaCa__Na_i, 3.0)) + (pow(var_INaCa__xmnao, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcai = 0.0036;
        const double var_INaCa__xmnai = 12.3;
        const double var_INaCa__yz3 = var_INaCa__xmcai * pow(var_INaCa__Na_o, 3.0) * (1.0 + pow(var_INaCa__Na_i / var_INaCa__xmnai, 3.0));
        const double var_INaCa__yz2 = pow(var_INaCa__xmnai, 3.0) * var_INaCa__Ca_o * (1.0 + (var_INaCa__csm / var_INaCa__xmcai));
        const double var_INaCa__zw8 = var_INaCa__yz1 + var_INaCa__yz2 + var_INaCa__yz3 + var_INaCa__yz4;
        const double var_INaCa__zw3 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o * exp(var_INaCa__V * 0.35 * var_INaCa__FonRT)) - (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__jNaCa = (var_INaCa__gNaCa * var_INaCa__aloss * var_INaCa__zw3) / (var_INaCa__zw4 * var_INaCa__zw8);
        const double var_INaCa__wca = var_cell__wca;
        const double var_INaCa__xiNaCa = var_INaCa__wca * var_INaCa__jNaCa;
        const double var_cell__xiNaCa = var_INaCa__xiNaCa;
        const double var_ICaL__Ca_o = var_Environment__Ca_o;
        const double var_ICaL__V = var_cell__V;
        const double var_ICaL__pca = 0.00054;
        const double var_ICaL__FonRT = var_Environment__FonRT;
        const double var_ICaL__F = var_Environment__F;
        const double var_ICaL__za = var_ICaL__V * 2.0 * var_ICaL__FonRT;
        const double var_ICaL__csm = var_Ca__csm;
        const double var_ICaL__rxa = (fabs(var_ICaL__za) < 0.001) ? ((4.0 * var_ICaL__pca * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (2.0 * var_ICaL__FonRT)) : ((4.0 * var_ICaL__pca * var_ICaL__V * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0));
        const double var_ICaL__gca = 182.0;
        const double var_ICaL__po = (((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2;
        const double var_ICaL__jca = var_ICaL__gca * var_ICaL__po * var_ICaL__rxa;
        const double var_ICaL__wca = var_cell__wca;
        const double var_ICaL__xica = 2.0 * var_ICaL__wca * var_ICaL__jca;
        const double var_cell__xica = var_ICaL__xica;
        const double var_IKr__gkr = 0.0125;
        const double var_IKr__K_o = var_Environment__K_o;
        const double var_IKr__V = var_cell__V;
        const double var_IKr__rg = 1.0 / (1.0 + exp((var_IKr__V + 33.0) / 22.4));
        const double var_IKr__ek = var_reversal_potentials__ek;
        const double var_IKr__xikr = var_IKr__gkr * sqrt(var_IKr__K_o / 5.4) * var_IKr__xr * var_IKr__rg * (var_IKr__V - var_IKr__ek);
        const double var_cell__xikr = var_IKr__xikr;
        const double var_IKs__Ca_i = var_Ca__Ca_i;
        const double var_IKs__gksx = 1.0 + (0.8 / (1.0 + pow(0.5 / var_IKs__Ca_i, 3.0)));
        const double var_IKs__gks = 0.1386;
        const double var_IKs__V = var_cell__V;
        const double var_reversal_potentials__prNaK = 0.01833;
        const double var_reversal_potentials__Na_o = var_Environment__Na_o;
        const double var_reversal_potentials__Na_i = var_Na__Na_i;
        const double var_reversal_potentials__eks = (1.0 / var_reversal_potentials__FonRT) * log((var_reversal_potentials__K_o + (var_reversal_potentials__prNaK * var_reversal_potentials__Na_o)) / (var_reversal_potentials__K_i + (var_reversal_potentials__prNaK * var_reversal_potentials__Na_i)));
        const double var_IKs__eks = var_reversal_potentials__eks;
        const double var_IKs__xiks = var_IKs__gks * var_IKs__gksx * var_IKs__xs1 * var_IKs__xs2 * (var_IKs__V - var_IKs__eks);
        const double var_cell__xiks = var_IKs__xiks;
        const double var_INa__V = var_cell__V;
        const double var_INa__gna = 12.0;
        const double var_reversal_potentials__ena = (1.0 / var_reversal_potentials__FonRT) * log(var_reversal_potentials__Na_o / var_reversal_potentials__Na_i);
        const double var_INa__ena = var_reversal_potentials__ena;
        const double var_INa__xina = var_INa__gna * var_INa__xh * var_INa__xj * var_INa__xm * var_INa__xm * var_INa__xm * (var_INa__V - var_INa__ena);
        const double var_cell__xina = var_INa__xina;

        /*
         * Currents in this model are given in microA/microF (as the results of conductance (mS/microF) times voltage (mV)).
         * Hence, no conversion factor is necessary because this term will be multiplied by Cm (in microF/cm2)
         * giving rise to mono/bidomain currents in the correct units.
         *
         * Note that the CellML file in the repository was coded up with currents as having "dimensionless" units and therefore Pycml raised a warning about the units of this.
         */
        return (var_cell__xik1+var_cell__xito+var_cell__xiNaK+var_cell__xiNaCa+var_cell__xica+var_cell__xina+var_cell__xikr+var_cell__xiks);
    }

    void BackwardEulerMahajanModel2008::ComputeResidual(double var_environment__time, const double rCurrentGuess[15], double rResidual[15])
    {
        std::vector<double>& rY = rGetStateVariables();
        double var_cell__V = rY[0];
        // Units: mV; Initial value: -87.169816169406
        double var_INa__xm = rY[1];
        // Units: dimensionless; Initial value: 0.001075453357
        double var_INa__xh = rY[2];
        // Units: dimensionless; Initial value: 0.990691306716
        double var_INa__xj = rY[3];
        // Units: dimensionless; Initial value: 0.993888937283

        double var_Ca__Ca_NSR = rCurrentGuess[0];
        double var_Ca__Ca_dyad = rCurrentGuess[1];
        double var_Ca__Ca_i = rCurrentGuess[2];
        double var_Ca__Ca_submem = rCurrentGuess[3];
        double var_Ca__tropi = rCurrentGuess[4];
        double var_Ca__trops = rCurrentGuess[5];
        double var_ICaL__c1 = rCurrentGuess[6];
        double var_ICaL__c2 = rCurrentGuess[7];
        double var_ICaL__xi1ba = rCurrentGuess[8];
        double var_ICaL__xi1ca = rCurrentGuess[9];
        double var_ICaL__xi2ba = rCurrentGuess[10];
        double var_ICaL__xi2ca = rCurrentGuess[11];
        double var_Irel__Ca_JSR = rCurrentGuess[12];
        double var_Irel__xir = rCurrentGuess[13];
        double var_Na__Na_i = rCurrentGuess[14];

        const double var_Environment__R = 8.314472;
        const double var_Environment__T = 308.0;
        const double var_Environment__F = 96.4853415;
        const double var_Environment__K_o = 5.4;
        const double var_Environment__Ca_o = 1.8;
        const double var_Environment__Na_o = 136.0;
        const double var_Environment__FonRT = var_Environment__F / (var_Environment__R * var_Environment__T);
        const double var_reversal_potentials__FonRT = var_Environment__FonRT;
        const double var_INaK__xkmko = 1.5;
        const double var_INaK__V = var_cell__V;
        const double var_INaK__FonRT = var_Environment__FonRT;
        const double var_INaK__Na_o = var_Environment__Na_o;
        const double var_INaK__sigma = (exp(var_INaK__Na_o / 67.3) - 1.0) / 7.0;
        const double var_INaK__fNaK = 1.0 / (1.0 + (0.1245 * exp((-0.1) * var_INaK__V * var_INaK__FonRT)) + (0.0365 * var_INaK__sigma * exp((-var_INaK__V) * var_INaK__FonRT)));
        const double var_INaK__Na_i = var_Na__Na_i;
        const double var_INaK__K_o = var_Environment__K_o;
        const double var_INaK__xkmnai = 12.0;
        const double var_INaK__gNaK = 1.5;
        const double var_INaK__xiNaK = (((var_INaK__gNaK * var_INaK__fNaK * var_INaK__Na_i) / (var_INaK__Na_i + var_INaK__xkmnai)) * var_INaK__K_o) / (var_INaK__K_o + var_INaK__xkmko);
        const double var_cell__wca = 8.0;
        const double var_INaCa__gNaCa = 0.84;
        const double var_INaCa__V = var_cell__V;
        const double var_INaCa__FonRT = var_Environment__FonRT;
        const double var_INaCa__zw4 = 1.0 + (0.2 * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__xkdna = 0.3;
        const double var_INaCa__Ca_submem = var_Ca__Ca_submem;
        const double var_INaCa__aloss = 1.0 / (1.0 + pow(var_INaCa__xkdna / var_INaCa__Ca_submem, 3.0));
        const double var_INaCa__Na_o = var_Environment__Na_o;
        const double var_INaCa__Na_i = var_Na__Na_i;
        const double var_INaCa__Ca_o = var_Environment__Ca_o;
        const double var_Ca__csm = var_Ca__Ca_submem / 1000.0;
        const double var_INaCa__csm = var_Ca__csm;
        const double var_INaCa__yz4 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o) + (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcao = 1.3;
        const double var_INaCa__xmnao = 87.5;
        const double var_INaCa__yz1 = (var_INaCa__xmcao * pow(var_INaCa__Na_i, 3.0)) + (pow(var_INaCa__xmnao, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcai = 0.0036;
        const double var_INaCa__xmnai = 12.3;
        const double var_INaCa__yz3 = var_INaCa__xmcai * pow(var_INaCa__Na_o, 3.0) * (1.0 + pow(var_INaCa__Na_i / var_INaCa__xmnai, 3.0));
        const double var_INaCa__yz2 = pow(var_INaCa__xmnai, 3.0) * var_INaCa__Ca_o * (1.0 + (var_INaCa__csm / var_INaCa__xmcai));
        const double var_INaCa__zw8 = var_INaCa__yz1 + var_INaCa__yz2 + var_INaCa__yz3 + var_INaCa__yz4;
        const double var_INaCa__zw3 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o * exp(var_INaCa__V * 0.35 * var_INaCa__FonRT)) - (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__jNaCa = (var_INaCa__gNaCa * var_INaCa__aloss * var_INaCa__zw3) / (var_INaCa__zw4 * var_INaCa__zw8);
        const double var_INaCa__wca = var_cell__wca;
        const double var_INaCa__xiNaCa = var_INaCa__wca * var_INaCa__jNaCa;
        const double var_ICaL__Ca_o = var_Environment__Ca_o;
        const double var_ICaL__V = var_cell__V;
        const double var_ICaL__pca = 0.00054;
        const double var_ICaL__FonRT = var_Environment__FonRT;
        const double var_ICaL__F = var_Environment__F;
        const double var_ICaL__za = var_ICaL__V * 2.0 * var_ICaL__FonRT;
        const double var_ICaL__csm = var_Ca__csm;
        const double var_ICaL__rxa = (fabs(var_ICaL__za) < 0.001) ? ((4.0 * var_ICaL__pca * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (2.0 * var_ICaL__FonRT)) : ((4.0 * var_ICaL__pca * var_ICaL__V * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0));
        const double var_ICaL__gca = 182.0;
        const double var_ICaL__po = (((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2;
        const double var_ICaL__jca = var_ICaL__gca * var_ICaL__po * var_ICaL__rxa;
        const double var_reversal_potentials__Na_o = var_Environment__Na_o;
        const double var_reversal_potentials__Na_i = var_Na__Na_i;
        const double var_INa__V = var_cell__V;
        const double var_INa__gna = 12.0;
        const double var_reversal_potentials__ena = (1.0 / var_reversal_potentials__FonRT) * log(var_reversal_potentials__Na_o / var_reversal_potentials__Na_i);
        const double var_INa__ena = var_reversal_potentials__ena;
        const double var_INa__xina = var_INa__gna * var_INa__xh * var_INa__xj * var_INa__xm * var_INa__xm * var_INa__xm * (var_INa__V - var_INa__ena);
        const double var_ICaL__Ca_dyad = var_Ca__Ca_dyad;
        const double var_ICaL__vth = 0.0;
        const double var_ICaL__s6 = 8.0;
        const double var_ICaL__poinf = 1.0 / (1.0 + exp((-(var_ICaL__V - var_ICaL__vth)) / var_ICaL__s6));
        const double var_ICaL__cat = 3.0;
        const double var_ICaL__fca = 1.0 / (1.0 + pow(var_ICaL__cat / var_ICaL__Ca_dyad, 3.0));
        const double var_ICaL__vx =  -40.0;
        const double var_ICaL__sx = 3.0;
        const double var_ICaL__vy =  -40.0;
        const double var_ICaL__sy = 4.0;
        const double var_ICaL__vyr =  -40.0;
        const double var_ICaL__syr = 11.32;
        const double var_ICaL__cpt = 6.09365;
        const double var_ICaL__taupo = 1.0;
        const double var_ICaL__alpha = var_ICaL__poinf / var_ICaL__taupo;
        const double var_ICaL__beta = (1.0 - var_ICaL__poinf) / var_ICaL__taupo;
        const double var_ICaL__k1 = 0.024168 * var_ICaL__fca;
        const double var_ICaL__k2 = 0.000103615;
        const double var_ICaL__k1t = 0.00413;
        const double var_ICaL__k2t = 0.00224;
        const double var_ICaL__poi = 1.0 / (1.0 + exp((-(var_ICaL__V - var_ICaL__vx)) / var_ICaL__sx));
        const double var_ICaL__tau3 = 3.0;
        const double var_ICaL__k3 = (1.0 - var_ICaL__poi) / var_ICaL__tau3;
        const double var_ICaL__k3t = var_ICaL__k3;
        const double var_ICaL__Ps = 1.0 / (1.0 + exp((-(var_ICaL__V - var_ICaL__vyr)) / var_ICaL__syr));
        const double var_ICaL__tca = 78.0329;
        const double var_ICaL__tau_ca = (var_ICaL__tca / (1.0 + pow(var_ICaL__Ca_dyad / var_ICaL__cpt, 4.0))) + 0.1;
        const double var_ICaL__recov = 10.0 + (4954.0 * exp(var_ICaL__V / 15.6));
        const double var_ICaL__Pr = 1.0 - (1.0 / (1.0 + exp((-(var_ICaL__V - var_ICaL__vy)) / var_ICaL__sy)));
        const double var_ICaL__tauca = ((var_ICaL__recov - var_ICaL__tau_ca) * var_ICaL__Pr) + var_ICaL__tau_ca;
        const double var_ICaL__k6 = (var_ICaL__fca * var_ICaL__Ps) / var_ICaL__tauca;
        const double var_ICaL__k5 = (1.0 - var_ICaL__Ps) / var_ICaL__tauca;
        const double var_ICaL__tauba = ((var_ICaL__recov - 450.0) * var_ICaL__Pr) + 450.0;
        const double var_ICaL__k6t = var_ICaL__Ps / var_ICaL__tauba;
        const double var_ICaL__k5t = (1.0 - var_ICaL__Ps) / var_ICaL__tauba;
        const double var_ICaL__k4 = (((((var_ICaL__k3 * var_ICaL__alpha) / var_ICaL__beta) * var_ICaL__k1) / var_ICaL__k2) * var_ICaL__k5) / var_ICaL__k6;
        const double var_ICaL__k4t = (((((var_ICaL__k3t * var_ICaL__alpha) / var_ICaL__beta) * var_ICaL__k1t) / var_ICaL__k2t) * var_ICaL__k5t) / var_ICaL__k6t;
        const double var_ICaL__r1 = 0.3;
        const double var_ICaL__r2 = 3.0;
        const double var_ICaL__s1 = 0.0182688 * var_ICaL__fca;
        const double var_ICaL__s1t = 0.00195;
        const double var_ICaL__s2 = (((var_ICaL__s1 * var_ICaL__r1) / var_ICaL__r2) * var_ICaL__k2) / var_ICaL__k1;
        const double var_ICaL__s2t = (((var_ICaL__s1t * var_ICaL__r1) / var_ICaL__r2) * var_ICaL__k2t) / var_ICaL__k1t;
        const double var_Irel__V = var_cell__V;
        const double var_Irel__Ca_NSR = var_Ca__Ca_NSR;
        const double var_Ileak_Iup_Ixfer__Ca_NSR = var_Ca__Ca_NSR;
        const double var_Ileak_Iup_Ixfer__gleak = 2.069e-05;
        const double var_Ileak_Iup_Ixfer__kj = 50.0;
        const double var_Ileak_Iup_Ixfer__Ca_i = var_Ca__Ca_i;
        const double var_Ileak_Iup_Ixfer__jleak = ((var_Ileak_Iup_Ixfer__gleak * var_Ileak_Iup_Ixfer__Ca_NSR * var_Ileak_Iup_Ixfer__Ca_NSR) / ((var_Ileak_Iup_Ixfer__Ca_NSR * var_Ileak_Iup_Ixfer__Ca_NSR) + (var_Ileak_Iup_Ixfer__kj * var_Ileak_Iup_Ixfer__kj))) * ((var_Ileak_Iup_Ixfer__Ca_NSR * 16.667) - var_Ileak_Iup_Ixfer__Ca_i);
        const double var_Ca__jleak = var_Ileak_Iup_Ixfer__jleak;
        const double var_Ileak_Iup_Ixfer__cup = 0.5;
        const double var_Ileak_Iup_Ixfer__vup = 0.4;
        const double var_Ileak_Iup_Ixfer__jup = (var_Ileak_Iup_Ixfer__vup * var_Ileak_Iup_Ixfer__Ca_i * var_Ileak_Iup_Ixfer__Ca_i) / ((var_Ileak_Iup_Ixfer__Ca_i * var_Ileak_Iup_Ixfer__Ca_i) + (var_Ileak_Iup_Ixfer__cup * var_Ileak_Iup_Ixfer__cup));
        const double var_Ca__jup = var_Ileak_Iup_Ixfer__jup;
        const double var_Ca__xir = var_Irel__xir;
        const double var_Ca__dCa_JSR = ((-var_Ca__xir) + var_Ca__jup) - var_Ca__jleak;
        const double var_Irel__dCa_JSR = var_Ca__dCa_JSR;
        const double var_Irel__po = var_ICaL__po;
        const double var_Irel__rxa = var_ICaL__rxa;
        const double var_Irel__cstar = 90.0;
        const double var_Irel__gryr = 2.58079;
        const double var_Irel__gbarsr = 26841.8;
        const double var_Irel__gdyad = 9000.0;
        const double var_Irel__ax = 0.3576;
        const double var_Irel__ay = 0.05;
        const double var_Irel__av = 11.3;
        const double var_Irel__bv = ((1.0 - var_Irel__av) * var_Irel__cstar) - 50.0;
        const double var_Irel__Qr0 = ((var_Irel__Ca_JSR > 50.0) && (var_Irel__Ca_JSR < var_Irel__cstar)) ? ((var_Irel__Ca_JSR - 50.0) / 1.0) : (var_Irel__Ca_JSR >= var_Irel__cstar) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0;
        const double var_Irel__Qr = (var_Irel__Ca_NSR * var_Irel__Qr0) / var_Irel__cstar;
        const double var_Irel__sparkV = exp((-var_Irel__ay) * (var_Irel__V + 30.0)) / (1.0 + exp((-var_Irel__ay) * (var_Irel__V + 30.0)));
        const double var_Irel__spark_rate = (var_Irel__gryr / 1.0) * var_Irel__po * fabs(var_Irel__rxa) * var_Irel__sparkV;
        const double var_Irel__taua = 100.0;
        const double var_Irel__taur = 30.0;
        const double var_Irel__xirp = (((var_Irel__po * var_Irel__Qr * fabs(var_Irel__rxa) * var_Irel__gbarsr) / 1.0) * exp((-var_Irel__ax) * (var_Irel__V + 30.0))) / (1.0 + exp((-var_Irel__ax) * (var_Irel__V + 30.0)));
        const double var_Irel__xicap = var_Irel__po * var_Irel__gdyad * fabs(var_Irel__rxa);
        const double var_Irel__xiryr = var_Irel__xirp + var_Irel__xicap;
        const double var_Na__wca = var_cell__wca;
        const double var_Na__xina = var_INa__xina;
        const double var_Na__xiNaK = var_INaK__xiNaK;
        const double var_Na__xiNaCa = var_INaCa__xiNaCa;
        const double var_Ca__xiryr = var_Irel__xiryr;
        const double var_Ca__bcal = 24.0;
        const double var_Ca__xkcal = 7.0;
        const double var_Ca__srmax = 47.0;
        const double var_Ca__srkd = 0.6;
        const double var_Ca__bmem = 15.0;
        const double var_Ca__kmem = 0.3;
        const double var_Ca__bsar = 42.0;
        const double var_Ca__ksar = 13.0;
        const double var_Ca__bpxs = (var_Ca__bcal * var_Ca__xkcal) / ((var_Ca__xkcal + var_Ca__Ca_submem) * (var_Ca__xkcal + var_Ca__Ca_submem));
        const double var_Ca__spxs = (var_Ca__srmax * var_Ca__srkd) / ((var_Ca__srkd + var_Ca__Ca_submem) * (var_Ca__srkd + var_Ca__Ca_submem));
        const double var_Ca__mempxs = (var_Ca__bmem * var_Ca__kmem) / ((var_Ca__kmem + var_Ca__Ca_submem) * (var_Ca__kmem + var_Ca__Ca_submem));
        const double var_Ca__sarpxs = (var_Ca__bsar * var_Ca__ksar) / ((var_Ca__ksar + var_Ca__Ca_submem) * (var_Ca__ksar + var_Ca__Ca_submem));
        const double var_Ca__dcsib = 1.0 / (1.0 + var_Ca__bpxs + var_Ca__spxs + var_Ca__mempxs + var_Ca__sarpxs);
        const double var_Ca__bpxi = (var_Ca__bcal * var_Ca__xkcal) / ((var_Ca__xkcal + var_Ca__Ca_i) * (var_Ca__xkcal + var_Ca__Ca_i));
        const double var_Ca__spxi = (var_Ca__srmax * var_Ca__srkd) / ((var_Ca__srkd + var_Ca__Ca_i) * (var_Ca__srkd + var_Ca__Ca_i));
        const double var_Ca__mempxi = (var_Ca__bmem * var_Ca__kmem) / ((var_Ca__kmem + var_Ca__Ca_i) * (var_Ca__kmem + var_Ca__Ca_i));
        const double var_Ca__sarpxi = (var_Ca__bsar * var_Ca__ksar) / ((var_Ca__ksar + var_Ca__Ca_i) * (var_Ca__ksar + var_Ca__Ca_i));
        const double var_Ca__dciib = 1.0 / (1.0 + var_Ca__bpxi + var_Ca__spxi + var_Ca__mempxi + var_Ca__sarpxi);
        const double var_Ca__xkon = 0.0327;
        const double var_Ca__xkoff = 0.0196;
        const double var_Ca__btrop = 70.0;
        const double var_Ca__xbi = (var_Ca__xkon * var_Ca__Ca_i * (var_Ca__btrop - var_Ca__tropi)) - (var_Ca__xkoff * var_Ca__tropi);
        const double var_Ca__xbs = (var_Ca__xkon * var_Ca__Ca_submem * (var_Ca__btrop - var_Ca__trops)) - (var_Ca__xkoff * var_Ca__trops);
        const double var_Ca__taud = 4.0;
        const double var_Ca__taups = 0.5;
        const double var_Ca__jd = (var_Ca__Ca_submem - var_Ca__Ca_i) / var_Ca__taud;
        const double var_Ca__jNaCa = var_INaCa__jNaCa;
        const double var_Ca__jca = var_ICaL__jca;
        const double d_dt_ICaL__c1 = ((var_ICaL__alpha * var_ICaL__c2) + (var_ICaL__k2 * var_ICaL__xi1ca) + (var_ICaL__k2t * var_ICaL__xi1ba) + (var_ICaL__r2 * var_ICaL__po)) - ((var_ICaL__beta + var_ICaL__r1 + var_ICaL__k1t + var_ICaL__k1) * var_ICaL__c1);
        const double d_dt_ICaL__c2 = ((var_ICaL__beta * var_ICaL__c1) + (var_ICaL__k5 * var_ICaL__xi2ca) + (var_ICaL__k5t * var_ICaL__xi2ba)) - ((var_ICaL__k6 + var_ICaL__k6t + var_ICaL__alpha) * var_ICaL__c2);
        const double d_dt_ICaL__xi1ca = ((var_ICaL__k1 * var_ICaL__c1) + (var_ICaL__k4 * var_ICaL__xi2ca) + (var_ICaL__s1 * var_ICaL__po)) - ((var_ICaL__k3 + var_ICaL__k2 + var_ICaL__s2) * var_ICaL__xi1ca);
        const double d_dt_ICaL__xi1ba = ((var_ICaL__k1t * var_ICaL__c1) + (var_ICaL__k4t * var_ICaL__xi2ba) + (var_ICaL__s1t * var_ICaL__po)) - ((var_ICaL__k3t + var_ICaL__k2t + var_ICaL__s2t) * var_ICaL__xi1ba);
        const double d_dt_ICaL__xi2ca = ((var_ICaL__k3 * var_ICaL__xi1ca) + (var_ICaL__k6 * var_ICaL__c2)) - ((var_ICaL__k5 + var_ICaL__k4) * var_ICaL__xi2ca);
        const double d_dt_ICaL__xi2ba = ((var_ICaL__k3t * var_ICaL__xi1ba) + (var_ICaL__k6t * var_ICaL__c2)) - ((var_ICaL__k5t + var_ICaL__k4t) * var_ICaL__xi2ba);
        const double d_dt_Irel__Ca_JSR = (var_Irel__Ca_NSR - var_Irel__Ca_JSR) / var_Irel__taua;
        const double d_dt_Irel__xir = (var_Irel__spark_rate * var_Irel__Qr) - ((var_Irel__xir * (1.0 - ((var_Irel__taur * var_Irel__dCa_JSR) / var_Irel__Ca_NSR))) / var_Irel__taur);
        const double d_dt_Na__Na_i = (-(var_Na__xina + (3.0 * var_Na__xiNaK) + (3.0 * var_Na__xiNaCa))) / (var_Na__wca * 1000.0);
        const double d_dt_Ca__Ca_dyad = var_Ca__xiryr - ((var_Ca__Ca_dyad - var_Ca__Ca_submem) / var_Ca__taups);
        const double d_dt_Ca__Ca_submem = var_Ca__dcsib * ((50.0 * (((var_Ca__xir - var_Ca__jd) - var_Ca__jca) + var_Ca__jNaCa)) - var_Ca__xbs);
        const double d_dt_Ca__Ca_i = var_Ca__dciib * (((var_Ca__jd - var_Ca__jup) + var_Ca__jleak) - var_Ca__xbi);
        const double d_dt_Ca__Ca_NSR = var_Ca__dCa_JSR;
        const double d_dt_Ca__tropi = var_Ca__xbi;
        const double d_dt_Ca__trops = var_Ca__xbs;

        rResidual[6] = rCurrentGuess[6] - rY[4] - mDt*1.0*d_dt_ICaL__c1;
        rResidual[7] = rCurrentGuess[7] - rY[5] - mDt*1.0*d_dt_ICaL__c2;
        rResidual[9] = rCurrentGuess[9] - rY[6] - mDt*1.0*d_dt_ICaL__xi1ca;
        rResidual[8] = rCurrentGuess[8] - rY[7] - mDt*1.0*d_dt_ICaL__xi1ba;
        rResidual[11] = rCurrentGuess[11] - rY[8] - mDt*1.0*d_dt_ICaL__xi2ca;
        rResidual[10] = rCurrentGuess[10] - rY[9] - mDt*1.0*d_dt_ICaL__xi2ba;
        rResidual[12] = rCurrentGuess[12] - rY[17] - mDt*1.0*d_dt_Irel__Ca_JSR;
        rResidual[13] = rCurrentGuess[13] - rY[18] - mDt*1.0*d_dt_Irel__xir;
        rResidual[14] = rCurrentGuess[14] - rY[19] - mDt*1.0*d_dt_Na__Na_i;
        rResidual[1] = rCurrentGuess[1] - rY[20] - mDt*1.0*d_dt_Ca__Ca_dyad;
        rResidual[3] = rCurrentGuess[3] - rY[21] - mDt*1.0*d_dt_Ca__Ca_submem;
        rResidual[2] = rCurrentGuess[2] - rY[22] - mDt*1.0*d_dt_Ca__Ca_i;
        rResidual[0] = rCurrentGuess[0] - rY[23] - mDt*1.0*d_dt_Ca__Ca_NSR;
        rResidual[4] = rCurrentGuess[4] - rY[24] - mDt*1.0*d_dt_Ca__tropi;
        rResidual[5] = rCurrentGuess[5] - rY[25] - mDt*1.0*d_dt_Ca__trops;
    }

    void BackwardEulerMahajanModel2008::ComputeJacobian(double var_environment__time, const double rCurrentGuess[15], double rJacobian[15][15])
    {
        std::vector<double>& rY = rGetStateVariables();
        double var_cell__V = rY[0];
        // Units: mV; Initial value: -87.169816169406
        double var_INa__xm = rY[1];
        // Units: dimensionless; Initial value: 0.001075453357
        double var_INa__xh = rY[2];
        // Units: dimensionless; Initial value: 0.990691306716
        double var_INa__xj = rY[3];
        // Units: dimensionless; Initial value: 0.993888937283

        double var_Ca__Ca_NSR = rCurrentGuess[0];
        double var_Ca__Ca_dyad = rCurrentGuess[1];
        double var_Ca__Ca_i = rCurrentGuess[2];
        double var_Ca__Ca_submem = rCurrentGuess[3];
        double var_Ca__tropi = rCurrentGuess[4];
        double var_Ca__trops = rCurrentGuess[5];
        double var_ICaL__c1 = rCurrentGuess[6];
        double var_ICaL__c2 = rCurrentGuess[7];
        double var_ICaL__xi1ba = rCurrentGuess[8];
        double var_ICaL__xi1ca = rCurrentGuess[9];
        double var_ICaL__xi2ba = rCurrentGuess[10];
        double var_ICaL__xi2ca = rCurrentGuess[11];
        double var_Irel__Ca_JSR = rCurrentGuess[12];
        double var_Irel__xir = rCurrentGuess[13];
        double var_Na__Na_i = rCurrentGuess[14];

        const double dt = 1.0 * mDt;

        const double var_Environment__R = 8.314472;
        const double var_Environment__T = 308.0;
        const double var_Environment__F = 96.4853415;
        const double var_Environment__Ca_o = 1.8;
        const double var_Environment__FonRT = var_Environment__F / (var_Environment__R * var_Environment__T);
        const double var_Ca__csm = var_Ca__Ca_submem / 1000.0;
        const double var_ICaL__Ca_o = var_Environment__Ca_o;
        const double var_ICaL__V = var_cell__V;
        const double var_ICaL__pca = 0.00054;
        const double var_ICaL__FonRT = var_Environment__FonRT;
        const double var_ICaL__F = var_Environment__F;
        const double var_ICaL__za = var_ICaL__V * 2.0 * var_ICaL__FonRT;
        const double var_ICaL__csm = var_Ca__csm;
        const double var_Irel__cstar = 90.0;
        const double var_Irel__av = 11.3;
        const double var_Irel__bv = ((1.0 - var_Irel__av) * var_Irel__cstar) - 50.0;

        rJacobian[0][0] = 1.0 - (dt * ((((((-4.138e-05) * var_Ca__Ca_NSR) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i)) + (((4.138e-05 * pow(var_Ca__Ca_NSR, 3.0)) / pow(pow(var_Ca__Ca_NSR, 2.0) + 2500.0, 2.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i))) - ((0.00034484023 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0))));
        rJacobian[0][1] = 0.0;
        rJacobian[0][2] = (-dt) * ((((0.8 * var_Ca__Ca_i) / (pow(var_Ca__Ca_i, 2.0) + 0.25)) - ((0.8 * pow(var_Ca__Ca_i, 3.0)) / pow(pow(var_Ca__Ca_i, 2.0) + 0.25, 2.0))) + ((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)));
        rJacobian[0][3] = 0.0;
        rJacobian[0][4] = 0.0;
        rJacobian[0][5] = 0.0;
        rJacobian[0][6] = 0.0;
        rJacobian[0][7] = 0.0;
        rJacobian[0][8] = 0.0;
        rJacobian[0][9] = 0.0;
        rJacobian[0][10] = 0.0;
        rJacobian[0][11] = 0.0;
        rJacobian[0][12] = 0.0;
        rJacobian[0][13] = dt;
        rJacobian[0][14] = 0.0;
        rJacobian[1][0] = ((((((-298.2422222) * dt) * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728));
        rJacobian[1][1] = 1.0 + (2.0 * dt);
        rJacobian[1][2] = 0.0;
        rJacobian[1][3] = (-2.0) * dt;
        rJacobian[1][4] = 0.0;
        rJacobian[1][5] = 0.0;
        rJacobian[1][6] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][7] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][8] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][9] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][10] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][11] = (-dt) * (((((((-298.2422222) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728))) - (9000.0 * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))));
        rJacobian[1][12] = (((((((-298.2422222) * dt) * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? 1.0 : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? var_Irel__av : 0.0))) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.3576) * var_cell__V) - 10.728)) / (1.0 + exp(((-0.3576) * var_cell__V) - 10.728));
        rJacobian[1][13] = 0.0;
        rJacobian[1][14] = 0.0;
        rJacobian[2][0] = ((-dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_i, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_i, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_i, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_i, 2.0)))) * (((((4.138e-05 * var_Ca__Ca_NSR) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i)) - (((4.138e-05 * pow(var_Ca__Ca_NSR, 3.0)) / pow(pow(var_Ca__Ca_NSR, 2.0) + 2500.0, 2.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i))) + ((0.00034484023 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)));
        rJacobian[2][1] = 0.0;
        rJacobian[2][2] = (1.0 + (((dt / pow((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_i, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_i, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_i, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_i, 2.0)), 2.0)) * (((((((1.0 / 4.0) * var_Ca__Ca_submem) - ((1.0 / 4.0) * var_Ca__Ca_i)) - ((0.4 * pow(var_Ca__Ca_i, 2.0)) / (pow(var_Ca__Ca_i, 2.0) + 0.25))) + (((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i))) - ((0.0327 * var_Ca__Ca_i) * (70.0 - var_Ca__tropi))) + (0.0196 * var_Ca__tropi))) * (((((-336.0) / pow(7.0 + var_Ca__Ca_i, 3.0)) - (56.4 / pow(0.6 + var_Ca__Ca_i, 3.0))) - (9.0 / pow(0.3 + var_Ca__Ca_i, 3.0))) - (1092.0 / pow(13.0 + var_Ca__Ca_i, 3.0))))) - ((dt / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_i, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_i, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_i, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_i, 2.0)))) * (((((-2.539) - ((0.8 * var_Ca__Ca_i) / (pow(var_Ca__Ca_i, 2.0) + 0.25))) + ((0.8 * pow(var_Ca__Ca_i, 3.0)) / pow(pow(var_Ca__Ca_i, 2.0) + 0.25, 2.0))) - ((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0))) + (0.0327 * var_Ca__tropi)));
        rJacobian[2][3] = (((-1.0) / 4.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_i, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_i, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_i, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_i, 2.0)));
        rJacobian[2][4] = ((-dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_i, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_i, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_i, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_i, 2.0)))) * ((0.0327 * var_Ca__Ca_i) + 0.0196);
        rJacobian[2][5] = 0.0;
        rJacobian[2][6] = 0.0;
        rJacobian[2][7] = 0.0;
        rJacobian[2][8] = 0.0;
        rJacobian[2][9] = 0.0;
        rJacobian[2][10] = 0.0;
        rJacobian[2][11] = 0.0;
        rJacobian[2][12] = 0.0;
        rJacobian[2][13] = 0.0;
        rJacobian[2][14] = 0.0;
        rJacobian[3][0] = 0.0;
        rJacobian[3][1] = 0.0;
        rJacobian[3][2] = (((-25.0) / 2.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)));
        rJacobian[3][3] = (1.0 + (((dt / pow((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)), 2.0)) * (((((((50.0 * var_Irel__xir) - ((25.0 / 2.0) * var_Ca__Ca_submem)) + ((25.0 / 2.0) * var_Ca__Ca_i)) - ((9100.0 * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0))))) + ((((42.0 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022))) - ((0.0327 * var_Ca__Ca_submem) * (70.0 - var_Ca__trops))) + (0.0196 * var_Ca__trops))) * (((((-336.0) / pow(7.0 + var_Ca__Ca_submem, 3.0)) - (56.4 / pow(0.6 + var_Ca__Ca_submem, 3.0))) - (9.0 / pow(0.3 + var_Ca__Ca_submem, 3.0))) - (1092.0 / pow(13.0 + var_Ca__Ca_submem, 3.0))))) - ((dt / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((((-14.789) + (((((3.402 / pow(1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)), 2.0)) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022)) / pow(var_Ca__Ca_submem, 4.0))) - ((((105649.152 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * exp((-0.02449002923) * var_cell__V)) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022))) - ((((172864.0778 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / pow(((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022, 2.0))) + (0.0327 * var_Ca__trops)));
        rJacobian[3][4] = 0.0;
        rJacobian[3][5] = ((-dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * ((0.0327 * var_Ca__Ca_submem) + 0.0196);
        rJacobian[3][6] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][7] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][8] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][9] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][10] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][11] = (((-9100.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * (((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)));
        rJacobian[3][12] = 0.0;
        rJacobian[3][13] = ((-50.0) * dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)));
        rJacobian[3][14] = ((-dt) / ((((1.0 + (168.0 / pow(7.0 + var_Ca__Ca_submem, 2.0))) + (28.2 / pow(0.6 + var_Ca__Ca_submem, 2.0))) + (4.5 / pow(0.3 + var_Ca__Ca_submem, 2.0))) + (546.0 / pow(13.0 + var_Ca__Ca_submem, 2.0)))) * ((((((226.8 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * pow(var_Na__Na_i, 2.0)) * exp(0.01318693882 * var_cell__V)) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022)) - (((((1003.760877 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / pow(((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022, 2.0)) * pow(var_Na__Na_i, 2.0)));
        rJacobian[4][0] = 0.0;
        rJacobian[4][1] = 0.0;
        rJacobian[4][2] = (-dt) * (2.289 - (0.0327 * var_Ca__tropi));
        rJacobian[4][3] = 0.0;
        rJacobian[4][4] = 1.0 - (dt * (((-0.0327) * var_Ca__Ca_i) - 0.0196));
        rJacobian[4][5] = 0.0;
        rJacobian[4][6] = 0.0;
        rJacobian[4][7] = 0.0;
        rJacobian[4][8] = 0.0;
        rJacobian[4][9] = 0.0;
        rJacobian[4][10] = 0.0;
        rJacobian[4][11] = 0.0;
        rJacobian[4][12] = 0.0;
        rJacobian[4][13] = 0.0;
        rJacobian[4][14] = 0.0;
        rJacobian[5][0] = 0.0;
        rJacobian[5][1] = 0.0;
        rJacobian[5][2] = 0.0;
        rJacobian[5][3] = (-dt) * (2.289 - (0.0327 * var_Ca__trops));
        rJacobian[5][4] = 0.0;
        rJacobian[5][5] = 1.0 - (dt * (((-0.0327) * var_Ca__Ca_submem) - 0.0196));
        rJacobian[5][6] = 0.0;
        rJacobian[5][7] = 0.0;
        rJacobian[5][8] = 0.0;
        rJacobian[5][9] = 0.0;
        rJacobian[5][10] = 0.0;
        rJacobian[5][11] = 0.0;
        rJacobian[5][12] = 0.0;
        rJacobian[5][13] = 0.0;
        rJacobian[5][14] = 0.0;
        rJacobian[6][0] = 0.0;
        rJacobian[6][1] = (((1.957608 * dt) / pow(1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)), 2.0)) / pow(var_Ca__Ca_dyad, 4.0)) * var_ICaL__c1;
        rJacobian[6][2] = 0.0;
        rJacobian[6][3] = 0.0;
        rJacobian[6][4] = 0.0;
        rJacobian[6][5] = 0.0;
        rJacobian[6][6] = 1.0 - (dt * (((-4.30413) + (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V)))) - (0.024168 / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0))))));
        rJacobian[6][7] = (-dt) * ((1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))) - 3.0);
        rJacobian[6][8] = 2.99776 * dt;
        rJacobian[6][9] = 2.999896385 * dt;
        rJacobian[6][10] = 3.0 * dt;
        rJacobian[6][11] = 3.0 * dt;
        rJacobian[6][12] = 0.0;
        rJacobian[6][13] = 0.0;
        rJacobian[6][14] = 0.0;
        rJacobian[7][0] = 0.0;
        rJacobian[7][1] = (-dt) * (((((-(1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) / pow(((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1, 2.0)) * var_ICaL__xi2ca) * ((((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) - ((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)))) - ((((((81.0 / pow(1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)), 2.0)) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1)) / pow(var_Ca__Ca_dyad, 4.0)) - ((((1.0 / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / pow(((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1, 2.0)) * ((((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) - ((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0))))) * var_ICaL__c2));
        rJacobian[7][2] = 0.0;
        rJacobian[7][3] = 0.0;
        rJacobian[7][4] = 0.0;
        rJacobian[7][5] = 0.0;
        rJacobian[7][6] = (-dt) * (1.0 - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))));
        rJacobian[7][7] = 1.0 - (dt * ((((((-1.0) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1)) - ((1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / ((((-440.0) + (4954.0 * exp(0.0641025641 * var_cell__V))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + 450.0))) - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V)))));
        rJacobian[7][8] = 0.0;
        rJacobian[7][9] = 0.0;
        rJacobian[7][10] = ((-dt) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) / ((((-440.0) + (4954.0 * exp(0.0641025641 * var_cell__V))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + 450.0);
        rJacobian[7][11] = ((-dt) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1);
        rJacobian[7][12] = 0.0;
        rJacobian[7][13] = 0.0;
        rJacobian[7][14] = 0.0;
        rJacobian[8][0] = 0.0;
        rJacobian[8][1] = 0.0;
        rJacobian[8][2] = 0.0;
        rJacobian[8][3] = 0.0;
        rJacobian[8][4] = 0.0;
        rJacobian[8][5] = 0.0;
        rJacobian[8][6] = (-0.00218) * dt;
        rJacobian[8][7] = 0.00195 * dt;
        rJacobian[8][8] = 1.0 - (dt * ((-0.337629096) + ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0))))));
        rJacobian[8][9] = 0.00195 * dt;
        rJacobian[8][10] = (-dt) * ((((((1.84375 * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))))) / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))) / (1.0 - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))))) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) * (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) - 0.00195);
        rJacobian[8][11] = 0.00195 * dt;
        rJacobian[8][12] = 0.0;
        rJacobian[8][13] = 0.0;
        rJacobian[8][14] = 0.0;
        rJacobian[9][0] = 0.0;
        rJacobian[9][1] = (-dt) * ((((1.957608 / pow(1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)), 2.0)) / pow(var_Ca__Ca_dyad, 4.0)) * var_ICaL__c1) + (((1.4797728 / pow(1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)), 2.0)) * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) / pow(var_Ca__Ca_dyad, 4.0)));
        rJacobian[9][2] = 0.0;
        rJacobian[9][3] = 0.0;
        rJacobian[9][4] = 0.0;
        rJacobian[9][5] = 0.0;
        rJacobian[9][6] = ((-0.0058992) * dt) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)));
        rJacobian[9][7] = (0.0182688 * dt) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)));
        rJacobian[9][8] = (0.0182688 * dt) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)));
        rJacobian[9][9] = 1.0 - (dt * ((((-0.0182688) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))) - 0.3334447806) + ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0))))));
        rJacobian[9][10] = (0.0182688 * dt) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)));
        rJacobian[9][11] = (-dt) * ((((((233.2480818 * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))))) / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))) / (1.0 - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))))) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) * (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) - (0.0182688 / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))));
        rJacobian[9][12] = 0.0;
        rJacobian[9][13] = 0.0;
        rJacobian[9][14] = 0.0;
        rJacobian[10][0] = 0.0;
        rJacobian[10][1] = 0.0;
        rJacobian[10][2] = 0.0;
        rJacobian[10][3] = 0.0;
        rJacobian[10][4] = 0.0;
        rJacobian[10][5] = 0.0;
        rJacobian[10][6] = 0.0;
        rJacobian[10][7] = ((-dt) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / ((((-440.0) + (4954.0 * exp(0.0641025641 * var_cell__V))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + 450.0);
        rJacobian[10][8] = (-dt) * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))));
        rJacobian[10][9] = 0.0;
        rJacobian[10][10] = 1.0 - (dt * (((-(1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) / ((((-440.0) + (4954.0 * exp(0.0641025641 * var_cell__V))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + 450.0)) - (((((1.84375 * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))))) / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))) / (1.0 - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))))) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) * (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904)))));
        rJacobian[10][11] = 0.0;
        rJacobian[10][12] = 0.0;
        rJacobian[10][13] = 0.0;
        rJacobian[10][14] = 0.0;
        rJacobian[11][0] = 0.0;
        rJacobian[11][1] = (-dt) * (((((((81.0 / pow(1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)), 2.0)) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1)) * var_ICaL__c2) / pow(var_Ca__Ca_dyad, 4.0)) - (((((1.0 / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / pow(((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1, 2.0)) * var_ICaL__c2) * ((((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) - ((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0))))) + ((((1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904)))) / pow(((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1, 2.0)) * var_ICaL__xi2ca) * ((((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) - ((0.2263746095 / pow(1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0)), 2.0)) * pow(var_Ca__Ca_dyad, 3.0)))));
        rJacobian[11][2] = 0.0;
        rJacobian[11][3] = 0.0;
        rJacobian[11][4] = 0.0;
        rJacobian[11][5] = 0.0;
        rJacobian[11][6] = 0.0;
        rJacobian[11][7] = (((-dt) / (1.0 + (27.0 / pow(var_Ca__Ca_dyad, 3.0)))) / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1);
        rJacobian[11][8] = 0.0;
        rJacobian[11][9] = (-dt) * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))));
        rJacobian[11][10] = 0.0;
        rJacobian[11][11] = 1.0 - (dt * (((-(1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) / (((((9.9 + (4954.0 * exp(0.0641025641 * var_cell__V))) - (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) * (1.0 - (1.0 / (1.0 + exp((((-1.0) / 4.0) * var_cell__V) - 10.0))))) + (78.0329 / (1.0 + (0.0007252537376 * pow(var_Ca__Ca_dyad, 4.0))))) + 0.1)) - (((((233.2480818 * ((1.0 / 3.0) - ((1.0 / 3.0) / (1.0 + exp((((-1.0) / 3.0) * var_cell__V) - (40.0 / 3.0)))))) / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))) / (1.0 - (1.0 / (1.0 + exp(((-1.0) / 8.0) * var_cell__V))))) * (1.0 - (1.0 / (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904))))) * (1.0 + exp(((-0.08833922261) * var_cell__V) - 3.533568904)))));
        rJacobian[11][12] = 0.0;
        rJacobian[11][13] = 0.0;
        rJacobian[11][14] = 0.0;
        rJacobian[12][0] = ((-1.0) / 100.0) * dt;
        rJacobian[12][1] = 0.0;
        rJacobian[12][2] = 0.0;
        rJacobian[12][3] = 0.0;
        rJacobian[12][4] = 0.0;
        rJacobian[12][5] = 0.0;
        rJacobian[12][6] = 0.0;
        rJacobian[12][7] = 0.0;
        rJacobian[12][8] = 0.0;
        rJacobian[12][9] = 0.0;
        rJacobian[12][10] = 0.0;
        rJacobian[12][11] = 0.0;
        rJacobian[12][12] = 1.0 + ((1.0 / 100.0) * dt);
        rJacobian[12][13] = 0.0;
        rJacobian[12][14] = 0.0;
        rJacobian[13][0] = (-dt) * ((((((0.02867544444 * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0))) - (((1.0 / 30.0) * var_Irel__xir) * ((((-30.0) * ((((((-4.138e-05) * var_Ca__Ca_NSR) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i)) + (((4.138e-05 * pow(var_Ca__Ca_NSR, 3.0)) / pow(pow(var_Ca__Ca_NSR, 2.0) + 2500.0, 2.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i))) - ((0.00034484023 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)))) / var_Ca__Ca_NSR) + ((30.0 * (((-var_Irel__xir) + ((0.4 * pow(var_Ca__Ca_i, 2.0)) / (pow(var_Ca__Ca_i, 2.0) + 0.25))) - (((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i)))) / pow(var_Ca__Ca_NSR, 2.0)))));
        rJacobian[13][1] = 0.0;
        rJacobian[13][2] = (((-dt) * var_Irel__xir) * ((((0.8 * var_Ca__Ca_i) / (pow(var_Ca__Ca_i, 2.0) + 0.25)) - ((0.8 * pow(var_Ca__Ca_i, 3.0)) / pow(pow(var_Ca__Ca_i, 2.0) + 0.25, 2.0))) + ((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)))) / var_Ca__Ca_NSR;
        rJacobian[13][3] = 0.0;
        rJacobian[13][4] = 0.0;
        rJacobian[13][5] = 0.0;
        rJacobian[13][6] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][7] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][8] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][9] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][10] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][11] = (((((0.02867544444 * dt) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? (var_Irel__Ca_JSR - 50.0) : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? ((var_Irel__av * var_Irel__Ca_JSR) + var_Irel__bv) : 0.0));
        rJacobian[13][12] = (((((((-0.02867544444) * dt) * ((((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2)) * fabs(((0.0753539361 * fabs(var_cell__V)) < 0.001) ? (((2.0 * var_ICaL__pca) * var_ICaL__F) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) : ((((((4.0 * var_ICaL__pca) * var_ICaL__V) * var_ICaL__F) * var_ICaL__FonRT) * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0)))) * exp(((-0.05) * var_cell__V) - 1.5)) / (1.0 + exp(((-0.05) * var_cell__V) - 1.5))) * var_Ca__Ca_NSR) * (((50.0 < var_Irel__Ca_JSR) && (var_Irel__Ca_JSR < 90.0)) ? 1.0 : ((var_Irel__cstar <= var_Irel__Ca_JSR) ? var_Irel__av : 0.0));
        rJacobian[13][13] = 1.0 - (dt * ((((-1.0) / 30.0) + ((((-var_Irel__xir) + ((0.4 * pow(var_Ca__Ca_i, 2.0)) / (pow(var_Ca__Ca_i, 2.0) + 0.25))) - (((2.069e-05 * pow(var_Ca__Ca_NSR, 2.0)) / (pow(var_Ca__Ca_NSR, 2.0) + 2500.0)) * ((16.667 * var_Ca__Ca_NSR) - var_Ca__Ca_i))) / var_Ca__Ca_NSR)) - (var_Irel__xir / var_Ca__Ca_NSR)));
        rJacobian[13][14] = 0.0;
        rJacobian[14][0] = 0.0;
        rJacobian[14][1] = 0.0;
        rJacobian[14][2] = 0.0;
        rJacobian[14][3] = (-dt) * ((((((((-0.00020412) / pow(1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)), 2.0)) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022)) / pow(var_Ca__Ca_submem, 4.0)) + ((((6.33894912 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * exp((-0.02449002923) * var_cell__V)) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022))) + ((((10.37184466 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / pow(((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022, 2.0)));
        rJacobian[14][4] = 0.0;
        rJacobian[14][5] = 0.0;
        rJacobian[14][6] = 0.0;
        rJacobian[14][7] = 0.0;
        rJacobian[14][8] = 0.0;
        rJacobian[14][9] = 0.0;
        rJacobian[14][10] = 0.0;
        rJacobian[14][11] = 0.0;
        rJacobian[14][12] = 0.0;
        rJacobian[14][13] = 0.0;
        rJacobian[14][14] = 1.0 - (dt * (((((((((-0.03981212072) * var_INa__xh) * var_INa__xj) * pow(var_INa__xm, 3.0)) / var_Na__Na_i) - ((0.0004402173911 / ((1.0 + (0.1245 * exp((-0.003767696805) * var_cell__V))) + (0.03412424596 * exp((-0.03767696805) * var_cell__V)))) / (var_Na__Na_i + 12.0))) + (((0.0004402173911 / ((1.0 + (0.1245 * exp((-0.003767696805) * var_cell__V))) + (0.03412424596 * exp((-0.03767696805) * var_cell__V)))) * var_Na__Na_i) / pow(var_Na__Na_i + 12.0, 2.0))) - (((((0.013608 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * pow(var_Na__Na_i, 2.0)) * exp(0.01318693882 * var_cell__V)) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / (((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022))) + (((((0.06022565262 / (1.0 + (0.027 / pow(var_Ca__Ca_submem, 3.0)))) * (((1.8 * pow(var_Na__Na_i, 3.0)) * exp(0.01318693882 * var_cell__V)) - (((314432.0 / 125.0) * var_Ca__Ca_submem) * exp((-0.02449002923) * var_cell__V)))) / (1.0 + (0.2 * exp((-0.02449002923) * var_cell__V)))) / pow(((7.966356167 * pow(var_Na__Na_i, 3.0)) + (4115.811375 * var_Ca__Ca_submem)) + 12405.2022, 2.0)) * pow(var_Na__Na_i, 2.0))));
    }

    void BackwardEulerMahajanModel2008::UpdateTransmembranePotential(double var_Environment__time)
    {
        // Time units: ms
        var_Environment__time *= 1.0;
        std::vector<double>& rY = rGetStateVariables();
        double var_cell__V = rY[0];
        // Units: mV; Initial value: -87.169816169406
        double var_INa__xm = rY[1];
        // Units: dimensionless; Initial value: 0.001075453357
        double var_INa__xh = rY[2];
        // Units: dimensionless; Initial value: 0.990691306716
        double var_INa__xj = rY[3];
        // Units: dimensionless; Initial value: 0.993888937283
        double var_ICaL__c1 = rY[4];
        // Units: dimensionless; Initial value: 0.000018211252
        double var_ICaL__c2 = rY[5];
        // Units: dimensionless; Initial value: 0.979322592773
        double var_ICaL__xi1ca = rY[6];
        // Units: dimensionless; Initial value: 0.001208153482
        double var_ICaL__xi1ba = rY[7];
        // Units: dimensionless; Initial value: 0.000033616596
        double var_ICaL__xi2ca = rY[8];
        // Units: dimensionless; Initial value: 0.004173008466
        double var_ICaL__xi2ba = rY[9];
        // Units: dimensionless; Initial value: 0.015242594688
        double var_IKr__xr = rY[10];
        // Units: dimensionless; Initial value: 0.007074239331
        double var_IKs__xs1 = rY[11];
        // Units: dimensionless; Initial value: 0.048267587131
        double var_IKs__xs2 = rY[12];
        // Units: dimensionless; Initial value: 0.105468807033
        double var_Ito__xtos = rY[13];
        // Units: dimensionless; Initial value: 0.00364776906
        double var_Ito__ytos = rY[14];
        // Units: dimensionless; Initial value: 0.174403618112
        double var_Ito__xtof = rY[15];
        // Units: dimensionless; Initial value: 0.003643592594
        double var_Ito__ytof = rY[16];
        // Units: dimensionless; Initial value: 0.993331326442
        double var_Na__Na_i = rY[19];
        // Units: mM; Initial value: 11.441712311614
        double var_Ca__Ca_submem = rY[21];
        // Units: uM; Initial value: 0.226941113355
        double var_Ca__Ca_i = rY[22];
        // Units: uM; Initial value: 0.256752008084

        const double var_Environment__R = 8.314472;
        const double var_Environment__T = 308.0;
        const double var_Environment__F = 96.4853415;
        const double var_Environment__K_o = 5.4;
        const double var_Environment__Ca_o = 1.8;
        const double var_Environment__Na_o = 136.0;
        const double var_Environment__FonRT = var_Environment__F / (var_Environment__R * var_Environment__T);
//        const double var_cell__time = var_Environment__time;
        const double var_reversal_potentials__K_i = 140.0;
        const double var_reversal_potentials__FonRT = var_Environment__FonRT;
        const double var_reversal_potentials__K_o = var_Environment__K_o;
        const double var_reversal_potentials__ek = (1.0 / var_reversal_potentials__FonRT) * log(var_reversal_potentials__K_o / var_reversal_potentials__K_i);
        const double var_IK1__ek = var_reversal_potentials__ek;
        const double var_IK1__V = var_cell__V;
        const double var_IK1__aki = 1.02 / (1.0 + exp(0.2385 * ((var_IK1__V - var_IK1__ek) - 59.215)));
        const double var_IK1__bki = ((0.49124 * exp(0.08032 * ((var_IK1__V - var_IK1__ek) + 5.476))) + (1.0 * exp(0.06175 * ((var_IK1__V - var_IK1__ek) - 594.31)))) / (1.0 + exp((-0.5143) * ((var_IK1__V - var_IK1__ek) + 4.753)));
        const double var_IK1__xkin = var_IK1__aki / (var_IK1__aki + var_IK1__bki);
        const double var_IK1__K_o = var_Environment__K_o;
        const double var_IK1__gkix = 0.3;
        const double var_IK1__xik1 = var_IK1__gkix * sqrt(var_IK1__K_o / 5.4) * var_IK1__xkin * (var_IK1__V - var_IK1__ek);
        const double var_cell__xik1 = var_IK1__xik1;
        const double var_Ito__ek = var_reversal_potentials__ek;
        const double var_Ito__gtos = 0.04;
        const double var_Ito__V = var_cell__V;
        const double var_Ito__rt2 = (var_Ito__V + 33.5) / 10.0;
        const double var_Ito__rs_inf = 1.0 / (1.0 + exp(var_Ito__rt2));
        const double var_Ito__xitos = var_Ito__gtos * var_Ito__xtos * (var_Ito__ytos + (0.5 * var_Ito__rs_inf)) * (var_Ito__V - var_Ito__ek);
        const double var_Ito__gtof = 0.11;
        const double var_Ito__xitof = var_Ito__gtof * var_Ito__xtof * var_Ito__ytof * (var_Ito__V - var_Ito__ek);
        const double var_Ito__xito = var_Ito__xitos + var_Ito__xitof;
        const double var_cell__xito = var_Ito__xito;
        const double var_INaK__xkmko = 1.5;
        const double var_INaK__V = var_cell__V;
        const double var_INaK__FonRT = var_Environment__FonRT;
        const double var_INaK__Na_o = var_Environment__Na_o;
        const double var_INaK__sigma = (exp(var_INaK__Na_o / 67.3) - 1.0) / 7.0;
        const double var_INaK__fNaK = 1.0 / (1.0 + (0.1245 * exp((-0.1) * var_INaK__V * var_INaK__FonRT)) + (0.0365 * var_INaK__sigma * exp((-var_INaK__V) * var_INaK__FonRT)));
        const double var_INaK__Na_i = var_Na__Na_i;
        const double var_INaK__K_o = var_Environment__K_o;
        const double var_INaK__xkmnai = 12.0;
        const double var_INaK__gNaK = 1.5;
        const double var_INaK__xiNaK = (((var_INaK__gNaK * var_INaK__fNaK * var_INaK__Na_i) / (var_INaK__Na_i + var_INaK__xkmnai)) * var_INaK__K_o) / (var_INaK__K_o + var_INaK__xkmko);
        const double var_cell__xiNaK = var_INaK__xiNaK;
        const double var_cell__wca = 8.0;
        const double var_INaCa__gNaCa = 0.84;
        const double var_INaCa__V = var_cell__V;
        const double var_INaCa__FonRT = var_Environment__FonRT;
        const double var_INaCa__zw4 = 1.0 + (0.2 * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__xkdna = 0.3;
        const double var_INaCa__Ca_submem = var_Ca__Ca_submem;
        const double var_INaCa__aloss = 1.0 / (1.0 + pow(var_INaCa__xkdna / var_INaCa__Ca_submem, 3.0));
        const double var_INaCa__Na_o = var_Environment__Na_o;
        const double var_INaCa__Na_i = var_Na__Na_i;
        const double var_INaCa__Ca_o = var_Environment__Ca_o;
        const double var_Ca__csm = var_Ca__Ca_submem / 1000.0;
        const double var_INaCa__csm = var_Ca__csm;
        const double var_INaCa__yz4 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o) + (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcao = 1.3;
        const double var_INaCa__xmnao = 87.5;
        const double var_INaCa__yz1 = (var_INaCa__xmcao * pow(var_INaCa__Na_i, 3.0)) + (pow(var_INaCa__xmnao, 3.0) * var_INaCa__csm);
        const double var_INaCa__xmcai = 0.0036;
        const double var_INaCa__xmnai = 12.3;
        const double var_INaCa__yz3 = var_INaCa__xmcai * pow(var_INaCa__Na_o, 3.0) * (1.0 + pow(var_INaCa__Na_i / var_INaCa__xmnai, 3.0));
        const double var_INaCa__yz2 = pow(var_INaCa__xmnai, 3.0) * var_INaCa__Ca_o * (1.0 + (var_INaCa__csm / var_INaCa__xmcai));
        const double var_INaCa__zw8 = var_INaCa__yz1 + var_INaCa__yz2 + var_INaCa__yz3 + var_INaCa__yz4;
        const double var_INaCa__zw3 = (pow(var_INaCa__Na_i, 3.0) * var_INaCa__Ca_o * exp(var_INaCa__V * 0.35 * var_INaCa__FonRT)) - (pow(var_INaCa__Na_o, 3.0) * var_INaCa__csm * exp(var_INaCa__V * (0.35 - 1.0) * var_INaCa__FonRT));
        const double var_INaCa__jNaCa = (var_INaCa__gNaCa * var_INaCa__aloss * var_INaCa__zw3) / (var_INaCa__zw4 * var_INaCa__zw8);
        const double var_INaCa__wca = var_cell__wca;
        const double var_INaCa__xiNaCa = var_INaCa__wca * var_INaCa__jNaCa;
        const double var_cell__xiNaCa = var_INaCa__xiNaCa;
        const double var_ICaL__Ca_o = var_Environment__Ca_o;
        const double var_ICaL__V = var_cell__V;
        const double var_ICaL__pca = 0.00054;
        const double var_ICaL__FonRT = var_Environment__FonRT;
        const double var_ICaL__F = var_Environment__F;
        const double var_ICaL__za = var_ICaL__V * 2.0 * var_ICaL__FonRT;
        const double var_ICaL__csm = var_Ca__csm;
        const double var_ICaL__rxa = (fabs(var_ICaL__za) < 0.001) ? ((4.0 * var_ICaL__pca * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (2.0 * var_ICaL__FonRT)) : ((4.0 * var_ICaL__pca * var_ICaL__V * var_ICaL__F * var_ICaL__FonRT * ((var_ICaL__csm * exp(var_ICaL__za)) - (0.341 * var_ICaL__Ca_o))) / (exp(var_ICaL__za) - 1.0));
        const double var_ICaL__gca = 182.0;
        const double var_ICaL__po = (((((1.0 - var_ICaL__xi1ca) - var_ICaL__xi2ca) - var_ICaL__xi1ba) - var_ICaL__xi2ba) - var_ICaL__c1) - var_ICaL__c2;
        const double var_ICaL__jca = var_ICaL__gca * var_ICaL__po * var_ICaL__rxa;
        const double var_ICaL__wca = var_cell__wca;
        const double var_ICaL__xica = 2.0 * var_ICaL__wca * var_ICaL__jca;
        const double var_cell__xica = var_ICaL__xica;
        const double var_IKr__gkr = 0.0125;
        const double var_IKr__K_o = var_Environment__K_o;
        const double var_IKr__V = var_cell__V;
        const double var_IKr__rg = 1.0 / (1.0 + exp((var_IKr__V + 33.0) / 22.4));
        const double var_IKr__ek = var_reversal_potentials__ek;
        const double var_IKr__xikr = var_IKr__gkr * sqrt(var_IKr__K_o / 5.4) * var_IKr__xr * var_IKr__rg * (var_IKr__V - var_IKr__ek);
        const double var_cell__xikr = var_IKr__xikr;
        const double var_IKs__Ca_i = var_Ca__Ca_i;
        const double var_IKs__gksx = 1.0 + (0.8 / (1.0 + pow(0.5 / var_IKs__Ca_i, 3.0)));
        const double var_IKs__gks = 0.1386;
        const double var_IKs__V = var_cell__V;
        const double var_reversal_potentials__prNaK = 0.01833;
        const double var_reversal_potentials__Na_o = var_Environment__Na_o;
        const double var_reversal_potentials__Na_i = var_Na__Na_i;
        const double var_reversal_potentials__eks = (1.0 / var_reversal_potentials__FonRT) * log((var_reversal_potentials__K_o + (var_reversal_potentials__prNaK * var_reversal_potentials__Na_o)) / (var_reversal_potentials__K_i + (var_reversal_potentials__prNaK * var_reversal_potentials__Na_i)));
        const double var_IKs__eks = var_reversal_potentials__eks;
        const double var_IKs__xiks = var_IKs__gks * var_IKs__gksx * var_IKs__xs1 * var_IKs__xs2 * (var_IKs__V - var_IKs__eks);
        const double var_cell__xiks = var_IKs__xiks;
//        const double var_cell__stim_offset = 0.0;
//        const double var_cell__stim_period = 400.0;
//        const double var_cell__past = floor(var_cell__time / var_cell__stim_period) * var_cell__stim_period;
//        const double var_cell__stim_amplitude =  -15.0;
//        const double var_cell__stim_duration = 3.0;
        double var_cell__i_Stim = GetStimulus((1.0/1)*var_Environment__time);
        const double var_INa__V = var_cell__V;
        const double var_INa__gna = 12.0;
        const double var_reversal_potentials__ena = (1.0 / var_reversal_potentials__FonRT) * log(var_reversal_potentials__Na_o / var_reversal_potentials__Na_i);
        const double var_INa__ena = var_reversal_potentials__ena;
        const double var_INa__xina = var_INa__gna * var_INa__xh * var_INa__xj * var_INa__xm * var_INa__xm * var_INa__xm * (var_INa__V - var_INa__ena);
        const double var_cell__xina = var_INa__xina;
        const double var_cell__Itotal = -(var_cell__xina + var_cell__xik1 + var_cell__xikr + var_cell__xiks + var_cell__xito + var_cell__xiNaCa + var_cell__xica + var_cell__xiNaK + var_cell__i_Stim);

        double d_dt_cell__V;

        d_dt_cell__V = var_cell__Itotal;

        rY[0] += mDt * 1.0*d_dt_cell__V;
    }

    void BackwardEulerMahajanModel2008::ComputeOneStepExceptVoltage(double var_Environment__time)
    {
        // Time units: ms
        var_Environment__time *= 1.0;
        std::vector<double>& rY = rGetStateVariables();
        double var_cell__V = rY[0];
        // Units: mV; Initial value: -87.169816169406

        const double var_Ito__V = var_cell__V;
        const double var_Ito__rt2 = (var_Ito__V + 33.5) / 10.0;
        const double var_IKr__V = var_cell__V;
        const double var_IKs__V = var_cell__V;
        const double var_INa__V = var_cell__V;
        const double var_INa__am = (fabs(var_INa__V + 47.13) > 0.001) ? ((0.32 * 1.0 * (var_INa__V + 47.13)) / (1.0 - exp((-0.1) * (var_INa__V + 47.13)))) : 3.2;
        const double var_INa__bm = 0.08 * exp((-var_INa__V) / 11.0);
        const double var_INa__ah = (var_INa__V < (-40.0)) ? (0.135 * exp((80.0 + var_INa__V) / (-6.8))) : 0.0;
        const double var_INa__bh = (var_INa__V < (-40.0)) ? ((3.56 * exp(0.079 * var_INa__V)) + (310000.0 * exp(0.35 * var_INa__V))) : (1.0 / (0.13 * (1.0 + exp((var_INa__V + 10.66) / (-11.1)))));
        const double var_INa__aj = (var_INa__V < (-40.0)) ? (((((-127140.0) * exp(0.2444 * var_INa__V)) - (3.474e-05 * exp((-0.04391) * var_INa__V))) * 1.0 * (var_INa__V + 37.78)) / (1.0 + exp(0.311 * (var_INa__V + 79.23)))) : 0.0;
        const double var_INa__bj = (var_INa__V < (-40.0)) ? ((0.1212 * exp((-0.01052) * var_INa__V)) / (1.0 + exp((-0.1378) * (var_INa__V + 40.14)))) : ((0.3 * exp((-2.535e-07) * var_INa__V)) / (1.0 + exp((-0.1) * (var_INa__V + 32.0))));
        const double var_IKr__xkrv1 = (fabs(var_IKr__V + 7.0) > 0.001) ? ((0.00138 * 1.0 * (var_IKr__V + 7.0)) / (1.0 - exp((-0.123) * (var_IKr__V + 7.0)))) : (0.00138 / 0.123);
        const double var_IKr__xkrv2 = (fabs(var_IKr__V + 10.0) > 0.001) ? ((0.00061 * 1.0 * (var_IKr__V + 10.0)) / (exp(0.145 * (var_IKr__V + 10.0)) - 1.0)) : (0.00061 / 0.145);
        const double var_IKr__taukr = 1.0 / (var_IKr__xkrv1 + var_IKr__xkrv2);
        const double var_IKr__xkrinf = 1.0 / (1.0 + exp((-(var_IKr__V + 50.0)) / 7.5));
        const double var_IKs__xs1ss = 1.0 / (1.0 + exp((-(var_IKs__V - 1.5)) / 16.7));
        const double var_IKs__xs2ss = var_IKs__xs1ss;
        const double var_IKs__tauxs1 = (fabs(var_IKs__V + 30.0) < (0.001 / 0.0687)) ? (1.0 / ((7.19e-05 / 0.148) + (0.000131 / 0.0687))) : (1.0 / (((7.19e-05 * (var_IKs__V + 30.0)) / (1.0 - exp((-0.148) * (var_IKs__V + 30.0)))) + ((0.000131 * (var_IKs__V + 30.0)) / (exp(0.0687 * (var_IKs__V + 30.0)) - 1.0))));
        const double var_IKs__tauxs2 = 4.0 * var_IKs__tauxs1;
        const double var_Ito__rt1 = (-(var_Ito__V + 3.0)) / 15.0;
        const double var_Ito__rt3 = (var_Ito__V + 60.0) / 10.0;
        const double var_Ito__rt4 = (((-var_Ito__V) / 30.0) * var_Ito__V) / 30.0;
        const double var_Ito__rt5 = (var_Ito__V + 33.5) / 10.0;
        const double var_Ito__xtos_inf = 1.0 / (1.0 + exp(var_Ito__rt1));
        const double var_Ito__ytos_inf = 1.0 / (1.0 + exp(var_Ito__rt2));
        const double var_Ito__xtof_inf = var_Ito__xtos_inf;
        const double var_Ito__ytof_inf = var_Ito__ytos_inf;
        const double var_Ito__txs = (9.0 / (1.0 + exp(-var_Ito__rt1))) + 0.5;
        const double var_Ito__tys = (3000.0 / (1.0 + exp(var_Ito__rt3))) + 30.0;
        const double var_Ito__txf = (3.5 * exp(var_Ito__rt4)) + 1.5;
        const double var_Ito__tyf = (20.0 / (1.0 + exp(var_Ito__rt5))) + 20.0;

        const double _g_0 = var_IKr__xkrinf / var_IKr__taukr;
        const double _h_0 = (-1.0) / var_IKr__taukr;
        const double _g_1 = var_IKs__xs1ss / var_IKs__tauxs1;
        const double _h_1 = (-1.0) / var_IKs__tauxs1;
        const double _g_2 = var_IKs__xs2ss / var_IKs__tauxs2;
        const double _h_2 = (-1.0) / var_IKs__tauxs2;
        const double _g_3 = var_INa__ah * 1.0;
        const double _h_3 = (var_INa__ah * (-1.0)) - (var_INa__bh * 1.0);
        const double _g_4 = var_INa__aj * 1.0;
        const double _h_4 = (var_INa__aj * (-1.0)) - (var_INa__bj * 1.0);
        const double _g_5 = var_INa__am * 1.0;
        const double _h_5 = (var_INa__am * (-1.0)) - (var_INa__bm * 1.0);
        const double _g_6 = var_Ito__xtof_inf / var_Ito__txf;
        const double _h_6 = (-1.0) / var_Ito__txf;
        const double _g_7 = var_Ito__xtos_inf / var_Ito__txs;
        const double _h_7 = (-1.0) / var_Ito__txs;
        const double _g_8 = var_Ito__ytof_inf / var_Ito__tyf;
        const double _h_8 = (-1.0) / var_Ito__tyf;
        const double _g_9 = var_Ito__ytos_inf / var_Ito__tys;
        const double _h_9 = (-1.0) / var_Ito__tys;

        const double dt = mDt*1.0;
        rY[1] = (rY[1] + _g_5*dt) / (1 - _h_5*dt);
        rY[2] = (rY[2] + _g_3*dt) / (1 - _h_3*dt);
        rY[3] = (rY[3] + _g_4*dt) / (1 - _h_4*dt);
        rY[10] = (rY[10] + _g_0*dt) / (1 - _h_0*dt);
        rY[11] = (rY[11] + _g_1*dt) / (1 - _h_1*dt);
        rY[12] = (rY[12] + _g_2*dt) / (1 - _h_2*dt);
        rY[13] = (rY[13] + _g_7*dt) / (1 - _h_7*dt);
        rY[14] = (rY[14] + _g_9*dt) / (1 - _h_9*dt);
        rY[15] = (rY[15] + _g_6*dt) / (1 - _h_6*dt);
        rY[16] = (rY[16] + _g_8*dt) / (1 - _h_8*dt);

        double _guess[15] = {rY[23],rY[20],rY[22],rY[21],rY[24],rY[25],rY[4],rY[5],rY[7],rY[6],rY[9],rY[8],rY[17],rY[18],rY[19]};
        CardiacNewtonSolver<15> *_solver = CardiacNewtonSolver<15>::Instance();
        _solver->Solve(*this, var_Environment__time, _guess);
        rY[23] = _guess[0];
        rY[20] = _guess[1];
        rY[22] = _guess[2];
        rY[21] = _guess[3];
        rY[24] = _guess[4];
        rY[25] = _guess[5];
        rY[4] = _guess[6];
        rY[5] = _guess[7];
        rY[7] = _guess[8];
        rY[6] = _guess[9];
        rY[9] = _guess[10];
        rY[8] = _guess[11];
        rY[17] = _guess[12];
        rY[18] = _guess[13];
        rY[19] = _guess[14];
    }

template<>
void OdeSystemInformation<BackwardEulerMahajanModel2008>::Initialise(void)
{
    // Time units: ms
    //
    this->mVariableNames.push_back("V");
    this->mVariableUnits.push_back("mV");
    this->mInitialConditions.push_back(-87.169816169406);

    this->mVariableNames.push_back("xm");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.001075453357);

    this->mVariableNames.push_back("xh");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.990691306716);

    this->mVariableNames.push_back("xj");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.993888937283);

    this->mVariableNames.push_back("c1");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.000018211252);

    this->mVariableNames.push_back("c2");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.979322592773);

    this->mVariableNames.push_back("xi1ca");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.001208153482);

    this->mVariableNames.push_back("xi1ba");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.000033616596);

    this->mVariableNames.push_back("xi2ca");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.004173008466);

    this->mVariableNames.push_back("xi2ba");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.015242594688);

    this->mVariableNames.push_back("xr");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.007074239331);

    this->mVariableNames.push_back("xs1");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.048267587131);

    this->mVariableNames.push_back("xs2");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.105468807033);

    this->mVariableNames.push_back("xtos");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.00364776906);

    this->mVariableNames.push_back("ytos");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.174403618112);

    this->mVariableNames.push_back("xtof");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.003643592594);

    this->mVariableNames.push_back("ytof");
    this->mVariableUnits.push_back("dimensionless");
    this->mInitialConditions.push_back(0.993331326442);

    this->mVariableNames.push_back("Ca_JSR");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(97.505463697266);

    this->mVariableNames.push_back("xir");
    this->mVariableUnits.push_back("uM_per_ms");
    this->mInitialConditions.push_back(0.006679257264);

    this->mVariableNames.push_back("Na_i");
    this->mVariableUnits.push_back("mM");
    this->mInitialConditions.push_back(11.441712311614);

    this->mVariableNames.push_back("Ca_dyad");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(1.716573130685);

    this->mVariableNames.push_back("Ca_submem");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(0.226941113355);

    this->mVariableNames.push_back("Ca_i");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(0.256752008084);

    this->mVariableNames.push_back("Ca_NSR");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(104.450004990523);

    this->mVariableNames.push_back("tropi");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(22.171689894953);

    this->mVariableNames.push_back("trops");
    this->mVariableUnits.push_back("uM");
    this->mInitialConditions.push_back(19.864701949854);

    this->mInitialised = true;
}


// Serialization for Boost >= 1.36
#include "SerializationExportWrapperForCpp.hpp"
CHASTE_CLASS_EXPORT(BackwardEulerMahajanModel2008)