Chaste Commit::1fd4e48e3990e67db148bc1bc4cf6991a0049d0c
CorriasBuistSMCModified.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36
37#include <cmath>
38#include <cassert>
39#include <memory>
40#include "Exception.hpp"
41#include "OdeSystemInformation.hpp"
42#include "CorriasBuistSMCModified.hpp"
43#include "HeartConfig.hpp"
44
45
46 CorriasBuistSMCModified::CorriasBuistSMCModified(boost::shared_ptr<AbstractIvpOdeSolver> pSolver, boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
47 : AbstractCardiacCell(
48 pSolver,
49 14,
50 0,
51 pIntracellularStimulus)
52 {
53 mpSystemInfo = OdeSystemInformation<CorriasBuistSMCModified>::Instance();
54 mScaleFactorCarbonMonoxide = 1.0; // initialise to 1 --> no effect
55 mFakeIccStimulusPresent = true;//by default we use the fake ICC stimulus
56
57 /* parameters */
58 Cm = 77.0*1e-6;// 77 pF --> microF
59
60 Asurf_in_cm_square = Cm / HeartConfig::Instance()->GetCapacitance();
61 Asurf = Asurf_in_cm_square / 0.01;//cm2 --> mm2 /*cell surface area (mm2)*/
62
63 VolCell = 3.5e-6; /*cell volume (mm3)*/
64 hCa = 2.014e-4; /*conc for half inactivation of fCa */
65 sCa = 1.31e-05; /*slope factor for inactivation of fCa */
66
67 /* concentrations */
68 Ki = 164.0; /*intra K conc (mM)*/
69 Nai = 10.0; /*intra Na conc (mM)*/
70 ACh = 1e-05; /*acetylcholine conc (mM)*/
71 CaiRest = 0.9e-04; /*baseline Ca conc (mM)*/
72
73 /* maximum conductances*/
74 gLVA_max = 2.33766E-05; /*max conductance of ILVA*/ // (0.18 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
75 gCaL_max = 0.008441558; /*max conductance of ICaL*/ // (65.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
76 gBK_max = 0.005935065; /*max conductance of IBK)*/ // (45.7 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
77 gKb_max = 1.87013E-06; /*max conductance of IKb*/ // (0.0144 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
78 gKA_max = 0.001168831; /*max conductance of IKA*/ // (9.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
79 gKr_max = 0.004545455; /*max conductance of IKr*/ // (35.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
80 gNa_max = 0.00038961; /*max conductance of INa*/ // (3.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
81 gnsCC_max = 0.006493506; /*max conductance of InsCC*/ // (50.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
82 gcouple = 1.3*1e-6 / Asurf; /* coupling conductance bewteen fake ICC and SMC*/ // 1.3 nS * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
83
84 JCaExt_max = 0.31705; /*max flux of CaSR (mM/ms)*/
85
86 /* Temperature corrections */
87 Q10Ca = 2.1; /*(dim)*/
88 Q10K = 1.5; /*(dim)*/ //1.365
89 Q10Na = 2.45; /*(dim)*/
90 Texp = 297.0; /*(degK)*/
91
92 Ca_o = 2.5; // mM
93 K_o = 7.0; // mM
94 Na_o = 137.0; // mM
95
96 /* Nernst parameters */
97 R = 8314.4; // pJ/nmol/K
98 T = 310.0; // degK
99 F = 96484.6; // nC/nmol
100 FoRT = 0.03743; // 1/mV
101 RToF = 26.7137; // mV
102
103 T_correct_Ca = pow(Q10Ca,(T-Texp)/10.0);/*temperature correction for Ca (dim)*/
104 T_correct_K = pow(Q10K,(T-Texp)/10.0); /*temperature correction for K (dim)*/
105 T_correct_Na = pow(Q10Na,(T-Texp)/10.0);/*temperature correction for Na (dim)*/
106 T_correct_gBK = gBK_max + 1.1*(T-Texp)*1e-6/Asurf; /*temperature correction for gBK*/ // (nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
107
108 /* Nernst potentials */
109 EK = RToF*log(K_o/Ki); /*Nernst potential for K (mV)*/
110 ENa = RToF*log(Na_o/Nai); /*Nernst potential for Na (mV)*/
111 EnsCC = -28.0; /*Nernst potential for nsCC (mV)*/
112
113 Init();
114
115 }
116
120
123
124 void CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor(double scaleFactor)
125 {
126 mScaleFactorCarbonMonoxide = scaleFactor;
127 }
128
133
138
143
144 double CorriasBuistSMCModified::GetIIonic(const std::vector<double>* pStateVariables)
145 {
146 if (!pStateVariables) pStateVariables = &rGetStateVariables();
147 const std::vector<double>& rY = *pStateVariables;
148
149 double ECa = 0.5*RToF*log(Ca_o/rY[13]);
150
151 /* inward sodium current */
152 double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);
153
154 /* L-type calcium current */
155 double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);
156
157 /* low voltage activated (T-type) calcium current */
158 double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);
159
160 /* large conductance calcium activated potassium current */
161 double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
162 double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);
163
164 /* delayed rectifier potassium current */
165 double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);
166
167 /* A-type potassium current */
168 double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);
169
170 /* background (leakage) potassium current */
171 double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);
172
173 /* non-specific cation current */
174 double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
175 double rACh_nsCC = 1.0/(1.0+(0.01/ACh));
176 double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);
177
178 /* phenomenological calcium extrusion current */
179 double JCaExt = JCaExt_max*pow(rY[13],1.34);
180
181 //i_ionic_in microA/mm2
182 double i_ionic = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);
183
184 assert(!std::isnan(i_ionic));
188 return i_ionic / 0.01;
189 }
190
191 void CorriasBuistSMCModified::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
192 {
193
194 // index [0] = -69.8; // Vm (mV)
195 // index [1] = 5.0e-6; // d_CaL (dim)
196 // index [2] = 0.953; // f_CaL (dim)
197 // index [3] = 1.0; // fCa_CaL (dim)
198 // index [4] = 0.0202; // d_LVA (dim)
199 // index [5] = 1.0; // f_LVA (dim)
200 // index [6] = 0.0086; // d_Na (dim)
201 // index [7] = 0.061; // f_Na (dim)
202 // index [8] = 0.096; // m_nsCC (dim)
203 // index [9] = 1.92e-4; // xr1 (dim)
204 // index [10] = 0.812; // xr2 (dim)
205 // index [11] = 0.00415; // xa1 (dim)
206 // index [12] = 0.716; // xa2 (dim)
207 // index [13] = 0.9e-04; // Cai (mM)
208
209 double ECa = 0.5*RToF*log(Ca_o/rY[13]);
210
211 /* inward sodium current */
212 double inf_d_Na = 1.0/(1.0+exp(-(rY[0]+47.0)/4.8));
213 double tau_d_Na = (0.44-0.017*rY[0])*T_correct_Na;
214 double inf_f_Na = 1.0/(1.0+exp((rY[0]+78.0)/3.0));
215 double tau_f_Na = (5.5-0.25*rY[0])*T_correct_Na;
216
217 double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);
218
219 /* L-type calcium current */
220 double inf_d_CaL = 1.0/(1.0+exp(-(rY[0]+17.0)/4.3));
221 double tau_d_CaL = 0.47*T_correct_Ca;
222
223 double inf_f_CaL = 1.0/(1.0+exp((rY[0]+43.0)/8.9));
224 double tau_f_CaL = 86.0*T_correct_Ca;
225
226 double inf_fCa_CaL = 1.0-(1.0/(1.0+exp(-((rY[13]-CaiRest)-hCa)/sCa)));
227 double tau_fCa_CaL = 2.0*T_correct_Ca;
228 double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);
229
230 /* low voltage activated (T-type) calcium current */
231 double inf_d_LVA = 1.0/(1.0+exp(-(rY[0]+27.5)/10.9));
232 double tau_d_LVA = 3.0*T_correct_Ca;
233
234 double inf_f_LVA = 1.0/(1.0+exp((rY[0]+15.8)/7.0));
235 double tau_f_LVA = 7.58*exp(rY[0]*0.00817)*T_correct_Ca;
236
237 double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);
238
239 /* large conductance calcium activated potassium current */
240 double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
241 double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);
242
243 /* delayed rectifier potassium current */
244 double inf_xr1 = 1.0/(1.0+exp(-(rY[0]+27.0)/5.0));
245 double tau_xr1 = 80.0*T_correct_K;
246
247 double inf_xr2 = 0.2+0.8/(1.0+exp((rY[0]+58.0)/10.0));
248 double tau_xr2 = (-707.0+1481.0*exp((rY[0]+36.0)/95.0))*T_correct_K;
249
250 double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);
251
252 /* A-type potassium current */
253 double inf_xa1 = 1.0/(1.0+exp(-(rY[0]+26.5)/7.9));
254 double tau_xa1 = (31.8+175.0*exp(-0.5*pow(((rY[0]+44.4)/22.3),2.0)))*T_correct_K;
255
256 double inf_xa2 = 0.1+0.9/(1.0+exp((rY[0]+65.0)/6.2));
257 double tau_xa2 = 90.0*T_correct_K;
258
259 double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);
260
261 /* background (leakage) potassium current */
262 double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);
263
264 /* non-specific cation current */
265 double inf_m_nsCC = 1.0/(1.0+exp(-(rY[0]+25.0)/20.0));
266 double tau_m_nsCC = 150.0/(1.0+exp(-(rY[0]+66.0)/26.0));
267 double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
268 double rACh_nsCC = 1.0/(1.0+(0.01/ACh));
269
270 double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);
271
272 /* phenomenological calcium extrusion current */
273 double JCaExt = JCaExt_max*pow(rY[13],1.34);
274
275
276 double t_ICCplateau = 7582.0; // time_units
277 double V_decay = 37.25; // voltage_units
278 double t_ICCpeak = 98.0; // time_units
279
280 double period = 20000.0; // time_units
281 double stim_start = ((time > (period * 1.0)) && (time <= (period * 2.0))) ? (period * 1.0) : ((time > (period * 2.0)) && (time <= (period * 3.0))) ? (period * 2.0) : ((time > (period * 3.0)) && (time <= (period * 4.0))) ? (period * 3.0) : ((time > (period * 4.0)) && (time <= (period * 5.0))) ? (period * 4.0) : 0.0; // time_units
282 double local_time = time - (stim_start + t_ICCpeak); // time_units
283 double t_ICC_stimulus = 10000.0; // time_units
284 double delta_VICC = 59.0; // voltage_units
285
286 double i_stim;
287 //see whether we are running this in isolaation (and we need the fake ICC stimulus) or coupled to a real ICC model
288 if (mFakeIccStimulusPresent)
289 {
290 //for single cell simulations where we want the fake ICC stimulus in
291 i_stim = (local_time < t_ICCpeak) ? (gcouple * delta_VICC) : ((local_time >= t_ICCpeak) && (local_time <= t_ICCplateau)) ? (gcouple * delta_VICC * (1.0 / (1.0 + exp((local_time - 8000.0) / 1000.0)))) : ((local_time > t_ICCplateau) && (local_time < t_ICC_stimulus)) ? (gcouple * V_decay * (1.0 / (1.0 + exp((local_time - 8000.0) / 150.0)))) : 0.0; // current_units
292 }
293 else
294 {
295 i_stim = GetStimulus(time);//for tissue simulations with current injected into SMC
296 }
297
298 /* membrane potential */
299 double Iion = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);
300
301 double voltage_derivative;
302 if (mSetVoltageDerivativeToZero)
303 {
304 voltage_derivative = 0.0;
305 }
306 else
307 {
308 voltage_derivative = (-1.0 / 0.01) * (-i_stim + Iion);//microA/mm2---> microA/cm2
309 //std::cout<<rY[0]<<std::endl;
310 assert(!std::isnan(voltage_derivative));
311 }
312
313 rDY[0] = voltage_derivative;/* Vm */
314 rDY[1] = (inf_d_CaL - rY[1])/tau_d_CaL;
315 rDY[2] = (inf_f_CaL - rY[2])/tau_f_CaL;
316 rDY[3] = (inf_fCa_CaL - rY[3])/tau_fCa_CaL;
317 rDY[4] = (inf_d_LVA - rY[4])/tau_d_LVA;
318 rDY[5] = (inf_f_LVA - rY[5])/tau_f_LVA;
319 rDY[6] = (inf_d_Na - rY[6])/tau_d_Na;
320 rDY[7] = (inf_f_Na - rY[7])/tau_f_Na;
321 rDY[8] = (inf_m_nsCC - rY[8])/tau_m_nsCC;
322 rDY[9] = (inf_xr1 - rY[9])/tau_xr1;
323 rDY[10] = (inf_xr2 - rY[10])/tau_xr2;
324 rDY[11] = (inf_xa1 - rY[11])/tau_xa1;
325 rDY[12] = (inf_xa2 - rY[12])/tau_xa2;
326 rDY[13] = (-(ICaL+ILVA)*Asurf/(2.0*F*VolCell)-JCaExt); /* intracellular calcium *1000 M-> mM; /1000 F units*/
327 }
328
329template<>
330void OdeSystemInformation<CorriasBuistSMCModified>::Initialise(void)
331{
332 // Time units: time_units
333 //
334 this->mSystemName = "SMC_model_Martincode";
335
336 this->mVariableNames.push_back("Vm_SM");
337 this->mVariableUnits.push_back("mV");
338 this->mInitialConditions.push_back(-69.8); // Vm (mV)
339
340 this->mVariableNames.push_back("d_CaL");
341 this->mVariableUnits.push_back("dim");
342 this->mInitialConditions.push_back(5.0e-6); // d_CaL (dim));
343
344 this->mVariableNames.push_back("f_CaL");
345 this->mVariableUnits.push_back("dim");
346 this->mInitialConditions.push_back(0.953); // f_CaL (dim)
347
348 this->mVariableNames.push_back("fCa_CaL");
349 this->mVariableUnits.push_back("dim");
350 this->mInitialConditions.push_back(1.0); // fCa_CaL (dim)
351
352 this->mVariableNames.push_back("d_LVA");
353 this->mVariableUnits.push_back("dim");
354 this->mInitialConditions.push_back(0.0202); // d_LVA (dim)
355
356 this->mVariableNames.push_back("f_LVA");
357 this->mVariableUnits.push_back("dim");
358 this->mInitialConditions.push_back(1.0); // f_LVA (dim)
359
360 this->mVariableNames.push_back("d_Na");
361 this->mVariableUnits.push_back("dim");
362 this->mInitialConditions.push_back(0.0086); // d_Na (dim)
363
364 this->mVariableNames.push_back("f_Na");
365 this->mVariableUnits.push_back("dim");
366 this->mInitialConditions.push_back(0.061); // f_Na (dim)
367
368 this->mVariableNames.push_back("m_nsCC");
369 this->mVariableUnits.push_back("dim");
370 this->mInitialConditions.push_back(0.096); // m_nsCC (dim)
371
372 this->mVariableNames.push_back("xr1");
373 this->mVariableUnits.push_back("dim");
374 this->mInitialConditions.push_back(1.92e-4); // xr1 (dim)
375
376 this->mVariableNames.push_back("xr2");
377 this->mVariableUnits.push_back("dim");
378 this->mInitialConditions.push_back(0.812); // xr2 (dim)
379
380 this->mVariableNames.push_back("xa1");
381 this->mVariableUnits.push_back("dim");
382 this->mInitialConditions.push_back(0.00415); // xa1 (dim)
383
384 this->mVariableNames.push_back("xa2");
385 this->mVariableUnits.push_back("dim");
386 this->mInitialConditions.push_back(0.716); // xa2 (dim)
387
388 this->mVariableNames.push_back("Cai");
389 this->mVariableUnits.push_back("mM");
390 this->mInitialConditions.push_back(0.9e-04); // Cai (mM)
391
392 this->mInitialised = true;
393}
394
395
396// Serialization for Boost >= 1.36
397#include "SerializationExportWrapperForCpp.hpp"
398CHASTE_CLASS_EXPORT(CorriasBuistSMCModified)
CHASTE_CLASS_EXPORT
#define CHASTE_CLASS_EXPORT(T)
Definition SerializationExportWrapper.hpp:347
AbstractUntemplatedParameterisedSystem::mpSystemInfo
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
Definition AbstractUntemplatedParameterisedSystem.hpp:301
CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor
double SetCarbonMonoxideScaleFactor()
CorriasBuistSMCModified::GetIIonic
double GetIIonic(const std::vector< double > *pStateVariables=NULL)
Definition CorriasBuistSMCModified.cpp:144
CorriasBuistSMCModified::CorriasBuistSMCModified
CorriasBuistSMCModified(boost::shared_ptr< AbstractIvpOdeSolver > pSolver, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
Definition CorriasBuistSMCModified.cpp:46
CorriasBuistSMCModified::SetFakeIccStimulusPresent
void SetFakeIccStimulusPresent(bool present)
Definition CorriasBuistSMCModified.cpp:129
CorriasBuistSMCModified::EvaluateYDerivatives
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
Definition CorriasBuistSMCModified.cpp:191
HeartConfig::GetCapacitance
double GetCapacitance() const
Definition HeartConfig.cpp:1575
HeartConfig::Instance
static HeartConfig * Instance()
Definition HeartConfig.cpp:208
OdeSystemInformation::Instance
static boost::shared_ptr< OdeSystemInformation< ODE_SYSTEM > > Instance()
Definition OdeSystemInformation.hpp:114