Nektar++
CourtemancheRamirezNattel98.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CourtemancheRamirezNattel.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Courtemanche-Ramirez-Nattel ionic atrial cell model.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iostream>
36 #include <string>
37 
38 #include <LibUtilities/BasicUtils/VmathArray.hpp>
39 #include <CardiacEPSolver/CellModels/CourtemancheRamirezNattel98.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  std::string CourtemancheRamirezNattel98::className
46  = GetCellModelFactory().RegisterCreatorFunction(
47  "CourtemancheRamirezNattel98",
48  CourtemancheRamirezNattel98::create,
49  "Ionic model of human atrial cell electrophysiology.");
50 
51  // Register cell model variants
52  std::string CourtemancheRamirezNattel98::lookupIds[2] = {
53  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
54  "Original", CourtemancheRamirezNattel98::eOriginal),
55  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
56  "AF", CourtemancheRamirezNattel98::eAF)
57  };
58 
59  // Register default variant
60  std::string CourtemancheRamirezNattel98::def =
61  LibUtilities::SessionReader::RegisterDefaultSolverInfo(
62  "CellModelVariant", "Original");
63 
64  /**
65  *
66  */
67  CourtemancheRamirezNattel98::CourtemancheRamirezNattel98(
68  const LibUtilities::SessionReaderSharedPtr& pSession,
69  const MultiRegions::ExpListSharedPtr& pField)
70  : CellModel(pSession, pField)
71  {
72  model_variant = pSession->GetSolverInfoAsEnum<
73  CourtemancheRamirezNattel98::Variants>("CellModelVariant");
74 
75  C_m = 100; // picoF
76  g_Na = 7.8; // nanoS_per_picoF
77  g_K1 = 0.09; // nanoS_per_picoF
78  g_Kr = 0.029411765;
79  g_Ks = 0.12941176;
80  g_b_Na = 0.0006744375;
81  g_b_Ca = 0.001131;
82  R = 8.3143;
83  T = 310.0;
84  F = 96.4867;
85  Na_o = 140.0; // millimolar
86  K_o = 5.4; // millimolar
87  sigma = 1.0/7.0*(exp(Na_o/67.3)-1);
88  K_i = 1.5;
89  K_m_Na_i = 10.0;
90  I_Na_K_max = 0.59933874;
91  I_NaCa_max = 1600.0;
92  gamma = 0.35;
93  Ca_o = 1.8;
94  K_m_Na = 87.5;
95  K_m_Ca = 1.38;
96  K_sat = 0.1;
97  I_p_Ca_max = 0.275;
98  Trpn_max = 0.07;
99  Km_Trpn = 0.0005;
100  Cmdn_max = 0.05;
101  Csqn_max = 10.0;
102  Km_Cmdn = 0.00238;
103  Km_Csqn = 0.8;
104  NSR_I_up_max = 0.005;
105  NSR_I_Ca_max = 15.0;
106  NSR_K_up = 0.00092;
107  JSR_K_rel = 30.0;
108  JSR_V_cell = 20100.0;
109  JSR_V_rel = 0.0048 * JSR_V_cell;
110  JSR_V_up = 0.0552 * JSR_V_cell;
111  tau_tr = 180.0;
112  K_Q10 = 3.0;
113  V_i = 0.68*JSR_V_cell;
114 
115  switch (model_variant) {
116  case eOriginal:
117  g_to = 0.1652; // nanoS_per_picoF
118  g_Kur_scaling = 1.0;
119  g_Ca_L = 0.12375;
120  break;
121  case eAF:
122  g_to = 0.0826; // nanoS_per_picoF
123  g_Kur_scaling = 0.5;
124  g_Ca_L = 0.037125;
125  break;
126  }
127 
128  m_nvar = 21;
129 
130  // List gates and concentrations
131  m_gates.push_back(1);
132  m_gates.push_back(2);
133  m_gates.push_back(3);
134  m_gates.push_back(4);
135  m_gates.push_back(5);
136  m_gates.push_back(6);
137  m_gates.push_back(7);
138  m_gates.push_back(8);
139  m_gates.push_back(9);
140  m_gates.push_back(10);
141  m_gates.push_back(11);
142  m_gates.push_back(12);
143  m_gates.push_back(13);
144  m_gates.push_back(14);
145  m_gates.push_back(15);
146  m_concentrations.push_back(16);
147  m_concentrations.push_back(17);
148  m_concentrations.push_back(18);
149  m_concentrations.push_back(19);
150  m_concentrations.push_back(20);
151  }
152 
153 
154 
155  /**
156  *
157  */
158  CourtemancheRamirezNattel98::~CourtemancheRamirezNattel98()
159  {
160 
161  }
162 
163 
164 
165  void CourtemancheRamirezNattel98::v_Update(
166  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
167  Array<OneD, Array<OneD, NekDouble> >&outarray,
168  const NekDouble time)
169  {
170  ASSERTL0(inarray.get() != outarray.get(),
171  "Must have different arrays for input and output.");
172 
173  // Variables
174  // 0 V membrane potential
175  // 2 m fast sodium current m gate
176  // 3 h fast sodium current h gate
177  // 4 j fast sodium current j gate
178  // 5 o_a transient outward potassium o_a gate
179  // 6 o_i transient outward potassium o_i gate
180  // 7 u_a ultra-rapid delayed rectifier K current gate
181  // 8 u_i ultra-rapid delayed rectifier K current gate
182  // 9 x_r rapid delayed rectifier K current gate
183  // 10 x_s slow delayed rectifier K current gate
184  // 11 d L_type calcium gate
185  // 12 f L-type calcium gate
186  // 13 f_Ca L-type calcium gate
187  // 14 u Ca release u gate
188  // 15 v Ca release v gate
189  // 16 w Ca release w gate
190  // 17 Na_i Sodium
191  // 18 Ca_i Calcium
192  // 19 K_i Potassium
193  // 20 Ca_rel Calcium Rel
194  // 21 Ca_up Calcium up
195  int n = m_nq;
196  int i = 0;
197  NekDouble alpha, beta;
198  Vmath::Zero(n, outarray[0], 1);
199 
200  Array<OneD, NekDouble> &tmp = outarray[11];
201  Array<OneD, NekDouble> &tmp2 = outarray[12];
202 
203  // E_Na
204  Array<OneD, NekDouble> &tmp_E_na = outarray[14];
205  Vmath::Sdiv(n, Na_o, inarray[16], 1, tmp_E_na, 1);
206  Vmath::Vlog(n, tmp_E_na, 1, tmp_E_na, 1);
207  Vmath::Smul(n, R*T/F, tmp_E_na, 1, tmp_E_na, 1);
208 
209  // Sodium I_Na
210  Array<OneD, NekDouble> &tmp_I_Na = outarray[15];
211  Vmath::Vsub(n, inarray[0], 1, tmp_E_na, 1, tmp_I_Na, 1);
212  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
213  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
214  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
215  Vmath::Vmul(n, inarray[2], 1, tmp_I_Na, 1, tmp_I_Na, 1);
216  Vmath::Vmul(n, inarray[3], 1, tmp_I_Na, 1, tmp_I_Na, 1);
217  Vmath::Smul(n, C_m*g_Na, tmp_I_Na, 1, tmp_I_Na, 1);
218  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na, 1, outarray[0], 1);
219  Vmath::Smul(n, -1.0, tmp_I_Na, 1, outarray[16], 1);
220 
221  // Background current, sodium
222  Array<OneD, NekDouble> &tmp_I_b_Na = outarray[15];
223  Vmath::Vsub(n, inarray[0], 1, tmp_E_na, 1, tmp_I_b_Na, 1);
224  Vmath::Smul(n, C_m*g_b_Na, tmp_I_b_Na, 1, tmp_I_b_Na, 1);
225  Vmath::Vsub(n, outarray[0], 1, tmp_I_b_Na, 1, outarray[0], 1);
226  Vmath::Vsub(n, outarray[16], 1, tmp_I_b_Na, 1, outarray[16], 1);
227 
228  // V - E_K
229  Array<OneD, NekDouble> &tmp_V_E_k = outarray[14];
230  Vmath::Sdiv(n, K_o, inarray[18], 1, tmp_V_E_k, 1);
231  Vmath::Vlog(n, tmp_V_E_k, 1, tmp_V_E_k, 1);
232  Vmath::Smul(n, R*T/F, tmp_V_E_k, 1, tmp_V_E_k, 1);
233  Vmath::Vsub(n, inarray[0], 1, tmp_V_E_k, 1, tmp_V_E_k, 1);
234 
235  // Potassium I_K1
236  Array<OneD, NekDouble> &tmp_I_K1 = outarray[15];
237  Vmath::Sadd(n, 80.0, inarray[0], 1, tmp_I_K1, 1);
238  Vmath::Smul(n, 0.07, tmp_I_K1, 1, tmp_I_K1, 1);
239  Vmath::Vexp(n, tmp_I_K1, 1, tmp_I_K1, 1);
240  Vmath::Sadd(n, 1.0, tmp_I_K1, 1, tmp_I_K1, 1);
241  Vmath::Vdiv(n, tmp_V_E_k, 1, tmp_I_K1, 1, tmp_I_K1, 1);
242  Vmath::Smul(n, C_m*g_K1, tmp_I_K1, 1, tmp_I_K1, 1);
243  Vmath::Vsub(n, outarray[0], 1, tmp_I_K1, 1, outarray[0], 1);
244  Vmath::Smul(n, -1.0, tmp_I_K1, 1, outarray[18], 1);
245 
246  // Transient Outward K+ current
247  Array<OneD, NekDouble> &tmp_I_to = outarray[15];
248  Vmath::Vmul(n, inarray[5], 1, tmp_V_E_k, 1, tmp_I_to, 1);
249  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
250  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
251  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
252  Vmath::Smul(n, C_m*g_to, tmp_I_to, 1, tmp_I_to, 1);
253  Vmath::Vsub(n, outarray[0], 1, tmp_I_to, 1, outarray[0], 1);
254  Vmath::Vsub(n, outarray[18], 1, tmp_I_to, 1, outarray[18], 1);
255 
256  // Ultrarapid Delayed rectifier K+ current
257  Array<OneD, NekDouble> &tmp_I_kur = outarray[15];
258  Vmath::Sadd(n, -15.0, inarray[0], 1, tmp_I_kur, 1);
259  Vmath::Smul(n, -1.0/13.0, tmp_I_kur, 1, tmp_I_kur, 1);
260  Vmath::Vexp(n, tmp_I_kur, 1, tmp_I_kur, 1);
261  Vmath::Sadd(n, 1.0, tmp_I_kur, 1, tmp_I_kur, 1);
262  Vmath::Sdiv(n, 0.05, tmp_I_kur, 1, tmp_I_kur, 1);
263  Vmath::Sadd(n, 0.005, tmp_I_kur, 1, tmp_I_kur, 1);
264  Vmath::Vmul(n, tmp_V_E_k, 1, tmp_I_kur, 1, tmp_I_kur, 1);
265  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
266  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
267  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
268  Vmath::Vmul(n, inarray[7], 1, tmp_I_kur, 1, tmp_I_kur, 1);
269  Vmath::Smul(n, C_m*g_Kur_scaling, tmp_I_kur, 1, tmp_I_kur, 1);
270  Vmath::Vsub(n, outarray[0], 1, tmp_I_kur, 1, outarray[0], 1);
271  Vmath::Vsub(n, outarray[18], 1, tmp_I_kur, 1, outarray[18], 1);
272 
273  // Rapid delayed outward rectifier K+ current
274  Array<OneD, NekDouble> &tmp_I_Kr = outarray[15];
275  Vmath::Sadd(n, 15.0, inarray[0], 1, tmp_I_Kr, 1);
276  Vmath::Smul(n, 1.0/22.4, tmp_I_Kr, 1, tmp_I_Kr, 1);
277  Vmath::Vexp(n, tmp_I_Kr, 1, tmp_I_Kr, 1);
278  Vmath::Sadd(n, 1.0, tmp_I_Kr, 1, tmp_I_Kr, 1);
279  Vmath::Vdiv(n, tmp_V_E_k, 1, tmp_I_Kr, 1, tmp_I_Kr, 1);
280  Vmath::Vmul(n, inarray[8], 1, tmp_I_Kr, 1, tmp_I_Kr, 1);
281  Vmath::Smul(n, C_m*g_Kr, tmp_I_Kr, 1, tmp_I_Kr, 1);
282  Vmath::Vsub(n, outarray[0], 1, tmp_I_Kr, 1, outarray[0], 1);
283  Vmath::Vsub(n, outarray[18], 1, tmp_I_Kr, 1, outarray[18], 1);
284 
285  // Slow delayed outward rectifier K+ Current
286  Array<OneD, NekDouble> &tmp_I_Ks = outarray[15];
287  Vmath::Vmul(n, inarray[9], 1, tmp_V_E_k, 1, tmp_I_Ks, 1);
288  Vmath::Vmul(n, inarray[9], 1, tmp_I_Ks, 1, tmp_I_Ks, 1);
289  Vmath::Smul(n, C_m*g_Ks, tmp_I_Ks, 1, tmp_I_Ks, 1);
290  Vmath::Vsub(n, outarray[0], 1, tmp_I_Ks, 1, outarray[0], 1);
291  Vmath::Vsub(n, outarray[18], 1, tmp_I_Ks, 1, outarray[18], 1);
292 
293  // Background current, calcium
294  Array<OneD, NekDouble> &tmp_I_b_Ca = outarray[1];
295  Vmath::Sdiv(n, Ca_o, inarray[17], 1, tmp_I_b_Ca, 1);
296  Vmath::Vlog(n, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
297  Vmath::Smul(n, 0.5*R*T/F, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
298  Vmath::Vsub(n, inarray[0], 1, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
299  Vmath::Smul(n, C_m*g_b_Ca, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
300  Vmath::Vsub(n, outarray[0], 1, tmp_I_b_Ca, 1, outarray[0], 1);
301 
302  // L-Type Ca2+ current
303  Array<OneD, NekDouble> &tmp_I_Ca_L = outarray[2];
304  Vmath::Sadd(n, -65.0, inarray[0], 1, tmp_I_Ca_L, 1);
305  Vmath::Vmul(n, inarray[10], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
306  Vmath::Vmul(n, inarray[11], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
307  Vmath::Vmul(n, inarray[12], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
308  Vmath::Smul(n, C_m*g_Ca_L, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
309  Vmath::Vsub(n, outarray[0], 1, tmp_I_Ca_L, 1, outarray[0], 1);
310 
311  // Na-K Pump Current
312  Array<OneD, NekDouble> &tmp_f_Na_k = outarray[14];
313  Vmath::Smul(n, -F/R/T, inarray[0], 1, tmp_f_Na_k, 1);
314  Vmath::Vexp(n, tmp_f_Na_k, 1, tmp, 1);
315  Vmath::Smul(n, 0.0365*sigma, tmp, 1, tmp, 1);
316  Vmath::Smul(n, -0.1*F/R/T, inarray[0], 1, tmp_f_Na_k, 1);
317  Vmath::Vexp(n, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
318  Vmath::Smul(n, 0.1245, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
319  Vmath::Vadd(n, tmp_f_Na_k, 1, tmp, 1, tmp_f_Na_k, 1);
320  Vmath::Sadd(n, 1.0, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
321 
322  Array<OneD, NekDouble> &tmp_I_Na_K = outarray[15];
323  Vmath::Sdiv(n, K_m_Na_i, inarray[16], 1, tmp_I_Na_K, 1);
324  Vmath::Vpow(n, tmp_I_Na_K, 1, 1.5, tmp_I_Na_K, 1);
325  Vmath::Sadd(n, 1.0, tmp_I_Na_K, 1, tmp_I_Na_K, 1);
326  Vmath::Vmul(n, tmp_f_Na_k, 1, tmp_I_Na_K, 1, tmp_I_Na_K, 1);
327  Vmath::Sdiv(n, C_m*I_Na_K_max*K_o/(K_o+K_i), tmp_I_Na_K, 1, tmp_I_Na_K, 1);
328  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na_K, 1, outarray[0], 1);
329  Vmath::Svtvp(n, -3.0, tmp_I_Na_K, 1, outarray[16], 1, outarray[16], 1);
330  Vmath::Svtvp(n, 2.0, tmp_I_Na_K, 1, outarray[18], 1, outarray[18], 1);
331 
332  // Na-Ca exchanger current
333  Array<OneD, NekDouble> &tmp_I_Na_Ca = outarray[3];
334  Vmath::Smul(n, (gamma-1)*F/R/T, inarray[0], 1, tmp, 1);
335  Vmath::Vexp(n, tmp, 1, tmp, 1);
336  Vmath::Smul(n, K_sat, tmp, 1, tmp_I_Na_Ca, 1);
337  Vmath::Sadd(n, 1.0, tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
338  Vmath::Smul(n, (K_m_Na*K_m_Na*K_m_Na + Na_o*Na_o*Na_o)*(K_m_Ca + Ca_o), tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
339 
340  Vmath::Smul(n, Na_o*Na_o*Na_o, tmp, 1, tmp2, 1);
341  Vmath::Vmul(n, tmp2, 1, inarray[17], 1, tmp2, 1);
342  Vmath::Smul(n, gamma*F/R/T, inarray[0], 1, tmp, 1);
343  Vmath::Vexp(n, tmp, 1, tmp, 1);
344  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
345  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
346  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
347  Vmath::Svtvm(n, Ca_o, tmp, 1, tmp2, 1, tmp, 1);
348  Vmath::Smul(n, C_m*I_NaCa_max, tmp, 1, tmp, 1);
349  Vmath::Vdiv(n, tmp, 1, tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
350  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na_Ca, 1, outarray[0], 1);
351  Vmath::Svtvp(n, -3.0, tmp_I_Na_Ca, 1, outarray[16], 1, outarray[16], 1);
352 
353  // Calcium Pump current
354  Array<OneD, NekDouble> &tmp_I_p_Ca = outarray[4];
355  Vmath::Sadd(n, 0.0005, inarray[17], 1, tmp_I_p_Ca, 1);
356  Vmath::Vdiv(n, inarray[17], 1, tmp_I_p_Ca, 1, tmp_I_p_Ca, 1);
357  Vmath::Smul(n, C_m*I_p_Ca_max, tmp_I_p_Ca, 1, tmp_I_p_Ca, 1);
358  Vmath::Vsub(n, outarray[0], 1, tmp_I_p_Ca, 1, outarray[0], 1);
359 
360  // Scale currents by capacitance
361  Vmath::Smul(n, 1.0/C_m, outarray[0], 1, outarray[0], 1);
362 
363  // Scale sodium and potassium by FV_i
364  Vmath::Smul(n, 1.0/F/V_i, outarray[16], 1, outarray[16], 1);
365  Vmath::Smul(n, 1.0/F/V_i, outarray[18], 1, outarray[18], 1);
366 
367  // I_tr
368  Array<OneD, NekDouble> &tmp_I_tr = outarray[5];
369  Vmath::Vsub(n, inarray[20], 1, inarray[19], 1, tmp_I_tr, 1);
370  Vmath::Smul(n, 1.0/tau_tr, tmp_I_tr, 1, tmp_I_tr, 1);
371 
372  // I_up_leak
373  Array<OneD, NekDouble> &tmp_I_up_leak = outarray[6];
374  Vmath::Smul(n, NSR_I_up_max/NSR_I_Ca_max, inarray[20], 1, tmp_I_up_leak, 1);
375 
376  // I_up
377  Array<OneD, NekDouble> &tmp_I_up = outarray[7];
378  Vmath::Sdiv(n, NSR_K_up, inarray[17], 1, tmp_I_up, 1);
379  Vmath::Sadd(n, 1.0, tmp_I_up, 1, tmp_I_up, 1);
380  Vmath::Sdiv(n, NSR_I_up_max, tmp_I_up, 1, tmp_I_up, 1);
381 
382  // I_rel
383  Array<OneD, NekDouble> &tmp_I_rel = outarray[8];
384  Vmath::Vsub(n, inarray[19], 1, inarray[17], 1, tmp_I_rel, 1);
385  Vmath::Vmul(n, tmp_I_rel, 1, inarray[13], 1, tmp_I_rel, 1);
386  Vmath::Vmul(n, tmp_I_rel, 1, inarray[13], 1, tmp_I_rel, 1);
387  Vmath::Vmul(n, tmp_I_rel, 1, inarray[14], 1, tmp_I_rel, 1);
388  Vmath::Vmul(n, tmp_I_rel, 1, inarray[15], 1, tmp_I_rel, 1);
389  Vmath::Smul(n, JSR_K_rel, tmp_I_rel, 1, tmp_I_rel, 1);
390 
391  // B1
392  Array<OneD, NekDouble> &tmp_B1 = outarray[9];
393  Vmath::Svtvm(n, 2.0, tmp_I_Na_Ca, 1, tmp_I_p_Ca, 1, tmp_B1, 1);
394  Vmath::Vsub(n, tmp_B1, 1, tmp_I_Ca_L, 1, tmp_B1, 1);
395  Vmath::Vsub(n, tmp_B1, 1, tmp_I_b_Ca, 1, tmp_B1, 1);
396  Vmath::Smul(n, 0.5/F, tmp_B1, 1, tmp_B1, 1);
397  Vmath::Svtvp(n, JSR_V_up, tmp_I_up_leak, 1, tmp_B1, 1, tmp_B1, 1);
398  Vmath::Svtvp(n, -JSR_V_up, tmp_I_up, 1, tmp_B1, 1, tmp_B1, 1);
399  Vmath::Svtvp(n, JSR_V_rel, tmp_I_rel, 1, tmp_B1, 1, tmp_B1, 1);
400  Vmath::Smul(n, 1.0/V_i, tmp_B1, 1, tmp_B1, 1);
401 
402  // B2
403  Array<OneD, NekDouble> &tmp_B2 = outarray[10];
404  Vmath::Sadd(n, Km_Cmdn, inarray[17], 1, tmp_B2, 1);
405  Vmath::Vmul(n, tmp_B2, 1, tmp_B2, 1, tmp_B2, 1);
406  Vmath::Sdiv(n, Cmdn_max*Km_Cmdn, tmp_B2, 1, tmp_B2, 1);
407  Vmath::Sadd(n, Km_Trpn, inarray[17], 1, tmp, 1);
408  Vmath::Vmul(n, tmp, 1, tmp, 1, tmp, 1);
409  Vmath::Sdiv(n, Trpn_max*Km_Trpn, tmp, 1, tmp, 1);
410  Vmath::Vadd(n, tmp, 1, tmp_B2, 1, tmp_B2, 1);
411  Vmath::Sadd(n, 1.0, tmp_B2, 1, tmp_B2, 1);
412 
413  // Calcium concentration (18)
414  Vmath::Vdiv(n, tmp_B1, 1, tmp_B2, 1, outarray[17], 1);
415 
416  // Calcium up (21)
417  Vmath::Vsub(n, tmp_I_up, 1, tmp_I_up_leak, 1, outarray[20], 1);
418  Vmath::Svtvp(n, -JSR_V_rel/JSR_V_up, tmp_I_tr, 1, outarray[20], 1, outarray[20], 1);
419 
420  // Calcium rel (20)
421  Vmath::Vsub(n, tmp_I_tr, 1, tmp_I_rel, 1, tmp, 1);
422  Vmath::Sadd(n, Km_Csqn, inarray[19], 1, outarray[19], 1);
423  Vmath::Vmul(n, outarray[19], 1, outarray[19], 1, outarray[19], 1);
424  Vmath::Sdiv(n, Csqn_max*Km_Csqn, outarray[19], 1, outarray[19], 1);
425  Vmath::Sadd(n, 1.0, outarray[19], 1, outarray[19], 1);
426  Vmath::Vdiv(n, tmp, 1, outarray[19], 1, outarray[19], 1);
427 
428  // Process gating variables
429  const NekDouble * v;
430  const NekDouble * x;
431  NekDouble * x_tau;
432  NekDouble * x_new;
433  // m
434  for (i = 0, v = &inarray[0][0], x = &inarray[1][0], x_new = &outarray[1][0], x_tau = &m_gates_tau[0][0];
435  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
436  {
437  alpha = (*v == (-47.13)) ? 3.2 : (0.32*(*v+47.13))/(1.0-exp((-0.1)*(*v + 47.13)));
438  beta = 0.08*exp(-(*v)/11.0);
439  *x_tau = 1.0/(alpha + beta);
440  *x_new = alpha*(*x_tau);
441  }
442  // h
443  for (i = 0, v = &inarray[0][0], x = &inarray[2][0], x_new = &outarray[2][0], x_tau = &m_gates_tau[1][0];
444  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
445  {
446  alpha = (*v >= -40.0) ? 0.0 : 0.135*exp(-((*v)+80.0)/6.8);
447  beta = (*v >= -40.0) ? 1.0/(0.13*(1.0+exp(-(*v + 10.66)/11.1)))
448  : 3.56*exp(0.079*(*v))+310000.0*exp(0.35*(*v));
449  *x_tau = 1.0/(alpha + beta);
450  *x_new = alpha*(*x_tau);
451  }
452  // j
453  for (i = 0, v = &inarray[0][0], x = &inarray[3][0], x_new = &outarray[3][0], x_tau = &m_gates_tau[2][0];
454  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
455  {
456  alpha = (*v >= -40.0) ? 0.0
457  : (-127140.0*exp(0.2444*(*v))-3.474e-05*exp(-0.04391*(*v)))*(((*v)+37.78)/(1.0+exp(0.311*((*v)+79.23))));
458  beta = (*v >= -40.0) ? (0.3*exp(-2.535e-07*(*v))/(1.0+exp(-0.1*(*v+32.0))))
459  : 0.1212*exp(-0.01052*(*v))/(1.0+exp(-0.1378*(*v+40.14)));
460  *x_tau = 1.0/(alpha + beta);
461  *x_new = alpha*(*x_tau);
462  }
463  // oa
464  for (i = 0, v = &inarray[0][0], x = &inarray[4][0], x_new = &outarray[4][0], x_tau = &m_gates_tau[3][0];
465  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
466  {
467  alpha = 0.65/(exp(-(*v+10.0)/8.5) + exp(-(*v-30.0)/59.0));
468  beta = 0.65/(2.5 + exp((*v+82.0)/17.0));
469  *x_tau = 1.0/K_Q10/(alpha + beta);
470  *x_new = (1.0/(1.0+exp(-(*v+20.47)/17.54)));
471  }
472  // oi
473  for (i = 0, v = &inarray[0][0], x = &inarray[5][0], x_new = &outarray[5][0], x_tau = &m_gates_tau[4][0];
474  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
475  {
476  alpha = 1.0/(18.53 + exp((*v+113.7)/10.95));
477  beta = 1.0/(35.56 + exp(-(*v+1.26)/7.44));
478  *x_tau = 1.0/K_Q10/(alpha + beta);
479  *x_new = (1.0/(1.0+exp((*v+43.1)/5.3)));
480  }
481  // ua
482  for (i = 0, v = &inarray[0][0], x = &inarray[6][0], x_new = &outarray[6][0], x_tau = &m_gates_tau[5][0];
483  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
484  {
485  alpha = 0.65/(exp(-(*v+10.0)/8.5)+exp(-(*v-30.0)/59.0));
486  beta = 0.65/(2.5+exp((*v+82.0)/17.0));
487  *x_tau = 1.0/K_Q10/(alpha + beta);
488  *x_new = 1.0/(1.0+exp(-(*v+30.3)/9.6));
489  }
490  // ui
491  for (i = 0, v = &inarray[0][0], x = &inarray[7][0], x_new = &outarray[7][0], x_tau = &m_gates_tau[6][0];
492  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
493  {
494  alpha = 1.0/(21.0 + exp(-(*v-185.0)/28.0));
495  beta = exp((*v-158.0)/16.0);
496  *x_tau = 1.0/K_Q10/(alpha + beta);
497  *x_new = 1.0/(1.0+exp((*v-99.45)/27.48));
498  }
499  // xr
500  for (i = 0, v = &inarray[0][0], x = &inarray[8][0], x_new = &outarray[8][0], x_tau = &m_gates_tau[7][0];
501  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
502  {
503  alpha = 0.0003*(*v+14.1)/(1-exp(-(*v+14.1)/5.0));
504  beta = 7.3898e-5*(*v-3.3328)/(exp((*v-3.3328)/5.1237)-1.0);
505  *x_tau = 1.0/(alpha + beta);
506  *x_new = 1.0/(1+exp(-(*v+14.1)/6.5));
507  }
508  // xs
509  for (i = 0, v = &inarray[0][0], x = &inarray[9][0], x_new = &outarray[9][0], x_tau = &m_gates_tau[8][0];
510  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
511  {
512  alpha = 4e-5*(*v-19.9)/(1.0-exp(-(*v-19.9)/17.0));
513  beta = 3.5e-5*(*v-19.9)/(exp((*v-19.9)/9.0)-1.0);
514  *x_tau = 0.5/(alpha + beta);
515  *x_new = 1.0/sqrt(1.0+exp(-(*v-19.9)/12.7));
516  }
517  // d
518  for (i = 0, v = &inarray[0][0], x = &inarray[10][0], x_new = &outarray[10][0], x_tau = &m_gates_tau[9][0];
519  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
520  {
521  *x_tau = (1-exp(-(*v+10.0)/6.24))/(0.035*(*v+10.0)*(1+exp(-(*v+10.0)/6.24)));
522  *x_new = 1.0/(1.0 + exp(-(*v+10)/8.0));
523  }
524  // f
525  for (i = 0, v = &inarray[0][0], x = &inarray[11][0], x_new = &outarray[11][0], x_tau = &m_gates_tau[10][0];
526  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
527  {
528  //alpha = 1.0/(1.0 + exp((*v+28.0)/6.9));
529  *x_tau = 9.0/(0.0197*exp(-0.0337*0.0337*(*v+10.0)*(*v+10.0))+0.02);
530  *x_new = exp((-(*v + 28.0)) / 6.9) / (1.0 + exp((-(*v + 28.0)) / 6.9));
531  }
532  // f_Ca
533  for (i = 0, v = &inarray[0][0], x = &inarray[12][0], x_new = &outarray[12][0], x_tau = &m_gates_tau[11][0];
534  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
535  {
536  *x_tau = 2.0;
537  *x_new = 1.0/(1.0+inarray[17][i]/0.00035);
538  }
539 
540  Array<OneD, NekDouble> &tmp_Fn = outarray[15];
541  Vmath::Svtsvtp(n, 0.5*5e-13/F, tmp_I_Ca_L, 1, -0.2*5e-13/F, tmp_I_Na_Ca, 1, tmp_Fn, 1);
542  Vmath::Svtvm(n, 1e-12*JSR_V_rel, tmp_I_rel, 1, tmp_Fn, 1, tmp_Fn, 1);
543 
544  // u
545  for (i = 0, v = &tmp_Fn[0], x = &inarray[13][0], x_new = &outarray[13][0], x_tau = &m_gates_tau[12][0];
546  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
547  {
548  *x_tau = 8.0;
549  *x_new = 1.0/(1.0 + exp(-(*v - 3.4175e-13)/1.367e-15));
550  }
551  // v
552  for (i = 0, v = &tmp_Fn[0], x = &inarray[14][0], x_new = &outarray[14][0], x_tau = &m_gates_tau[13][0];
553  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
554  {
555  *x_tau = 1.91 + 2.09/(1.0+exp(-(*v - 3.4175e-13)/13.67e-16));
556  *x_new = 1.0 - 1.0/(1.0 + exp(-(*v - 6.835e-14)/13.67e-16));
557  }
558  // w
559  for (i = 0, v = &inarray[0][0], x = &inarray[15][0], x_new = &outarray[15][0], x_tau = &m_gates_tau[14][0];
560  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
561  {
562  *x_tau = 6.0*(1.0-exp(-(*v-7.9)/5.0))/(1.0+0.3*exp(-(*v-7.9)/5.0))/(*v-7.9);
563  *x_new = 1.0 - 1.0/(1.0 + exp(-(*v - 40.0)/17.0));
564  }
565 
566  }
567 
568 
569  /**
570  *
571  */
572  void CourtemancheRamirezNattel98::v_GenerateSummary(SummaryList& s)
573  {
574  SolverUtils::AddSummaryItem(s, "Cell model","CourtemancheRamirezNattel98");
575  SolverUtils::AddSummaryItem(s, "Cell model var.", lookupIds[model_variant]);
576  }
577 
578 
579  void CourtemancheRamirezNattel98::v_SetInitialConditions()
580  {
581  Vmath::Fill(m_nq, -81.0, m_cellSol[0], 1);
582  Vmath::Fill(m_nq, 2.908e-03, m_cellSol[1], 1);
583  Vmath::Fill(m_nq, 9.649e-01, m_cellSol[2], 1);
584  Vmath::Fill(m_nq, 9.775e-01, m_cellSol[3], 1);
585  Vmath::Fill(m_nq, 3.043e-02, m_cellSol[4], 1);
586  Vmath::Fill(m_nq, 9.992e-01, m_cellSol[5], 1);
587  Vmath::Fill(m_nq, 4.966e-03, m_cellSol[6], 1);
588  Vmath::Fill(m_nq, 9.986e-01, m_cellSol[7], 1);
589  Vmath::Fill(m_nq, 3.296e-05, m_cellSol[8], 1);
590  Vmath::Fill(m_nq, 1.869e-02, m_cellSol[9], 1);
591  Vmath::Fill(m_nq, 1.367e-04, m_cellSol[10], 1);
592  Vmath::Fill(m_nq, 9.996e-01, m_cellSol[11], 1);
593  Vmath::Fill(m_nq, 7.755e-01, m_cellSol[12], 1);
594  Vmath::Fill(m_nq, 2.35e-112, m_cellSol[13], 1);
595  Vmath::Fill(m_nq, 1.0, m_cellSol[14], 1);
596  Vmath::Fill(m_nq, 0.9992, m_cellSol[15], 1);
597  Vmath::Fill(m_nq, 1.117e+01, m_cellSol[16], 1);
598  Vmath::Fill(m_nq, 1.013e-04, m_cellSol[17], 1);
599  Vmath::Fill(m_nq, 1.39e+02, m_cellSol[18], 1);
600  Vmath::Fill(m_nq, 1.488, m_cellSol[19], 1);
601  Vmath::Fill(m_nq, 1.488, m_cellSol[20], 1);
602  }
603 
604  std::string CourtemancheRamirezNattel98::v_GetCellVarName(unsigned int idx)
605  {
606  switch (idx)
607  {
608  case 0: return "u";
609  case 1: return "m";
610  case 2: return "h";
611  case 3: return "j";
612  case 4: return "o_a";
613  case 5: return "o_i";
614  case 6: return "u_a";
615  case 7: return "u_i";
616  case 8: return "x_r";
617  case 9: return "x_s";
618  case 10: return "d";
619  case 11: return "f";
620  case 12: return "f_Ca";
621  case 13: return "U";
622  case 14: return "V";
623  case 15: return "W";
624  case 16: return "Na_i";
625  case 17: return "Ca_i";
626  case 18: return "K_i";
627  case 19: return "Ca_rel";
628  case 20: return "Ca_up";
629  default: return "unknown";
630  }
631  }
632 
633 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::CellModel
Cell model base class.
Definition: CellModel.h:65
Nektar::CellModel::m_cellSol
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:125
Nektar::CellModel::m_nq
int m_nq
Number of physical points.
Definition: CellModel.h:116
Nektar::CellModel::m_concentrations
std::vector< int > m_concentrations
Indices of cell model variables which are concentrations.
Definition: CellModel.h:138
Nektar::CellModel::m_gates
std::vector< int > m_gates
Indices of cell model variables which are gates.
Definition: CellModel.h:140
Nektar::CellModel::m_gates_tau
Array< OneD, Array< OneD, NekDouble > > m_gates_tau
Storage for gate tau values.
Definition: CellModel.h:142
Nektar::CellModel::m_nvar
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:118
Nektar::CourtemancheRamirezNattel98::v_Update
virtual void v_Update(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Computes the reaction terms $f(u,v)$ and $g(u,v)$.
Definition: CourtemancheRamirezNattel98.cpp:165
Nektar::CourtemancheRamirezNattel98::v_GetCellVarName
virtual std::string v_GetCellVarName(unsigned int idx)
Definition: CourtemancheRamirezNattel98.cpp:604
Nektar::CourtemancheRamirezNattel98::v_GenerateSummary
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
Definition: CourtemancheRamirezNattel98.cpp:572
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::GetCellModelFactory
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:46
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Svtsvtp
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:691
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:192
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:104
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:116
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:565
Vmath::Svtvm
void Svtvm(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x - y
Definition: Vmath.cpp:602
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Sdiv
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:291
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:341
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372
Vmath::Vpow
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
Definition: Vmath.hpp:127
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:26 for Nektar++ by doxygen 1.9.1