Chaste  Release::2018.1
CorriasBuistSMCModified.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
117  CorriasBuistSMCModified::~CorriasBuistSMCModified()
118  {
119  }
120 
121  void CorriasBuistSMCModified::VerifyStateVariables()
122  {}
123 
124  void CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor(double scaleFactor)
125  {
126  mScaleFactorCarbonMonoxide = scaleFactor;
127  }
128 
129  void CorriasBuistSMCModified::SetFakeIccStimulusPresent(bool present)
130  {
131  mFakeIccStimulusPresent = present;
132  }
133 
134  bool CorriasBuistSMCModified::GetFakeIccStimulusPresent()
135  {
136  return mFakeIccStimulusPresent;
137  }
138 
139  double CorriasBuistSMCModified::GetCarbonMonoxideScaleFactor()
140  {
141  return mScaleFactorCarbonMonoxide;
142  }
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 
329 template<>
330 void 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"
398 CHASTE_CLASS_EXPORT(CorriasBuistSMCModified)
OdeSystemInformation::Instance
static boost::shared_ptr< OdeSystemInformation< ODE_SYSTEM > > Instance()
Definition: OdeSystemInformation.hpp:114
CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor
double SetCarbonMonoxideScaleFactor()
CorriasBuistSMCModified::CorriasBuistSMCModified
CorriasBuistSMCModified(boost::shared_ptr< AbstractIvpOdeSolver > pSolver, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
Definition: CorriasBuistSMCModified.cpp:46
CorriasBuistSMCModified::EvaluateYDerivatives
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
Definition: CorriasBuistSMCModified.cpp:191
CorriasBuistSMCModified::SetFakeIccStimulusPresent
void SetFakeIccStimulusPresent(bool present)
Definition: CorriasBuistSMCModified.cpp:129
CorriasBuistSMCModified::GetIIonic
double GetIIonic(const std::vector< double > *pStateVariables=NULL)
Definition: CorriasBuistSMCModified.cpp:144
AbstractUntemplatedParameterisedSystem::mpSystemInfo
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
Definition: AbstractUntemplatedParameterisedSystem.hpp:301
HeartConfig::GetCapacitance
double GetCapacitance() const
Definition: HeartConfig.cpp:1575
CHASTE_CLASS_EXPORT
#define CHASTE_CLASS_EXPORT(T)
Definition: SerializationExportWrapper.hpp:347
HeartConfig::Instance
static HeartConfig * Instance()
Definition: HeartConfig.cpp:208