GCC Code Coverage Report

Directory: ./
File: Solver/PaciSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 438 0.0%
Branches: 0 238 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: F. Raphel
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/PaciSolver.hpp"
21
22
23 //************************************
24 //Constants
25 //************************************
26
27 static constexpr double R = 8.314472; // [J/mol/K] Gas constant
28 static constexpr double F = 96485.3415; // [C/mol] Faraday constant
29 static constexpr double T = 310.0; // [K] Temperature
30 static constexpr double Q = 2.3;
31 static constexpr double Lo = 0.025; //[dimensionless]
32 static constexpr double Nao = 151.; // [mM]
33 static constexpr double Ko = 5.4; // [mM]
34 static constexpr double Cao = 1.8; // [mM]
35 static constexpr double Ki = 150.; // [mM]//150 doc---- 140 Matlab
36 static constexpr double gCaL = 8.635702e-5; //[m3/F.s]
37 static constexpr double gKr = 29.8667; //[S/F]
38 static constexpr double gKs = 2.041;
39 static constexpr double gf = 30.10312;
40 static constexpr double arel = 16.464; //[mM/s]
41 static constexpr double brel = 0.25;
42 static constexpr double crel = 8.232;
43 static constexpr double Vleak = 4.4444e-4;
44 static constexpr double gpCa = 0.4125; //[A/F]
45 static constexpr double gbNa = 0.9; //[S/F]
46 static constexpr double gbCa = 0.69264;
47 static constexpr double Bufc = 0.25; //[mM]
48 static constexpr double Bufsr = 10.;
49 static constexpr double Kbufc = 0.001; //[mM]
50 static constexpr double Kbufsr = 0.3;
51 static constexpr double Kup = 0.00025; //[mM]
52 static constexpr double KpCa = 0.0005;
53 static constexpr double Pkna = 0.03;
54 static constexpr double Ksat = 0.1;
55 static constexpr double KmCa = 1.38; //[mM]
56 static constexpr double KmNai = 87.5;
57 static constexpr double alpha = 2.8571432; //[dimensionless]
58 static constexpr double gamma_value = 0.35;
59 static constexpr double KmNa = 40.; //[mM]
60 static constexpr double Kmk = 1.;
61
62 #ifdef FELISCE_WITH_SUNDIALS
63
64 #define Ith(v,i) NV_Ith_S(v,i-1) /* Ith numbers components 1..NEQ */
65 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */
66
67 static constexpr double Y0 = RCONST(0.0); /* initial y components */
68 static constexpr double Y1 = RCONST(1.0);
69 static constexpr double Y075 = RCONST(0.75);
70 static constexpr double Y01 = RCONST(0.1);
71 static constexpr double Ycai = RCONST(0.0002);
72 static constexpr double Ycasr = RCONST(0.3);
73 static constexpr double Ynai = RCONST(10.);
74 static constexpr double Y0bis = RCONST(1e-15); /* initial y components */
75
76 static constexpr double YmORd = RCONST(0.0074621);//minfORd
77 static constexpr double YjORd = RCONST( 0.692477);//jinfORd
78 static constexpr double Yhslow = RCONST(0.692574);//hslow
79 static constexpr double Yhfast = RCONST(0.692591);//hfast
80 static constexpr double Yhl = RCONST(0.496116);//hl
81 static constexpr double Yml = RCONST(0.000194015);//ml
82
83 //static constexpr double RTOL = RCONST(1.0e-6); /* scalar relative tolerance ----- 1e-6: doc Sundials */
84 //static constexpr double RTOL = RCONST(1.0e-4); //Nouveau test: avant: spontane
85
86 //static constexpr double RTOL = RCONST(1.0e-1); //Non sponta ok!
87 //static constexpr double RTOL = RCONST(1.0e-4);
88 static constexpr double RTOL = RCONST(1.0e0); //Non sponta ok!
89
90
91
92 //static constexpr double ATOL4 = RCONST(1.0e-4); /* std::vector absolute tolerance components */
93 static constexpr double ATOL4 = RCONST(1.0e-4); /* std::vector absolute tolerance components */
94 static constexpr double ATOL6 = RCONST(1.0e-6);
95 static constexpr double ATOL8 = RCONST(1.0e-8);
96 static constexpr double ATOL10 = RCONST(1.0e-10);
97 static constexpr double ATOL12 = RCONST(1.0e-12);
98 static constexpr double ATOL14 = RCONST(1.0e-14);
99
100 //static constexpr double RTOL = RCONST(1.0e-4); /* scalar relative tolerance 1e-4: 0.01% */
101 //static constexpr double ATOL = RCONST(1.0e-4); /* std::vector absolute tolerance components */
102
103 #endif
104
105 namespace felisce {
106
107 /////////////////
108 // Constructor //
109 /////////////////
110
111 ✗ PaciSolver::PaciSolver(FelisceTransient::Pointer fstransient):
112 ✗ IonicSolver(fstransient) {
113 ✗ m_nComp = 42;
114 }
115
116 /////////////////
117 /////////////////
118
119
120 ////////////////
121 // Destructor //
122 ////////////////
123
124 ✗ PaciSolver::~PaciSolver() = default;
125
126 ////////////////
127 ////////////////
128
129
130 //////////////////////////////////////////////////////////////
131 // ComputeRHS: useless in the non Spontaneous model (cvode) //
132 //////////////////////////////////////////////////////////////
133
134 ✗ void PaciSolver::computeRHS() {
135
136 #ifndef FELISCE_WITH_SUNDIALS
137 //Only spontaneous model
138 ✗ double& dt = m_fstransient->timeStep;
139
140 ✗ m_bdf.computeRHSTime(dt, m_vecRHS);
141
142
143 felInt pos;
144 double value_uExtrap,value_Cai,value_df1,value_fCa,value_g;
145
146
147 ✗ std::vector<double> values_RHS;
148 ✗ std::vector<double> valuem_solEDO;
149
150
151 ✗ values_RHS.resize(42);
152 ✗ valuem_solEDO.resize(42);
153
154
155 ✗ values_RHS.clear();
156 ✗ valuem_solEDO.clear();
157
158 ✗ int& numIt = m_fstransient->iteration;
159
160 ✗ for (felInt i = 0; i < m_size; i++) {
161 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
162 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);
163
164 //Initialization
165 ✗ m_vecSolEDO[4].getValues(1,&pos,&value_fCa);
166 ✗ m_vecSolEDO[13].getValues(1,&pos,&value_g);
167 ✗ m_vecSolEDO[14].getValues(1,&pos,&value_Cai);
168 ✗ m_vecSolEDO[17].getValues(1,&pos,&value_df1);
169
170 ✗ for(int iEDO=0; iEDO<42; iEDO++) {
171 ✗ m_vecSolEDO[iEDO].getValues(1,&pos,&valuem_solEDO[iEDO]);//value_RHS = _RHS(i)
172 }
173
174 ✗ if(numIt == 1) {
175 ✗ value_fCa = 1.;
176 ✗ value_g = 1.;
177 ✗ value_Cai = 0.0002;
178 ✗ value_df1 = 0.;
179 }
180
181 double hinf,tauh;
182 ✗ PaciSolver::hgate(hinf,tauh,value_uExtrap);
183
184 double jinf,tauj;
185 ✗ PaciSolver::jgate(jinf,tauj,value_uExtrap);
186
187 double minf,taum;
188 ✗ PaciSolver::mgate(minf,taum,value_uExtrap);
189
190 double dinf,taud;
191 ✗ PaciSolver::dgate(dinf,taud,value_uExtrap);
192
193 double fCainf,taufCa,constfCa;
194 ✗ PaciSolver::fcagate(fCainf,taufCa,constfCa,value_uExtrap,value_Cai,value_fCa);
195
196 double f1inf,tauf1;
197 ✗ PaciSolver::f1gate(f1inf,tauf1,value_uExtrap,value_df1,value_Cai);
198
199 double f2inf,tauf2;
200 ✗ PaciSolver::f2gate(f2inf,tauf2,value_uExtrap);
201
202 double rinf,taur;
203 ✗ PaciSolver::rgate(rinf,taur,value_uExtrap);
204
205 double qinf,tauq;
206 ✗ PaciSolver::qgate(qinf,tauq,value_uExtrap);
207
208 double xr1inf,tauxr1;
209 ✗ PaciSolver::xr1gate(xr1inf,tauxr1,value_uExtrap);
210
211 double xr2inf,tauxr2;
212 ✗ PaciSolver::xr2gate(xr2inf,tauxr2,value_uExtrap);
213
214 double xsinf,tauxs;
215 ✗ PaciSolver::xsgate(xsinf,tauxs,value_uExtrap);
216
217 double xfinf,tauf;
218 ✗ PaciSolver::xfgate(xfinf,tauf,value_uExtrap);
219
220 double ginf,taug,constg;
221 ✗ PaciSolver::ggate(ginf,taug,constg,value_uExtrap,value_Cai,value_g);
222
223
224 //! For non spontaneous AP
225 double minfORd,taumORd;
226 ✗ PaciSolver::mgateORd(minfORd, taumORd,value_uExtrap);
227
228 double jinfORd,taujORd;
229 ✗ PaciSolver::jgateORd(jinfORd, taujORd,value_uExtrap);
230
231 double hslowinf,tauhslow;
232 ✗ PaciSolver::hslowgateORd(hslowinf, tauhslow,value_uExtrap);
233
234 double hfastinf,tauhfast;
235 ✗ PaciSolver::hfastgateORd(hfastinf, tauhfast,value_uExtrap);
236
237 double hlinf,tauhl;
238 ✗ PaciSolver::hlgateORd(hlinf, tauhl,value_uExtrap);
239
240 double mlinf,tauml;
241 ✗ PaciSolver::mlgateORd(mlinf, tauml,value_uExtrap);
242
243
244 //RHS:
245 ✗ values_RHS[0]=hinf/tauh;
246 ✗ values_RHS[1]=jinf/tauj;
247 ✗ values_RHS[2]=minf/taum;
248 ✗ values_RHS[3]=dinf/taud;
249 ✗ values_RHS[4]=constfCa*fCainf/taufCa;
250 ✗ values_RHS[5]=f1inf/tauf1;
251 ✗ values_RHS[6]=f2inf/tauf2;
252 ✗ values_RHS[7]=rinf/taur;
253 ✗ values_RHS[8]=qinf/tauq;
254 ✗ values_RHS[9]=xr1inf/tauxr1;
255 ✗ values_RHS[10]=xr2inf/tauxr2;
256 ✗ values_RHS[11]=xsinf/tauxs;
257 ✗ values_RHS[12]=xfinf/tauf;
258 ✗ values_RHS[13]=constg*ginf/taug;
259 ✗ values_RHS[14]=0.;
260 ✗ values_RHS[15]=0.;
261 ✗ values_RHS[16]=0.;
262 ✗ values_RHS[17]=0.;
263 ✗ values_RHS[18]=0.;
264 ✗ values_RHS[19]=0.;
265 ✗ values_RHS[20]=0.;
266 ✗ values_RHS[21]=0.;
267 ✗ values_RHS[22]=0.;
268 ✗ values_RHS[23]=0.;
269 ✗ values_RHS[24]=0.;
270 ✗ values_RHS[25]=0.;
271 ✗ values_RHS[26]=0.;
272 ✗ values_RHS[27]=0.;
273 ✗ values_RHS[28]=0.;
274 ✗ values_RHS[29]=0.;
275 ✗ values_RHS[30]=0.;
276 ✗ values_RHS[31]=0.;
277 ✗ values_RHS[32]=0.;
278 ✗ values_RHS[33]=0.;
279 ✗ values_RHS[34]=0.;
280 ✗ values_RHS[35]=0.;
281
282
283 //! For non spontaneous AP: ORd
284
285 ✗ values_RHS[36]=minfORd/taumORd;
286 ✗ values_RHS[37]=jinfORd/taujORd;
287 ✗ values_RHS[38]=hslowinf/tauhslow;
288 ✗ values_RHS[39]=hfastinf/tauhfast;
289 ✗ values_RHS[40]=hlinf/tauhl;
290 ✗ values_RHS[41]=mlinf/tauml;
291
292 ✗ for (std::size_t index = 0 ; index < 42 ; ++index)
293 ✗ m_vecRHS[index].setValue(pos,values_RHS[index], ADD_VALUES);
294
295
296 }
297
298
299 ✗ for(int index =0 ; index < 42 ; ++index) {
300 ✗ m_vecRHS[index].assembly();
301 }
302
303 #endif
304
305 }
306
307
308 ✗ void PaciSolver::solveEDO() {
309
310
311 //################//
312 //Ionic parameters//
313 //################//
314
315 ✗ const double Cm = 98.7109;
316 ✗ const double Vc = 8800.;
317 ✗ const double Vsr = 583.73;
318 ✗ const double gNa = 3.6712302e3;
319 ✗ const double PNaK = 1.841424;
320 ✗ const double KNaCa = 4900.;
321 ✗ const double Vmaxup = 0.56064;
322
323
324 //##################
325 //##################
326
327
328
329 double value_uExtrap;
330 ✗ std::vector<double> values_RHS;
331 ✗ std::vector<double> valuem_solEDO;
332
333 ✗ values_RHS.resize(42);
334 ✗ valuem_solEDO.resize(42);
335
336 ✗ values_RHS.clear();
337 ✗ valuem_solEDO.clear();
338
339 ✗ double& dt = m_fstransient->timeStep;
340 ✗ double& coeffDeriv = m_bdf.coeffDeriv0();
341
342
343 double betaEDOh,betaEDOj,betaEDOm,betaEDOd,betaEDOfca,betaEDOf1,betaEDOf2,betaEDOr,betaEDOq,betaEDOxr1,betaEDOxr2,betaEDOxs,betaEDOxf,betaEDOg,betaEDOcai,betaEDOcasr,betaEDOnai;
344 double Cai,fca,df1,g,Casr,Nai,f1;
345
346 ✗ int& numIt = m_fstransient->iteration;
347 felInt pos;
348
349 #ifdef FELISCE_WITH_SUNDIALS
350
351 if(numIt == 1) { //resizes
352 M_userData.M_u_extrap.resize(m_size);
353 M_userData.value_df1.resize(m_size);
354 M_userData.vINaCa.resize(m_size);
355 M_userData.vIPCa.resize(m_size);
356 M_userData.vIbCa.resize(m_size);
357 M_userData.vENa.resize(m_size);
358 M_userData.vEKs.resize(m_size);
359 M_userData.vECa.resize(m_size);
360 M_userData.vINa.resize(m_size);
361 M_userData.vICaL.resize(m_size);
362 M_userData.vINaK.resize(m_size);
363 M_userData.vIrel.resize(m_size);
364 M_userData.vIup.resize(m_size);
365 M_userData.vIleak.resize(m_size);
366 M_userData.vIbNa.resize(m_size);
367 cvode_mem.resize(m_size);
368 initvalue.resize(m_size);
369 abstol.resize(m_size);
370 }
371
372 #endif
373
374 ✗ for (felInt i = 0; i < m_size; i++) {
375
376 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
377 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
378
379
380
381 #ifdef FELISCE_WITH_SUNDIALS
382
383 double time = m_fstransient->time;
384 int flag;
385
386 M_userData.M_u_extrap[i] = value_uExtrap;
387 M_userData.pos = i;
388
389 if(numIt == 1) { //1st iteration --> initializations for each Dof
390
391
392 M_userData.value_df1[i] = 0.;
393
394 cvode_mem[i] = NULL;
395 initvalue[i] = N_VNew_Serial(23);
396 abstol[i] = N_VNew_Serial(23);
397
398
399 Ith(initvalue[i],1) = Y075;//h
400 Ith(initvalue[i],2) = Y075;//j
401 Ith(initvalue[i],3) = Y0bis;//m
402 Ith(initvalue[i],4) = Y0bis;//d
403 Ith(initvalue[i],5) = Y1;//fca
404 Ith(initvalue[i],6) = Y1;//f1
405 Ith(initvalue[i],7) = Y1;//f2
406 Ith(initvalue[i],8) = Y0bis;//r
407 Ith(initvalue[i],9) = Y1;//q
408 Ith(initvalue[i],10) = Y0bis;//xr1
409 Ith(initvalue[i],11) = Y1;//xr2
410 Ith(initvalue[i],12) = Y0bis;//xs
411 Ith(initvalue[i],13) = Y01;//xf
412 Ith(initvalue[i],14) = Y1;//g
413 Ith(initvalue[i],15) = Ycai;//Cai
414 Ith(initvalue[i],16) = Ycasr;//Casr
415 Ith(initvalue[i],17) = Ynai;//Nai
416
417
418 Ith(initvalue[i],18) = YmORd;//minfORd
419 Ith(initvalue[i],19) = YjORd;//jinfORd
420 Ith(initvalue[i],20) = Yhslow;//hslow
421 Ith(initvalue[i],21) = Yhfast;//hfast
422 Ith(initvalue[i],22) = Yhl;//hl
423 Ith(initvalue[i],23) = Yml;//ml
424
425
426 Ith(abstol[i],1) = ATOL4;
427 Ith(abstol[i],2) = ATOL4;
428 Ith(abstol[i],3) = ATOL4;
429 Ith(abstol[i],4) = ATOL4;
430 Ith(abstol[i],5) = ATOL4;
431 Ith(abstol[i],6) = ATOL4;
432 Ith(abstol[i],7) = ATOL4;
433 Ith(abstol[i],8) = ATOL4;
434 Ith(abstol[i],9) = ATOL4;
435 Ith(abstol[i],10) = ATOL4;
436 Ith(abstol[i],11) = ATOL4;
437 Ith(abstol[i],12) = ATOL4;
438 Ith(abstol[i],13) = ATOL4;
439 Ith(abstol[i],14) = ATOL4;
440 Ith(abstol[i],15) = ATOL4;
441 Ith(abstol[i],16) = ATOL4;
442 Ith(abstol[i],17) = ATOL4;
443
444 Ith(abstol[i],18) = ATOL4;
445 Ith(abstol[i],19) = ATOL4;
446 Ith(abstol[i],20) = ATOL4;
447 Ith(abstol[i],21) = ATOL4;
448 Ith(abstol[i],22) = ATOL4;
449 Ith(abstol[i],23) = ATOL4;
450
451 cvode_mem[i] = CVodeCreate(CV_BDF,CV_NEWTON);//CV_BDF or CV_ADAMS and CV_NEWTON or CV_FUNCTIONAL || (CV_BDF,CV_NEWTON) for stiff pb
452
453 datatest = &M_userData;
454 flag = CVodeSetUserData(cvode_mem[i],datatest);
455 flag = CVodeInit(cvode_mem[i],M_f,0.0,initvalue[i]);
456
457
458 //Tolerances************
459 //
460 flag = CVodeSVtolerances(cvode_mem[i],RTOL,abstol[i]);
461 //flag = CVodeWFtolerances(cvode_mem, ewt);
462
463 //**********************
464
465
466 //Optional inputs********************
467 //
468 //flag = CVSpgmr(cvode_mem, PREC_LEFT, 0);
469 //flag = CVSpilsSetGSType(cvode_mem, MODIFIED_GS);
470 //flag = CVodeRootInit(cvode_mem,2,g);//Optional
471 //
472 //************************************
473
474 flag = CVDense(cvode_mem[i],23);//--->if ORd? else 17?
475
476
477 //Optional inputs******************
478 //
479
480 //flag = CVodeSetMaxNumSteps(cvode_mem[i],500);//default: 500
481 //flag = CVodeSetMaxNumSteps(cvode_mem[i],1000);//ok!
482 flag = CVodeSetMaxNumSteps(cvode_mem[i],1000);//default: 500
483
484
485 flag = CVodeSetMaxStep(cvode_mem[i],0.1);//if 0: inf
486
487 //---------------
488 //
489 //For more informations, print flag
490 //
491 // If flag = 0: SUCCESS
492 // If flag = 1: TSTOP_RETURN
493 //....see cvode.h
494
495
496 realtype tout;
497 int flag;
498
499
500 flag = CVode(cvode_mem[i], time, initvalue[i], &tout, CV_NORMAL);
501
502
503 if ( flag != CV_SUCCESS ) {
504 std::cout << "CVODE solver failed " << std::endl;
505 exit(1);
506 }
507
508 } else {
509
510 double df1;
511 m_vecSolEDO[17].getValues(1,&pos,&df1);
512 M_userData.value_df1[i] = df1;
513 double f1;
514 m_vecSolEDO[5].getValues(1,&pos,&f1);
515
516 realtype tout;
517 int flag;
518
519
520 flag = CVode(cvode_mem[i], time, initvalue[i], &tout, CV_NORMAL);
521
522 if ( flag != CV_SUCCESS ) {
523 std::cout << "CVODE solver failed " << std::endl;
524 exit(1);
525 }
526 }
527
528 for(int j=0; j<17; j++) {
529 valuem_solEDO[j] = Ith (initvalue[i],j+1);
530 }
531
532 valuem_solEDO[17] = (valuem_solEDO[5]-f1)/dt;
533 valuem_solEDO[18] = M_userData.vINaCa[i];
534 valuem_solEDO[19] = M_userData.vIPCa[i];
535 valuem_solEDO[20] = M_userData.vIbCa[i];
536 valuem_solEDO[21] = M_userData.vENa[i];
537 valuem_solEDO[22] = M_userData.vEKs[i];
538 valuem_solEDO[23] = M_userData.vECa[i];
539 valuem_solEDO[24] = M_userData.vINa[i];
540 valuem_solEDO[25] = M_userData.vICaL[i];
541 valuem_solEDO[31] = M_userData.vINaK[i];
542 valuem_solEDO[32] = M_userData.vIrel[i];
543 valuem_solEDO[33] = M_userData.vIup[i];
544 valuem_solEDO[34] = M_userData.vIleak[i];
545 valuem_solEDO[35] = M_userData.vIbNa[i];
546
547 for(int j=0; j<6; j++) {
548 valuem_solEDO[36+j] = Ith (initvalue[i],j+18);
549 }
550
551
552
553 #else
554
555 ✗ for(int iEDO=0; iEDO<36; iEDO++)
556 ✗ m_vecRHS[iEDO].getValues(1,&pos,&values_RHS[iEDO]);//value_RHS
557
558
559 ✗ m_vecSolEDO[4].getValues(1,&pos,&fca);
560 ✗ m_vecSolEDO[5].getValues(1,&pos,&f1);
561 ✗ m_vecSolEDO[13].getValues(1,&pos,&g);
562 ✗ m_vecSolEDO[14].getValues(1,&pos,&Cai);
563 ✗ m_vecSolEDO[15].getValues(1,&pos,&Casr);
564 ✗ m_vecSolEDO[16].getValues(1,&pos,&Nai);
565 ✗ m_vecSolEDO[17].getValues(1,&pos,&df1);
566
567
568 ✗ if(numIt==1) {
569 ✗ fca = 1.;
570 ✗ g = 1.;
571 ✗ Cai = 0.0002;
572 ✗ Casr = 0.3;
573 ✗ Nai = 10.;
574 ✗ f1 = 1.;
575 ✗ df1 = 0.;
576 }
577
578 double hinf,tauh;
579 ✗ PaciSolver::hgate(hinf,tauh,value_uExtrap);
580
581 double jinf,tauj;
582 ✗ PaciSolver::jgate(jinf,tauj,value_uExtrap);
583
584 double minf,taum;
585 ✗ PaciSolver::mgate(minf,taum,value_uExtrap);
586
587 double dinf,taud;
588 ✗ PaciSolver::dgate(dinf,taud,value_uExtrap);
589
590 double fCainf,taufCa,constfCa;
591 ✗ PaciSolver::fcagate(fCainf,taufCa,constfCa,value_uExtrap,Cai,fca);
592 double f1inf,tauf1;
593 ✗ PaciSolver::f1gate(f1inf,tauf1,value_uExtrap,df1,Cai);
594
595 double f2inf,tauf2;
596 ✗ PaciSolver::f2gate(f2inf,tauf2,value_uExtrap);
597
598 double rinf,taur;
599 ✗ PaciSolver::rgate(rinf,taur,value_uExtrap);
600
601 double qinf,tauq;
602 ✗ PaciSolver::qgate(qinf,tauq,value_uExtrap);
603
604 double xr1inf,tauxr1;
605 ✗ PaciSolver::xr1gate(xr1inf,tauxr1,value_uExtrap);
606
607 double xr2inf,tauxr2;
608 ✗ PaciSolver::xr2gate(xr2inf,tauxr2,value_uExtrap);
609
610 double xsinf,tauxs;
611 ✗ PaciSolver::xsgate(xsinf,tauxs,value_uExtrap);
612
613 double xfinf,tauf;
614 ✗ PaciSolver::xfgate(xfinf,tauf,value_uExtrap);
615
616 double ginf,taug,constg;
617 ✗ PaciSolver::ggate(ginf,taug,constg,value_uExtrap,Cai,g);
618
619
620
621 //! For non spontaneous AP
622 double minfORd,taumORd;
623 ✗ PaciSolver::mgateORd(minfORd, taumORd,value_uExtrap);
624
625 double jinfORd,taujORd;
626 ✗ PaciSolver::jgateORd(jinfORd, taujORd,value_uExtrap);
627
628 double hslowinf,tauhslow;
629 ✗ PaciSolver::hslowgateORd(hslowinf, tauhslow,value_uExtrap);
630
631 double hfastinf,tauhfast;
632 ✗ PaciSolver::hfastgateORd(hfastinf, tauhfast,value_uExtrap);
633
634 double hlinf,tauhl;
635 ✗ PaciSolver::hlgateORd(hlinf, tauhl,value_uExtrap);
636
637 double mlinf,tauml;
638 ✗ PaciSolver::mlgateORd(mlinf, tauml,value_uExtrap);
639
640
641
642
643 //h
644 ✗ betaEDOh = coeffDeriv/dt+1./tauh;
645 ✗ valuem_solEDO[0] = values_RHS[0]/betaEDOh;
646
647 //j
648 ✗ betaEDOj = coeffDeriv/dt+1./tauj;
649 ✗ valuem_solEDO[1] = values_RHS[1]/betaEDOj;
650
651 //m
652 ✗ betaEDOm = coeffDeriv/dt+1./taum;
653 ✗ valuem_solEDO[2] = values_RHS[2]/betaEDOm;
654
655 //d
656 ✗ betaEDOd = coeffDeriv/dt+1./taud;
657 ✗ valuem_solEDO[3] = values_RHS[3]/betaEDOd;
658
659 //fCa
660 ✗ betaEDOfca = coeffDeriv/dt+constfCa/taufCa;
661 ✗ valuem_solEDO[4] = values_RHS[4]/betaEDOfca;
662
663 //f1
664 ✗ betaEDOf1 = coeffDeriv/dt+1./tauf1;
665 ✗ valuem_solEDO[5] = values_RHS[5]/betaEDOf1;
666
667 //f2
668 ✗ betaEDOf2 = coeffDeriv/dt+1./tauf2;
669 ✗ valuem_solEDO[6] = values_RHS[6]/betaEDOf2;
670
671 //r
672 ✗ betaEDOr = coeffDeriv/dt+1./taur;
673 ✗ valuem_solEDO[7] = values_RHS[7]/betaEDOr;
674
675 //q
676 ✗ betaEDOq = coeffDeriv/dt+1./tauq;
677 ✗ valuem_solEDO[8] = values_RHS[8]/betaEDOq;
678
679 //xr1
680 ✗ betaEDOxr1 = coeffDeriv/dt+1./tauxr1;
681 ✗ valuem_solEDO[9] = values_RHS[9]/betaEDOxr1;
682
683 //xr2
684 ✗ betaEDOxr2 = coeffDeriv/dt+1./tauxr2;
685 ✗ valuem_solEDO[10] = values_RHS[10]/betaEDOxr2;
686
687 //xs
688 ✗ betaEDOxs = coeffDeriv/dt+1./tauxs;
689 ✗ valuem_solEDO[11] = values_RHS[11]/betaEDOxs;
690
691 //xf
692 ✗ betaEDOxf = coeffDeriv/dt+1./tauf;
693 ✗ valuem_solEDO[12] = values_RHS[12]/betaEDOxf;
694
695 //g
696 ✗ betaEDOg = coeffDeriv/dt+constg/taug;
697 ✗ valuem_solEDO[13] = values_RHS[13]/betaEDOg;
698
699
700
701 //! For non spontaneous AP
702
703 //mORd
704 ✗ const double betaEDOmORd = coeffDeriv/dt+1./taumORd;
705 ✗ valuem_solEDO[18] = values_RHS[18]/betaEDOmORd;
706
707 //jORd
708 ✗ const double betaEDOjORd = coeffDeriv/dt+1./taujORd;
709 ✗ valuem_solEDO[19] = values_RHS[19]/betaEDOjORd;
710
711 //hslow
712 ✗ const double betaEDOhslow = coeffDeriv/dt+1./tauhslow;
713 ✗ valuem_solEDO[20] = values_RHS[20]/betaEDOhslow;
714
715 //hfast
716 ✗ const double betaEDOhfast = coeffDeriv/dt+1./tauhfast;
717 ✗ valuem_solEDO[21] = values_RHS[21]/betaEDOhfast;
718
719 //hl
720 ✗ const double betaEDOhl = coeffDeriv/dt+constfCa/tauhl;
721 ✗ valuem_solEDO[22] = values_RHS[22]/betaEDOhl;
722
723 //ml
724 ✗ const double betaEDOml = coeffDeriv/dt+1./tauml;
725 ✗ valuem_solEDO[23] = values_RHS[23]/betaEDOml;
726
727
728
729 //Calcul des courants pour les concentrations
730 ✗ double value_iup = Vmaxup/(1.+std::pow(Kup/Cai,2));
731 ✗ double value_ileak = Vleak*(Casr-Cai);
732
733
734 ✗ double value_irel = (crel+(arel*std::pow(Casr,2))/(std::pow(brel,2)+std::pow(Casr,2)))*valuem_solEDO[3]*valuem_solEDO[13];
735 ✗ value_irel = value_irel*0.0411;//ventricular
736 //Atrial
737 //value_irel = value_irel*0.0556;
738
739
740
741 //iNaK
742 ✗ const double num1 = PNaK*Ko*Nai/(Ko+Kmk);
743 ✗ const double num2 = num1/(Nai+KmNa);
744 ✗ double value_iNaK = num2/(1.+0.1245*std::exp(-0.1*value_uExtrap*F/(R*T))+0.0353*std::exp(-value_uExtrap*F/(R*T)));
745
746 //iNaCa
747 ✗ const double prodnum1 = std::pow(Nai,3)*Cao*std::exp(gamma_value*value_uExtrap*F/(R*T));
748 ✗ const double prodnum2 = alpha*std::pow(Nao,3)*Cai*std::exp((gamma_value-1.)*value_uExtrap*F/(R*T));
749 ✗ const double numINaCa = KNaCa*(prodnum1-prodnum2);
750 ✗ const double prodden1 = std::pow(KmNai,3)+std::pow(Nao,3);
751 ✗ const double prodden2 = KmCa+Cao;
752 ✗ const double prodden3 = 1.+Ksat*std::exp((gamma_value-1.)*value_uExtrap*F/(R*T));
753 ✗ const double denINaCa = prodden1*prodden2*prodden3;
754 ✗ double value_iNaCa = numINaCa/denINaCa;
755
756 ✗ valuem_solEDO[18] = value_iNaCa;
757
758 //iCaL
759 ✗ const double term1 = (gCaL*4.*value_uExtrap*F*F)/(R*T);
760 ✗ const double term2 = (Cai*std::exp((2.*value_uExtrap*F)/(R*T))-0.341*Cao)/(std::expm1((2.*value_uExtrap*F)/(R*T)));
761
762 ✗ double value_iCaL = term1*term2*valuem_solEDO[3]*valuem_solEDO[5]*valuem_solEDO[6]*valuem_solEDO[4];
763
764
765 //ipCa
766 ✗ double value_iPCa = gpCa*Cai/(Cai+KpCa);
767
768 ✗ valuem_solEDO[19] = value_iPCa;
769
770
771 ✗ const double ENa = (R*T/F)*std::log(Nao/Nai);
772 ✗ const double ECa = 0.5*(R*T/F)*std::log(Cao/Cai);
773
774
775
776 //iNa
777 ✗ const double m3 = valuem_solEDO[2]*valuem_solEDO[2]*valuem_solEDO[2];
778 ✗ double value_iNa=gNa*m3*valuem_solEDO[0]*valuem_solEDO[1]*(value_uExtrap-ENa);
779
780 //ibNa-ibCa
781 ✗ double value_ibNa=gbNa*(value_uExtrap-ENa);
782 ✗ double value_ibCa=gbCa*(value_uExtrap-ECa);
783
784 ✗ valuem_solEDO[20] = value_ibCa;
785
786 ✗ double Caibufc = 1./(1.+(Bufc*Kbufc)/std::pow(Cai+Kbufc,2));
787 ✗ double Casrbufsr = 1./(1.+(Bufsr*Kbufsr)/std::pow(Casr+Kbufsr,2));
788
789
790 //Cai
791 ✗ const double valCai = Caibufc*(value_ileak-value_iup+value_irel-Cm*(value_iCaL+value_ibCa+value_iPCa-2.*value_iNaCa)/(2.*Vc*F));
792
793 //Casr
794 ✗ const double valCasr = (Casrbufsr*Vc)*(value_iup-(value_irel+value_ileak))/Vsr;
795
796 //Nai
797 ✗ const double valNai = -Cm*(value_iNa+value_ibNa+3.*value_iNaK+3.*value_iNaCa)/(F*Vc);
798
799
800 //Cai
801 ✗ betaEDOcai = coeffDeriv/dt;
802 ✗ values_RHS[14] += valCai;
803 ✗ valuem_solEDO[14] = values_RHS[14]/betaEDOcai;
804
805 //Casr
806 ✗ betaEDOcasr = coeffDeriv/dt;
807 ✗ values_RHS[15] += valCasr;
808 ✗ valuem_solEDO[15] = values_RHS[15]/betaEDOcasr;
809
810 //Nai
811 ✗ betaEDOnai = coeffDeriv/dt;
812 ✗ values_RHS[16] += valNai;
813 ✗ valuem_solEDO[16] = values_RHS[16]/betaEDOnai;
814
815
816 //Calcul de df1
817 ✗ valuem_solEDO[17] = (valuem_solEDO[5]-f1)/dt;
818
819
820 //Calcul de ENa EKs et ECa
821 ✗ valuem_solEDO[21] = (R*T/F)*std::log(Nao/Nai);
822 ✗ valuem_solEDO[22] = (R*T/F)*std::log((Ko+Pkna*Nao)/(Ki+Pkna*Nai));
823 ✗ valuem_solEDO[23] = 0.5*(R*T/F)*std::log(Cao/Cai);
824
825
826
827 #endif
828
829
830 ✗ for(int iEDO=0; iEDO<17; iEDO++) {
831 ✗ m_vecSolEDO[iEDO].setValue(pos,valuem_solEDO[iEDO], INSERT_VALUES);
832 }
833
834 ✗ m_vecSolEDO[17].setValue(pos,valuem_solEDO[17], INSERT_VALUES);
835 ✗ m_vecSolEDO[18].setValue(pos,valuem_solEDO[18], INSERT_VALUES);
836 ✗ m_vecSolEDO[19].setValue(pos,valuem_solEDO[19], INSERT_VALUES);
837 ✗ m_vecSolEDO[20].setValue(pos,valuem_solEDO[20], INSERT_VALUES);
838
839
840 //ENa EKs ECa
841 ✗ m_vecSolEDO[21].setValue(pos,valuem_solEDO[21], INSERT_VALUES);
842 ✗ m_vecSolEDO[22].setValue(pos,valuem_solEDO[22], INSERT_VALUES);
843 ✗ m_vecSolEDO[23].setValue(pos,valuem_solEDO[23], INSERT_VALUES);
844
845
846
847 //Currennts Iup,Ileak,Irel,INaK,INaCa,ICaL,IpCa,INa,IbNa,IbCa
848
849 #ifdef FELISCE_WITH_SUNDIALS
850 m_vecSolEDO[24].setValue(pos,valuem_solEDO[24], INSERT_VALUES);
851 m_vecSolEDO[25].setValue(pos,valuem_solEDO[25], INSERT_VALUES);
852 m_vecSolEDO[31].setValue(pos,valuem_solEDO[31], INSERT_VALUES);
853 m_vecSolEDO[32].setValue(pos,valuem_solEDO[32], INSERT_VALUES);
854 m_vecSolEDO[33].setValue(pos,valuem_solEDO[33], INSERT_VALUES);
855 m_vecSolEDO[34].setValue(pos,valuem_solEDO[34], INSERT_VALUES);
856 m_vecSolEDO[35].setValue(pos,valuem_solEDO[35], INSERT_VALUES);
857
858 //! For non spontaneous AP
859 m_vecSolEDO[36].setValue(pos,valuem_solEDO[18], INSERT_VALUES);
860 m_vecSolEDO[37].setValue(pos,valuem_solEDO[19], INSERT_VALUES);
861 m_vecSolEDO[38].setValue(pos,valuem_solEDO[20], INSERT_VALUES);
862 m_vecSolEDO[39].setValue(pos,valuem_solEDO[21], INSERT_VALUES);
863 m_vecSolEDO[40].setValue(pos,valuem_solEDO[22], INSERT_VALUES);
864 m_vecSolEDO[41].setValue(pos,valuem_solEDO[23], INSERT_VALUES);
865
866 #else
867 ✗ m_vecSolEDO[24].setValue(pos,value_iNa, INSERT_VALUES);
868 ✗ m_vecSolEDO[25].setValue(pos,value_iCaL, INSERT_VALUES);
869 ✗ m_vecSolEDO[31].setValue(pos,value_iNaK, INSERT_VALUES);
870 ✗ m_vecSolEDO[32].setValue(pos,value_irel, INSERT_VALUES);
871 ✗ m_vecSolEDO[33].setValue(pos,value_iup, INSERT_VALUES);
872 ✗ m_vecSolEDO[34].setValue(pos,value_ileak, INSERT_VALUES);
873 ✗ m_vecSolEDO[35].setValue(pos,value_ibNa, INSERT_VALUES);
874 #endif
875
876 }
877
878 ✗ for (std::size_t index = 0 ; index < 42; ++index)
879 ✗ m_vecSolEDO[index].assembly();
880
881
882 }
883
884
885
886
887
888 ✗ void PaciSolver::computeIon() {
889
890 //################//
891 //Ionic parameters//
892 //################//
893
894 ✗ const double gNa = 3.6712302e3;
895 ✗ const double gto = 29.9038;
896 ✗ const double gK1 = 28.1492;
897 ✗ const double PNaK = 1.841424;
898 ✗ const double KNaCa = 4900.;
899
900 //##################
901 //##################
902
903 //*******************************************************//
904 //Declaration des variables: courants et concentrations
905 //*******************************************************//
906
907 double value_ion;
908 double value_uExtrap;
909
910 double value_iCaL;
911 double value_iPCa;
912 double value_iNaCa;
913 double value_iNaK;
914 double value_ik1;
915 double value_ito;
916 double value_ikr;
917 double value_iks;
918 double value_if;
919 double value_iNa;
920
921 //gates:
922 ✗ std::vector<double> valuem_solEDO;
923 ✗ valuem_solEDO.resize(42);
924
925
926 ✗ valuem_solEDO.clear();
927
928 felInt pos;
929
930 //*****************************************//
931 //******** Calculation ******************//
932 //*****************************************//
933
934 ✗ for(felInt i = 0 ; i < m_size; i++) {
935
936
937 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
938 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
939
940 ✗ for(int iEDO=0; iEDO<42; iEDO++) {
941 ✗ m_vecSolEDO[iEDO].getValues(1,&pos,&valuem_solEDO[iEDO]);
942 }
943
944
945 //******************************************
946 //Reversal potential: Nernst equation
947 //
948 ✗ const double Ef = -17e-3; //[mV]
949 ✗ const double EK = (R*T/F)*std::log(Ko/Ki);
950 ✗ const double ENa = valuem_solEDO[21];
951 ✗ const double EKs = valuem_solEDO[22];
952 ✗ const double ECa = valuem_solEDO[23];
953
954 //******************************************
955 //******************************************
956
957
958
959 //*************
960 //Na+ current
961 //
962 ✗ const double m3sponta = valuem_solEDO[2]*valuem_solEDO[2]*valuem_solEDO[2];
963 //-------------------------------------
964 //value_iNa=gNa*m3*valuem_solEDO[0]*valuem_solEDO[1]*(value_uExtrap-ENa);
965 //-------------------------------------
966
967 //! For non spontaneous AP
968
969 ✗ const double hORd = 0.99*valuem_solEDO[39]+0.01*valuem_solEDO[38];
970 ✗ const double m3 = valuem_solEDO[36]*valuem_solEDO[36]*valuem_solEDO[36];
971 ✗ const double iNafast = 75.*(value_uExtrap-ENa)*m3*valuem_solEDO[37]*hORd;
972 ✗ const double iNalate = 0.0075*(value_uExtrap-ENa)*valuem_solEDO[41]*valuem_solEDO[40];
973
974 ✗ if(FelisceParam::instance().spontaneous == false) {
975 ✗ value_iNa = iNafast+iNalate;
976 } else {
977 ✗ value_iNa=gNa*m3sponta*valuem_solEDO[0]*valuem_solEDO[1]*(value_uExtrap-ENa);
978 }
979
980
981 //*********************
982 //L-Type Ca2+ current
983 //
984 ✗ const double term1 = (gCaL*4.*value_uExtrap*F*F)/(R*T);
985 ✗ const double term2 = (valuem_solEDO[14]*std::exp((2.*value_uExtrap*F)/(R*T))-0.341*Cao)/(std::expm1((2.*value_uExtrap*F)/(R*T)));
986 //------------------------------------------
987 ✗ value_iCaL = term1*term2*valuem_solEDO[3]*valuem_solEDO[5]*valuem_solEDO[6]*valuem_solEDO[4];
988 //------------------------------------------
989
990
991
992 //***************************
993 //Transient outward current
994 //---------------------------------
995 ✗ value_ito = gto*valuem_solEDO[7]*valuem_solEDO[8]*(value_uExtrap-EK);
996 //---------------------------------
997
998 //************************************
999 //Rapid delayed rectifier K+ current
1000 //--------------------------------------------------
1001 ✗ value_ikr = gKr*std::sqrt(Ko/5.4)*valuem_solEDO[9]*valuem_solEDO[10]*(value_uExtrap-EK);
1002 //--------------------------------------------------
1003
1004
1005 //***********************************
1006 //Slow delayed rectifier K+ current
1007 //--------------------------------------------------------------------------
1008 ✗ value_iks = gKs*std::pow(valuem_solEDO[11],2)*(1.+0.6/(1.+std::pow(3.8e-5/valuem_solEDO[14],1.4)))*(value_uExtrap-EKs);
1009 //--------------------------------------------------------------------------
1010
1011
1012 //*****************************
1013 //Inward rectifier K+ current
1014 //
1015 ✗ const double alphak1 = 3.91/(1.+std::exp(1000.*0.5942*(value_uExtrap-EK-200e-3)));
1016 ✗ const double betak1 = (-1.509*std::exp(1000.*0.0002*(value_uExtrap-EK+100e-3))+std::exp(1000.*0.5886*(value_uExtrap-EK-10e-3)))/(1.+std::exp(1000.*0.4547*(value_uExtrap-EK)));
1017 ✗ const double xK1inf = alphak1/(alphak1+betak1);
1018 //-------------------------------------------------
1019 ✗ value_ik1 = gK1*xK1inf*std::sqrt(Ko/5.4)*(value_uExtrap-EK);
1020 //-------------------------------------------------
1021
1022
1023 //*******************************************
1024 //Hyperpolarization activated funny current
1025 //------------------------------
1026 ✗ value_if = gf*valuem_solEDO[12]*(value_uExtrap-Ef);
1027 //------------------------------
1028
1029
1030
1031 //********************
1032 //Na+K+ pump current
1033 //
1034 ✗ const double num1 = PNaK*Ko*valuem_solEDO[16]/(Ko+Kmk);
1035 ✗ const double num2 = num1/(valuem_solEDO[16]+KmNa);
1036 //--------------------------------------------------------------------------------
1037 ✗ value_iNaK = num2/(1.+0.1245*std::exp(-0.1*value_uExtrap*F/(R*T))+0.0353*std::exp(-value_uExtrap*F/(R*T)));
1038 //--------------------------------------------------------------------------------
1039
1040
1041
1042 //***************************
1043 //Na+Ca2+ exchanger current
1044 //
1045 ✗ const double prodnum1 = std::pow(valuem_solEDO[16],3)*Cao*std::exp(gamma_value*value_uExtrap*F/(R*T));
1046 ✗ const double prodnum2 = alpha*std::pow(Nao,3)*valuem_solEDO[14]*std::exp((gamma_value-1.)*value_uExtrap*F/(R*T));
1047 ✗ const double numINaCa = KNaCa*(prodnum1-prodnum2);
1048 ✗ const double prodden1 = std::pow(KmNai,3)+std::pow(Nao,3);
1049 ✗ const double prodden2 = KmCa+Cao;
1050 ✗ const double prodden3 = 1.+Ksat*std::exp((gamma_value-1.)*value_uExtrap*F/(R*T));
1051 ✗ const double denINaCa = prodden1*prodden2*prodden3;
1052 //-------------------------------
1053 ✗ value_iNaCa = numINaCa/denINaCa;
1054 //-------------------------------
1055
1056
1057
1058 //*******************
1059 //Ca2+ pump current
1060 //--------------------------------------
1061 ✗ value_iPCa = gpCa*valuem_solEDO[14]/(valuem_solEDO[14]+KpCa);
1062 //--------------------------------------
1063
1064
1065
1066 //*********************
1067 //Background currents
1068 //--------------------------------------
1069 ✗ double value_ibNa=gbNa*(value_uExtrap-ENa);
1070 ✗ double value_ibCa=gbCa*(value_uExtrap-ECa);
1071 //--------------------------------------
1072
1073
1074 //********************
1075 //Membrane potential
1076 //********************
1077
1078 ✗ value_iNaCa = valuem_solEDO[18];
1079 ✗ value_iPCa = valuem_solEDO[19];
1080 ✗ value_ibCa = valuem_solEDO[20];
1081
1082 //Currents Ito,IKr,IKs,IK1,If
1083 ✗ m_vecSolEDO[26].setValue(pos,value_ito, INSERT_VALUES);
1084 ✗ m_vecSolEDO[27].setValue(pos,value_ikr, INSERT_VALUES);
1085 ✗ m_vecSolEDO[28].setValue(pos,value_iks, INSERT_VALUES);
1086 ✗ m_vecSolEDO[29].setValue(pos,value_ik1, INSERT_VALUES);
1087 ✗ m_vecSolEDO[30].setValue(pos,value_if, INSERT_VALUES);
1088
1089 ✗ value_ion = -(value_ik1+value_ito+value_ikr+value_iks+value_iCaL+value_iNaK+value_iNa+value_iNaCa+value_iPCa+value_if+value_ibNa+value_ibCa);
1090 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
1091
1092 }
1093
1094 ✗ m_ion.assembly();
1095
1096 ✗ m_vecSolEDO[26].assembly();
1097 ✗ m_vecSolEDO[27].assembly();
1098 ✗ m_vecSolEDO[28].assembly();
1099 ✗ m_vecSolEDO[29].assembly();
1100 ✗ m_vecSolEDO[30].assembly();
1101
1102
1103 //********************
1104 //********************
1105
1106
1107 }
1108
1109 //****************************************************************************************************
1110 //****************************************************************************************************
1111
1112 //Gates functions
1113
1114 //****************************************************************************************************
1115 //****************************************************************************************************
1116
1117
1118 ✗ void PaciSolver::hgate(double& hinf,double& tauh, double u) {
1119
1120 ✗ hinf = 1./std::sqrt(1.+std::exp(1000.*(u+72.1e-3)/5.7));
1121 double alphah,betah;
1122
1123 ✗ if(u<-40e-3) {
1124 ✗ alphah = 0.057*std::exp(-1000.*(u+80e-3)/6.8);
1125 ✗ betah = 2.7*std::exp(0.079*u*1000.)+(3.1e5)*std::exp(0.3485*u*1000.);
1126 ✗ tauh = 1.5/(alphah+betah);
1127 } else {
1128 ✗ alphah = 0.;
1129 ✗ betah = 0.77/(0.13*(1.+std::exp(-1000.*(u+10.66e-3)/11.1)));
1130 ✗ tauh = 2.542;
1131 }
1132 }
1133
1134 ✗ void PaciSolver::jgate(double& jinf,double& tauj, double u) {
1135
1136 ✗ jinf = 1./std::sqrt(1.+std::exp(1000.*(u+72.1e-3)/5.7));
1137 double alphaj,betaj;
1138
1139 ✗ if(u<-40e-3) {
1140 ✗ const double numalphaj = 1000.*(u+37.78e-3)*(-25428.*std::exp(1000.*0.2444*u)-(6.948e-6)*std::exp(-0.04391*u*1000.));
1141 ✗ const double denalphaj = 1./(1.+std::exp(0.311*(u+79.23e-3)*1000.));
1142
1143 ✗ alphaj = numalphaj/denalphaj;
1144 ✗ betaj = (0.02424*std::exp(-0.01052*u*1000.))/(1.+std::exp(-0.1378*(u+40.14e-3)*1000.));
1145 } else {
1146 ✗ alphaj = 0.;
1147 ✗ betaj = (0.6*std::exp(0.057*u*1000.))/(1.+std::exp(-0.1*(u+32.e-3)*1000.));
1148 }
1149 ✗ tauj = 7./(alphaj+betaj);
1150 }
1151
1152 ✗ void PaciSolver::mgate(double& minf,double& taum, double u) {
1153
1154 ✗ const double minfradicand = 1.+std::exp(1000.*(-u-34.1e-3)/5.9);
1155 ✗ minf = 1./std::pow(minfradicand,1./3.);
1156 ✗ const double alpham = 1./(1.+std::exp(1000.*(-u-60e-3)/5.));
1157 ✗ const double betam = (0.1/(1.+std::exp(1000.*(u+35e-3)/5.)))+(0.1/(1.+std::exp(1000.*(u-50e-3)/200.)));
1158 ✗ taum = alpham*betam;
1159 }
1160
1161 ✗ void PaciSolver::dgate(double& dinf,double& taud, double u) {
1162
1163 ✗ dinf = 1./(1.+std::exp(-1000.*(u+9.1e-3)/7.));
1164 //Atrial
1165 //dinf = 1./(1.+std::exp(-1000.*(value_uExtrap+5.986e-3)/7.));
1166 ✗ const double alphad = (1.4/(1.+std::exp(1000.*(-u-35e-3)/13.)))+0.25;
1167 ✗ const double betad = 1.4/(1.+std::exp(1000.*(u+5e-3)/5.));
1168 ✗ const double gammad = 1./(1.+std::exp(1000.*(-u+50e-3)/20.));
1169 ✗ taud = (alphad*betad)+gammad;
1170 }
1171
1172 ✗ void PaciSolver::fcagate(double& fcainf,double& taufca, double& constfca, double u, double Cai, double fCa) {
1173
1174 ✗ const double alphafCa = 1./(1.+std::pow(Cai/0.0006,8));
1175 ✗ const double betafCa = 0.1/(1.+std::exp((Cai-0.0009)/0.0001));
1176 ✗ const double gammafCa = 0.3/(1.+std::exp((Cai-0.00075)/0.0008));
1177 ✗ fcainf = (alphafCa+betafCa+gammafCa)/1.3156;
1178 ✗ taufca = 2.;//[ms]
1179 ✗ constfca = 1.;
1180 ✗ if(fcainf>fCa && u>-60e-3) {
1181 ✗ constfca = 0.;
1182 }
1183 }
1184
1185 ✗ void PaciSolver::f1gate(double& f1inf,double& tauf1, double u, double df1, double Cai) {
1186
1187 ✗ f1inf = 1./(1.+std::exp(1000.*(u+26e-3)/3.));//ventricular
1188 //Atrial
1189 //f1inf = 1./(1.+std::exp(-1000.*(value_uExtrap+25.266e-3)/3.));
1190
1191 ✗ const double sum1 = 1102.5*std::exp(-std::pow(std::pow(1000.*(u+27e-3),2.)/15.,2.));
1192 ✗ const double sum2 = 200./(1.+std::exp(1000.*(-u+13e-3)/10.));
1193 ✗ const double sum3 = 180./(1.+std::exp(1000.*(u+30e-3)/10.));
1194 ✗ tauf1 = sum1+sum2+sum3+20.;
1195 ✗ if(df1>0.) {
1196 ✗ tauf1 = tauf1*(1.+1433*(Cai-(50e-6)));
1197 }
1198 }
1199
1200 ✗ void PaciSolver::f2gate(double& f2inf,double& tauf2, double u) {
1201
1202 ✗ const double sum12 = 600.*std::exp(-std::pow(1000.*(u+25e-3),2.)/170.);
1203 ✗ const double sum22 = 31./(1.+std::exp(1000.*(-u+25e-3)/10.));
1204 ✗ const double sum32 = 16./(1.+std::exp(1000.*(u+30e-3)/10.));
1205 ✗ tauf2 = sum12+sum22+sum32;
1206 ✗ f2inf = (0.67/(1.+std::exp(1000.*(u+35e-3)/4.)))+0.33;
1207 }
1208
1209 ✗ void PaciSolver::rgate(double& rinf,double& taur, double u) {
1210 ✗ rinf = 1./(1.+std::exp(1000.*(-u+22.3e-3)/18.75));
1211 ✗ taur = (14.40516/(1.037*std::exp(1000.*0.09*(u+30.61e-3))+0.369*std::exp(-0.12*(u+23.84e-3)*1000.)))+2.75352;
1212 }
1213
1214 ✗ void PaciSolver::qgate(double& qinf,double& tauq, double u) {
1215 ✗ qinf = 1./(1.+std::exp(1000.*(u+53e-3)/13.));
1216 ✗ tauq = (39.102/(0.57*std::exp(-0.08*(u+44e-3)*1000.)+0.065*std::exp(1000.*0.1*(u+45.93e-3))))+6.06;
1217 }
1218
1219 ✗ void PaciSolver::xr1gate(double& xr1inf,double& tauxr1, double u) {
1220 ✗ const double V05 = (-((R*T)/(F*Q))*std::log(std::pow(1.+Cao/2.6,4)/(Lo*std::pow(1.+Cao/0.58,4)))-0.019);
1221 ✗ xr1inf = 1./(1.+std::exp(1000.*(V05-u)/4.9));
1222 ✗ const double alphaxr1 = 450./(1.+std::exp(1000.*(-u-45e-3)/10.));
1223 ✗ const double betaxr1 = 6./(1.+std::exp(1000.*(u+30e-3)/11.5));
1224 ✗ tauxr1 = alphaxr1*betaxr1;
1225 }
1226
1227 ✗ void PaciSolver::xr2gate(double& xr2inf,double& tauxr2, double u) {
1228 ✗ xr2inf = 1./(1.+std::exp(1000.*(u+88e-3)/50.));
1229 ✗ const double alphaxr2 = 3./(1.+std::exp(1000.*(-u-60e-3)/20.));
1230 ✗ const double betaxr2 = 1.12/(1.+std::exp(1000.*(u-60e-3)/20.));
1231 ✗ tauxr2 = alphaxr2*betaxr2;
1232 }
1233
1234 ✗ void PaciSolver::xsgate(double& xsinf,double& tauxs, double u) {
1235 ✗ xsinf = 1./(1.+std::exp(1000.*(-u-20e-3)/16.));
1236 ✗ const double alphaxs = 1100./std::sqrt((1.+std::exp(1000.*(-u-10e-3)/6.)));
1237 ✗ const double betaxs = 1./(1.+std::exp(1000.*(u-60e-3)/20.));
1238 ✗ tauxs = alphaxs*betaxs;
1239 }
1240
1241 ✗ void PaciSolver::xfgate(double& xfinf,double& tauf, double u) {
1242 ✗ xfinf = 1./(1.+std::exp(1000.*(u+77.85e-3)/5.));
1243 ✗ tauf = 1900./(1.+std::exp(1000.*(u+15e-3)/10.));
1244 }
1245
1246 ✗ void PaciSolver::ggate(double& ginf,double& taug,double& constg, double u, double Cai, double g) {
1247
1248 ✗ if(Cai<=0.00035) {
1249 ✗ ginf = 1./(1.+std::pow(Cai/0.00035,6));
1250 } else {
1251 ✗ ginf = 1./(1.+std::pow(Cai/0.00035,16));
1252 }
1253
1254 ✗ taug = 2.; //[ms]
1255 ✗ constg = 1.;
1256
1257 ✗ if(ginf>g && u>-60e-3) {
1258 ✗ constg = 0.;
1259 }
1260 }
1261
1262 //For non spontaneous AP, only used for INa
1263
1264 ✗ void PaciSolver::mgateORd(double& minf, double& taum, double u) {
1265 ✗ minf = 1./(1.+std::exp(-1000.*(u+39.57e-3)/9.871));
1266 ✗ const double exp1 = 6.765*std::exp(1000.*(u+11.64e-3)/34.77);
1267 ✗ const double exp2 = 8.552*std::exp(-1000.*(u+77.42e-3)/5.955);
1268 ✗ taum = 1./(exp1+exp2);
1269 }
1270
1271 ✗ void PaciSolver::jgateORd(double& jinf, double& tauj, double u) {
1272 ✗ jinf = 1./(1.+std::exp(1000.*(u+82.9e-3)/6.086));
1273 ✗ const double exp1 = 0.02136*std::exp(-1000.*(u+100.6e-3)/8.281);
1274 ✗ const double exp2 = 0.3052*std::exp(1000.*(u+0.9941e-3)/38.45);
1275 ✗ tauj = 2.038+1./(exp1+exp2);
1276 }
1277
1278 ✗ void PaciSolver::hslowgateORd(double& hslowinf, double& tauhslow, double u) {
1279 ✗ hslowinf = 1./(1.+std::exp(1000.*(u+82.9e-3)/6.086));
1280 ✗ const double exp1 = 0.009764*std::exp(-1000.*(u+17.95e-3)/28.05);
1281 ✗ const double exp2 = 0.3343*std::exp(1000.*(u+5.730e-3)/56.66);
1282 ✗ tauhslow = 1./(exp1+exp2);
1283 }
1284
1285 ✗ void PaciSolver::hfastgateORd(double& hfastinf, double& tauhfast, double u) {
1286 ✗ hfastinf = 1./(1.+std::exp(1000.*(u+82.9e-3)/6.086));
1287 ✗ const double exp1 = 1.432e-5*std::exp(-1000.*(u+1.196e-3)/6.285);
1288 ✗ const double exp2 = 6.149*std::exp(1000.*(u+0.5096e-3)/20.27);
1289 ✗ tauhfast = 1./(exp1+exp2);
1290
1291 }
1292
1293 ✗ void PaciSolver::hlgateORd(double& hlinf, double& tauhl, double u) {
1294 ✗ hlinf = 1./(1.+std::exp(1000.*(u+87.61e-3)/7.488));
1295 ✗ tauhl = 200.;
1296 }
1297
1298 ✗ void PaciSolver::mlgateORd(double& mlinf, double& tauml, double u) {
1299 ✗ mlinf = 1./(1.+std::exp(-1000.*(u+42.85e-3)/5.264));
1300 ✗ const double exp1 = 6.765*std::exp(1000.*(u+11.64e-3)/34.77);
1301 ✗ const double exp2 = 8.552*std::exp(-1000.*(u+77.42e-3)/5.955);
1302 ✗ tauml = 1./(exp1+exp2);
1303 }
1304
1305
1306
1307 //****************************************************************************************************
1308 //****************************************************************************************************
1309
1310 // CVODE: y' = f(y)
1311
1312 //****************************************************************************************************
1313 //****************************************************************************************************
1314
1315
1316 #ifdef FELISCE_WITH_SUNDIALS
1317
1318 //for each node
1319 int PaciSolver::M_f(realtype t, N_Vector y, N_Vector ydot, void *f_data) {
1320
1321 //-----------------------------------------------------//
1322 // RHS
1323 //-----------------------------------------------------//
1324
1325
1326 UserData* data = static_cast< UserData* >(f_data);
1327 const int pos = data->pos;
1328 const double UExtrap_forM_f = data->M_u_extrap[pos];
1329 const double value_df1 = data->value_df1[pos];
1330
1331 double hinf,tauh;
1332 PaciSolver::hgate(hinf,tauh,UExtrap_forM_f);
1333
1334 double jinf,tauj;
1335 PaciSolver::jgate(jinf,tauj,UExtrap_forM_f);
1336
1337 double minf,taum;
1338 PaciSolver::mgate(minf,taum,UExtrap_forM_f);
1339
1340 double dinf,taud;
1341 PaciSolver::dgate(dinf,taud,UExtrap_forM_f);
1342
1343 double fcainf,taufca,constfca;
1344 PaciSolver::fcagate(fcainf,taufca,constfca,UExtrap_forM_f,Ith(y,15),Ith(y,5));
1345
1346 double f1inf,tauf1;
1347 PaciSolver::f1gate(f1inf,tauf1,UExtrap_forM_f,value_df1,Ith(y,15));
1348
1349 double f2inf,tauf2;
1350 PaciSolver::f2gate(f2inf,tauf2,UExtrap_forM_f);
1351
1352 double rinf,taur;
1353 PaciSolver::rgate(rinf,taur,UExtrap_forM_f);
1354
1355 double qinf,tauq;
1356 PaciSolver::qgate(qinf,tauq,UExtrap_forM_f);
1357
1358 double xr1inf,tauxr1;
1359 PaciSolver::xr1gate(xr1inf,tauxr1,UExtrap_forM_f);
1360
1361 double xr2inf,tauxr2;
1362 PaciSolver::xr2gate(xr2inf,tauxr2,UExtrap_forM_f);
1363
1364 double xsinf,tauxs;
1365 PaciSolver::xsgate(xsinf,tauxs,UExtrap_forM_f);
1366
1367 double xfinf,tauf;
1368 PaciSolver::xfgate(xfinf,tauf,UExtrap_forM_f);
1369
1370 double ginf,taug,constg;
1371 PaciSolver::ggate(ginf,taug,constg,UExtrap_forM_f,Ith(y,15),Ith(y,14));
1372
1373
1374 //! Non spontaneous (ORd)
1375 double minfORd,taumORd;
1376 PaciSolver::mgateORd(minfORd, taumORd,UExtrap_forM_f);
1377
1378 double jinfORd,taujORd;
1379 PaciSolver::jgateORd(jinfORd, taujORd,UExtrap_forM_f);
1380
1381 double hslowinf,tauhslow;
1382 PaciSolver::hslowgateORd(hslowinf, tauhslow,UExtrap_forM_f);
1383
1384 double hfastinf,tauhfast;
1385 PaciSolver::hfastgateORd(hfastinf, tauhfast,UExtrap_forM_f);
1386
1387 double hlinf,tauhl;
1388 PaciSolver::hlgateORd(hlinf, tauhl,UExtrap_forM_f);
1389
1390 double mlinf,tauml;
1391 PaciSolver::mlgateORd(mlinf, tauml,UExtrap_forM_f);
1392
1393
1394
1395 //------------------------------------
1396 // Calculations for Cai, Casr and Nai
1397 //------------------------------------
1398
1399 const double Cm = 98.7109;
1400 const double Vc = 8800.;
1401 const double Vsr = 583.73;
1402 const double gNa = 3.6712302e3;
1403 const double PNaK = 1.841424;
1404 const double KNaCa = 4900.;
1405 const double Vmaxup = 0.56064;
1406
1407
1408 const double Ileak = Vleak*(Ith(y,16)-Ith(y,15));
1409 const double Iup = Vmaxup/(1.+std::pow(Kup/Ith(y,15),2));
1410 double Irel = (crel+(arel*std::pow(Ith(y,16),2))/(std::pow(brel,2)+std::pow(Ith(y,16),2)))*Ith(y,4)*Ith(y,14);
1411 Irel = Irel*0.0411;//ventricular
1412
1413 //Atrial
1414 //Irel = value_irel*0.0556;
1415
1416 //iNaK
1417 const double num1 = PNaK*Ko*Ith(y,17)/(Ko+Kmk);
1418 const double num2 = num1/(Ith(y,17)+KmNa);
1419 const double INaK = num2/(1.+0.1245*std::exp(-0.1*UExtrap_forM_f*F/(R*T))+0.0353*std::exp(-UExtrap_forM_f*F/(R*T)));
1420
1421
1422 //iNaCa
1423 const double prodnum1 = std::pow(Ith(y,17),3)*Cao*std::exp(gamma_value*UExtrap_forM_f*F/(R*T));
1424 const double prodnum2 = alpha*std::pow(Nao,3)*Ith(y,15)*std::exp((gamma_value-1.)*UExtrap_forM_f*F/(R*T));
1425 const double numINaCa = KNaCa*(prodnum1-prodnum2);
1426 const double prodden1 = std::pow(KmNai,3)+std::pow(Nao,3);
1427 const double prodden2 = KmCa+Cao;
1428 const double prodden3 = 1.+Ksat*std::exp((gamma_value-1.)*UExtrap_forM_f*F/(R*T));
1429 const double denINaCa = prodden1*prodden2*prodden3;
1430 const double INaCa = numINaCa/denINaCa;
1431
1432
1433 //iCaL
1434 const double term1 = (gCaL*4.*UExtrap_forM_f*F*F)/(R*T);
1435 const double term2 = (Ith(y,15)*std::exp((2.*UExtrap_forM_f*F)/(R*T))-0.341*Cao)/(std::expm1((2.*UExtrap_forM_f*F)/(R*T)));
1436 const double ICaL = term1*term2*Ith(y,4)*Ith(y,5)*Ith(y,6)*Ith(y,7);
1437
1438
1439 //ipCa
1440 const double IPCa = gpCa*Ith(y,15)/(Ith(y,15)+KpCa);
1441
1442
1443 const double ENa = (R*T/F)*std::log(Nao/Ith(y,17));
1444 const double EKs = (R*T/F)*std::log((Ko+Pkna*Nao)/(Ki+Pkna*Ith(y,17)));
1445 const double ECa = 0.5*(R*T/F)*std::log(Cao/Ith(y,15));
1446
1447
1448 //iNa
1449 //const double m3 = std::pow(Ith(y,3),3.);
1450 //const double INa=gNa*m3*Ith(y,1)*Ith(y,2)*(UExtrap_forM_f-ENa);
1451
1452 //! For non spontaneous AP
1453 const double hORd = 0.99*Ith(y,21)+0.01*Ith(y,20);
1454 const double m3 = std::pow(Ith(y,18),3.);
1455 const double iNafast = 75.*(UExtrap_forM_f-ENa)*m3*Ith(y,19)*hORd;
1456 const double iNalate = 0.0075*(UExtrap_forM_f-ENa)*Ith(y,23)*Ith(y,22);
1457
1458 const double INa = iNafast+iNalate;
1459
1460
1461 //ibNa-ibCa
1462 const double IbNa=gbNa*(UExtrap_forM_f-ENa);
1463 const double IbCa=gbCa*(UExtrap_forM_f-ECa);
1464
1465
1466 const double Caibufc = 1./(1.+(Bufc*Kbufc)/std::pow(Ith(y,15)+Kbufc,2));
1467 const double Casrbufsr = 1./(1.+(Bufsr*Kbufsr)/std::pow(Ith(y,16)+Kbufsr,2));
1468
1469 //Cai-Casr-Nai
1470 double cstteCai = Caibufc*(Ileak-Iup+Irel-Cm*(ICaL+IbCa+IPCa-2.*INaCa)/(2.*Vc*F));
1471 double cstteCasr = (Casrbufsr*Vc)*(Iup-(Irel+Ileak))/Vsr;
1472 double cstteNai = -Cm*(INa+IbNa+3.*INaK+3.*INaCa)/(F*Vc);
1473
1474
1475 //Add currents in data
1476 data->vINaCa[pos] = INaCa;
1477 data->vIPCa[pos] = IPCa;
1478 data->vIbCa[pos] = IbCa;
1479 data->vENa[pos] = ENa;
1480 data->vEKs[pos] = EKs;
1481 data->vECa[pos] = ECa;
1482 data->vINa[pos] = INa;
1483 data->vICaL[pos] = ICaL;
1484 data->vINaK[pos] = INaK;
1485 data->vIrel[pos] = Irel;
1486 data->vIup[pos] = Iup;
1487 data->vIleak[pos] = Ileak;
1488 data->vIbNa[pos] = IbNa;
1489
1490
1491 //-----------------------------------------------------//
1492 // Derivates
1493 //-----------------------------------------------------//
1494
1495 Ith(ydot,1) = (hinf-Ith(y,1))/tauh;//h
1496 Ith(ydot,2) = (jinf-Ith(y,2))/tauj;//j
1497 Ith(ydot,3) = (minf-Ith(y,3))/taum;//m
1498 Ith(ydot,4) = (dinf-Ith(y,4))/taud;//d
1499 Ith(ydot,5) = constfca*(fcainf-Ith(y,5))/taufca;//fca
1500 Ith(ydot,6) = (f1inf-Ith(y,6))/tauf1;//f1
1501 Ith(ydot,7) = (f2inf-Ith(y,7))/tauf2;//f2
1502 Ith(ydot,8) = (rinf-Ith(y,8))/taur;//r
1503 Ith(ydot,9) = (qinf-Ith(y,9))/tauq;//q
1504 Ith(ydot,10) = (xr1inf-Ith(y,10))/tauxr1;//xr1
1505 Ith(ydot,11) = (xr2inf-Ith(y,11))/tauxr2;//xr2
1506 Ith(ydot,12) = (xsinf-Ith(y,12))/tauxs;//xs
1507 Ith(ydot,13) = (xfinf-Ith(y,13))/tauf;//xf
1508 Ith(ydot,14) = constg*(ginf-Ith(y,14))/taug;//g
1509 Ith(ydot,15) = cstteCai;//Cai
1510 Ith(ydot,16) = cstteCasr;//Casr
1511 Ith(ydot,17) = cstteNai;//Nai
1512
1513 //! For non spontaneous AP
1514 Ith(ydot,18) = (minfORd-Ith(y,18))/taumORd;//mORd
1515 Ith(ydot,19) = (jinfORd-Ith(y,19))/taujORd;//jORd
1516 Ith(ydot,20) = (hslowinf-Ith(y,20))/tauhslow;//hslow
1517 Ith(ydot,21) = (hfastinf-Ith(y,21))/tauhfast;//hfast
1518 Ith(ydot,22) = (hlinf-Ith(y,22))/tauhl;//hl
1519 Ith(ydot,23) = (mlinf-Ith(y,23))/tauml;//ml
1520
1521 return(0);
1522 }
1523
1524
1525
1526
1527
1528 static int ewt(N_Vector y, N_Vector w, void *user_data) {
1529 int i;
1530 realtype yy, ww, rtol, atol[17];
1531
1532 rtol = RTOL;
1533 atol[0] = ATOL8;
1534 atol[1] = ATOL8;
1535 atol[2] = ATOL8;
1536 atol[3] = ATOL8;
1537 atol[4] = ATOL8;
1538 atol[5] = ATOL8;
1539 atol[6] = ATOL8;
1540 atol[7] = ATOL8;
1541 atol[8] = ATOL8;
1542 atol[9] = ATOL8;
1543 atol[10] = ATOL8;
1544 atol[11] = ATOL8;
1545 atol[12] = ATOL8;
1546 atol[13] = ATOL8;
1547 atol[14] = ATOL8;
1548 atol[15] = ATOL8;
1549 atol[16] = ATOL8;
1550 atol[17] = ATOL8;
1551
1552
1553 for (i=1; i<=17; i++) {
1554 yy = Ith(y,i);
1555 ww = rtol * ABS(yy) + atol[i-1];
1556 if (ww <= 0.0) return (-1);
1557 Ith(w,i) = 1.0/ww;
1558 }
1559
1560 return(0);
1561 }
1562
1563
1564
1565
1566
1567 void PaciSolver::PrintOutput(realtype t, realtype y1, realtype y2, realtype y3, realtype y4, realtype y5, realtype y6) {
1568 #if defined(SUNDIALS_EXTENDED_PRECISION)
1569 printf("At t = %0.4Le y =%14.6Le %14.6Le %14.6Le %14.6Le %14.6Le %14.6Le\n", t, y1, y2, y3, y4, y5, y6);
1570 #elif defined(SUNDIALS_DOUBLE_PRECISION)
1571 printf("At t = %0.4le y =%14.6le %14.6le %14.6le %14.6le %14.6le %14.6le\n", t, y1, y2, y3, y4, y5, y6);
1572 #else
1573 printf("At t = %0.4e y =%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", t, y1, y2, y3, y4, y5, y6);
1574 #endif
1575
1576 return;
1577 }
1578
1579 #endif
1580
1581
1582
1583
1584 }
1585
1586
1587
1588
1589 /*
1590 Stock std::cout and others informations
1591
1592
1593
1594 if(pos==700){
1595 std::cout << "Iion: " << value_ion << std::endl;
1596 }
1597 */
1598 /*
1599 if(pos==700){
1600 std::cout << " " << std::endl;
1601 std::cout << "ik1: " << value_ik1 << std::endl;
1602 std::cout << "ito: " << value_ito << std::endl;
1603 std::cout << "ikr: " << value_ikr << std::endl;
1604 std::cout << "iks: " << value_iks << std::endl;
1605 std::cout << "iCaL: " << value_iCaL << std::endl;
1606 std::cout << "iNaK: " << value_iNaK << std::endl;
1607 std::cout << "iNa: " << value_iNa << std::endl;
1608 std::cout << "iNaCa: " << value_iNaCa << std::endl;
1609 std::cout << "iPCa: " << value_iPCa << std::endl;
1610 std::cout << "if: " << value_if << std::endl;
1611 std::cout << "ibNa: " << value_ibNa << std::endl;
1612 std::cout << "ibCa: " << value_ibCa << std::endl;
1613 std::cout << " " << std::endl;
1614 std::cout << "Iion: " << value_ion << " sans Cm: " << value_ion*Cm << std::endl;
1615 std::cout << "u: " << value_uExtrap << std::endl;
1616 std::cout << " " << std::endl;
1617 }
1618 */
1619
1620 //const double coefCm = FelisceParam::instance().Cm;//Test pour recuperer la valeur Cm du data
1621 /*
1622 if(coefCm == 98.7109){
1623 //Ventricular
1624 Vc = 8800.;
1625 Vsr = 583.73;
1626 gNa = 3.6712302e3;
1627 gto = 29.9038;
1628 gk1 = 28.1492;
1629 PNaK = 1.841424;
1630 KNaCa = 4900.;
1631 Vmaxup = 0.56064;
1632 }
1633 else if(coefCm == 78.6671){
1634 //Atrial
1635 Vc = 7012.;
1636 Vsr = 465.20;
1637 gNa = 6.646185e3;
1638 gto = 59.8077;
1639 gk1 = 19.1925;
1640 PNaK = 1.4731392;
1641 KNaCa = 2450.;
1642 Vmaxup = 0.22;
1643 }
1644 */
1645
1646
1647 /*
1648 if(pos==700){
1649 std::cout << "avant: " << std::endl;
1650 PaciSolver::PrintOutput(time-dt, Ith(initvalue,1), Ith(initvalue,2), Ith(initvalue,3), Ith(initvalue,4), Ith(initvalue,5), Ith(initvalue,6));
1651 PaciSolver::PrintOutput(time-dt, Ith(initvalue,7), Ith(initvalue,8), Ith(initvalue,9), Ith(initvalue,10), Ith(initvalue,11), Ith(initvalue,12));
1652 PaciSolver::PrintOutput(time-dt, Ith(initvalue,13), Ith(initvalue,14), Ith(initvalue,15), Ith(initvalue,16), Ith(initvalue,17), Ith(initvalue,18));
1653 std::cout << " " << std::endl;
1654 }
1655 */
1656
1657
1658