// This is a program for parallel simulations of a model of cardiac tissue // // written by Christian Zemlin March 9, 2007 // Basic simulation parameters const double msSimulationTime=10.0; const double outputInterval=1.0; // ms const int beats= (int) (msSimulationTime/outputInterval); const int width = 200; const int height = 80; const int length = (width+2)*(height+2); const int numberAttempts=2; const double stimWidth=3 ; // sidelength of square to be stimulated double BCL=500.0; const bool EPgradient=false; // Vary EP properties from LA to PV? // Scan pMin const double pMinMin=0.0; const double pMinMax=0.01; const double pMinStep=0.05; // PV Electrophysiology const double iNaFactor=1.0; // to weaken iNa const double PViK1=0.58; const double PVito=0.75; const double PViCaL=0.7; const double PViKs=1.6; const double PViKr=1.5; const double PVKi=145.0; const double LAKi=145.0; // noise const double noiselevel=1.0; const double p0=0.0; // Diffusion const double Dmin=1.0; // minimum diffusion constant const double Dmax=4.0; // maximum diffusion constant const double kx=1.0; const double kymin=0.25; const double kymax=1.0; const double breakThroughX=width*4/5; // where (in x dir) trans. in D occurs const double transitionY=height*2/3+0.5; // where (in y dir) trans. occurs #include #include #include #include #include vector Dx; // diffusion in x dir (as function of space) vector Dy; // diffusion in y dir (as function of space) double randgauss( double min, double max, double sigma, double center); class cmMedium // defines a 1d org Courtmanche medium { public: cmMedium(); void Allocate(); // Allocate memory for tabulation arrays void Tabulate(); // Tabulate functions void ShowState(int pos); // print out all state variables void ComputeIndices(); // compute tabulaton indices void DumpBinary(); // dump state var content in bin form, needs outfile void Advance(); void Initialize(double pMin); void ComputeGatingVariables(); void Create_Cai_LookupTables(); void Create_zNai_LookupTables(); void Create_zCarel_LookupTables(); void CreateCurrentLookupTables(); void CreateGatingLookupTables(); void SetStableStartValues(); void ComputeCurrents(); void ComputeVmAndConcentrations(); void appendMovie(char* filename, int upperBound); void writePicture(char* filename, int upperBound); void writeState(char* filename); void readState(char* filename); void widen(int newHeight); void addToMovie2d(ofstream& movie); // state variables vector Vm; vector m; vector h; vector j; vector xr; vector xs; vector oa; vector oi; vector ua; vector ui; vector d; vector f; vector fCa; vector Cai; vector Caup; vector qCa; vector Carel; vector Nai; vector Ki; vector u; vector v; vector Cli; vector w; vector Fn; vector dv; vector iK1Factor; vector iCaLFactor; vector itoFactor; vector iKrFactor; vector iKsFactor; vector tstim; vector stimtime; //other accessible Variables double StimPeriod ; // stimulation period double deltat ; // time step int TimeSteps ; // number of time steps double t; // model time double IstimSize; // size of stimulus current int Step ; // current simulation step double tmod ; // model time modulo stimulation period double Itot ; // total current int Index ; // Vm index double RVmTab ; // 1/(tabulation step) for Vm double dtstim; // stimulation time double duration; // of the simulation double Istim ; // strength of stimulus current double stimDuration; int length; int pos; int geoHeight; // for 2D media: geometry (I use this 1D class for 2D) int geoWidth; // for 2D media: geometry bool dimensionsChanged; //to inform myStepAdi that arrays need to be resized private: int zCarel_size; int zfCass_size; int IndexMax ; double zNai_start ; double zNai_end ; double zNai_incr ; double zNai_invincr ; int zNai_size ; double zCarel_start ; double zCarel_end ; double zCarel_incr ; double zCarel_invincr ; double B1, one_over_B2; double Cm ; double gbCa_bar ; double gbNa_bar ; double gCaL_bar ; double gK1_bar ; double gKr_bar ; double gKs_bar ; double gNa_bar ; double gto_bar ; double gCaCl_bar; double INaCamax_bar ; double INaKmax_bar ; double IpCamax_bar ; double gbCa ; double gbNa ; double gCaL ; double gK1 ; double gKr ; double gKs ; double gNa ; double gto ; double gCaCl; double INaCamax ; double INaKmax ; double IpCamax ; double Cao ; double Caupmax ; double Cmdnmax ; double Csqnmax ; double dCaith ; double Clo; double Fc ; double CMgamma ; double Iupmax ; double Ko ; double KmCa ; double Kmcmdn ; double Kmcsqn ; double KmKo ; double KmNai ; double KmNao ; double Kmtrpn ; double KQ10 ; double krel ; double ksat ; double Kup ; double Nao ; double PRNaK ; double tau_fCa ; double tau_tr ; double Trpnmax ; double Vcell ; double Vi ; double Vrel ; double Vup ; double cmdnmax_times_kmcmdn ; double csqnmax_times_kmcsqn ; double one_over_2FcVi ; double one_over_Caupmax ; double one_over_FcVi ; double one_over_tau_tr ; double one_over_Vi ; double Fn_k1 ; double Fn_k2 ; double Ko_on_KoKmKo ; double RTon2F ; double RTonF ; double trpnmax_times_kmtrpn ; double Vrel_over_Vup ; double VmTab ; vector Ca_cmdn ; vector Ca_csqn ; vector Ca_trpn ; double dCa_cmdn_dt; double dCa_trpn_dt; double dCa_csqn_dt; int Cai_nM ; vector dss; double ECa ; double EK ; double ENa ; double ECl; vector expd; vector CMexpf; double expfCa ; vector exph; vector expj; vector expm; vector expoa; vector expoi; double expu ; vector expua; vector expui; double expv_191 ; double expv_400 ; vector expw; vector expxr; vector expxs; double fCass ; vector fss; double qCass; vector hss; double IbCa ; double IbK ; double IbNa ; double ICaL ; double IK1 ; double IKr ; double IKs ; double IKur ; double INa ; double INaCa ; double INaK ; double IpCa ; double ICaCl ; double aux; double expqCa; double Irel ; double Ito ; double Itr ; double Iup ; double Iupleak ; double IV ; vector jss; vector mss; vector oass; vector oiss; double SampleRate ; vector uass; vector uiss; double vss ; vector wss; vector xrss; vector xsss; vector zB2; vector zCarel; int zCarel_index ; vector zfCass; vector zfNaK; vector zgKur; vector zICaLCa1; vector zICaLCa2; vector zICaLK; vector zICaLNa; vector zIK1; vector zIKr; vector zINaCa1; vector zINaCa2; vector zIpCa; vector zIup; vector zIV; vector zNai; int zNai_index ; double sigma ; double microsecs ; double StimRate ; }; #define MCW MPI_COMM_WORLD void splitDomain(int* begin, int* end, int* slaveShare, int width, int numProc); void diffuse2dSimple(double* v, double* dv, double dt, int width, int height, int* begin, int* end, bool first); void readInput2d(double* v, char* filename, int width, int height); void initialize2d(double* v, int width, int height); void resetUpper(double* v, int width, int height); void resetSpecial(double* v, int width, int height); void binDump(char *filename); void binDumpPart(char *filename, int xmin, int xmax, int ymin, int ymax); void binRestore(char *filename); void output(ofstream& scrollFile, vector, int width, int height, int beats,bool first, int interrupt); void reset(int pos); void makeRefractory(int pos); const int MASTER = 0; const int TAG = 1; int myRank, numProc; int numSlaves; main (int argc, char** argv) { int i; int beat; int xpos, ypos; // ofstream outFile("CMmovie.var"); int attempt; int pos; //position on the ring // set up MPI MPI_Init(&argc, &argv); MPI_Comm_rank(MCW, &myRank); MPI_Comm_size(MCW, &numProc); MPI_Status status; numSlaves = numProc-1; char final[]="./in"; // define intervals on which the individual slaves work int* begin = new int[numProc]; int* end = new int[numProc]; int* slaveShare = new int[numProc]; splitDomain(begin, end, slaveShare, width, numSlaves); int rank; double pMin=0.0; cmMedium cm; cm.length=length; cm.Initialize(pMin); int steps = (int)(outputInterval/cm.deltat); if (myRank == MASTER) { int whichSlave; int frame; char filename[255]; ofstream scrollFile; // where the output for Scroll is written bool first; // is this the first frame of the current attempt // (there can be multiple attempts per program run ofstream DxFile; ofstream DyFile; char DxFilename[255]; char DyFilename[255]; int exCount; bool done; int interrupt; char dirName[255]; char command[255]; for (pMin= pMinMin; pMin-20) exCount++; } if ((exCount==-1) && (frame>10)) { cout << "We are done" << endl; done=true; interrupt=frame; } // let the slaves know if they should continue for (whichSlave = 1; whichSlave <= numSlaves; whichSlave++) { MPI_Send(&interrupt,1 , MPI_INT, whichSlave, TAG, MCW); } output(scrollFile, cm.Vm, width, height,beats,first,interrupt); first=false; } DxFile.close(); DyFile.close(); scrollFile.close(); } } } else // not master { int increment; char dummy; bool first=true; int interrupt; for (pMin=pMinMin; pMin 1) // transmit information at left border { MPI_Send(&(cm.Vm[0])+begin[myRank]*(height+2), (height+2), MPI_DOUBLE, myRank - 1, TAG, MCW); MPI_Recv(&(cm.Vm[0])+(begin[myRank]-1)*(height+2), (height+2), MPI_DOUBLE, myRank-1,TAG, MCW, &status); } if (myRank < numSlaves) // transmit information at right border { MPI_Send(&(cm.Vm[0])+(end[myRank]-1)*(height+2), (height+2), MPI_DOUBLE, myRank +1, TAG, MCW); MPI_Recv(&(cm.Vm[0])+end[myRank]*(height+2), (height+2), MPI_DOUBLE, myRank+1,TAG, MCW, &status); } */ // Uncomment the following section for periodic boundary conditions /* p if (myRank == 1) // transmit information at left border { MPI_Send(&(cm.Vm[0])+begin[myRank]*(height+2), (height+2), MPI_DOUBLE, numSlaves, TAG, MCW); MPI_Recv(&(cm.Vm[0])+(begin[myRank]-1)*(height+2), (height+2), MPI_DOUBLE, numSlaves,TAG, MCW, &status); } if (myRank == numSlaves) // transmit information at right border { MPI_Send(&(cm.Vm[0])+(end[myRank]-1)*(height+2), (height+2), MPI_DOUBLE, 1, TAG, MCW); MPI_Recv(&(cm.Vm[0])+end[myRank]*(height+2), (height+2), MPI_DOUBLE, 1,TAG, MCW, &status); } */ cm.t += cm.deltat; }// for increment .. } // for beat .. } } // if myRank == SLAVE // if (myRank==1) // cm.ShowState(height+2+2); MPI_Finalize(); } void splitDomain(int* begin, int* end, int* slaveShare, int width, int numSlaves) { int i; begin[1] = 1; for (i=1; i 0.001 ) tau = z2 * (1.0-z1) / (0.035*(VmBuf+10.0)); // [54] else tau = z2 * (1.0/6.24) / 0.035; expd[Index] = exp(-deltat/tau); dss[Index] = 1.0 / (1.0 + exp(-(VmBuf+10.0)/6.0)); // [54] // f tau = 400.0 / ( 1.0+4.5* exp(-0.0007*(VmBuf-9.0)*(VmBuf-9.0))); // [55] CMexpf[Index] = exp(-deltat/tau); fss[Index] = 1.0 / ( 1.0 + exp((VmBuf+24.6)/6.2)); // [55] // fCa expfCa = exp(-deltat/tau_fCa); // constant once deltat known [56] // *** Gating of iNa: // m if ( fabs(VmBuf+47.13) > 0.001) alpha = 0.32 * (VmBuf+47.13) / ( 1.0 - exp(-0.1*(VmBuf+47.13)) ); // [30] else alpha = 3.2; beta = 0.08 * exp(-VmBuf/11.0); // [30] z1 = alpha + beta; tau = 1.0 / z1; expm[Index] = exp(-deltat/tau); mss[Index] = alpha / z1; // h if ( VmBuf < -40.0 ) { alpha = 0.135 * exp(-(VmBuf+80.0)/6.8); beta = 3.56 * exp(0.079*VmBuf) + 3.1E+5 * exp(0.35*VmBuf); // [31] } else { alpha = 0.0; beta = 1.0 / ( 0.13 * ( 1.0 + exp(-(VmBuf+10.66)/11.1) ) ); } z1 = alpha + beta; tau = 1.0 / z1; exph[Index] = exp(-deltat/tau); hss[Index] = alpha / z1; // j if ( VmBuf < -40.0 ) { alpha = ( -1.2714E+5 * exp(0.2444*VmBuf) - 3.474E-5 * exp(-0.04391*VmBuf) ) * (VmBuf+37.78) / ( 1.0 + exp(0.311*(VmBuf+79.23)) ); // [32] beta = 0.1212 * exp(-0.01052*VmBuf) / ( 1.0 + exp(-0.1378*(VmBuf+40.14)) ); // [33] } else { alpha = 0.0; beta = 0.3 * exp(-2.535E-7*VmBuf) / ( 1.0 + exp(-0.1*(VmBuf+32.0)) ); } z1 = alpha + beta; tau = 1.0 / z1; expj[Index] = exp(-deltat/tau); jss[Index] = alpha / z1; // *** Gating of IKr and IKs: // xr if (( fabs ( VmBuf-248.0 ) < 0.005) || (fabs(VmBuf+163.0)<0.005)) { expxr[Index] = expxr[Index - 1]; // where functions ill-defined xrss[Index] = xrss[Index-1]; // ensure continuity } else { alpha = 0.04 * (VmBuf - 248.0) / (1.0 - exp(-(VmBuf-248.0)/28.0)); beta = 0.028 * (VmBuf + 163.0) / (exp((VmBuf+163.0)/21.0) - 1.0); tau = 1.0 / (alpha + beta); expxr[Index] = exp(-deltat / tau); xrss[Index] = 1.0 / (1.0 + exp(-(VmBuf + 7.654) / 5.377)); // [49] } // xs if ( fabs ( VmBuf +28.5 ) < 0.005) { expxs[Index] = expxs[Index - 1]; // where functions ill-defined xsss[Index] = xsss[Index-1]; // ensure continuity } else { alpha = 1.0E-5 * (VmBuf +28.5) / (1.0 - exp(-(VmBuf+28.5)/115.0)); beta = 2.3E-4 * (VmBuf +28.5) / (exp((VmBuf+28.5)/3.3) - 1.0); tau = 1.0 / (alpha + beta); expxs[Index] = exp(-deltat / tau); xsss[Index] = 1.0 / sqrt(1.0 + exp(-(VmBuf-13.0) / 12.0)); // [52] } // *** Gating of Ito: // oa alpha = 0.65 / (exp(-(VmBuf+18.0)/8.5) + exp(-(VmBuf-16.0)/59.0)); // [37] beta = 1.2 / (2.2 + exp((VmBuf+75.0)/18.0)); tau = 1.0 / ((alpha + beta)); // [38] expoa[Index] = exp(-deltat/tau); oass[Index] = pow(1.0 + exp(-(VmBuf+0.5)/10.5),-0.333333333); // [38] // oi alpha = 1.0 / (6.2 + exp((VmBuf + 105.2)/9.85)); // [39] beta = 1.0 / (7.54 + exp(-(VmBuf -8.89)/12.87)); tau = 1.0 / ((alpha + beta)); // [40] expoi[Index] = exp(-deltat/tau); oiss[Index] = 1.0 / (1.0 + exp((VmBuf+43.377)/6.45)); // [40] // *** Gating of IKur: // ua alpha = 1.47 / (exp(-(VmBuf+33.2)/30.63) + exp(-(VmBuf-27.6)/30.65)); // [43] beta = 0.42 / (exp((VmBuf+26.64)/2.49) + exp((VmBuf+44.41)/20.36)); tau = 1.0 / ((alpha + beta)); // [44] expua[Index] = exp(-deltat/tau); uass[Index] = pow(1.0 + exp(-(VmBuf+2.81)/9.49),-0.333333333); // [44] // ui alpha = 1.0 / (21.0 + exp(-(VmBuf-185.0)/28.0)); // [45] beta = exp((VmBuf-158.0)/16.0); tau = 1.0 / ((alpha + beta)); // [46] expui[Index] = exp(-deltat/tau); uiss[Index] = 1.0/(1.0 + exp((VmBuf-99.45)/27.48)); // [46] // *** Gating of Irel: // (Note that gating variables u and v cannot be in lookup tables because of Fn dependence) // w z1 = VmBuf - 7.9; if (fabs(z1) > 0.001) { z2 = exp(-z1/5.0); tau = 6.0 * (1.0 - z2) / ((1.0 + 0.3 * z2) * z1); // [67] } else tau = 0.923077; expw[Index] = exp(-deltat/tau); wss[Index] = 1.0 - 1.0 / (1.0 + exp(-(VmBuf-40.0)/17.0)); // [67] // u tau_u = 8.0; expu = exp(-deltat/tau_u); } } void cmMedium::SetStableStartValues() { int i; // allocate for state variables // initialize state variables for (i=0; i 6.835E-14) vss = 0.0; else vss = 1.0; if (Fn[pos] > 3.4175E-13) { u[pos] = 1.0 + (u[pos] - 1.0 ) * expu; v[pos] = vss + (v[pos] - vss ) * expv_400; } else { u[pos] = (u[pos] ) * expu; v[pos] = vss + (v[pos] - vss ) * expv_191; } w[pos] = wss[Index] + (w[pos] - wss[Index]) * expw[Index]; } // Subroutine to compute the currents of the model cell: void cmMedium::ComputeCurrents() { double Vm_minus_EK=Vm[pos]-EK;; // spare variable ECa = 13.35 * log(Cao / Cai[pos]); //[28] ENa = 26.71 * log(Nao / Nai[pos]); EK = 26.71 * log(Ko / Ki[pos]); ECl = -26.71 * log(Clo/ Cli[pos]); INa = iNaFactor*gNa * m[pos] * m[pos] * m[pos] * h[pos] * j[pos] * (Vm[pos] - ENa); Vm_minus_EK = Vm[pos] - EK; IK1 = iK1Factor[pos]*gK1 * Vm_minus_EK * zIK1[Index]; //[35] Ito = itoFactor[pos]*gto * oa[pos] * oa[pos] * oa[pos] * oi[pos] * Vm_minus_EK; //[36] IKur = zgKur[Index] * ua[pos] * ua[pos] * ua[pos] * ui[pos] * Vm_minus_EK; IKr = iKrFactor[pos]*xr[pos] * Vm_minus_EK * zIKr[Index]; IKs = iKsFactor[pos]*gKs * (xs[pos] * xs[pos]) * Vm_minus_EK; //[50] ICaL = iCaLFactor[pos]*gCaL * d[pos] * f[pos] * fCa[pos] * (Vm[pos] - 65.0); //[53] ICaCl = gCaCl * qCa[pos] *(Vm[pos]-ECl); // - because Cl is negative ion INaK = zfNaK[Index] * zNai[zNai_index]; //[57] INaCa = zINaCa1[Index] * (Nai[pos] * Nai[pos] * Nai[pos]) - Cai[pos] * zINaCa2[Index]; //[60] IbCa = gbCa * (Vm[pos] - ECa); //[61] IbK = 0.0; IbNa = gbNa * (Vm[pos] - ENa); //[62] IpCa = zIpCa[Cai_nM]; //[63] Irel = krel * u[pos] * u[pos] * v[pos] * w[pos] * (Carel[pos] - Cai[pos]); Itr = (Caup[pos] - Carel[pos]) * one_over_tau_tr; //[69] Iup = zIup[Cai_nM]; Iupleak = (Caup[pos] * Iupmax) * one_over_Caupmax; //[72] IV = IK1 + INaK + IpCa + INaCa + IbNa + IbCa; Itot = INa + IK1 + Ito + IKur + IKr + IKs + ICaL + IpCa + INaK + INaCa + IbNa + IbCa + ICaCl; //[19] if (t>=tstim[pos] && t<(tstim[pos]+deltat)) { cout << "start stimulating at " << pos << " " << t << " " << tstim[pos] << " " << Istim << endl; stimtime[pos] = 0; tstim[pos]=-1 ; if (t>2000) BCL=1000.0; } if((stimtime[pos]>=0) && (stimtime[pos]=100*connectionFraction) Dy[pos]=0.0; Ki[pos]=PVKi+(LAKi-PVKi)/(1+exp(- (kymax+(kymin-kymax)*1/(1+exp(-kx*(xpos-breakThroughX))))*(ypos-transitionY))); if ((ypos<2) && (xposheight/2) //iK1Factor[pos]=1.0; // if (ypos>height-5) // tstim[pos]=0.0; } } void cmMedium::Advance() { for (pos=0; pos v, int width, int height, int beats, bool first, int interrupt) { char vUC; int xpos, ypos; if (first) { int dimensions = 3; int frames=beats; int bytes=1; int framerate=100; // cout << "beats in there " << beats << endl; //cout << "movie has " << beats << " frames." << endl; scrollFile.write((char*) &dimensions, 4); scrollFile.write((char*) &width, 4); scrollFile.write((char*) &height, 4); scrollFile.write((char*) &beats, 4); scrollFile.write((char*) &bytes, 4); scrollFile.write((char*) &framerate, 4); first=false; } for (ypos=1; ypos<=height; ypos++) for (xpos=1; xpos<=width; xpos++) { vUC = (char) ((int) ((v[xpos*(height+2) + ypos]+ 100.0))); scrollFile.write(&vUC,1); } if (interrupt) { scrollFile.seekp(12); scrollFile.write((char*) &interrupt,4); } } double randgauss( double min, double max, double sigma, double center) { double random = (min + (max-min) * (double)rand()/RAND_MAX); //create random domain between [min,max] double tmp = (random-center)/sigma; double gauss= exp(-tmp*tmp/2); //gaussian formula return gauss; }