- Author:
- WeiweiAi <wai484@aucklanduni.ac.nz>
- Date:
- 2022-10-19 10:24:37+13:00
- Desc:
- Add readme to link the documentation on overleaf.
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/701/rawfile/6741283829975cd4f084159e7734dcbece29f7b1/Components/uscm.cpp
/*
Uterine Single Cell Model
Author: Craig Testrow
Group: Biological Physics, The University of Manchester
Supervisor: Prof. Henggui Zhang
161214
To compile:
g++ -lm -Wall uscm.cpp -o uscm
To run the compiled program:
./uscm
*/
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <fstream>
using namespace std;
//Global constants and associated variables
const double R = 8314; // Molar gas constant
const double Temperature = 308.; // Temp in Kelvin (35C)
const double Faraday = 96485.; // Faraday's constant mC mM-1 // PRETEST
const double Cell_Volume = 2.65e-8; // cm3
const double Cell_Area = 7.422e-5; // cm2
const double Membrane_Capacitance = 0.0014; // mF cm-2 to give a total capacitance of about 100 pF
const double zCa = 2.0; // Calcium valence
const double zCl = -1.0; // Chlorine valence
const double zNa = 1.0; // Sodium valence
const double zK = 1.0; // Potassium valence
const double dt = 0.01; // time step (10microsecs)
const int num_cells = 1; // total number of cells
int cur_cell = 0;
const double PK(1.), PNa(0.35); // Permeabilities (rates/ease of flow) of K and Na ions
//Extracellular ion concentrations
double Ko(5.4);
double Nao(140.);
double Cao(2.5);
double Clo(146.4);
double Mgo(0.5);
//Intracellular ion concentrations
double Cai[num_cells], Nai[num_cells], Ki[num_cells], Cli[num_cells];
//SR ion concentrations
double Ca_u[num_cells], Ca_r[num_cells];
//Membrane voltage
double V[num_cells];
//Gating variables
double Ih_y[num_cells], ICaL_d[num_cells], ICaL_f1[num_cells], ICaL_f2[num_cells], ICaT_b[num_cells], ICaT_g[num_cells], INa_m[num_cells], INa_h[num_cells];
double IK1_q[num_cells], IK1_r1[num_cells], IK1_r2[num_cells], IK2_p[num_cells], IK2_k1[num_cells], IK2_k2[num_cells], IKA_s[num_cells], IKA_x[num_cells];
double V_KCa_half[num_cells], IKCa_pf[num_cells], IKCa_ps[num_cells];
double ICl_c1[num_cells], ICl_c2[num_cells], ICl_n[num_cells];
//Ryanodine receptors
double R00[num_cells],R10[num_cells],R01[num_cells],R11[num_cells];
//Current totals
double Itot;
double ICa_tot, INa_tot, IK_tot, ICl_tot;
//Other parameters
double Isus = 0.6; //sustained background current
double Istim = 0, Iinitial = 0; //reinitialised in the main program if used
double InsCa, InsNa, InsK, Insleak; // Ins currents that need to be accessed outside of function
double IsocCa, IsocNa; // Isoc currents that need to be accessed outside of function
double Iup, Itr, Irel; // SR currents that need to be accessed outside of function
//Steady states (declared here to be accessable to class_vclamp and output)
double Ih_y_steady_state;
double ICaL_d_steady_state, ICaL_f_steady_state;
double ICaT_b_steady_state, ICaT_g_steady_state;
double INa_m_steady_state, INa_h_steady_state;
double IK1_q_steady_state, IK1_r_steady_state;
double IK2_p_steady_state, IK2_k_steady_state;
double IKA_s_steady_state, IKA_x_steady_state;
double IKCa_p_steady_state;
double ICl_c_steady_state,ICl_n_steady_state;
double Inscc_fMg_steady_state;
//Time constants (declared here to be accessable to class_vclamp and output)
double Ih_y_Tau;
double ICaL_d_Tau, ICaL_f1_Tau, ICaL_f2_Tau;
double ICaT_b_Tau, ICaT_g_Tau;
double INa_m_Tau, INa_h_Tau;
double IK1_q_Tau, IK1_r1_Tau, IK1_r2_Tau;
double IK2_p_Tau, IK2_k1_Tau, IK2_k2_Tau;
double IKA_s_Tau, IKA_x_Tau;
double IKCa_pf_Tau, IKCa_ps_Tau;
double ICl_c1_Tau, ICl_c2_Tau,ICl_n_Tau;
//Kinetic cross-bridge cycle
double M, Mp, AMp, AM;
double Frac_phos, Frac_stress;
//Force
double Fp, Fx, Fa, Fs, Ft;
double lc, ls, la, lx; //lengths of cell and force elements
/* Current magnitudes */
double Ih_mag(0.75);
double ICaL_mag(0.530*1.0);
double ICaT_mag(1.333*1.0);
double INa_mag(0.96*1);//0.75);
double Isoc_mag(1);
double IsocCa_mag(0.1);
double IsocNa_mag(1);
double IK1_mag(0.9*1.0);
double IK2_mag(0.4*1.0);
double IKA_mag(0.0075*1);
double IKCa_mag(10*1.0);
double IKleak_mag(0.08);
double Inscc_mag(1*1.0);
double InsCa_mag(1);
double InsNa_mag(1*1.0);
double InsK_mag(1);
double Insleak_mag(0.8*1);
double ICl_mag(1);
double ICapump_mag(5e0);
double INaK_mag; // see below in time loop
double INaCa_mag(5e0);
double INaKCl_mag(0.02);
double Iup_mag(1);
double Irel_mag(1);
double Itr_mag(1);
// Cell class
class cell
{
public:
// Constructor to initialise all current variables, etc
cell(void)
{
//Initial conditions for parameters
for (int i=0;i<num_cells;i++)
{
//Initialise gating variables
Ih_y[i]=0.00223178;
ICaL_d[i]=0.0125799;
ICaL_f1[i]=0.942219;
ICaL_f2[i]=0.942218;
ICaT_b[i]=0.169361;
ICaT_g[i]=0.0318277;
INa_m[i]=0.124663;
INa_h[i]=0.371776;
IK1_q[i]=0.07298;
IK1_r1[i]=0.159599;
IK1_r2[i]=0.129642;
IK2_p[i]=0.0782046;
IK2_k1[i]=0.995883;
IK2_k2[i]=0.95664;
IKA_s[i]=0.0366284;
IKA_x[i]=0.653819;
V_KCa_half[i]=-41.7 * log10(Cai[i]) - 128.2;
IKCa_pf[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
IKCa_ps[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
ICl_c1[i]=0.000864388;
ICl_c2[i]=0.000864392;
ICl_n[i]=0.445644;
//Initialise each element of concentration arrays
Cai[i]=8.97e-5;
Nai[i]=4.;
Ki[i]=140;
Cli[i]=130;
Ca_u[i]=0.0273854;
Ca_r[i]=0.0232185;
//Initialise each element of membrane voltage array
V[i]=-52.541;
//Initialise ryanodine receptors
R00[i]=6.99E-09;//0.25;
R01[i]=7.85E-12;//.25;
R10[i]=1.88E-11;//.25;
R11[i]=2.11E-14;//.25;
/* //Initialise gating variables
Ih_y[i]=0.00212991;
ICaL_d[i]=0.0133195;
ICaL_f1[i]=0.938986;
ICaL_f2[i]=0.938985;
ICaT_b[i]=0.179073;
ICaT_g[i]=0.0308496;
INa_m[i]=0.129659;
INa_h[i]=0.360711;
IK1_q[i]=0.0741461;
IK1_r1[i]=0.151067;
IK1_r2[i]=0.108987;
IK2_p[i]=0.0794842;
IK2_k1[i]=0.995353;
IK2_k2[i]=0.954862;
IKA_s[i]=0.0385646;
IKA_x[i]=0.644685;
V_KCa_half[i]=-41.7 * log10(Cai[i]) - 128.2;
IKCa_pf[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
IKCa_ps[i]=1/(1+exp(-(V[i]-V_KCa_half[i])/18.25));
ICl_c1[i]=0.00089031;
ICl_c2[i]=0.000890317;
ICl_n[i]=0.440717;
//Initialise each element of concentration arrays
Cai[i]=9.43e-5;
Nai[i]=4.;
Ki[i]=140;
Cli[i]=130;
Ca_u[i]=2.64648;
Ca_r[i]=2.24248;
//Initialise each element of membrane voltage array
V[i]=-52.1358;
//Initialise ryanodine receptors
R00[i]=1;//0.25;
R01[i]=0;//.25;
R10[i]=0;//.25;
R11[i]=0;//.25;
*/
//Initialise force elements
//These do not need to be arrays as they're never used beyond single cell
M = 0.97;
Mp = 0.01;
AMp = 0.01;
AM = 1-(M+Mp+AMp);
ls = 30.08;
la = 25.3;
}
}
//Current functions
/* ##############################################
* # #
* # Hyperpolarisation-activated Current - Ih #
* # #
* ############################################## */
double Ih ()
{
double gh = 0.0542;
double Eh = (R * Temperature / Faraday) * log((Ko + (PNa / PK) * Nao) / (Ki[cur_cell] + (PNa / PK) * Nai[cur_cell]));
// Activation
Ih_y_steady_state = 1/(1+exp((V[cur_cell]+105.39)/8.66));
if (V[cur_cell]>160)
{Ih_y_Tau = 1/(3.5e-6*(1-0.0487*V[cur_cell]+pow(-0.0487*V[cur_cell],2)/2
+pow(-0.0487*V[cur_cell],3)/6)+0.04*(1+0.0521*V[cur_cell]+pow(0.0521*V[cur_cell],2)/2+pow(0.0521*V[cur_cell],3)/6));}
else
{Ih_y_Tau = 1/(3.5e-6*exp(-0.0497*V[cur_cell])+0.04*exp(0.0521*V[cur_cell]));}
Ih_y[cur_cell] = Ih_y[cur_cell] + dt * ((Ih_y_steady_state - Ih_y[cur_cell]) / Ih_y_Tau);
return (gh * pow(Ih_y[cur_cell],1) * (V[cur_cell] - Eh));
}
/* ###################################
* # #
* # L-Type Calcium Current - ICaL #
* # #
* ################################### */
double ICaL (double nif)
{
double gCaL(0.6/* *0.4*/);//reduced 60% for estradiol
double ECaL = 45;
const double KmCaL(0.0006);
// Activation
ICaL_d_steady_state = 1/(1+exp(-(V[cur_cell]+22)/7));
ICaL_d_Tau = 2.29+5.7/(1+pow((V[cur_cell]+29.97)/9,2));
ICaL_d[cur_cell] = ICaL_d[cur_cell] + dt * ((ICaL_d_steady_state - ICaL_d[cur_cell]) / ICaL_d_Tau);
// Inactivation
ICaL_f_steady_state = 1/(1+exp((V[cur_cell]+33)/7));
ICaL_f1_Tau = 12;
ICaL_f2_Tau = 90.97-(90.97/((1+exp((V[cur_cell]+13.96)/45.38))*(1+exp(-(V[cur_cell]+9.5)/3.39))));
ICaL_f1[cur_cell] = ICaL_f1[cur_cell] + dt * ((ICaL_f_steady_state - ICaL_f1[cur_cell]) / ICaL_f1_Tau);
ICaL_f2[cur_cell] = ICaL_f2[cur_cell] + dt * ((ICaL_f_steady_state - ICaL_f2[cur_cell]) / ICaL_f2_Tau);
const double IC50_nif =1e-6;
double Hill=1/(1+pow((nif/IC50_nif),1));
gCaL = gCaL*Hill;
return (gCaL * pow(ICaL_d[cur_cell],2.0) * (1 / (1 + pow((Cai[cur_cell]/KmCaL),1))) * pow((0.8 * ICaL_f1[cur_cell] + 0.2 * ICaL_f2[cur_cell]),1))*(V[cur_cell] - ECaL);
}
/* ###################################
* # #
* # T-Type Calcium Current - ICaT #
# #
################################### */
double ICaT ()
{
double gCaT(0.0174/* *0.4*/);//reduced 60% for estradiol
double ECaT = 42;
// Activation
ICaT_b_steady_state = 1/(1+exp(-(V[cur_cell]+43.)/6.)); //best fit to data 14xxxx. This provides closest IV curve, despite being odd fit to activation
ICaT_b_Tau = (0.45+ 3.9/(1+((V[cur_cell]+66)/26)*((V[cur_cell]+66)/26)));
ICaT_b[cur_cell] = ICaT_b[cur_cell] + dt * ((ICaT_b_steady_state - ICaT_b[cur_cell]) / ICaT_b_Tau);
// Inactivation
ICaT_g_steady_state = 0.02+0.98/(1+exp((V[cur_cell]+72.98)/4.64));
ICaT_g_Tau = 150 - 150/((1+exp((V[cur_cell]-417.43)/203.18))*(1+exp(-(V[cur_cell]+61.11)/8.07)));
ICaT_g[cur_cell] = ICaT_g[cur_cell] + dt * ((ICaT_g_steady_state - ICaT_g[cur_cell]) / ICaT_g_Tau);
return (gCaT * pow(ICaT_b[cur_cell],2) * (pow(ICaT_g[cur_cell],1)) * (V[cur_cell] - ECaT));
}
/* ###############################
* # #
* # Fast Sodium Current - INa #
# #
############################### */
double INa ()
{
double gNa(0.12);
double ENa = (R * Temperature / Faraday) * log(Nao / Nai[cur_cell]);
// Activation
INa_m_steady_state = 1/(1+exp(-(V[cur_cell]+35.)/9));
INa_m_Tau = 0.25+7/(1+exp((V[cur_cell]+38)/10));
INa_m[cur_cell] = INa_m[cur_cell] + dt * ((INa_m_steady_state - INa_m[cur_cell]) / INa_m_Tau);
// Inactivation
INa_h_steady_state = 1/(1+exp((V[cur_cell]+57)/8.5));
INa_h_Tau = 0.9+1002.85/(1+((V[cur_cell]+47.5)/1.5)*((V[cur_cell]+47.5)/1.5));
INa_h[cur_cell] = INa_h[cur_cell] + dt * ((INa_h_steady_state - INa_h[cur_cell]) / INa_h_Tau);
return (gNa * pow(INa_m[cur_cell],3) * pow(INa_h[cur_cell],1) * (V[cur_cell] - ENa));
}
/* ########################################################
* # #
# Store-Operated Non-Selective Cation Current - Isoc #
# #
########################################################*/
double Isoc ()
{
double ECaL = 45;
double ECaT = 42;
double ECa = ECaL + ECaT;
double ENa = (R * Temperature / Faraday) * log(Nao / Nai[cur_cell]);
double gsocCa(8.3e-5);
double gsocNa(5.75e-4);
const double Kmsoc(0.0001);
double Psoc = 1 / (1 + Cai[cur_cell] / Kmsoc);
IsocNa = gsocNa * Psoc * (V[cur_cell] - ENa);
IsocCa = gsocCa * Psoc * (V[cur_cell] - ECa);
return (IsocCa + IsocNa);
}
/* ################################################
# #
# Delayed Rectifying Potassium Current - IK1 #
# #
################################################ */
double IK1 ()
{
double gK1(0.2275/* *0.7*/);//reduced 30% for estradiol
double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
// Activation
IK1_q_steady_state = 1/(1+exp(-(V[cur_cell]-7.7)/23.7));
IK1_q_Tau = (500/(1+((V[cur_cell]+60.71)/15.79)*((V[cur_cell]+60.71)/15.79)));
IK1_q[cur_cell] = IK1_q[cur_cell] + dt * ((IK1_q_steady_state - IK1_q[cur_cell]) / IK1_q_Tau);
// Inactivation
IK1_r_steady_state = 1/(1+exp((V[cur_cell]+63)/6.3));
IK1_r1_Tau = (5e3/(1+((V[cur_cell]+62.71)/35.86)*((V[cur_cell]+62.71)/35.86)));
IK1_r2_Tau = (3e4+2.2e5/(1+exp((V[cur_cell]+22)/4)));
IK1_r1[cur_cell] = IK1_r1[cur_cell] + dt * ((IK1_r_steady_state - IK1_r1[cur_cell]) / IK1_r1_Tau);
IK1_r2[cur_cell] = IK1_r2[cur_cell] + dt * ((IK1_r_steady_state - IK1_r2[cur_cell]) / IK1_r2_Tau);
return (gK1 * pow(IK1_q[cur_cell],1) * pow((0.38 * IK1_r1[cur_cell] + 0.62 * IK1_r2[cur_cell]),1) * (V[cur_cell] - EK));
}
/* ################################################
* # #
* # Delayed Rectifying Potassium Current - IK2 #
# #
################################################ */
double IK2 ()
{
double gK2(0.07875/* *0.5*/);//reduced 50% for estradiol
double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
// Activation
IK2_p_steady_state = 1/(1+exp(-(V[cur_cell]-4.2)/23));
IK2_p_Tau = 100/(1+((V[cur_cell]+64.1)/28.67)*((V[cur_cell]+64.1)/28.67));
IK2_p[cur_cell] = IK2_p[cur_cell] + dt * ((IK2_p_steady_state - IK2_p[cur_cell]) / IK2_p_Tau);
// Inactivation
IK2_k_steady_state = 1/(1+exp((V[cur_cell]+21.2)/5.7));
IK2_k1_Tau = 1e6-1e6/((1+exp((V[cur_cell]-315)/50))*(1+exp(-(V[cur_cell]+74.9)/8)));
IK2_k2_Tau = 2.5e6-2.5e6/((1+exp((V[cur_cell]-132.87)/25.4))*(1+exp(-(V[cur_cell]+24.92)/2.68)));
IK2_k1[cur_cell] = IK2_k1[cur_cell] + dt * ((IK2_k_steady_state - IK2_k1[cur_cell]) / IK2_k1_Tau);
IK2_k2[cur_cell] = IK2_k2[cur_cell] + dt * ((IK2_k_steady_state - IK2_k2[cur_cell]) / IK2_k2_Tau);
return (gK2 * IK2_p[cur_cell] * (0.75 * IK2_k1[cur_cell] + 0.25 * IK2_k2[cur_cell]) * (V[cur_cell] - EK));
}
/* ##############################################
* # #
* # A-Type Transient Potassium Current - IKA #
# #
############################################## */
double IKA ()
{
double gKA(0.189/* *0.4*/);//reduced 60% for estradiol
double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
// Activation
IKA_s_steady_state = 1/(1+exp(-(V[cur_cell]+27.79)/7.57));
IKA_s_Tau = 17/(1+((V[cur_cell]+20.52)/35)*((V[cur_cell]+20.52)/35));
IKA_s[cur_cell] = IKA_s[cur_cell] + dt * ((IKA_s_steady_state - IKA_s[cur_cell]) / IKA_s_Tau);
// Inactivation
IKA_x_steady_state = 0.3293+0.71165/(1+exp((V[cur_cell]+53.919)/7.8111));
IKA_x_Tau = 7.5+10/(1+pow((V[cur_cell]+34.18)/120,2));
IKA_x[cur_cell] = IKA_x[cur_cell] + dt * ((IKA_x_steady_state - IKA_x[cur_cell]) / IKA_x_Tau);
return (gKA * pow(IKA_s[cur_cell],1) * exp(pow(IKA_x[cur_cell],1)) * (V[cur_cell] - EK));
}
/* ################################################
* # #
* # Calcium-Activated Potassium Current - IKCa #
* # #
* ################################################ */
double IKCa ()
{
double PBKCa = 3.9e-13;
double NBKCa = 5.94e6;
double iKCa;
// Activation
V_KCa_half[cur_cell]=-41.7 * log10(Cai[cur_cell]) - 128.2;
IKCa_p_steady_state = 1 / (1 + exp(-(V[cur_cell] - V_KCa_half[cur_cell])/18.25));
IKCa_pf_Tau = 0.84;
IKCa_ps_Tau = 35.9;
// Need an approximation to avoid /zero at higher voltages. Based on taylor expansion of exponential
if(V[cur_cell]>160)
{
iKCa = 1e4 * PBKCa * V[cur_cell] * (Faraday * Faraday / (R * Temperature)) *
((Ko - Ki[cur_cell] * (1+(V[cur_cell] * Faraday / (R * Temperature))+pow(V[cur_cell] * Faraday / (R * Temperature),2)/2+pow(V[cur_cell] * Faraday / (R * Temperature),3)/6)) /
(1-(1+(V[cur_cell] * Faraday / (R * Temperature))+pow(V[cur_cell] * Faraday / (R * Temperature),2)/2+pow(V[cur_cell] * Faraday / (R * Temperature),3)/6)));
}
else
{
iKCa = 1e4 * PBKCa * V[cur_cell] * (Faraday * Faraday / (R * Temperature)) * ((Ko - Ki[cur_cell] * exp(V[cur_cell] * Faraday / (R * Temperature))) /
(1 - exp(V[cur_cell] * Faraday / (R * Temperature))));
}
IKCa_pf[cur_cell] = IKCa_pf[cur_cell] + dt * ((IKCa_p_steady_state - IKCa_pf[cur_cell]) / IKCa_pf_Tau);
IKCa_ps[cur_cell] = IKCa_ps[cur_cell] + dt * ((IKCa_p_steady_state - IKCa_ps[cur_cell]) / IKCa_ps_Tau);
double PKCa = 0.17 * IKCa_pf[cur_cell] + 0.83 * IKCa_ps[cur_cell];
return (Cell_Area * NBKCa * PKCa * iKCa / 100);
}
/* ###################################################
* # #
* # Background Potassium Leakage Current - IKleak #
* # #
* ###################################################*/
double IKleak ()
{
double gKleak(0.005);
double EK = (R * Temperature / Faraday) * log(Ko / Ki[cur_cell]);
return (gKleak * (V[cur_cell] - EK));
}
/* ##########################################
* # #
* # Non-Selective Cation Current - Inscc #
* # #
* ########################################## */
double Inscc ()
{
const double gns(0.0123), g_leak(0.009685);
const double KdMg(0.28);
double PNa_PCs(0.9), PK_PCs(1.3), PCa_PCs(0.89); // Permeability ratios for Inscc
double Ens = (R * Temperature / Faraday) *
log((PNa_PCs * Nao + PK_PCs * Ko + 4 * PCa_PCs * Cao * (Cao / (1 + exp((V[cur_cell] * Faraday) / (R * Temperature))))) / //should there be a permeability of calcium
(PNa_PCs * Nai[cur_cell] + PK_PCs * Ki[cur_cell] + 4 * PCa_PCs * Cai[cur_cell] * (Cai[cur_cell] / (1 + exp((V[cur_cell] * Faraday) / (R * Temperature)))))); //factor in here near PCa_PCs? if so, what is it?
Inscc_fMg_steady_state = 0.1 + 0.9 / ((1 + pow((Mgo / KdMg),1.3)));
double gnscc_Ca = (0.03 / (1 + pow((150 / (Cao + 1e-8)),2)))/0.000525;
double gnscc_Na = (0.03 / (1 + pow((150 / (Nao + 1e-8)),2)))/0.0123;
double gnscc_K = (0.03 / (1 + pow((150 / (Ko + 1e-8)),2)))/0.0123;
InsCa = InsCa_mag*0.15*Inscc_fMg_steady_state * gnscc_Ca * gns * (V[cur_cell]-Ens);
InsNa = InsNa_mag*Inscc_fMg_steady_state * gnscc_Na * gns * (V[cur_cell]-Ens);
InsK = InsK_mag*1.19*Inscc_fMg_steady_state * gnscc_K * gns * (V[cur_cell]-Ens);
Insleak = Insleak_mag*g_leak * (V[cur_cell] - Ens);
return (InsCa + InsNa + InsK + Insleak);
}
/* ##############################################
* # #
* # Calcium-Activated Chloride Current - ICl #
* # #
* ############################################## */
double ICl ()
{
double gCl(0.1875);
double ECl_pos = -30, ECl_neg = -5;
double KdCl(256.98e-3);
double gCl_value;
// Activation
ICl_c_steady_state = 1/(1+exp(-(V[cur_cell]-44.08)/13.7));
ICl_c1_Tau = 24.-16.88/(1+exp(-(V[cur_cell]-48.14)/13.));
ICl_c2_Tau = 32.+81.54/(1+exp(-(V[cur_cell]+43.51)/39.5));
ICl_c1[cur_cell] = ICl_c1[cur_cell] + dt * ((ICl_c_steady_state - ICl_c1[cur_cell]) / ICl_c1_Tau);
ICl_c2[cur_cell] = ICl_c2[cur_cell] + dt * ((ICl_c_steady_state - ICl_c2[cur_cell]) / ICl_c2_Tau);
// Inactivation
ICl_n_steady_state = 1 / (1 + exp((V[cur_cell] + 56.97) / 20.29));
ICl_n_Tau = 0.69;
ICl_n[cur_cell] = ICl_n[cur_cell] + dt * ((ICl_n_steady_state - ICl_n[cur_cell]) / ICl_n_Tau);
if (V[cur_cell] <= -20.0)
{ gCl_value = gCl * (19.26 + 365.8 * exp(0.0564 * V[cur_cell]));}
else
{ gCl_value = gCl * (0.752 + 0.0349 * exp(-0.414 * V[cur_cell]));}
if (V[cur_cell] > 0)
{
return gCl_value * (Cai[cur_cell] / (Cai[cur_cell] + KdCl)) * pow(ICl_c1[cur_cell],0.5) * (V[cur_cell] - ECl_pos);
}
else
{
return gCl_value * (Cai[cur_cell] / (Cai[cur_cell] + KdCl)) * pow(ICl_c2[cur_cell],1) * pow(ICl_n[cur_cell],1) * (V[cur_cell] - ECl_neg);
}
}
/* ####################################
* # #
* # Calcium-Pump Current - ICapump #
* # #
* #################################### */
double ICapump ()
{
const double JpCa(7e-7); //Max flux in mM/ms
const double KpCa(0.0005); //Half saturation of cai in mM
const double hpCa(2); //Hill coefficient
return ((zCa*Faraday*Cell_Volume)/(Membrane_Capacitance*Cell_Area))*(JpCa/(1+pow((KpCa/Cai[cur_cell]),hpCa)));
}
/* ###############################################
* # #
* # Sodium/Potassium Exchanger Current - INaK #
* # #
* ############################################### */
double INaK ()
{
const double gNaK(1.7);
const double KdK(2), KdNa(22);
const double nK(1.5), nNa(2); //Hill coefficients
double kNaK = 1 / (1 + pow((KdK / Ko),nK));
double nNaK = 1 / (1 + pow((KdNa / Nai[cur_cell]),nNa));
double fNaK = 1 / (1 + 0.125 * exp(-(0.1 * V[cur_cell] * Faraday) / (R * Temperature)) +
0.00219 * exp(Nao / 49.71) * exp(-(1.9 * V[cur_cell] * Faraday) / (R * Temperature)));
return (gNaK * fNaK * kNaK * nNaK);
}
/* ##############################################
* # #
* # Sodium/Calcium Exchanger Current - INaCa #
* # #
* ############################################## */
double INaCa ()
{
double JNaCa(3.5e-6); // Maximum flux in mM/ms
const double KmAllo(0.0003), KmNai(30), KmCai(0.007), KmNao(87.5), KmCao(1.3);
const double nallo(4);
const double NaCa_energy_barrier(0.35); // Gamma, energy barrier of Na/Ca exchanger in the electric field
const double ksat(0.27); // Saturation if INaCa at negative potentials
double VF_RT = (V[cur_cell]*Faraday)/(R*Temperature);
double barrier = exp(NaCa_energy_barrier * VF_RT);
double barrier_1 = exp((NaCa_energy_barrier-1) * VF_RT);
double Allo = 1 / (1 + pow((KmAllo / Cai[cur_cell]),nallo));
double E1 = JNaCa * pow(Nai[cur_cell],3) * Cao * barrier - JNaCa * pow(Nao,3) * Cai[cur_cell] * barrier_1;
double E2 = 1 + ksat * barrier_1;
double E3 = KmCao * pow(Nai[cur_cell],3) + pow(KmNao,3) * Cai[cur_cell] + pow(KmNai,3) * Cao * (1 + (Cai[cur_cell] / KmCai));
double E4 = Cao * pow(Nai[cur_cell],3) + pow(Nao,3) * Cai[cur_cell] + pow(Nao,3) * KmCai * (1 + (pow(Nai[cur_cell],3)/KmNai));
double delta_E = E1 / (E2 * (E3 + E4));
double INaCa_JNaCa = Allo * delta_E;
return ((INaCa_JNaCa * zCa * Cell_Volume * Faraday) / (Cell_Area * Membrane_Capacitance));
}
/* #############################################################
* # #
* # Sodium-Potassium-Chloride Co-Transport Current - INaKCl #
* # #
* ############################################################# */
double INaKCl ()
{
double LNaKCl(1.79e-8);
double RNaKCl(1); // Co-Transport regulation
return (-RNaKCl * zCl * Cell_Area * LNaKCl * R * Faraday * Temperature *
log(((Nao * Ko) / (Nai[cur_cell] * Ki[cur_cell])) *pow( (Clo / Cli[cur_cell]),2)));
}
/* ##############################
# #
# Kinetic Cross-Bridge Cycle #
# #
############################## */
void KCBC ()
{
double Wortmannin = 1; //1=off, 0=on
const double MLCK_max = 0.84;
//KCa_MLCK is a key component here due to it being relative to the Cai and then raised to a power
//Reducing it causes KCBC to be more sensistive to calcium levels
const double K_MLCK = 721.746e-6,KCa_MLCK = 6*180e-6;
const double pM=1, nM=8.7613;
double MLCK = Wortmannin * MLCK_max/(1+pow((Cai[cur_cell]/K_MLCK),pM));
double K1 = MLCK/(1+pow((KCa_MLCK/Cai[cur_cell]),nM));
double K6 = K1;
const double K2 = 0.4, K3 = 1.8, K4 = 0.1, K5 = 0.4, K7 = 0.045; //Yang 2003
M = M+dt*(-K1*M+K2*Mp+K7*AM);
Mp = Mp+dt*(K4*AMp+K1*M-(K2+K3)*Mp);
AMp = AMp+dt*(K3*Mp+K6*AM-(K4+K5)*AMp);
AM = AM+dt*(K5*AMp-(K7+K6)*AM);
Frac_phos = AMp + Mp;
Frac_stress = AMp + AM;
}
/* ####################
* # #
* # Force Production #
* # #
* #################### */
void force ()
{
lc=120; // cell length
double kp=0.1, kx1=12.5, kx2=8.8, ks=0.2, mu_s=0.01, lc0=40, ls0=30;
double lopt=100, fAMp=1.3, fAM=85.5, vx=5, beta=7.5, alpha_p=0.1, alpha_s=4.5;
KCBC();
Fp = kp*(exp(alpha_p*(lc-lc0)/lc0)-1);
Fx = (kx1*AMp+kx2*AM)*lx*exp(-beta*(pow(((la-lopt)/lopt),2)));
Fa = (fAMp*AMp*(vx+la)+fAM*AM*la)*exp(-beta*(pow(((-la-lopt)/lopt),2)));
Fs = Fa;
Ft = Fa+Fp;
// Following lengths are initialised in main program
ls = ls+dt*((1/mu_s)*(Fs-ks*(exp(alpha_s*(ls-ls0)/ls0)-1)));
la = la+dt*((Fa/exp(-beta*(pow(((la-lopt)/lopt),2)))-(fAMp*AMp*vx))/(fAM*AM+fAMp*AMp));
lx=(lc-la-ls);
}
/* ######################
* # #
* # Ryanodine receptor #
* # #
* ###################### */
void RR ()
{
const double Kr1=2500,Kr2=1.05,K_r1=0.0076,K_r2=0.084; //according to Yang 2003
R00[cur_cell]+=dt*(-(Kr1*Cai[cur_cell]+Kr2)*Cai[cur_cell])*R00[cur_cell];
R10[cur_cell]+=dt*(Kr1*Cai[cur_cell]*Cai[cur_cell]*R00[cur_cell]-(K_r1+Kr2*Cai[cur_cell])*R10[cur_cell]+K_r2*R11[cur_cell]);
R01[cur_cell]+=dt*(Kr2*Cai[cur_cell]*R00[cur_cell]+K_r1*R11[cur_cell]-(K_r2+Kr1*Cai[cur_cell]*Cai[cur_cell])*R01[cur_cell]);
R11[cur_cell]+=dt*(Kr2*Cai[cur_cell]*R10[cur_cell]-(K_r1+K_r2)*R11[cur_cell]+Kr1*Cai[cur_cell]*Cai[cur_cell]*R01[cur_cell]);
// Initialised above
}
/* ##################################################
* # #
* # Sarcoplasmic Reticulum control and SR currents #
* # #
* ################################################## */
void ISR ()
{
double Iup_max=6.68*0.1, Kmup=1e-3*10, vol_u=1.33e-9, vol_r=1.33e-11, tau_tr=180,Frel=0.2,Rleak=75.07e-4,tau_rel=0.015,K_CSQN=0.8,CSQN=15,conc_red=1e6;
RR();
Iup = Iup_max * (Cai[cur_cell]/(Cai[cur_cell]+Kmup));
Itr = ((Ca_u[cur_cell]-Ca_r[cur_cell])*zCa*Faraday*vol_u*conc_red)/tau_tr;
Irel = Frel*((R10[cur_cell]*R10[cur_cell]+Rleak)*(Ca_r[cur_cell]-Cai[cur_cell])*zCa*Faraday*vol_r*conc_red)/tau_rel;
Ca_u[cur_cell] += dt*((Iup-Itr)/(zCa*Faraday*vol_u*conc_red));
Ca_r[cur_cell] += dt*((Itr-Irel)/(zCa*Faraday*vol_r*conc_red))*pow(1+(CSQN*K_CSQN)/pow(K_CSQN+Ca_r[cur_cell],2),-1);
}
/* ##################
* # #
* # Calcium buffer #
* # #
* ################## */
double CaBuffer ()
{
double FBCa = 2.0e1, Scm = 0.3, Bf = 0.2, Kd = 2.6e-4, Kdb = 5.39810e-4;
double frac1 = (Scm*Kd)/pow(Kd+Cai[cur_cell],2);
double frac2 = (Bf*Kdb)/pow(Kdb+Cai[cur_cell],2);
return FBCa*pow(1+frac1+frac2,-1);
}
/* ########################
* # #
* # Total current - Itot #
* # #
* ######################## */
double Itotal(double nif)
{
return Ih_mag*Ih()+ICaL_mag*ICaL(nif)+ICaT_mag*ICaT()+INa_mag*INa()
+IK1_mag*IK1()+IK2_mag*IK2()+IKA_mag*IKA()+IKCa_mag*IKCa()
+Inscc_mag*Inscc()+ICl_mag*ICl()+ICapump_mag*ICapump()+INaK_mag*INaK()
+INaCa_mag*INaCa()+Isus+Isoc_mag*Isoc()+IKleak_mag*IKleak()
+INaKCl_mag*INaKCl();
}
/* #############################################
* # #
* # Intracellular calcium concentration - Cai #
* # #
* ############################################# */
double ConCai(double nif)
{
ISR();
ICa_tot = ICaL_mag*ICaL(nif) + ICaT_mag*ICaT() - 2*INaCa_mag*INaCa()
+ InsCa_mag*InsCa + IsocCa_mag*IsocCa
+ ICapump_mag*ICapump() + Iup_mag*Iup - Irel_mag*Irel;
Cai[cur_cell] +=dt * (-(((Cell_Area * Membrane_Capacitance)
/(zCa*Cell_Volume*Faraday))
*(ICa_tot)*CaBuffer()));
return Cai[cur_cell];
}
/* ########################################
* # #
* # Intracellular potassium current - Ki #
* # #
* ######################################## */
void IK()
{
ISR();
IK_tot = IK1_mag*IK1() + IK2_mag*IK2() + IKA_mag*IKA()
+ IKCa_mag*IKCa() - 2*InsK_mag*InsK
+ IKleak_mag*IKleak() + INaKCl_mag*INaKCl() - 2*INaK_mag*INaK();
}
}; // end of cell class
int main ()
{
// Generate output files
ofstream output_file1, output_file2, output_file3;
// Open the files
output_file1.open ("Output_AP.dat");
output_file2.open ("Output_sstau.dat");
output_file3.open ("Output_force.dat");
ofstream IC;
IC.open("IC.txt");
//Column titles in output files
// Currents
output_file1 <<
"time" <<"\t"<< "v" <<"\t"<< "Cai" <<"\t"<< "Itot" <<"\t"<< "ICa_tot" <<"\t"<< "IK_tot" <<"\t"<<
"Ih" <<"\t"<< "ICaL" <<"\t"<< "ICaT" <<"\t"<< "INa" <<"\t"<< "IK1" <<"\t"<< "IK2" <<"\t"<< "IKA" <<"\t"<<
"IKCa" <<"\t"<< "Inscc" <<"\t"<< "ICl" <<"\t"<< "ICapump" <<"\t"<< "INaK" <<"\t"<< "INaCa" <<"\t"<<
"Isus" <<"\t"<< "Isoc" <<"\t"<< "INaKCl" <<"\t"<< "IKleak" << "\t" <<
"InsCa" << "\t" << "InsNa" << "\t" << "InsK" <<"\t"<< "Insleak" << "\t"
"Iup" <<"\t"<< "Irel" <<"\t"<< "Ca_u" <<"\t"<< "Ca_r" <<"\t"<<"Itr"<< "\t"<<
"cur_cell" <<"\t"<< "Istim" <<"\t"<< "Ft" <<"\n";
/*SS and time const*/
output_file2 <<
"time" << "\t" << "v" << "\t" << "cai" << "\t" << "Ih_y_steady_state" << "\t" << "Ih_y_Tau" << "\t" <<
"ICaL_d_steady_state" << "\t" << "ICaL_d_Tau" << "\t" <<
"ICaL_f1_steady_state" << "\t" << "ICaL_f1_Tau" << "\t" << "ICaL_f2_Tau" << "\t" <<
"ICaT_b_steady_state" << "\t" << "ICaT_b_Tau" << "\t" << "ICaT_g_steady_state" << "\t" << "ICaT_g_Tau" << "\t" <<
"INa_m_steady_state" << "\t" << "INa_m_Tau" << "\t" << "INa_h_steady_state" << "\t" << "INa_h_Tau" << "\t" <<
"IK1_q_steady_state" << "\t" << "IK1_q_Tau" << "\t" << "IK1_r_steady_state" << "\t" << "IK1_r1_Tau" << "\t" << "IK1_r2_Tau" << "\t" <<
"IK2_p_steady_state" << "\t" << "IK2_p_Tau" << "\t" << "IK2_k_steady_state" << "\t" << "IK2_k1_Tau" << "\t" << "IK2_k2_Tau" << "\t" <<
"IKA_s_steady_state" << "\t" << "IKA_s_Tau" << "\t" << "IKA_x_steady_state" << "\t" << "IKA_x_Tau" << "\t" <<
"IKCa_p_steady_state" << "\t" << "IKCa_pf_Tau" << "\t" << "IKCa_ps_Tau" << "\t" <<
"Inscc_fMg_steady_state" << "\t" <<
"ICl_c_steady_state" << "\t" << "ICl_c1_Tau" << "\t" << "ICl_c2_Tau" << "\t" << "ICl_n_steady_state" << "\t" << "ICl_n_Tau" << "\n";
/*Force*/
output_file3 <<
"time" << "\t" << "v" << "\t" << "Cai" << "\t" <<
"M" << "\t" << "AM" << "\t" << "Mp" << "\t" << "AMp" << "\t" <<
"Frac_phos" << "\t" << "Frac_stress" << "\t" <<
"Ft" << "\t" << "Fs" << "\t" << "Fa" << "\t" << "Fx" << "\t" << "Fp" << "\t" <<
"Istim" << "\t" << "ICa_tot" << "\t" <<
"lc" << "\t" << "ls" << "\t" << "la" << "\t" << "lx" << "\n";
// Generate cells
cell myocyte[num_cells];
// Intialise time counter for outputs
int counter = 0;
// Set nifedipine conc
double nif = 1e-12;
// Set Stim condutions
// int pulse_time=5e3;
// int pulse_duration=100;
// int pulse_counter=0;
//Time loop (remember time is in ms, dt is 0.01 ms. So after 100 time loops Time == 1. E.g. Time == 1e3 == 1 second)
for (double Time=0; Time<=25e3; (Time += dt) && (counter++))
{
//if (counter==10e5) ICl_mag=ICl_mag*2;
// Check voltage and set INaK_mag
if (V[cur_cell]>0) {INaK_mag = 2;}
else {INaK_mag = 0;}
//----------------Stimuli---------------------
//Default stim to zero (or reset for loops)
Istim=0;
//Stimulus current (simple single pulse)
if ((Time > 5e3 && Time <= 15e3))
{
//if (Time>9850 && Time<=9860) Istim = -0.400; //for short pulses during AP change Istim to higher value
//else Istim = -0.400;
Istim = -0.400;
}//0.05 //Iclamp;
else {Istim = 0.;}
//Repetative stimulus
/* if (Time > 5e3 && pulse_counter<10 && (counter*dt)>pulse_time && (counter*dt)<=pulse_time+pulse_duration)
{
Istim = -0.400;
}//0.05 //Iclamp;
else {Istim = 0.;}
//if (counter%10000==0) cout << counter*dt << "\t" << pulse_time << "\t" << pulse_duration << endl;
if ((counter*dt)==pulse_time+pulse_duration)
{
pulse_time+=400;
pulse_counter++;
}
*/
//--------------------------------------------
//Calculate intracellular calcium
myocyte[cur_cell].ConCai(nif);
//Calculate net current through membrane
Itot = myocyte[cur_cell].Itotal(nif);
//Calculate forces
myocyte[cur_cell].force();
//Calculate total K current
myocyte[cur_cell].IK();
//------------Voltage protocols---------------
//Voltage protocol (forward Euler)
V[cur_cell] += (-(Itot + Istim) * dt)/Membrane_Capacitance;
//Voltage clamp pulses
/* if (Time > 5e3 && pulse_counter<10 && (counter*dt)>pulse_time && (counter*dt)<=pulse_time+pulse_duration)
{
V[cur_cell] = -30;
}
else V[cur_cell] = -60.;
if ((counter*dt)==pulse_time+pulse_duration)
{
pulse_time+=400;
pulse_counter++;
}
*/
//--------------------------------------------
//Output to terminal
if (counter%10000==0) // 100 for every ms, 1000 for every 10 ms
{
cout << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << endl;
}
//Output data to files
if (counter%100==0 /*&& Time>11e3*/) // 100 for every ms, 1000 for every 10 ms
{
//Currents
output_file1 <<
Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" <<
Itot << "\t" <<
ICa_tot << "\t" <<
IK_tot << "\t" <<
Ih_mag*myocyte[cur_cell].Ih() << "\t" <<
ICaL_mag*myocyte[cur_cell].ICaL(nif) << "\t" <<
ICaT_mag*myocyte[cur_cell].ICaT() << "\t" <<
INa_mag*myocyte[cur_cell].INa() << "\t" <<
IK1_mag*myocyte[cur_cell].IK1() << "\t" <<
IK2_mag*myocyte[cur_cell].IK2() << "\t" <<
IKA_mag*myocyte[cur_cell].IKA() << "\t" <<
IKCa_mag*myocyte[cur_cell].IKCa() << "\t" <<
Inscc_mag*myocyte[cur_cell].Inscc() << "\t" <<
ICl_mag*myocyte[cur_cell].ICl() << "\t" <<
ICapump_mag*myocyte[cur_cell].ICapump() << "\t" <<
INaK_mag*myocyte[cur_cell].INaK() << "\t" <<
INaCa_mag*myocyte[cur_cell].INaCa() << "\t" <<
Isus << "\t" <<
Isoc_mag*myocyte[cur_cell].Isoc() << "\t" <<
INaKCl_mag*myocyte[cur_cell].INaKCl() << "\t" <<
IKleak_mag*myocyte[cur_cell].IKleak() << "\t" <<
InsCa_mag*InsCa << "\t" <<
InsNa_mag*InsNa << "\t" <<
InsK_mag*InsK << "\t" <<
Insleak_mag*Insleak << "\t" <<
Iup_mag*Iup << "\t" <<
Irel_mag*Irel << "\t" <<
Ca_u[cur_cell] << "\t" <<
Ca_r[cur_cell] << "\t" <<
Itr_mag*Itr << "\t" <<
cur_cell << "\t" <<
Istim << "\t" <<
Ft << "\n";
/*SS and time const*/
output_file2 << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" << Ih_y_steady_state << "\t" << Ih_y_Tau << "\t" <<
ICaL_d_steady_state << "\t" << ICaL_d_Tau << "\t" <<
ICaL_f_steady_state << "\t" << ICaL_f1_Tau << "\t" << ICaL_f2_Tau << "\t" <<
ICaT_b_steady_state << "\t" << ICaT_b_Tau << "\t" << ICaT_g_steady_state << "\t" << ICaT_g_Tau << "\t" <<
INa_m_steady_state << "\t" << INa_m_Tau << "\t" << INa_h_steady_state << "\t" << INa_h_Tau << "\t" <<
IK1_q_steady_state << "\t" << IK1_q_Tau << "\t" << IK1_r_steady_state << "\t" << IK1_r1_Tau << "\t" << IK1_r2_Tau << "\t" <<
IK2_p_steady_state << "\t" << IK2_p_Tau << "\t" << IK2_k_steady_state << "\t" << IK2_k1_Tau << "\t" << IK2_k2_Tau << "\t" <<
IKA_s_steady_state << "\t" << IKA_s_Tau << "\t" << IKA_x_steady_state << "\t" << IKA_x_Tau << "\t" <<
IKCa_p_steady_state << "\t" << IKCa_pf_Tau << "\t" << IKCa_ps_Tau << "\t" <<
Inscc_fMg_steady_state << "\t" <<
ICl_c_steady_state << "\t" << ICl_c1_Tau << "\t" << ICl_c2_Tau << "\t" << ICl_n_steady_state << "\t" << ICl_n_Tau << "\n";
/*Force*/
output_file3 << Time << "\t" << V[cur_cell] << "\t" << Cai[cur_cell] << "\t" <<
M << "\t" << AM << "\t" << Mp << "\t" << AMp << "\t" <<
Frac_phos << "\t" << Frac_stress << "\t" <<
Ft << "\t" << Fs << "\t" << Fa << "\t" << Fx << "\t" << Fp<< "\t" <<
Istim << "\t" << ICa_tot << "\t" <<
lc << "\t" << ls << "\t" << la << "\t" << lx << "\n";
}
if (counter>1100e4 && counter<=1101e4)//Time>(Time-1e2))
{
IC << "Time" << "\t" << "Ih_y[0]" << "\t" << "ICaL_d[0]" << "\t" << "ICaL_f1[0]" << "\t" << "ICaL_f2[0]" << "\t" <<"ICaT_b[0]"
<< "\t" << "ICaT_g[0]" << "\t" << "INa_m[0]" << "\t" << "INa_h[0]" << "\t" << "IK1_q[0]"
<< "\t" << "IK1_r1[0]" << "\t" << "IK1_r2[0]" << "\t" << "IK2_p[0]" << "\t" << "IK2_k1[0]" << "\t" << "IK2_k2[0]"
<< "\t" << "IKA_s[0]" << "\t" << "IKA_x[0]"
<< "\t" << "ICl_c1[0]" << "\t" << "ICl_c2[0]" << "\t" << "ICl_n[0]" << "\t" << "V[0]"
<< "\t" << "Cai[0]" << "\t" << "Ca_u[0]" << "\t" << "Ca_r[0]" << "\t" << "R00[0]"
<< "\t" << "R01[0]" << "\t" << "R10[0]" << "\t" << "R11[0]" << "\t" << "V_KCa_half"
<< "\t" << "IKCa_pf[0]" << "\t" << "IKCa_ps[0]" << endl;
IC << Time << "\t" << Ih_y[cur_cell] << "\t" << ICaL_d[cur_cell] << "\t" << ICaL_f1[cur_cell] << "\t" << ICaL_f2[cur_cell] << "\t" <<ICaT_b[cur_cell]
<< "\t" << ICaT_g[cur_cell] << "\t" << INa_m[cur_cell] << "\t" << INa_h[cur_cell] << "\t" << IK1_q[cur_cell]
<< "\t" << IK1_r1[cur_cell] << "\t" << IK1_r2[cur_cell] << "\t" << IK2_p[cur_cell] << "\t" << IK2_k1[cur_cell] << "\t" << IK2_k2[cur_cell]
<< "\t" << IKA_s[cur_cell] << "\t" << IKA_x[cur_cell]
<< "\t" << ICl_c1[cur_cell] << "\t" << ICl_c2[cur_cell] << "\t" << ICl_n[cur_cell] << "\t" << V[cur_cell]
<< "\t" << Cai[cur_cell] << "\t" << Ca_u[cur_cell] << "\t" << Ca_r[cur_cell] << "\t" << R00[cur_cell]
<< "\t" << R01[cur_cell] << "\t" << R10[cur_cell] << "\t" << R11[cur_cell] << "\t" << V_KCa_half[cur_cell]
<< "\t" << IKCa_pf[cur_cell] << "\t" << IKCa_ps[cur_cell] << endl;
}
} //end of time loop
//cell A;
} // end of main
/*Notes:
* Elements cell class must contain:
* channel functions to work out current through membrane
* membrane voltage
* intracellular ion handling
*
* Elements outside of class:
* diffusion between cells
* number of cells etc
*
* */
