File contents
// heart10Doc.cpp : implementation of the CHeart10Doc class
//
#include "stdafx.h"
#include "heart10.h"
#include "math.h"
#include "heart10Doc.h"
#include "heart10View.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
/////////////////////////////////////////////////////////////////////////////
// CHeart10Doc
IMPLEMENT_DYNCREATE(CHeart10Doc, CDocument)
BEGIN_MESSAGE_MAP(CHeart10Doc, CDocument)
//{{AFX_MSG_MAP(CHeart10Doc)
// NOTE - the ClassWizard will add and remove mapping macros here.
// DO NOT EDIT what you see in these blocks of generated code!
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CHeart10Doc construction/destruction
CHeart10Doc::CHeart10Doc()
{
NewCell = TRUE;
bUserRate = TRUE;
iVariablePulse = 2; //which pulse is variable pulse?
}
CHeart10Doc::~CHeart10Doc()
{
}
BOOL CHeart10Doc::OnNewDocument()
{
if (!CDocument::OnNewDocument())
return FALSE;
NewCell=TRUE;
CString Title;
Title = "Virtual Heart from SNU Laboratory";
SetTitle(Title);
return TRUE;
}
/////////////////////////////////////////////////////////////////////////////
// CHeart10Doc serialization
void CHeart10Doc::Serialize(CArchive& ar)
{
if (ar.IsStoring())
{
// TODO: add storing code here
}
else
{
// TODO: add loading code here
}
}
/////////////////////////////////////////////////////////////////////////////
// CHeart10Doc diagnostics
#ifdef _DEBUG
void CHeart10Doc::AssertValid() const
{
CDocument::AssertValid();
}
void CHeart10Doc::Dump(CDumpContext& dc) const
{
CDocument::Dump(dc);
}
#endif //_DEBUG
/////////////////////////////////////////////////////////////////////////////
// CHeart10Doc commands
void CHeart10Doc::DefIonVar()
{
Ke = Ion[2];
Nae =Ion[0];
Cae =Ion[4];
Mgi =Ion[6];
}
void CHeart10Doc::INITIALSETTING()
{
//----------------------- Default Variable Values before User Defined Value ---------//
if (!bUserRate)
{
//----------------- [Permeability ratio] -----------------------//
RGKNA = 0.086; // PK/PNa of INa (voltage gated Na channel)
RGITONa = 0.09; // PNa/PK of Ito (transient outward K-channel)
RGIsusNa = 0.09; // PNa/PK of Sustained Outward Na component
RGIKrNa = 0.09; // PNa/PK of IKr
RGIKsNa = 0.04; // PNa/PK of IKs
RGILNSCK = 1.4; // PK/PNa of NSC (non-selective cation channel)
RGILNSCCA = 0.7; // PCa/PNa of NSC (non-selective cation channel)
PCaLNa = 0.0001; // PNa/PCa of ICaL (L-type Ca-channel)
PCaLK = 0.002; // PK/PCa of ICaL (L-type Ca-channel)
RGfNa = 0.6216; // PNa/PK of If (hyperpolarization activated current)
//----------------- [ratio of second component rate constants-IKr] ------//
FS = 0.4; // IKr channel activation II
//--------------------- [Ca dynamics] ------------------------------//
KmTRPN = 0.00077; // dissociation constant of troponin (KmTRPN)
KmCMDN= 0.00238; // dissociation constant of calmodulin (KmCMDN)
KmCSQN = 0.8; // dissociation constant of calsequestrin (KmCSQN)
TRPNt = 70.0e-3; // total concentration of of troponin in mM/L (TRPNt)
CMDNt = 50.0e-3; // total concentration of of calmodulin in mM/L (CMDNt)
CSQNt = 10.0; // total concentration of of calsequestrin in mM/L (CSQNt)
RVrel = 0.02; //0.004, ratio of Ca releasing store volume to cell volume (RVrel)
RVup = 0.05; // 0.18, ratio of Ca uptaking store volume to cell volume (RVup)
Caupmax = 1.0; // maximum Ca concentration in uptake site (Caupmax)
Caupinitial = 2.590; // initial Ca concentration in uptake pool (C[4][6][L])
Carelinitial = 9.422; // initial Ca concentration in release pool (C[5][6][L])
Taurel = 1.0; // time constant of release (Taurel)
Tautr = 2020.0; // time constant of transfer (Tautr)
TauUp = 0.05; // time constant of uptake (TauUp)
KmCa = 0.001; // Km value of activation of the Ca release from SR (KmCa)
SF_Arel = 0.01*2.0; // amplitude factor for Ca release (Arel)
SF_Atr = 2.0; // amplitude factor for Ca transfer (Atr)
SF_Aup = 2.0; // amplitude factor for Ca uptake (Aup)
//------------------------- [Na-K pump] --------------------------------//
KmNa = 11.0; // Km value for Na sensitivity of Na-K pump (KmNa)
Partition_v = 0.32; // Partition_v of Na-K pump
//------------------------- [Na-Ca Exchanger] ------------------------------//
KNAI = 670; // dissociation constant for Na at the inner side of the membrane (NaCa exchanger)
KNAO = 670000; // dissociation constant for Na at the outer side of the membrane (NaCa exchanger)
KCaI = 0.00138; // dissociation constant for Ca at the inner side of the membrane (NaCa exchanger)
KCaO = 1.38; // dissociation constant for Ca at the outer side of the membrane (NaCa exchanger)
BasalRate = 250.0; // basal rate constant for the Na-Ca exchanger (BasalRate)
//----------------------- [rectifying field of IK channel] -------------------//
RectifyingK = 50; // if minus, then inward rectifying
}
/*----------------------- Physical Dimension of Cell-------------------------------
Width = 20 micron, Length = 100 micron, Thickness = 8 micron
Surface Area = 2 * 100 * 20 + 2 * 100 * 8 + 2 * 20 * 8 = 5920 micron^2 = 5920 * 0.01 pF = 59 pF
Capacitance by Folding Surface Area = 2.0 * Surface Area = 2.0 * 59 pF = 118 pF
Cell Volume = 100 * 20 * 8 * 0.5 (um^3)
multiply 0.5 because osmoticaly active volume is smaller than the bulk volume
See also Fujioka et al, 1998. 3.14 * (30/2)^2 * 100 = 70650 micron^3
----------------------------------------------------------------------------------*/
Cap =SF[23] * 132.0;
Vi = 100.0 * 20.0 * 8.0 / 2.0 * pow(SF[23]*0.132 / 0.132,3.0 / 2.0);
Vtotal =Condition[2] * Vi; // Cell Volume + Pipette Volume
Vrel = RVrel * Vi; //0.02 original VOLUME OF SR RELEASING STORE
Vup = RVup * Vi; //0.05 original VOLUME OF SR UPTAKE STORE
//------------------------ Conductance of Ion Channel ----------------------------//
GIRK1[1] = SF[3] * 0.8 * 1800.0 / 20.0 * SF[23]; //IK1 component, adjusted to give 1.2 nA outward current at -40 mV
GNa[1] =SF[4]* 1.875 * 58.0 * SF[23]; //400 'sodium current
GKr[1] =SF[5]* 132 * 0.108 * 0.3 * SF[23]; //*0.4 0.099 original from Ono. Conductance factor for IK,rapid?
GITOK[1] =SF[6] * 1.5 * SF[23]; //Transient Outward K component, Cambell et al., J.Gen. Physiol,vol 101,571-601,1993
GIsusK[1] =SF[7] * 0.3 * SF[23]; //Ca dependent K, K component, Cambell et al., J.Gen. Physiol,vol 101,571-601,1993
PCaL[1] =SF[8] * 24 * 40 * SF[23]/2.0; //IcaLAmplitude '20 Ca Permeability factor for L-Ca channel
PCaT[1] =SF[9] * 66 * 20 * SF[23]; //T-type Ca channel amplitude
GIKsK[1] =SF[10] * 0.247 * 2.5 / 2 * SF[23]; //'IKs Amplitude factor , ZengJL et al, 1995, CR 77:140-152
GfK[1] = SF[11] * 1.621 * 7.4 * SF[23] / (1 + RGfNa);
GfNa[1] = SF[11] * 1.621 * 7.4 * SF[23] * RGfNa / (1 + RGfNa);
ACh[1] = SF[12] * SF[23];
GbK[1] = SF[18] * SF[23] / 3.0;
GSAC[1] = SF[19] * 0.01 * 275 * SF[23] / 2.0; //0.04 Ca component conductance
GILNSCNa[1] = 1.8 * SF[20] * 0.0561 * 1.5 * SF[23]; //leak Na component
RGILNSCCA = 1.5; // PCa/PNa of NSC (non-selective cation channel)
//------------------------- Conductance of Na-K pump --------------------------------//
INaKmax[1] = SF[22] * 132.0 * 2.88 * 1.3 * SF[23]; //1.19 'SAKAI ET AL. J.P. 490:51-62
//------------------------- Conductance of Na-Ca exchanger ------------------------------//
NumberofCareer = 132.0 * 100.0 * 2000.0 * 4.0 * SF[23]; //(pF*100)um^2 * DENSITY/um^2
INaCaMax[1] = SF[21] * NumberofCareer * Elementarycharge * 1000000000000.0 * BasalRate * 1; //pA
//--------------------- Ca dynamics ------------------------------//
Arel = SF_Arel * Vrel * F / Taurel; //Amplitude factor for Ca release
Atr = SF_Atr * Vi * F / Tautr; //Amplitude factor for Ca transfer
Aup = SF_Aup * Vi * F / TauUp; //Amplitude factor for Ca uptake
//Aleak = 0.2 * 450.0 / 2.0; //Amplitude factor for Ca leak
}
void CHeart10Doc::calculatevm()
{
for (int L = 1;L<=cellNo;L++)
{
if (m_AP == 0)
{
Vm[L] = Erev;
iterative = 200;
}
else
{
Vm[L] = Pulse[0];
iterative = 40;
}
qr = 0.0;
}
do
{
for (L = 1;L<=cellNo;L++)
{
if (fabs(d[L]) > limitfactor)
{
Vm[L] = Vm[L] + pow(-1.0, qr) * 0.0005 * qr * iterative;
}
}
CABUFFER(); //to calculate concentration of free Ca in cytosol from x(16,5).
for (L = 1;L<=cellNo;L++)
{
VNa[L] = RTF * log(Nae / C[2][6][L]); //REVERSAL POTENTIAL FOR Na
VK[L] = RTF * log(Ke / C[1][6][L]); //REVERSAL POTENTIAL FOR K
VCa[L] = RTF2 * log(Cae / Caf[L]); //REVERSAL POTENTIAL FOR Ca
}
constantfield_noma();
RateConstants(); //determine x(y,1) & x(y,2)
for (L = 1;L<=cellNo;L++)
{
for (int N = 3;N<=NOP;N++)
{
m[N][6][L] = m[N][1][L] / (m[N][1][L] + m[N][2][L]);
m[N][7][L] = m[N][6][L];
h[N][6][L] = h[N][1][L] / (h[N][1][L] + h[N][2][L]);
h[N][7][L] = h[N][6][L];
}
for (N = 1;N<=3;N++)
{
Org[N][6][L] = Org[N][1][L] / (Org[N][1][L] + Org[N][2][L]);
Org[N][7][L] = Org[N][6][L];
}
mm[5][6][L] = mm[5][1][L] / (mm[5][1][L] + mm[5][2][L]);
mm[5][7][L] = mm[5][6][L];
hh[4][6][L] = hh[4][1][L] / (hh[4][1][L] + hh[4][2][L]);
hh[4][7][L] = hh[4][6][L];
m[3][6][L] = m[3][1][L] / (m[3][1][L] + m[3][2][L] * pow(1 - h[3][3][L],3.0)); //Inward rectifier K current
m[3][7][L] = m[3][6][L];
}
DefSteps(FALSE);
currentamplitude();
for (L = 1;L<=cellNo;L++)
{
if (fabs(d[L]) > limitfactor)
{
Itot_vm[L] = Itot_vm[L] + Vm[L] / (0.001*Condition[1]);
}
if (m_AP == 0)
{
d[1] = Itot_vm[1];
limitfactor = 20;
}
else
{
d[1] = Pulse[0] - Vm[1] - Itot_vm[1] * (0.001*Condition[0]);
limitfactor = 0.1;
}
}
qr = qr + 1;
} while(!(fabs(d[1]) < limitfactor));
for (L = 1;L<=cellNo;L++)
{
for (int N = 3;N<=NOP;N++)
{
m[N][7][L] = m[N][6][L];
h[N][7][L] = h[N][6][L];
}
mm[5][7][L] = mm[5][6][L];
hh[4][7][L] = hh[4][6][L];
for (N = 1;N<=3;N++)
{
Org[N][7][L] = Org[N][6][L];
}
C[6][7][L] = Vm[L];
C[6][6][L] = C[6][7][L];
}
}
void CHeart10Doc::CalSteps()
{
DefIonVar();
MemoryCurrent();
CABUFFER(); //to calculate concentration of free Ca in cytosol from x(16,5).
for (int L = 1;L<=cellNo;L++)
{
VNa[L] = RTF * log(Nae / C[2][6][L]); //REVERSAL POTENTIAL FOR Na
VK[L] = RTF * log(Ke / C[1][6][L]); //REVERSAL POTENTIAL FOR K
VCa[L] = RTF2 * log(Cae / Caf[L]); //REVERSAL POTENTIAL FOR Ca
}
RateConstants(); //determine x(y,1) & x(y,2)
constantfield_noma();
DefSteps(m_AP);
currentamplitude();
differen();
RUNGEKUTTA(); //INTEGRATION of parameters
CurrentTime = CurrentTime + timeinterval;
}
void CHeart10Doc::dataseparate()
{
}
void CHeart10Doc::CABUFFER()
{
for (int L = 1;L<=cellNo;L++)
{
int N = 1;
do
{
Caf[L] = C[3][6][L] / (1.0 + (TRPNf[L] / KmTRPN) + (CMDNf[L] / KmCMDN));
TRPNf[L] = TRPNt / (1.0 + (Caf[L] / KmTRPN));
CMDNf[L] = CMDNt / (1.0 + (Caf[L] / KmCMDN));
N = N + 1;
} while(N<=10);
N = 1;
do
{
Cafrel[L] = C[5][6][L] / (1.0 + (CSQNf[L] / KmCSQN));
CSQNf[L] = CSQNt / (1.0 + (Cafrel[L] / KmCSQN));
N = N + 1;
} while(N<=10);
}
}
void CHeart10Doc::constantfield_noma()
{
for (int L = 1;L<=cellNo;L++)
{
float RectifyingFactor;
if (bRectify) //rectifing current equation= Vm[L]
{
RectifyingFactor = 50;
}
else
{
RectifyingFactor = 0;
}
RTF = R * (Condition[3]+273.15) / F;
RTF2 = R * (Condition[3]+273.15) / F / 2.0;
if (Vm[L] == RectifyingFactor) //******* constant field for Ca *********
{
CaCF[L] = 0.99906;
}
else
{
CaCF[L] = ((Vm[L] - RectifyingFactor) / RTF2) / (1 - exp(-(Vm[L] - RectifyingFactor) / RTF2));
}
CaCF[L] = CaCF[L] * (Caf[L] * exp(RectifyingFactor / RTF2) - Cae * exp(-(Vm[L] - RectifyingFactor) / RTF2));
if (Vm[L] == RectifyingFactor) //******* constant field for K *********
{
KCF[L] = 0.99906;
}
else
{
KCF[L] = ((Vm[L] - RectifyingFactor) / RTF) / (1.0 - exp(-(Vm[L] - RectifyingFactor) / RTF));
}
KCF[L] = KCF[L] * (C[1][6][L] * exp(RectifyingFactor / RTF) - Ke * exp(-(Vm[L] - RectifyingFactor) / RTF));
if (Vm[L] == RectifyingFactor) //******* constant field for Na *********
{
NaCF[L] = 0.99906;
}
else
{
NaCF[L] = ((Vm[L] - RectifyingFactor) / RTF) / (1.0 - exp(-(Vm[L] - RectifyingFactor) / RTF));
}
NaCF[L] = NaCF[L] * (C[2][6][L] * exp(RectifyingFactor / RTF) - Nae * exp(-(Vm[L] - RectifyingFactor) / RTF));
}
}
void CHeart10Doc::RateConstants()
{
if (bUserRate == TRUE)
{
RateConstantsUser();
}
else
{
RateConstantsDefault(); //RateConstantsUser();
}
for (int L=1;L<=cellNo;L++)
{
//********** Ih current******************************
m[11][1][L] = exp(-0.0220741 * (Vm[L] + 386.9)); //wilders
m[11][2][L] = exp(0.052 * (Vm[L] - 73.08)); //wilders
m[11][3][L] = m[11][1][L] / (m[11][1][L] + m[11][2][L]);
m[11][4][L] = 1.0 / (m[11][1][L] + m[11][2][L]);
h[11][1][L] = exp(-0.0220741 * (Vm[L] + 386.9)); //wilders
h[11][2][L] = exp(0.052 * (Vm[L] - 73.08)); //wilders
h[11][3][L] = h[11][1][L] / (h[11][1][L] + h[11][2][L]);
h[11][4][L] = 1.0 / (h[11][1][L] + h[11][2][L]);
//***********Ach-activated K current****************
if (ACh[L] == 0.0)
{
m[12][1][L] = 0.165 * 0.001 * 1.5; //spontaneous opening /msec"
}
else
{
m[12][1][L] = 12.32 / (1.0 + 4.2 / ACh[L]) * 0.001;
}
m[12][2][L] = 10.0 * exp(0.0133 * (Vm[L] + 40.0)) * 0.001; //beta of IKACh
h[12][1][L] = exp(-0.0220741 * (Vm[L] + 386.9)); //wilders
h[12][2][L] = exp(0.052 * (Vm[L] - 73.08)); //wilders
h[12][3][L] = h[12][1][L] / (h[12][1][L] + h[12][2][L]);
h[12][4][L] = 1.0 / (h[12][1][L] + h[12][2][L]);
//************ Ca transfer in the Ca store***************************
if (Vm[L] == -64.0)
{
Org[2][1][L] = 0.0200047;
}
else
{
Org[2][1][L] = 0.000625 * 8.0 * (Vm[L] + 64.0) / (exp(0.25 * (Vm[L] + 64.0)) - 1.0);
}
Org[2][2][L] = 0.005 * 8.0 / (1.0 + exp(-0.25 * (Vm[L] + 64.0)));
//------------- Ca uptake in the network SR ---------------------------//
float k1SERCA = 0.01; // Kyoto-model, 2003
float k2SERCA = 1.0; // Kyoto-model, 2003
float k3SERCA = 1.0/(1.0+0.1/4.663); // Kyoto-model, 2003
float k4SERCA = 0.01; // Kyoto-model, 2003
float KmCai = 0.08; // mM
float KmCao = 0.0008; // mM
float pE1Ca = 1.0/(1+KmCai/C[4][6][L]);
float pE2Ca = 1.0/(1+KmCao/Caf[L]);
float pE1 = 1 - pE1Ca;
float pE2 = 1 - pE2Ca;
Org[3][1][L] = k2SERCA*pE2Ca+k4SERCA*pE2;
Org[3][2][L] = k1SERCA*pE1Ca+k3SERCA*pE1;
//************ NACA EXCHANGE PINGPON ****************
PNAO = 1.0 / (1.0 + ((1.0 + Cae / KCaO) / (pow(Nae,3.0) / KNAO))); //'PROBABILITY OF Na BINDING
PNAI = 1.0 / (1.0 + ((1.0 + Caf[L] / KCaI) / (pow(C[2][6][L],3.0) / KNAI))); //'PROBABILITY OF Na BINDING
PCAO = 1.0 / (1.0 + ((1.0 + pow(Nae,3.0) / KNAO) / (Cae / KCaO))); //'PROBABILITY OF Na BINDING
PCAI = 1.0 / (1.0 + ((1.0 + pow(C[2][6][L],3.0) / KNAI) / (Caf[L] / KCaI))); //'PROBABILITY OF Na BINDING
k12 = PCAO;
k21 = PCAI;
K30[L] = exp(Partition_v * Vm[L] / RTF) * PNAI;
K03[L] = exp((Partition_v - 1.0) * Vm[L] / RTF) * PNAO;
Org[1][2][L] = K30[L] + k21; //RATE FROM PXI
Org[1][1][L] = K03[L] + k12; //RATE TO PXI
/*------------------ RyR and SR Ca2+ kinetics -------------------------------//
RyR channel
(1) C -> O
kco = 280000*[Ca2+]i^2 - 150*ICaL*P(openCaL)
(2) O -> U
kou = 0.08/(1+0.36/[Ca]rel)
(3) U <-> C
kuc = 0.000377*[Ca]rel^2
kcu = 0.000849
I(RyR) = P(RyR)*([Ca]rel - [Ca]i)*P(openRyR)
P(RyR) = 62000
//----------------------------------------------------------------------------*/
kco[L] = 280000.0*pow(Caf[L], 2.0) - 150.0*0.0676*CaCF[L]*m[8][6][L]*h[8][6][L];
kou[L] = 0.08/(1+0.36/Cafrel[L]);
kuc[L] = 0.000377*pow(Cafrel[L], 2.0);
kcu[L] = 0.000849;
}
}
void CHeart10Doc::currentamplitude()
{
float RKCF;
float RKCF2;
float Tune;
for (int L = 1;L<=cellNo;L++)
{
//------------------- SODIUM CURRENT ----------------------------------------------------------//
ii[4][1][L] = pow(m[4][6][L],3.0) * h[4][6][L] * hh[4][6][L] * GNa[L] * RGKNA * KCF[L];
ii[4][2][L] = pow(m[4][6][L],3.0) * h[4][6][L] * hh[4][6][L] * GNa[L] * NaCF[L]; //(E - VNa)
ii[4][4][L] = ii[4][1][L] + ii[4][2][L];
//************* Capacitive transient****
if (timeintervalb == 0)
{
Icap[L] = Cap * C[6][5][L];
}
else
{
Icap[L] = Cap * (C[6][6][L] - C[6][7][L]) / timeintervalb;
}
//******** Ito **********
ii[6][1][L] = (m[6][6][L] * h[6][6][L]) * GITOK[L] * KCF[L];
ii[6][2][L] = (m[6][6][L] * h[6][6][L]) * GITOK[L] * RGITONa * NaCF[L];
ii[6][4][L] = ii[6][1][L] + ii[6][2][L];
//******** Isus **********
ii[7][1][L] = m[7][6][L] * h[7][6][L] * GIsusK[L] * KCF[L];
ii[7][2][L] = m[7][6][L] * h[7][6][L] * GIsusK[L] * RGIsusNa * NaCF[L];
ii[7][4][L] = ii[7][1][L] + ii[7][2][L];
//******* L-type Ca channel ***********
ICaLCainf = PCaL[L]* CaCF[L];
ICaLKinf = PCaLK * PCaL[L] * KCF[L];
ICaLNainf = PCaLNa * PCaL[L] * NaCF[L];
ii[8][3][L] = ICaLCainf * m[8][6][L] * h[8][6][L];
ii[8][1][L] = ICaLKinf * m[8][6][L] * h[8][6][L];
ii[8][2][L] = ICaLNainf * m[8][6][L] * h[8][6][L];
ii[8][4][L] = ii[8][1][L] + ii[8][2][L] + ii[8][3][L];
//************ ICaT ****************
ICaTInf = PCaT[L] * CaCF[L];
ii[9][4][L] = m[9][6][L] * h[9][6][L]*ICaTInf;
//*************** IKr ***************
FF = 1-FS;
PA = FF * m[5][6][L] + FS * mm[5][6][L];
PI = h[5][6][L];
ii[5][1][L] = GKr[L] * PA * PI * (Vm[L] - VK[L]); //0.099
ii[5][2][L] = GKr[L] * PA * PI * RGIKrNa * (Vm[L] - VNa[L]); //0.099
ii[5][4][L] = ii[5][1][L] + ii[5][2][L]; //0.099
//*********** IRK1 **************
ii[3][4][L] = GIRK1[L] * (Vm[L] - VK[L]) * pow((Ke / 5.4),0.62) * m[3][6][L] * (1-h[3][3][L]);
//************** IKs *****************
ii[10][1][L] = m[10][6][L] * GIKsK[L] * KCF[L];
ii[10][2][L] = m[10][6][L] * RGIKsNa * GIKsK[L] * NaCF[L];
ii[10][4][L] = ii[10][1][L] + ii[10][2][L];
//******** Iha **********
ii[11][1][L] = GfK[L] * (Vm[L] - VK[L]) * pow(m[11][6][L],2.0); //Wilders
ii[11][2][L] = GfNa[L] * (Vm[L] - VNa[L]) * pow(m[11][6][L],2.0); //Wilders
ii[11][4][L] = ii[11][1][L] + ii[11][2][L];
//********* iKACh *******
ii[12][4][L] = 0.27 / 2.0 * Cap * m[12][6][L] * (Vm[L] - VK[L]) / (1.0 + exp((Vm[L] + 20.0) / 20.0));
//*********** exchangers ***************************
float RectifyingFactor = 50;
ii[21][4][L] = INaCaMax[L] * (Org[1][6][L]* K30[L] - (1.0 - Org[1][6][L]) * K03[L]);
ii[22][4][L] = (1.0 / (1.0 + pow(KmNa / C[2][6][L],1.36))) * (1.0 - pow((Vm[L] - RectifyingFactor) / 125.0 / 2.0 ,2.0)) * INaKmax[L];
//************ background K currents ******************
if (Vm[L] == RectifyingK) //constant field for background K channel
{
RKCF = 0.99906;
}
else
{
RKCF = ((Vm[L] - RectifyingK) / RTF) / (1.0 - exp(-(Vm[L] - RectifyingK) / RTF));
}
if (RectifyingK == 0)
{
Tune = 1.0;
}
else
{
Tune = pow(fabs(RectifyingK), 0.1*fabs(RectifyingK)/RectifyingK);
}
RKCF2 = RKCF * (C[1][6][L] - Ke * exp(-Vm[L] / RTF));
ii[18][4][L] = GbK[L] * RKCF2 * Tune ; //Conductance X constant field X amplitude tuning
/*------------ Stretch-Activated Channel ------------------------------------*/
ii[19][2][L] = 1.0 * GSAC[L] * NaCF[L]; //(E - VCa)
ii[19][1][L] = 1.32 * GSAC[L] * KCF[L]; //(E - VCa)
ii[19][3][L] = 0.7 * GSAC[L] * CaCF[L]; //(E - VCa)
ii[19][4][L] = ii[19][2][L]+ii[19][1][L]+ii[19][3][L]; //(E - VCa)
/*------------ NONSELECTIVE CATION currents ----------------------------------
------------------------------------------------------------------------------*/
ii[20][2][L] = GILNSCNa[L] * NaCF[L] * 1 / (1 + RGILNSCK + RGILNSCCA) ;
ii[20][1][L] = GILNSCNa[L] * KCF[L] * RGILNSCK / (1 + RGILNSCK + RGILNSCCA);
ii[20][3][L] = GILNSCNa[L] * CaCF[L] * RGILNSCCA / (1 + RGILNSCK + RGILNSCCA);
ii[20][4][L] = ii[20][1][L] + ii[20][2][L] + ii[20][3][L];
/*------------------ RyR and SR Ca2+ kinetics -------------------------------//
RyR channel
(1) C -> O
kco = 280000*[Ca2+]i^2 - 150*ICaL*P(openCaL)
(2) O -> U
kou = 0.08/(1+0.36/[Ca]rel)
(3) U <-> C
kuc = 0.000377*[Ca]rel^2
kcu = 0.000849
I(RyR) = P(RyR)*([Ca]rel - [Ca]i)*P(openRyR)
P(RyR) = 62000
//----------------------------------------------------------------------------*/
Irel[L] = 280000.0*50.0/132.0*(Cafrel[L] - Caf[L])*RyROpen[L];
Atr = 386*50.0/132.0;
Itr[L] = Atr * (C[4][6][L] - Cafrel[L]);
Ileak[L] = 0.1*459.0*50.0/132.0*(C[4][6][L] - Caf[L]); // leak component
float k1SERCA = 0.01; // Kyoto-model, 2003
float k2SERCA = 1.0; // Kyoto-model, 2003
float KmCai = 0.08; // mM
float KmCao = 0.0008; // mM
float pE1Ca = 1.0/(1+KmCai/C[4][6][L]);
float pE2Ca = 1.0/(1+KmCao/Caf[L]);
float pE1 = 1 - pE1Ca;
float pE2 = 1 - pE2Ca;
Iup[L] = 42500.0 * 50.0/132.0 * (k2SERCA*pE2Ca*(1-Org[3][6][L]) - k1SERCA*pE1Ca*Org[3][6][L]); // Iup[L] = Aup * Caf[L] * (1.0 - (C[4][6][L] / Caupmax));
//******************* net ION flux to estimate cellular concentration *****
ICanet[L] = ii[19][3][L] + ii[20][3][L] + ii[5][2][L] + ii[8][3][L] + ii[9][4][L] - 2.0 * ii[21][4][L] - Irel[L] - Ileak[L] + Iup[L];
IKnet[L] = ii[19][1][L] + ii[20][1][L] + ii[8][1][L] + ii[5][1][L] + ii[6][1][L]+ ii[7][1][L] + ii[10][1][L] + ii[4][1][L] + ii[3][4][L] - 2.0 * ii[22][4][L] + ii[11][1][L] + ii[12][4][L] + ii[18][4][L];
INanet[L] = ii[19][2][L] + ii[20][2][L] + ii[8][2][L] + ii[6][2][L] + ii[7][2][L] + ii[4][2][L] + ii[10][2][L] + 3.0 * ii[21][4][L] + 3.0 * ii[22][4][L] + ii[11][2][L];
//******************* net current to calculate dV/dt ***********************
Itot_vm[L] = ii[3][4][L] + ii[4][4][L] + ii[5][4][L] + ii[6][4][L] + ii[7][4][L] + ii[8][4][L] + ii[9][4][L] + ii[10][4][L] + ii[11][4][L] + ii[12][4][L] + ii[18][4][L] + ii[19][4][L] + ii[20][4][L] + ii[21][4][L] + ii[22][4][L];
float tempTot;
if (m_AP == 0)
{
Itot[L] = Itot_vm[L] + Iinj[L] + Vm[L] / (0.001*Condition[1]);
tempTot = fabs(Iinj[L]) + fabs(Vm[L] / (0.001*Condition[1])) + fabs(ii[3][4][L]) + fabs(ii[4][4][L]) + fabs(ii[5][4][L]) + fabs(ii[6][4][L]) + fabs(ii[7][4][L]) + fabs(ii[8][4][L]) + fabs(ii[9][4][L]) + fabs(ii[10][4][L]) + fabs(ii[11][4][L]) + fabs(ii[12][4][L]) + fabs(ii[18][4][L]) + fabs(ii[19][4][L]) + fabs(ii[20][4][L]) + fabs(ii[21][4][L]) + fabs(ii[22][4][L]);
}
else
{
Itot[L] = Itot_vm[L] + Icap[L] + Vm[L] / (0.001*Condition[1]);
tempTot = fabs(Icap[L]) + fabs(Vm[L] / (0.001*Condition[1])) + fabs(ii[3][4][L]) + fabs(ii[4][4][L]) + fabs(ii[5][4][L]) + fabs(ii[6][4][L]) + fabs(ii[7][4][L]) + fabs(ii[8][4][L]) + fabs(ii[9][4][L]) + fabs(ii[10][4][L]) + fabs(ii[11][4][L]) + fabs(ii[12][4][L]) + fabs(ii[18][4][L]) + fabs(ii[19][4][L]) + fabs(ii[20][4][L]) + fabs(ii[21][4][L]) + fabs(ii[22][4][L]);
}
Itot_up = (tempTot+Itot[L])*0.5;
Itot_down = (-tempTot+Itot[L])*0.5;
}
}
void CHeart10Doc::RateConstantsDefault()
{
for (int L = 1;L<=cellNo;L++)
{
//************* IRK1 CURRENT ******************
h[3][1][L] = 12.75*exp(0.035*(Vm[L]-VK[L]-10))/(1+exp(0.015*(Vm[L]-VK[L]-140)));
h[3][2][L] = 3*exp(-0.048*(Vm[L]-VK[L]-10))*(1+exp(0.064*(Vm[L]-VK[L]-38)))/(1+exp(0.03*(Vm[L]-VK[L]-70)));
h[3][3][L] = h[3][1][L] / (h[3][1][L] + h[3][2][L]);
h[3][4][L] = 1 / (h[3][1][L] + h[3][2][L]);
m[3][1][L] = 1.0/(8000*exp((Vm[L] - VK[L] - 97)/8.5)+7*exp((Vm[L] - VK[L] - 97)/300));
m[3][2][L] = 1.0/(0.00014*exp((Vm[L] - VK[L] - 97)/(-9.1))+0.2*exp((Vm[L] - VK[L] - 97)/(-500)));
m[3][3][L] = m[3][1][L] / (m[3][1][L] + m[3][2][L]);
m[3][4][L] = 1 / (m[3][1][L] + m[3][2][L]);
//*************** INa (voltage gated Na current) from Luo and Rudy, circ. res, vol 74(6),1071-1113 ***************
m[4][3][L] = 1.0 / (1.0 + exp((Vm[L] + 52.0) / (-5.27))); //steady-state open of m-gate, 1999.4.27 modified from Luo and Rudy
m[4][4][L] = 0.02 + 1.0 / (0.4 * exp((Vm[L] + 68.0) / (-2.9)) + 0.09 * exp((Vm[L] + 55.0) / 4.7) + 8.0); //tau of m-gate, 1999.4.27 modified
m[4][1][L] = m[4][3][L] / m[4][4][L]; //alfa of m-gate
m[4][2][L] = (1.0 - m[4][3][L]) / m[4][4][L]; //beta of m-gate
if (Vm[L] >= -40.0)
{
h[4][1][L] = 0.0; //alpha h
h[4][2][L] = 1.0 / (0.13 * (1.0 + exp((Vm[L] + 10.66) / (-11.1)))); //beta h
h[4][3][L] = h[4][1][L] / (h[4][1][L] + h[4][2][L]);
h[4][4][L] = 1.0 / (h[4][1][L] + h[4][2][L]);
hh[4][1][L] = 0.0; //alpha j
hh[4][2][L] = 0.3 * exp(-0.0000002535 * Vm[L]) / (1.0 + exp(-0.1 * (Vm[L] + 32.0))); //beta j
hh[4][3][L] = hh[4][1][L] / (hh[4][1][L] + hh[4][2][L]);
hh[4][4][L] = 1.0 / (hh[4][1][L] + hh[4][2][L]);
}
else
{
h[4][1][L] = 0.135 * exp((80.0 + Vm[L]) / (-6.8)); //alpha h
h[4][2][L] = 3.56 * exp(0.079 * Vm[L]) + 310000.0 * exp(0.35 * Vm[L]);
h[4][3][L] = h[4][1][L] / (h[4][1][L] + h[4][2][L]);
h[4][4][L] = 1.0 / (h[4][1][L] + h[4][2][L]);
hh[4][1][L] = (0.000012714 * exp(0.2444 * Vm[L]) - 0.00003474 * exp(-0.04391 * Vm[L])) * (Vm[L] + 37.78) / (1.0 + exp(0.311 * (Vm[L] + 79.23)));
hh[4][2][L] = 0.1212 * exp(-0.01052 * Vm[L]) / (1.0 + exp(-0.1378 * (Vm[L] + 40.14)));
hh[4][3][L] = hh[4][1][L] / (hh[4][1][L] + hh[4][2][L]);
hh[4][4][L] = 1.0 / (hh[4][1][L] + hh[4][2][L]);
}
//------------------ Ik RAPID -------------------------------------------//
Painfinite = 1.0 / (1.0 + exp((-23.2 - Vm[L]) / 10.6));
tauaf = 1000.0 / (37.2 * exp(Vm[L] / 15.9) + 0.96 * exp(-Vm[L] / 22.5));
tauas = 1000.0 / (4.2 * exp(Vm[L] / 17.0) + 0.15 * exp(-Vm[L] / 21.6));
m[5][1][L] = Painfinite / tauaf;
m[5][2][L] = (1.0 - (1.0 / (1.0 + exp((-23.2 - Vm[L]) / 10.6)))) / (1000.0 / (37.2 * exp(Vm[L] / 15.9) + 0.96 * exp(-Vm[L] / 22.5)));
m[5][3][L] = (1.0 - Painfinite) / tauaf;
m[5][4][L] = 1.0 / (m[5][1][L] + m[5][2][L]);
mm[5][1][L] = Painfinite / tauas;
mm[5][2][L] = (1.0 - Painfinite) / tauas;
mm[5][3][L] = mm[5][1][L] / (mm[5][1][L] + mm[5][2][L]);
mm[5][4][L] = 1.0 / (mm[5][1][L] + mm[5][2][L]);
Piinfinite = 1.0 / (1.0 + exp((Vm[L] + 28.6) / 17.1));
taui = 2.0;
h[5][1][L] = Piinfinite / taui;
h[5][2][L] = (1 - Piinfinite) / taui;
h[5][3][L] = h[5][1][L] / (h[5][1][L] + h[5][2][L]);
h[5][4][L] = 1.0 / (h[5][1][L] + h[5][2][L]);
/* ------------------------------- Ito CURRENT -----------------------------------
Apkon M, Nerbonne JM (1991) Characterization of two distinct depolarization-activated
K+ currents in isolated adult rat ventricular myocytes. J Gen Physiol 97: 973-1011
1. Reversal potential: -74.2 mV(peak), -66 mV (plateau)
-----------------------------------------------------------------------------------*/
m[6][1][L] = 1 / (1.4 * exp(Vm[L]/(-14241)) + 0.35 * exp(Vm[L]/(-15)));
m[6][2][L] = 1 / (22.7 * exp(Vm[L]/37.9) + 1775 * exp(Vm[L]/8.4));
m[6][3][L] = m[6][1][L] / (m[6][1][L] + m[6][2][L]);
m[6][4][L] = 1/ (m[6][1][L] + m[6][2][L]);
h[6][1][L] = 1 / (2.5 * exp(Vm[L]/3.2) + 2195 * exp(Vm[L]/23.6));
h[6][2][L] = 1 / (40.4 * exp(Vm[L]/(-76022)) + 0.29 * exp(Vm[L]/(-5.4)));
h[6][3][L] = h[6][1][L] / (h[6][1][L] + h[6][2][L]);
h[6][4][L] = 1 / (h[6][1][L] + h[6][2][L]);
//**************Sustained Outward K+ current (Human atrium) **************
m[7][1][L] = 1 / (0.0002 * exp(Vm[L]/(-7.2)) + 1.37 * exp(Vm[L]/(-47.8)));
m[7][2][L] = 1 / (1346 * exp(Vm[L]/(16.0)) + 1600 * exp(Vm[L]/(11.3)));
m[7][3][L] = m[7][1][L] / (m[7][1][L] + m[7][2][L]);
m[7][4][L] = 1/ (m[7][1][L] + m[7][2][L]);
h[7][1][L] = 1 / (6613 * exp(Vm[L]/135222) + 143392 * exp(Vm[L]/12.1));
h[7][2][L] = 1 / (4801 * exp(Vm[L]/(-125)) + 119 * exp(Vm[L]/(-10.3)));
h[7][3][L] = h[7][1][L] / (h[7][1][L] + h[7][2][L]);
h[7][4][L] = 1 / (h[7][1][L] + h[7][2][L]);
//----------------- ICA L CURRENT ------------------------------//
if (Vm[L]== -30.0)
{
m[8][1][L] = 0.02792845;
}
else
{
m[8][1][L] = 0.01045 * (Vm[L] + 30.0) / (1.0 - exp(-(Vm[L] + 30.0) / 2.5)); //+ .03125 * E / (1 - EXP(-E / 4.8)) 'alpha, activation parameter for L-type Ca channel
}
if (Vm[L] == 0.0)
{
m[8][2][L] = 0.01052639;
}
else
{
m[8][2][L] = 0.00421 * Vm[L] / (exp(Vm[L] / 2.5) - 1.0); //beta
}
m[8][3][L] = m[8][1][L] / (m[8][1][L] + m[8][2][L]);
m[8][4][L] = 1.0 / (m[8][1][L] + m[8][2][L]);
if (Vm[L] == -34.0)
{
h[8][1][L] = 0.00200021;
}
else
{
h[8][1][L] = 0.000355 * (Vm[L] + 34.0) / (exp((Vm[L] + 34.0) / 5.633) - 1.0);
}
h[8][2][L] = 0.0007 * (Vm[L] + 64.0) / (1.0 + exp(-(Vm[L] + 44.0) / 4.16)); //beta
h[8][3][L] = h[8][1][L] / (h[8][1][L] + h[8][2][L]);
h[8][4][L] = 1.0 / (h[8][1][L] + h[8][2][L]);
//************ ICA T CURRENT *******************************
DTinfinite = 1.0 / (1.0 + exp(-(Vm[L] + 23.0) / 6.1));
tauDT = 0.6 + 5.4 / (1.0 + exp(0.03 * (Vm[L] + 100.0)));
m[9][1][L] = DTinfinite / tauDT;
m[9][2][L] = (1.0 - DTinfinite) / tauDT;
m[9][3][L] = m[9][1][L] / (m[9][1][L] + m[9][2][L]);
m[9][4][L] = 1.0 / (m[9][1][L] + m[9][2][L]);
FTinfinite = 1.0 / (1 + exp((Vm[L] + 75) / 6.6));
tauFT = 1.0 + 40.0 / (1 + exp(0.08 * (Vm[L] + 65.0)));
h[9][1][L] = FTinfinite / tauFT;
h[9][2][L] = (1.0 - FTinfinite) / tauFT;
h[9][3][L] = h[9][1][L] / (h[9][1][L] + h[9][2][L]);
h[9][4][L] = 1.0 / (h[9][1][L] + h[9][2][L]);
//************* IKs CURRENT *****************
XsINF = 1.0 / (1.0 + exp(-(Vm[L] - 1.5) / 16.7));
if (Vm[L] == -30.0)
{
TAUxs = 2058.4;
}
else
{
TAUxs = 0.0000719 * (Vm[L] + 30.0) / (1.0 - exp(-0.148 * (Vm[L] + 30.0)));
TAUxs = TAUxs + 0.000131 * (Vm[L] + 30.0) / ((exp(0.0687 * (Vm[L] + 30.0)) - 1));
TAUxs = 1.0 * 2.0 / TAUxs;
}
m[10][1][L] = XsINF / TAUxs;
m[10][2][L] = 1.0 / TAUxs - m[10][1][L];
m[10][3][L] = m[10][1][L] / (m[10][1][L] + m[10][2][L]);
m[10][4][L] = 1.0 / (m[10][1][L] + m[10][2][L]);
h[10][1][L] = XsINF / TAUxs;
h[10][2][L] = 1.0 / TAUxs - h[10][1][L];
h[10][3][L] = h[10][1][L] / (h[10][1][L] + h[10][2][L]);
h[10][4][L] = 1.0 / (h[10][1][L] + h[10][2][L]);
}
}
void CHeart10Doc::RateConstantsUser()
{
for (int L = 1;L<= cellNo;L++)
{
for (int N = 3;N<=NOP-2;N++)
{
m[N][3][L] = 1 / (1 + exp(-(Vm[L] - RateM[N][1]) / RateM[N][2])); //voltage dependence of gate
m[N][4][L] = (RateM[N][3] - RateM[N][4]) * exp(-pow(Vm[L] - RateM[N][5], 2) / RateM[N][6]) + RateM[N][4]; //time constants of gate
m[N][1][L] = m[N][3][L] / (m[N][4][L] + pow(10,-20));
m[N][2][L] = (1 - m[N][3][L]) / (m[N][4][L] + pow(10,-20));
h[N][3][L] = 1 / (1 + exp(-(Vm[L] - RateH[N][1]) / RateH[N][2])); //voltage dependence of gate
h[N][4][L] = (RateH[N][3] - RateH[N][4]) * exp(-pow(Vm[L] - RateH[N][5], 2) / RateH[N][6]) + RateH[N][4]; //time constants of gate
h[N][1][L] = h[N][3][L] / (h[N][4][L] + pow(10,-20));
h[N][2][L] = (1 - h[N][3][L]) / (h[N][4][L] + pow(10,-20));
}
h[3][3][L] = 1 / (1 + exp(-((Vm[L] - VK[L]) - RateH[3][1]) / (RateH[3][2] + pow(10,-20)))); //'voltage dependence of gate
h[3][4][L] = (RateH[3][3] - RateH[3][4]) * exp(-pow((Vm[L] - VK[L]) - RateH[3][5], 2) / (RateH[3][6] + 1E-20)) + RateH[3][4]; //time constants of gate
h[3][1][L] = Mgi * h[3][3][L] / (h[3][4][L] + pow(10,-20));
h[3][2][L] = (1 - h[3][3][L]) / (h[3][4][L] + pow(10,-20));
h[3][3][L] = h[3][1][L] / (h[3][1][L] + h[3][2][L]);
h[3][4][L] = 1 / (h[3][1][L] + h[3][2][L]);
m[3][3][L] = 0.5 / (1 + exp(-((Vm[L] - VK[L]) - RateM[3][1]) / (RateM[3][2] + pow(10,-20)))); //'voltage dependence of gate
m[3][4][L] = (RateM[3][3] - RateM[3][4]) * exp(-pow((Vm[L]- VK[L]) - RateM[3][5], 2) / (RateM[3][6] + 1E-20)) + RateM[3][4]; //time constants of gate
m[3][1][L] = m[3][3][L] / (m[3][4][L] + pow(10,-20));
m[3][2][L] = (1 - m[3][3][L]) / (m[3][4][L] + pow(10,-20));
hh[4][3][L] = 1 / (1 + exp(-(Vm[L] - RateHH[4][1]) / RateHH[4][2])); //voltage dependence of gate
hh[4][4][L] = (RateHH[4][3] - RateHH[4][4]) * exp(-pow(Vm[L] - RateHH[4][5], 2) / RateHH[4][6]) + RateHH[4][4]; //time constants of gate
hh[4][1][L] = hh[4][3][L] / (hh[4][4][L] + pow(10,-20));
hh[4][2][L] = (1 - hh[4][3][L]) / (hh[4][4][L] + pow(10,-20));
mm[5][3][L] = 1 / (1 + exp(-(Vm[L] - RateMM[5][1]) / RateMM[5][2])); //voltage dependence of gate
mm[5][4][L] = (RateMM[5][3] - RateMM[5][4]) * exp(-pow(Vm[L] - RateMM[5][5], 2) / RateMM[5][6]) + RateMM[5][4]; //time constants of gate
mm[5][1][L] = mm[5][3][L] / (mm[5][4][L] + pow(10,-20));
mm[5][2][L] = (1 - mm[5][3][L]) / (mm[5][4][L] + pow(10,-20));
}
}
void CHeart10Doc::MemoryCurrent()
{
for (int L = 1;L<=cellNo;L++)
{
for (int N = 3;N<=NOP;N++)
{
ii[N][5][L] = ii[N][4][L];
}
for (N = 18;N<=22;N++)
{
ii[N][5][L] = ii[N][4][L];
}
Itotb[L] = Itot[L];
Icapb[L] = Icap[L];
Cafb[L] = Caf[L];
Cafrelb[L] = Cafrel[L];
Itot_downb = Itot_down;
Itot_upb = Itot_up;
}
}
void CHeart10Doc::RUNGEKUTTA()
{
//-------- temporary variables for Runge-Kutta integration ----------------------//
float TempNext_1[11], TempNext_2[11], TempNext_3[11], TempNext_4[11];
float TempRyROpen, TempRyRClose;
for (int L = 1;L<=cellNo;L++)
{
for (int N = 3;N<=NOP;N++) //dy*t
{
m[N][7][L] = m[N][6][L];
h[N][7][L] = h[N][6][L];
m[N][6][L] = m[N][7][L] + m[N][5][L] / (m[N][1][L] + m[N][2][L]) * (1 - exp(-(timeinterval) * (m[N][1][L] + m[N][2][L])));
h[N][6][L] = h[N][7][L] + h[N][5][L] / (h[N][1][L] + h[N][2][L]) * (1 - exp(-(timeinterval) * (h[N][1][L] + h[N][2][L])));
}
mm[5][7][L] = mm[5][6][L];
hh[4][7][L] = hh[4][6][L];
mm[5][6][L] = mm[5][7][L] + mm[5][5][L] / (mm[5][1][L] + mm[5][2][L]) * (1 - exp(-(timeinterval) * (mm[5][1][L] + mm[5][2][L])));
hh[4][6][L] = hh[4][7][L] + hh[4][5][L] / (hh[4][1][L] + hh[4][2][L]) * (1 - exp(-(timeinterval) * (hh[4][1][L] + hh[4][2][L])));
for (N = 1;N<=3;N++) //dy*t
{
Org[N][7][L] = Org[N][6][L];
Org[N][6][L] = Org[N][7][L] + Org[N][5][L] / (Org[N][1][L] + Org[N][2][L]) * (1 - exp(-(timeinterval) * (Org[N][1][L] + Org[N][2][L])));
}
//------------- Runge-Kutta Integration ----------------------------------------//
//--------------1st integration
for (N = 1;N<=5;N++) //dy*t
{
C[N][7][L] = C[N][6][L];
TempNext_1[N] = C[N][5][L] * timeinterval;
C[N][6][L] = C[N][7][L] + 0.5 * TempNext_1[N]; //limit positive number of concentration
}
C[6][7][L] = C[6][6][L];
TempNext_1[6] = C[6][5][L] * timeinterval;
C[6][6][L] = C[6][7][L] + 0.5 * TempNext_1[6];
//---------- RyR release kinetics --------------------//
TempRyROpen = RyROpen[L];
TempRyRClose = RyRClose[L];
TempNext_1[7] = dRyROpen_dt[L] * timeinterval;
RyROpen[L] = TempRyROpen + 0.5 * TempNext_1[7];
TempNext_1[8] = dRyRClose_dt[L] * timeinterval;
RyRClose[L] = TempRyRClose + 0.5 * TempNext_1[8];
differen();
//--------------2nd integration
for (N = 1;N<=5;N++) //dy*t
{
TempNext_2[N] = C[N][5][L] * timeinterval;
C[N][6][L] = C[N][7][L] + 0.5 * TempNext_2[N]; //limit positive number of concentration
}
TempNext_2[6] = C[6][5][L] * timeinterval;
C[6][6][L] = C[6][7][L] + 0.5 * TempNext_2[6];
//---------- RyR release kinetics --------------------//
TempNext_2[7] = dRyROpen_dt[L] * timeinterval;
RyROpen[L] = TempRyROpen + 0.5 * TempNext_2[7];
TempNext_2[8] = dRyRClose_dt[L] * timeinterval;
RyRClose[L] = TempRyRClose + 0.5 * TempNext_2[8];
differen();
//--------------3rd integration
for (N = 1;N<=5;N++) //dy*t
{
TempNext_3[N] = C[N][5][L] * timeinterval;
C[N][6][L] = C[N][7][L] + 0.5 * TempNext_3[N]; //limit positive number of concentration
}
TempNext_3[6] = C[6][5][L] * timeinterval;
C[6][6][L] = C[6][7][L] + 0.5 * TempNext_3[6];
//---------- RyR release kinetics --------------------//
TempNext_3[7] = dRyROpen_dt[L] * timeinterval;
RyROpen[L] = TempRyROpen + 0.5 * TempNext_3[7];
TempNext_3[8] = dRyRClose_dt[L] * timeinterval;
RyRClose[L] = TempRyRClose + 0.5 * TempNext_3[8];
differen();
//--------------4th integration
for (N = 1;N<=5;N++) //dy*t
{
TempNext_4[N] = C[N][5][L] * timeinterval;
C[N][6][L] = C[N][7][L] + (TempNext_1[N] + 2.0*TempNext_2[N] + 2.0*TempNext_3[N] + TempNext_4[N])/6.0; //limit positive number of concentration
}
TempNext_4[6] = C[6][5][L] * timeinterval;
C[6][6][L] = C[6][7][L] + (TempNext_1[6] + 2.0*TempNext_2[6] + 2.0*TempNext_3[6] + TempNext_4[6])/6.0;
Vm[L] = C[6][6][L];
//---------- RyR release kinetics --------------------//
TempNext_4[7] = dRyROpen_dt[L] * timeinterval;
RyROpen[L] = TempRyROpen + (TempNext_1[7] + 2.0*TempNext_2[7] + 2.0*TempNext_3[7] + TempNext_4[7])/6.0;
TempNext_4[8] = dRyRClose_dt[L] * timeinterval;
RyRClose[L] = TempRyRClose + (TempNext_1[8] + 2.0*TempNext_2[8] + 2.0*TempNext_3[8] + TempNext_4[8])/6.0;
}
}
void CHeart10Doc::differen()
{
//----------- Determine dy/dt of Hodgkin-Huxley gating --------------------//
for (int L = 1;L<=cellNo;L++)
{
for (int N = 3;N<=NOP;N++) //NumberOfParameters
{
m[N][5][L] = m[N][1][L] * (1.0 - m[N][6][L]) - m[N][2][L] * m[N][6][L];
h[N][5][L] = h[N][1][L] * (1.0 - h[N][6][L]) - h[N][2][L] * h[N][6][L];
}
mm[5][5][L] = mm[5][1][L] * (1.0 - mm[5][6][L]) - mm[5][2][L] * mm[5][6][L];
hh[4][5][L] = hh[4][1][L] * (1.0 - hh[4][6][L]) - hh[4][2][L] * hh[4][6][L];
m[3][5][L] = m[3][1][L] * (1.0 - m[3][6][L]) - m[3][2][L] * m[3][6][L] * pow((1.0 - h[3][3][L]),3.0);
for (N=1;N<=3;N++) //NumberOfParameters
{
Org[N][5][L] = Org[N][1][L] * (1.0 - Org[N][6][L]) - Org[N][2][L] * Org[N][6][L];
}
/*------------------ RyR and SR Ca2+ kinetics -------------------------------//
RyR channel
(1) C -> O
kco = 280000*[Ca2+]i^2 - 150*ICaL*P(openCaL)
(2) O -> U
kou = 0.08/(1+0.36/[Ca]rel)
(3) U <-> C
kuc = 0.000377*[Ca]rel^2
kcu = 0.000849
I(RyR) = P(RyR)*([Ca]rel - [Ca]i)*P(openRyR)
P(RyR) = 62000
//----------------------------------------------------------------------------*/
dRyROpen_dt[L] = kco[L]*RyRClose[L] - kou[L]*RyROpen[L];
dRyRClose_dt[L] = kuc[L]*(1-RyRClose[L]-RyROpen[L]) - (kco[L]+kcu[L])*RyRClose[L];
C[1][5][L] = -IKnet[L] / Vtotal / F;
C[2][5][L] = -INanet[L] / Vtotal / F;
C[3][5][L] = -ICanet[L] / 2.0 / Vtotal / F; //calculate ion concentration changes, d[I]/dt
C[4][5][L] = (Iup[L] - Itr[L] - Ileak[L]) / 2.0 / Vup / F;
C[5][5][L] = (Itr[L] - Irel[L]) / 2.0 / Vrel / F;
//******** Determine dy/dt of Membrane Potential*******
if (m_AP == 0)
{
C[6][5][L] = -1.0 / Cap * Itot[L];
}
else
{
C[6][5][L] = 1.0 / Cap * ((E - C[6][6][L]) / (0.001*Condition[0]) - Itot[L] + Icap[L]);
}
}
}
void CHeart10Doc::PrepareNextPulse()
{
if (iVariablePulse == 1) //if step pulse is Pulse I
{
Pulse[2] = Pulse[2] + Pulse[13];
Pulse[3] = Pulse[3] + Pulse[13];
Pulse[4] = Pulse[4] + Pulse[14];
Pulse[7] = Pulse[7] - Pulse[14];
}
else if (iVariablePulse == 2) //if step pulse is Pulse I
{
Pulse[5] = Pulse[5] + Pulse[13];
Pulse[6] = Pulse[6] + Pulse[13];
Pulse[7] = Pulse[7] + Pulse[14];
Pulse[10] = Pulse[10] - Pulse[14];
}
else if (iVariablePulse == 3) //if step pulse is Pulse I
{
Pulse[8] = Pulse[8] + Pulse[13];
Pulse[9] = Pulse[9] + Pulse[13];
Pulse[10] = Pulse[10] + Pulse[14];
Pulse[12] = Pulse[12] - Pulse[14];
}
CurrentTime = 0;
dataseparate();
episodeN = episodeN + 1;
}
void CHeart10Doc::SetTimeInterval()
{
if (m_FixedTimeInterval<=0)
{
if (m_AP == 0) //current clamp mode
{
if (fabs(Iinj[1])>0)
{
timeinterval =0.03;
}
else
{
timeinterval = 0.2 / (1.0 + exp(-(((-fabs(TotalTime - 700.0) + fabs(TotalTime + 700.0)) / 2.0) - 500.0) / 200.0)) * (0.1 + 1.9 / (1 + exp(fabs((-fabs(C[6][5][1] - 0.1) + fabs(C[6][5][1] + 0.1)) / 2.0) / 0.12)));
}
}
else //voltage clamp mode
{
timeinterval = 0.2 / (1.0 + exp(-(((-fabs(TotalTime - 2500.0) + fabs(TotalTime + 2500.0)) / 2.0) - 500.0) / 130.0)) * (0.1 + 1.9 / (1 + exp(fabs((-fabs(C[6][5][1] - 0.1) + fabs(C[6][5][1] + 0.1)) / 2.0)/ 0.12)));
}
}
else
{
timeinterval=m_FixedTimeInterval;
}
timeintervalb = timeinterval;
}
void CHeart10Doc::DefInitialVarPreset()
{
for(int L=1;L<=cellNo;L++)
{
//------------------- Initialize intracellular ion concentration ----------------------------//
C[1][6][L] = Ion[3]; //intracellular K concentration
C[2][6][L] = Ion[1]; //intracellular Na concentration
C[3][6][L] = Ion[5]; //cytosolic total Ca concentration
C[4][6][L] = Caupinitial; //Ca concentration in Ca uptake pool
C[5][6][L] = Carelinitial; //Ca concentration of Ca release pool
TRPNf[L] = TRPNt; //FREE CONCENTRATION OF TROPONIN, mM/L necessary to start calculation "buffer"
CMDNf[L] = CMDNt; //FREE CONCENTRATION OF CALMODULIN, mM/L
//------ Initialize diferential value --------------//
int N;
for (N=3;N<=NOP;N++) //dy*t
{
m[N][5][L] = 0;
}
mm[2][5][L] = 0;
hh[1][5][L] = 0;
for (N=1;N<=6;N++) //dy*t
{
C[N][5][L] = 0;
}
for (N=1;N<=3;N++) //dy*t
{
Org[N][5][L] = 0;
}
//------ Initialize NaCa exchanger value --------------//
Org[1][6][L] = 0.9459136; //INACA PXI
//------ Initialize RyR release kinetics ----------------//
dRyROpen_dt[L] = 0;
RyRClose[L] = 1;
RyROpen[L] = 0;
dRyRClose_dt[L] = 0;
// ---- Steady-state values --------------------------//
C[1][6][L] = 133.300308228; //intracellular K concentration
C[2][6][L] = 5.757790565; //intracellular Na concentration
C[3][6][L] = 0.00019059; //cytosolic total Ca concentration
C[4][6][L] = 0.76; //Ca concentration in Ca uptake pool
C[5][6][L] = 5.59; //Ca concentration of Ca release pool
TRPNf[L] = 0.06984663; //FREE CONCENTRATION OF TROPONIN, mM/L necessary to start calculation "buffer"
CMDNf[L] = 0.049964506; //FREE CONCENTRATION OF CALMODULIN, mM/L
Org[1][6][L] = 0.982397199;
RyRClose[L] = 0.000002070;
RyROpen[L] = 0.0975252;
}
//----------------------- Setting for Calculation -----------------------------------//
timeinterval = 0.01; //0.01 msec
timeintervalb = 0.01;
CurrentTime = 0.0; //TIME temporally in msec for Runge kutta
//------------------- Initialize Difference value to calculate Vm by iterative method ------------------//
d[1]=1;
d[2]=0;
//----------------------- Initial Voltage Clamp Setting -----------------------------//
E = Pulse[0];
RTF = R * (Condition[3]+273.15) / F;
RTF2 = R * (Condition[3]+273.15) / F / 2.0;
Ke = Ion[2];
Nae =Ion[0];
Cae =Ion[4];
Mgi =Ion[6];
Erev = RTF * log((Ke + Nae * 0.002 + Cae * 4 * 0.002) / (C[1][6][1] + C[2][6][1]*0.005 + C[3][6][1] * 4 * 0.005));
}
void CHeart10Doc::DefSteps(int iAP)
{
E0 = E;
float StimulusAmplitude;
float TimeElapseWithinPulse; //Elapsed time within pulse
TimeElapseWithinPulse=(CurrentTime / Pulse[18] - int(CurrentTime / Pulse[18]))*Pulse[18];
if ((TimeElapseWithinPulse >= 0.0) && (TimeElapseWithinPulse <= Pulse[1]))
StimulusAmplitude = Pulse[0];
if ((TimeElapseWithinPulse > Pulse[1]) && (TimeElapseWithinPulse <= (Pulse[4] + Pulse[1])))
StimulusAmplitude = (Pulse[3] - Pulse[2]) * (TimeElapseWithinPulse - Pulse[1]) / Pulse[4] + Pulse[2];
if ((TimeElapseWithinPulse > (Pulse[1] + Pulse[4])) && (TimeElapseWithinPulse <= (Pulse[1] + Pulse[4] + Pulse[7])))
StimulusAmplitude = (Pulse[6] - Pulse[5]) * (TimeElapseWithinPulse - Pulse[1] - Pulse[4]) / Pulse[7] + Pulse[5];
if ((TimeElapseWithinPulse > (Pulse[1] + Pulse[4] + Pulse[7])) && (TimeElapseWithinPulse <= (Pulse[1] + Pulse[4] + Pulse[7] + Pulse[10])))
StimulusAmplitude = (Pulse[9] - Pulse[8]) * (TimeElapseWithinPulse - Pulse[1] - Pulse[4] - Pulse[7]) / Pulse[10] + Pulse[8];
if ((TimeElapseWithinPulse > (Pulse[1] + Pulse[4] + Pulse[7] + Pulse[10])) && (TimeElapseWithinPulse <= (Pulse[1] + Pulse[4] + Pulse[7] + Pulse[10] + Pulse[12])))
StimulusAmplitude = Pulse[11];
if (TimeElapseWithinPulse > (Pulse[1] + Pulse[4] + Pulse[7] + Pulse[10] + Pulse[12]))
StimulusAmplitude=Pulse[0];
if (CurrentTime>(Pulse[18]*Pulse[17]))
StimulusAmplitude=Pulse[0];
if (E == E0) {}
else
{
//LabelCurrentTestPulse.Caption = E & " mV";
}
if (iAP==0)
{
Iinj[1]=StimulusAmplitude*1000.0;
}
else
{
E=StimulusAmplitude;
}
}
Click here to get the file