Commit e7a2496b authored by Thomas Kummer's avatar Thomas Kummer
Browse files

Example on how to generate fibers in biventricular geometry

Further coding in Ohara-Rudy model
parent 221aeb20
No preview for this file type
This diff is collapsed.
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
<plist version="1.0">
<dict>
<key>IDESourceControlProjectFavoriteDictionaryKey</key>
<false/>
<key>IDESourceControlProjectIdentifier</key>
<string>1DF34D31-C038-4001-9458-B53DAF11E65A</string>
<key>IDESourceControlProjectName</key>
<string>lifev-em</string>
<key>IDESourceControlProjectOriginsDictionary</key>
<dict>
<key>D38DB74F40EF7AAED3633D413489ADC23CCADBF5</key>
<string>https://gitlab.com/simioftheclouds/LifeV_EM.git</string>
</dict>
<key>IDESourceControlProjectPath</key>
<string>lifev-em.xcworkspace</string>
<key>IDESourceControlProjectRelativeInstallPathDictionary</key>
<dict>
<key>D38DB74F40EF7AAED3633D413489ADC23CCADBF5</key>
<string>..</string>
</dict>
<key>IDESourceControlProjectURL</key>
<string>https://gitlab.com/simioftheclouds/LifeV_EM.git</string>
<key>IDESourceControlProjectVersion</key>
<integer>111</integer>
<key>IDESourceControlProjectWCCIdentifier</key>
<string>D38DB74F40EF7AAED3633D413489ADC23CCADBF5</string>
<key>IDESourceControlProjectWCConfigurations</key>
<array>
<dict>
<key>IDESourceControlRepositoryExtensionIdentifierKey</key>
<string>public.vcs.git</string>
<key>IDESourceControlWCCIdentifierKey</key>
<string>D38DB74F40EF7AAED3633D413489ADC23CCADBF5</string>
<key>IDESourceControlWCCName</key>
<string>lifev-em</string>
</dict>
</array>
</dict>
</plist>
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
<plist version="1.0">
<dict>
<key>BuildLocationStyle</key>
<string>UseAppPreferences</string>
<key>CustomBuildLocationType</key>
<string>RelativeToDerivedData</string>
<key>DerivedDataLocationStyle</key>
<string>Default</string>
<key>IssueFilterStyle</key>
<string>ShowActiveSchemeOnly</string>
<key>LiveSourceIssuesEnabled</key>
<true/>
<key>SnapshotAutomaticallyBeforeSignificantChanges</key>
<true/>
<key>SnapshotLocationStyle</key>
<string>Default</string>
</dict>
</plist>
<?xml version="1.0" encoding="UTF-8"?>
<Bucket
type = "0"
version = "2.0">
</Bucket>
No preview for this file type
......@@ -270,28 +270,54 @@ void IonicOharaRudy::computeRhs ( const std::vector<Real>& v,
Real xs1 = v[35];
Real xs2 = v[35];
Real xk1 = v[37];
Real Jrelnp = v[37];
Real Jrelnp = v[38];
Real Jrelp = v[39];
Real CaMKt = v[40];
revpots(nai, ki);
//V
rhs[0] = - Itot (V, m , h , j , d, f, X, Ca);
//m
rhs[1] = dm (V, m);
//h
rhs[2] = dh (V, h);
//j
rhs[3] = dj (V, j);
//d
rhs[4] = dd (V, d);
//f
rhs[5] = df (V, f);
//X
rhs[6] = dX (V, X);
//Ca
rhs[7] = dCa (V, d, f, Ca);
revpots ( nai, ki );
updateConstants ( nass, nai, kss, ki, cass, cai, Jrelnp, Jrelp );
rhs[0] =
rhs[1] = dnai ();
rhs[2] = dnass ();
rhs[3] = dki ();
rhs[4] = dkss ();
rhs[5] = dcai ( cai );
rhs[6] = dcass ( cass );
rhs[7] = dcansr ();
rhs[8] = dcajsr ();
rhs[9] =
rhs[10] =
rhs[11] =
rhs[12] =
rhs[13] =
rhs[14] =
rhs[15] =
rhs[16] =
rhs[17] =
rhs[18] =
rhs[19] =
rhs[20] =
rhs[21] =
rhs[22] =
rhs[23] =
rhs[24] =
rhs[25] =
rhs[26] =
rhs[27] =
rhs[28] =
rhs[29] =
rhs[30] =
rhs[31] =
rhs[32] =
rhs[33] =
rhs[34] =
rhs[35] =
rhs[36] =
rhs[37] =
rhs[38] = dJrelnp ( Jrelnp, cajsr );
rhs[39] = dJrelp ( Jrelp, cajsr );
rhs[40] = dCaMKt ( CaMKt, cass );
}
void IonicOharaRudy::computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt )
......@@ -368,45 +394,64 @@ Real IonicOharaRudy::dCaMKt ( Real CaMKt,
return (aCaMK*CaMKb*(CaMKb+CaMKt)-bCaMK*CaMKt);
}
// updateConstants
void IonicOharaRudy::updateConstants( Real nass,
Real nai,
Real kss,
Real ki,
Real cass,
Real cai,
Real Jrelnp,
Real Jrelp)
{
M_JdiffNa=(nass-nai)/2.0;
M_JdiffK=(kss-ki)/2.0;
M_Jdiff=(cass-cai)/0.2;
double fJrelp=(1.0/(1.0+M_KmCaMK/M_CaMKa));
M_Jrel=(1.0-fJrelp)*Jrelnp+fJrelp*Jrelp;
double Jupnp=0.004375*cai/(cai+0.00092);
double Jupp=2.75*0.004375*cai/(cai+0.00092-0.00017);
//if (celltype==1) // ?
//{
// Jupnp*=1.3;
// Jupp*=1.3;
//}
double fJupp=(1.0/(1.0+M_KmCaMK/M_CaMKa));
M_Jleak=0.0039375*cansr/15.0;
M_Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-M_Jleak;
M_Jtr=(cansr-cajsr)/100.0;
}
// dJrelnp
Real IonicOharaRudy::dJrelnp ( Real Jrelnp,
Real cajsr )
{
double bt=4.75;
double a_rel=0.5*bt;
double Jrel_inf=a_rel*(-M_ICaL)/(1.0+pow(1.5/cajsr,8.0));
if (celltype==2) // ?
{
Jrel_inf*=1.7;
}
Jrel_inf*=1.7; // ? if celltype ==2 ?
double tau_rel=bt/(1.0+0.0123/cajsr);
if (tau_rel<0.005)
{
tau_rel=0.005;
}
Jrelnp=Jrel_inf-(Jrel_inf-Jrelnp)*exp(-dt/tau_rel); // ?
return (Jrel_inf - Jrelnp) / tau_rel; // Jrelnp=Jrel_inf-(Jrel_inf-Jrelnp)*exp(-dt/tau_rel);
}
// dJrelp
Real IonicOharaRudy::dJrelp ( Real Jrelp,
Real cajsr )
{
double btp=1.25*bt;
double a_relp=0.5*btp;
double Jrel_infp=a_relp*(-M_ICaL)/(1.0+pow(1.5/cajsr,8.0));
if (celltype==2) // ?
{
Jrel_infp*=1.7;
}
Jrel_infp*=1.7; // ? if celltype ==2 ?
double tau_relp=btp/(1.0+0.0123/cajsr);
if (tau_relp<0.005)
......@@ -414,73 +459,451 @@ Real IonicOharaRudy::dCaMKt ( Real CaMKt,
tau_relp=0.005;
}
Jrelp=Jrel_infp-(Jrel_infp-Jrelp)*exp(-dt/tau_relp); // ?
return (Jrel_infp - Jrelp) / tau_relp; //Jrelp=Jrel_infp-(Jrel_infp-Jrelp)*exp(-dt/tau_relp);
}
Real IonicOharaRudy::dnai ()
{
return (-(M_INa+M_INaL+3.0*M_INaCa_i+3.0*M_INaK+M_INab)*M_Acap/(M_F*M_vmyo)+M_JdiffNa*M_vss/M_vmyo);
}
Real IonicOharaRudy::dnass ()
{
return (-(M_ICaNa+3.0*M_INaCa_ss)*M_Acap/(M_F*M_vss)-M_JdiffNa);
}
Real IonicOharaRudy::dki ()
{
return (-(M_Ito+M_IKr+M_IKs+M_IK1+M_IKb+M_Ist-2.0*M_INaK)*M_Acap/(M_F*M_vmyo)+M_JdiffK*M_vss/M_vmyo);
}
Real IonicOharaRudy::dkss ()
{
return (-(M_ICaK)*M_Acap/(M_F*M_vss)-M_JdiffK);
}
// updateConstants
double fJrelp=(1.0/(1.0+KmCaMK/CaMKa));
Jrel=(1.0-fJrelp)*Jrelnp+fJrelp*Jrelp; // ?
Real IonicOharaRudy::dcai ( Real cai )
{
double Bcai;
//if (celltype==1)
//{
// Bcai=1.0/(1.0+1.3*cmdnmax*kmcmdn/pow(kmcmdn+cai,2.0)+trpnmax*kmtrpn/pow(kmtrpn+cai,2.0));
//}
//else
//{
Bcai=1.0/(1.0+M_cmdnmax*M_kmcmdn/pow(M_kmcmdn+cai,2.0)+M_trpnmax*M_kmtrpn/pow(M_kmtrpn+cai,2.0));
//}
return (Bcai*(-(M_IpCa+M_ICab-2.0*M_INaCa_i)*M_Acap/(2.0*M_F*M_vmyo)-M_Jup*M_vnsr/M_vmyo+M_Jdiff*M_vss/M_vmyo));
}
Real IonicOharaRudy::dcass ( Real cass )
{
double Bcass=1.0/(1.0+M_BSRmax*M_KmBSR/pow(M_KmBSR+cass,2.0)+M_BSLmax*M_KmBSL/pow(M_KmBSL+cass,2.0));
return (Bcass*(-(M_ICaL-2.0*M_INaCa_ss)*M_Acap/(2.0*M_F*M_vss)+M_Jrel*M_vjsr/M_vss-M_Jdiff));
}
double Jupnp=0.004375*cai/(cai+0.00092);
double Jupp=2.75*0.004375*cai/(cai+0.00092-0.00017);
if (celltype==1) // ?
{
Jupnp*=1.3;
Jupp*=1.3;
}
Real dcansr IonicOharaRudy::dcansr ()
{
return (M_Jup-M_Jtr*M_vjsr/M_vnsr);
}
double fJupp=(1.0/(1.0+KmCaMK/CaMKa));
M_Jleak=0.0039375*cansr/15.0;
Jup=(1.0-fJupp)*Jupnp+fJupp*Jupp-M_Jleak;
Jtr=(cansr-cajsr)/100.0;
Real dcajsr IonicOharaRudy::dcajsr ( Real cajsr )
{
double Bcajsr=1.0/(1.0+M_csqnmax*M_kmcsqn/pow(M_kmcsqn+cajsr,2.0));
return (Bcajsr*(M_Jtr-M_Jrel));
}
}
// dnai
nai+=dt*(-(INa+INaL+3.0*INaCa_i+3.0*INaK+INab)*Acap/(F*vmyo)+JdiffNa*vss/vmyo);
// dnass
nass+=dt*(-(ICaNa+3.0*INaCa_ss)*Acap/(F*vss)-JdiffNa);
// dki
ki+=dt*(-(Ito+IKr+IKs+IK1+IKb+Ist-2.0*INaK)*Acap/(F*vmyo)+JdiffK*vss/vmyo);
// dkss
kss+=dt*(-(ICaK)*Acap/(F*vss)-JdiffK);
// update (CaMKt, cass)
M_CaMKb=M_CaMKo*(1.0-CaMKt)/(1.0+M_KmCaM/cass);
M_CaMKa=M_CaMKb+CaMKt;
// dcai
double Bcai;
// where to put them?
double vffrt=M_v*M_F*M_F/(M_R*M_T);
double vfrt=V*F/(M_R*M_T);
Real dm ( m )
{
double mss=1.0/(1.0+exp((-(V+39.57))/9.871));
double tm=1.0/(6.765*exp((V+11.64)/34.77)+8.552*exp(-(V+77.42)/5.955));
return ( mss - m ) / tm; //m=mss-(mss-m)*exp(-dt/tm);
}
double hss=1.0/(1+exp((v+82.90)/6.086));
double thf=1.0/(1.432e-5*exp(-(v+1.196)/6.285)+6.149*exp((v+0.5096)/20.27));
double ths=1.0/(0.009794*exp(-(v+17.95)/28.05)+0.3343*exp((v+5.730)/56.66));
double Ahf=0.99;
double Ahs=1.0-Ahf;
hf=hss-(hss-hf)*exp(-dt/thf);
hs=hss-(hss-hs)*exp(-dt/ths);
double h=Ahf*hf+Ahs*hs;
double jss=hss;
double tj=2.038+1.0/(0.02136*exp(-(v+100.6)/8.281)+0.3052*exp((v+0.9941)/38.45));
j=jss-(jss-j)*exp(-dt/tj);
double hssp=1.0/(1+exp((v+89.1)/6.086));
double thsp=3.0*ths;
hsp=hssp-(hssp-hsp)*exp(-dt/thsp);
double hp=Ahf*hf+Ahs*hsp;
double tjp=1.46*tj;
jp=jss-(jss-jp)*exp(-dt/tjp);
double GNa=75;
double fINap=(1.0/(1.0+KmCaMK/CaMKa));
INa=GNa*(v-ENa)*m*m*m*((1.0-fINap)*h*j+fINap*hp*jp);
double mLss=1.0/(1.0+exp((-(v+42.85))/5.264));
double tmL=tm;
mL=mLss-(mLss-mL)*exp(-dt/tmL);
double hLss=1.0/(1.0+exp((v+87.61)/7.488));
double thL=200.0;
hL=hLss-(hLss-hL)*exp(-dt/thL);
double hLssp=1.0/(1.0+exp((v+93.81)/7.488));
double thLp=3.0*thL;
hLp=hLssp-(hLssp-hLp)*exp(-dt/thLp);
double GNaL=0.0075;
if (celltype==1)
{
Bcai=1.0/(1.0+1.3*cmdnmax*kmcmdn/pow(kmcmdn+cai,2.0)+trpnmax*kmtrpn/pow(kmtrpn+cai,2.0));
GNaL*=0.6;
}
double fINaLp=(1.0/(1.0+KmCaMK/CaMKa));
INaL=GNaL*(v-ENa)*mL*((1.0-fINaLp)*hL+fINaLp*hLp);
double ass=1.0/(1.0+exp((-(v-14.34))/14.82));
double ta=1.0515/(1.0/(1.2089*(1.0+exp(-(v-18.4099)/29.3814)))+3.5/(1.0+exp((v+100.0)/29.3814)));
a=ass-(ass-a)*exp(-dt/ta);
double iss=1.0/(1.0+exp((v+43.94)/5.711));
double delta_epi;
if (celltype==1)
{
delta_epi=1.0-(0.95/(1.0+exp((v+70.0)/5.0)));
}
else
{
Bcai=1.0/(1.0+cmdnmax*kmcmdn/pow(kmcmdn+cai,2.0)+trpnmax*kmtrpn/pow(kmtrpn+cai,2.0));
delta_epi=1.0;
}
cai+=dt*(Bcai*(-(IpCa+ICab-2.0*INaCa_i)*Acap/(2.0*F*vmyo)-Jup*vnsr/vmyo+Jdiff*vss/vmyo));
// dcass
Real dcass
double Bcass=1.0/(1.0+M_BSRmax*M_KmBSR/pow(M_KmBSR+cass,2.0)+BSLmax*KmBSL/pow(KmBSL+cass,2.0));
cass+=dt*(Bcass*(-(ICaL-2.0*INaCa_ss)*Acap/(2.0*F*vss)+Jrel*vjsr/vss-Jdiff));
double tiF=4.562+1/(0.3933*exp((-(v+100.0))/100.0)+0.08004*exp((v+50.0)/16.59));
double tiS=23.62+1/(0.001416*exp((-(v+96.52))/59.05)+1.780e-8*exp((v+114.1)/8.079));
tiF*=delta_epi;
tiS*=delta_epi;
double AiF=1.0/(1.0+exp((v-213.6)/151.2));
double AiS=1.0-AiF;
iF=iss-(iss-iF)*exp(-dt/tiF);
iS=iss-(iss-iS)*exp(-dt/tiS);
double i=AiF*iF+AiS*iS;
double assp=1.0/(1.0+exp((-(v-24.34))/14.82));
ap=assp-(assp-ap)*exp(-dt/ta);
double dti_develop=1.354+1.0e-4/(exp((v-167.4)/15.89)+exp(-(v-12.23)/0.2154));
double dti_recover=1.0-0.5/(1.0+exp((v+70.0)/20.0));
double tiFp=dti_develop*dti_recover*tiF;
double tiSp=dti_develop*dti_recover*tiS;
iFp=iss-(iss-iFp)*exp(-dt/tiFp);
iSp=iss-(iss-iSp)*exp(-dt/tiSp);
double ip=AiF*iFp+AiS*iSp;
double Gto=0.02;
if (celltype==1)
{
Gto*=4.0;
}
if (celltype==2)
{
Gto*=4.0;
}
double fItop=(1.0/(1.0+KmCaMK/CaMKa));
Ito=Gto*(v-EK)*((1.0-fItop)*a*i+fItop*ap*ip);
double dss=1.0/(1.0+exp((-(v+3.940))/4.230));
double td=0.6+1.0/(exp(-0.05*(v+6.0))+exp(0.09*(v+14.0)));
d=dss-(dss-d)*exp(-dt/td);
double fss=1.0/(1.0+exp((v+19.58)/3.696));
double tff=7.0+1.0/(0.0045*exp(-(v+20.0)/10.0)+0.0045*exp((v+20.0)/10.0));
double tfs=1000.0+1.0/(0.000035*exp(-(v+5.0)/4.0)+0.000035*exp((v+5.0)/6.0));
double Aff=0.6;
double Afs=1.0-Aff;
ff=fss-(fss-ff)*exp(-dt/tff);
fs=fss-(fss-fs)*exp(-dt/tfs);
double f=Aff*ff+Afs*fs;
double fcass=fss;
double tfcaf=7.0+1.0/(0.04*exp(-(v-4.0)/7.0)+0.04*exp((v-4.0)/7.0));
double tfcas=100.0+1.0/(0.00012*exp(-v/3.0)+0.00012*exp(v/7.0));
double Afcaf=0.3+0.6/(1.0+exp((v-10.0)/10.0));
double Afcas=1.0-Afcaf;
fcaf=fcass-(fcass-fcaf)*exp(-dt/tfcaf);
fcas=fcass-(fcass-fcas)*exp(-dt/tfcas);
double fca=Afcaf*fcaf+Afcas*fcas;
double tjca=75.0;
jca=fcass-(fcass-jca)*exp(-dt/tjca);
double tffp=2.5*tff;
ffp=fss-(fss-ffp)*exp(-dt/tffp);
double fp=Aff*ffp+Afs*fs;
double tfcafp=2.5*tfcaf;
fcafp=fcass-(fcass-fcafp)*exp(-dt/tfcafp);
double fcap=Afcaf*fcafp+Afcas*fcas;
double Kmn=0.002;
double k2n=1000.0;
double km2n=jca*1.0;
double anca=1.0/(k2n/km2n+pow(1.0+Kmn/cass,4.0));
nca=anca*k2n/km2n-(anca*k2n/km2n-nca)*exp(-km2n*dt);
double PhiCaL=4.0*vffrt*(cass*exp(2.0*vfrt)-0.341*cao)/(exp(2.0*vfrt)-1.0);
double PhiCaNa=1.0*vffrt*(0.75*nass*exp(1.0*vfrt)-0.75*nao)/(exp(1.0*vfrt)-1.0);
double PhiCaK=1.0*vffrt*(0.75*kss*exp(1.0*vfrt)-0.75*ko)/(exp(1.0*vfrt)-1.0);
double zca=2.0;
double PCa=0.0001;
if (celltype==1)
{
PCa*=1.2;
}
if (celltype==2)
{
PCa*=2.5;
}
double PCap=1.1*PCa;
double PCaNa=0.00125*PCa;
double PCaK=3.574e-4*PCa;
double PCaNap=0.00125*PCap;
double PCaKp=3.574e-4*PCap;
double fICaLp=(1.0/(1.0+KmCaMK/CaMKa));
ICaL=(1.0-fICaLp)*PCa*PhiCaL*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCap*PhiCaL*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaNa=(1.0-fICaLp)*PCaNa*PhiCaNa*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaNap*PhiCaNa*d*(fp*(1.0-nca)+jca*fcap*nca);
ICaK=(1.0-fICaLp)*PCaK*PhiCaK*d*(f*(1.0-nca)+jca*fca*nca)+fICaLp*PCaKp*PhiCaK*d*(fp*(1.0-nca)+jca*fcap*nca);
double xrss=1.0/(1.0+exp((-(v+8.337))/6.789));
double txrf=12.98+1.0/(0.3652*exp((v-31.66)/3.869)+4.123e-5*exp((-(v-47.78))/20.38));
double txrs=1.865+1.0/(0.06629*exp((v-34.70)/7.355)+1.128e-5*exp((-(v-29.74))/25.94));
double Axrf=1.0/(1.0+exp((v+54.81)/38.21));
double Axrs=1.0-Axrf;
xrf=xrss-(xrss-xrf)*exp(-dt/txrf);
xrs=xrss-(xrss-xrs)*exp(-dt/txrs);
double xr=Axrf*xrf+Axrs*xrs;
double rkr=1.0/(1.0+exp((v+55.0)/75.0))*1.0/(1.0+exp((v-10.0)/30.0));
double GKr=0.046;
if (celltype==1)
{
GKr*=1.3;
}
if (celltype==2)
{
GKr*=0.8;
}
IKr=GKr*sqrt(ko/5.4)*xr*rkr*(v-EK);
double xs1ss=1.0/(1.0+exp((-(v+11.60))/8.932));
double txs1=817.3+1.0/(2.326e-4*exp((v+48.28)/17.80)+0.001292*exp((-(v+210.0))/230.0));
xs1=xs1ss-(xs1ss-xs1)*exp(-dt/txs1);
double xs2ss=xs1ss;
double txs2=1.0/(0.01*exp((v-50.0)/20.0)+0.0193*exp((-(v+66.54))/31.0));
xs2=xs2ss-(xs2ss-xs2)*exp(-dt/txs2);
double KsCa=1.0+0.6/(1.0+pow(3.8e-5/cai,1.4));
double GKs=0.0034;
if (celltype==1)
{
GKs*=1.4;
}
IKs=GKs*KsCa*xs1*xs2*(v-EKs);
double xk1ss=1.0/(1.0+exp(-(v+2.5538*ko+144.59)/(1.5692*ko+3.8115)));
double txk1=122.2/(exp((-(v+127.2))/20.36)+exp((v+236.8)/69.33));
xk1=xk1ss-(xk1ss-xk1)*exp(-dt/txk1);
double rk1=1.0/(1.0+exp((v+105.8-2.6*ko)/9.493));
double GK1=0.1908;
if (celltype==1)
{
GK1*=1.2;
}
if (celltype==2)
{
GK1*=1.3;
}
IK1=GK1*sqrt(ko)*rk1*xk1*(v-EK);
double kna1=15.0;
double kna2=5.0;
double kna3=88.12;