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 <CardiacEPSolver/CellModels/CourtemancheRamirezNattel98.h>
39 #include <LibUtilities/BasicUtils/VmathArray.hpp>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 std::string CourtemancheRamirezNattel98::className =
46  GetCellModelFactory().RegisterCreatorFunction(
47  "CourtemancheRamirezNattel98", CourtemancheRamirezNattel98::create,
48  "Ionic model of human atrial cell electrophysiology.");
49 
50 // Register cell model variants
51 std::string CourtemancheRamirezNattel98::lookupIds[2] = {
52  LibUtilities::SessionReader::RegisterEnumValue(
53  "CellModelVariant", "Original", CourtemancheRamirezNattel98::eOriginal),
54  LibUtilities::SessionReader::RegisterEnumValue(
55  "CellModelVariant", "AF", CourtemancheRamirezNattel98::eAF)};
56 
57 // Register default variant
58 std::string CourtemancheRamirezNattel98::def =
59  LibUtilities::SessionReader::RegisterDefaultSolverInfo("CellModelVariant",
60  "Original");
61 
62 /**
63  *
64  */
65 CourtemancheRamirezNattel98::CourtemancheRamirezNattel98(
66  const LibUtilities::SessionReaderSharedPtr &pSession,
67  const MultiRegions::ExpListSharedPtr &pField)
68  : CellModel(pSession, pField)
69 {
70  model_variant =
71  pSession->GetSolverInfoAsEnum<CourtemancheRamirezNattel98::Variants>(
72  "CellModelVariant");
73 
74  C_m = 100; // picoF
75  g_Na = 7.8; // nanoS_per_picoF
76  g_K1 = 0.09; // nanoS_per_picoF
77  g_Kr = 0.029411765;
78  g_Ks = 0.12941176;
79  g_b_Na = 0.0006744375;
80  g_b_Ca = 0.001131;
81  R = 8.3143;
82  T = 310.0;
83  F = 96.4867;
84  Na_o = 140.0; // millimolar
85  K_o = 5.4; // millimolar
86  sigma = 1.0 / 7.0 * (exp(Na_o / 67.3) - 1);
87  K_i = 1.5;
88  K_m_Na_i = 10.0;
89  I_Na_K_max = 0.59933874;
90  I_NaCa_max = 1600.0;
91  gamma = 0.35;
92  Ca_o = 1.8;
93  K_m_Na = 87.5;
94  K_m_Ca = 1.38;
95  K_sat = 0.1;
96  I_p_Ca_max = 0.275;
97  Trpn_max = 0.07;
98  Km_Trpn = 0.0005;
99  Cmdn_max = 0.05;
100  Csqn_max = 10.0;
101  Km_Cmdn = 0.00238;
102  Km_Csqn = 0.8;
103  NSR_I_up_max = 0.005;
104  NSR_I_Ca_max = 15.0;
105  NSR_K_up = 0.00092;
106  JSR_K_rel = 30.0;
107  JSR_V_cell = 20100.0;
108  JSR_V_rel = 0.0048 * JSR_V_cell;
109  JSR_V_up = 0.0552 * JSR_V_cell;
110  tau_tr = 180.0;
111  K_Q10 = 3.0;
112  V_i = 0.68 * JSR_V_cell;
113 
114  switch (model_variant)
115  {
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 CourtemancheRamirezNattel98::~CourtemancheRamirezNattel98()
157 {
158 }
159 
160 void CourtemancheRamirezNattel98::v_Update(
161  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
162  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
163 {
164  ASSERTL0(inarray.get() != outarray.get(),
165  "Must have different arrays for input and output.");
166 
167  // Variables
168  // 0 V membrane potential
169  // 2 m fast sodium current m gate
170  // 3 h fast sodium current h gate
171  // 4 j fast sodium current j gate
172  // 5 o_a transient outward potassium o_a gate
173  // 6 o_i transient outward potassium o_i gate
174  // 7 u_a ultra-rapid delayed rectifier K current gate
175  // 8 u_i ultra-rapid delayed rectifier K current gate
176  // 9 x_r rapid delayed rectifier K current gate
177  // 10 x_s slow delayed rectifier K current gate
178  // 11 d L_type calcium gate
179  // 12 f L-type calcium gate
180  // 13 f_Ca L-type calcium gate
181  // 14 u Ca release u gate
182  // 15 v Ca release v gate
183  // 16 w Ca release w gate
184  // 17 Na_i Sodium
185  // 18 Ca_i Calcium
186  // 19 K_i Potassium
187  // 20 Ca_rel Calcium Rel
188  // 21 Ca_up Calcium up
189  int n = m_nq;
190  int i = 0;
191  NekDouble alpha, beta;
192  Vmath::Zero(n, outarray[0], 1);
193 
194  Array<OneD, NekDouble> &tmp = outarray[11];
195  Array<OneD, NekDouble> &tmp2 = outarray[12];
196 
197  // E_Na
198  Array<OneD, NekDouble> &tmp_E_na = outarray[14];
199  Vmath::Sdiv(n, Na_o, inarray[16], 1, tmp_E_na, 1);
200  Vmath::Vlog(n, tmp_E_na, 1, tmp_E_na, 1);
201  Vmath::Smul(n, R * T / F, tmp_E_na, 1, tmp_E_na, 1);
202 
203  // Sodium I_Na
204  Array<OneD, NekDouble> &tmp_I_Na = outarray[15];
205  Vmath::Vsub(n, inarray[0], 1, tmp_E_na, 1, tmp_I_Na, 1);
206  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
207  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
208  Vmath::Vmul(n, inarray[1], 1, tmp_I_Na, 1, tmp_I_Na, 1);
209  Vmath::Vmul(n, inarray[2], 1, tmp_I_Na, 1, tmp_I_Na, 1);
210  Vmath::Vmul(n, inarray[3], 1, tmp_I_Na, 1, tmp_I_Na, 1);
211  Vmath::Smul(n, C_m * g_Na, tmp_I_Na, 1, tmp_I_Na, 1);
212  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na, 1, outarray[0], 1);
213  Vmath::Smul(n, -1.0, tmp_I_Na, 1, outarray[16], 1);
214 
215  // Background current, sodium
216  Array<OneD, NekDouble> &tmp_I_b_Na = outarray[15];
217  Vmath::Vsub(n, inarray[0], 1, tmp_E_na, 1, tmp_I_b_Na, 1);
218  Vmath::Smul(n, C_m * g_b_Na, tmp_I_b_Na, 1, tmp_I_b_Na, 1);
219  Vmath::Vsub(n, outarray[0], 1, tmp_I_b_Na, 1, outarray[0], 1);
220  Vmath::Vsub(n, outarray[16], 1, tmp_I_b_Na, 1, outarray[16], 1);
221 
222  // V - E_K
223  Array<OneD, NekDouble> &tmp_V_E_k = outarray[14];
224  Vmath::Sdiv(n, K_o, inarray[18], 1, tmp_V_E_k, 1);
225  Vmath::Vlog(n, tmp_V_E_k, 1, tmp_V_E_k, 1);
226  Vmath::Smul(n, R * T / F, tmp_V_E_k, 1, tmp_V_E_k, 1);
227  Vmath::Vsub(n, inarray[0], 1, tmp_V_E_k, 1, tmp_V_E_k, 1);
228 
229  // Potassium I_K1
230  Array<OneD, NekDouble> &tmp_I_K1 = outarray[15];
231  Vmath::Sadd(n, 80.0, inarray[0], 1, tmp_I_K1, 1);
232  Vmath::Smul(n, 0.07, tmp_I_K1, 1, tmp_I_K1, 1);
233  Vmath::Vexp(n, tmp_I_K1, 1, tmp_I_K1, 1);
234  Vmath::Sadd(n, 1.0, tmp_I_K1, 1, tmp_I_K1, 1);
235  Vmath::Vdiv(n, tmp_V_E_k, 1, tmp_I_K1, 1, tmp_I_K1, 1);
236  Vmath::Smul(n, C_m * g_K1, tmp_I_K1, 1, tmp_I_K1, 1);
237  Vmath::Vsub(n, outarray[0], 1, tmp_I_K1, 1, outarray[0], 1);
238  Vmath::Smul(n, -1.0, tmp_I_K1, 1, outarray[18], 1);
239 
240  // Transient Outward K+ current
241  Array<OneD, NekDouble> &tmp_I_to = outarray[15];
242  Vmath::Vmul(n, inarray[5], 1, tmp_V_E_k, 1, tmp_I_to, 1);
243  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
244  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
245  Vmath::Vmul(n, inarray[4], 1, tmp_I_to, 1, tmp_I_to, 1);
246  Vmath::Smul(n, C_m * g_to, tmp_I_to, 1, tmp_I_to, 1);
247  Vmath::Vsub(n, outarray[0], 1, tmp_I_to, 1, outarray[0], 1);
248  Vmath::Vsub(n, outarray[18], 1, tmp_I_to, 1, outarray[18], 1);
249 
250  // Ultrarapid Delayed rectifier K+ current
251  Array<OneD, NekDouble> &tmp_I_kur = outarray[15];
252  Vmath::Sadd(n, -15.0, inarray[0], 1, tmp_I_kur, 1);
253  Vmath::Smul(n, -1.0 / 13.0, tmp_I_kur, 1, tmp_I_kur, 1);
254  Vmath::Vexp(n, tmp_I_kur, 1, tmp_I_kur, 1);
255  Vmath::Sadd(n, 1.0, tmp_I_kur, 1, tmp_I_kur, 1);
256  Vmath::Sdiv(n, 0.05, tmp_I_kur, 1, tmp_I_kur, 1);
257  Vmath::Sadd(n, 0.005, tmp_I_kur, 1, tmp_I_kur, 1);
258  Vmath::Vmul(n, tmp_V_E_k, 1, tmp_I_kur, 1, tmp_I_kur, 1);
259  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
260  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
261  Vmath::Vmul(n, inarray[6], 1, tmp_I_kur, 1, tmp_I_kur, 1);
262  Vmath::Vmul(n, inarray[7], 1, tmp_I_kur, 1, tmp_I_kur, 1);
263  Vmath::Smul(n, C_m * g_Kur_scaling, tmp_I_kur, 1, tmp_I_kur, 1);
264  Vmath::Vsub(n, outarray[0], 1, tmp_I_kur, 1, outarray[0], 1);
265  Vmath::Vsub(n, outarray[18], 1, tmp_I_kur, 1, outarray[18], 1);
266 
267  // Rapid delayed outward rectifier K+ current
268  Array<OneD, NekDouble> &tmp_I_Kr = outarray[15];
269  Vmath::Sadd(n, 15.0, inarray[0], 1, tmp_I_Kr, 1);
270  Vmath::Smul(n, 1.0 / 22.4, tmp_I_Kr, 1, tmp_I_Kr, 1);
271  Vmath::Vexp(n, tmp_I_Kr, 1, tmp_I_Kr, 1);
272  Vmath::Sadd(n, 1.0, tmp_I_Kr, 1, tmp_I_Kr, 1);
273  Vmath::Vdiv(n, tmp_V_E_k, 1, tmp_I_Kr, 1, tmp_I_Kr, 1);
274  Vmath::Vmul(n, inarray[8], 1, tmp_I_Kr, 1, tmp_I_Kr, 1);
275  Vmath::Smul(n, C_m * g_Kr, tmp_I_Kr, 1, tmp_I_Kr, 1);
276  Vmath::Vsub(n, outarray[0], 1, tmp_I_Kr, 1, outarray[0], 1);
277  Vmath::Vsub(n, outarray[18], 1, tmp_I_Kr, 1, outarray[18], 1);
278 
279  // Slow delayed outward rectifier K+ Current
280  Array<OneD, NekDouble> &tmp_I_Ks = outarray[15];
281  Vmath::Vmul(n, inarray[9], 1, tmp_V_E_k, 1, tmp_I_Ks, 1);
282  Vmath::Vmul(n, inarray[9], 1, tmp_I_Ks, 1, tmp_I_Ks, 1);
283  Vmath::Smul(n, C_m * g_Ks, tmp_I_Ks, 1, tmp_I_Ks, 1);
284  Vmath::Vsub(n, outarray[0], 1, tmp_I_Ks, 1, outarray[0], 1);
285  Vmath::Vsub(n, outarray[18], 1, tmp_I_Ks, 1, outarray[18], 1);
286 
287  // Background current, calcium
288  Array<OneD, NekDouble> &tmp_I_b_Ca = outarray[1];
289  Vmath::Sdiv(n, Ca_o, inarray[17], 1, tmp_I_b_Ca, 1);
290  Vmath::Vlog(n, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
291  Vmath::Smul(n, 0.5 * R * T / F, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
292  Vmath::Vsub(n, inarray[0], 1, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
293  Vmath::Smul(n, C_m * g_b_Ca, tmp_I_b_Ca, 1, tmp_I_b_Ca, 1);
294  Vmath::Vsub(n, outarray[0], 1, tmp_I_b_Ca, 1, outarray[0], 1);
295 
296  // L-Type Ca2+ current
297  Array<OneD, NekDouble> &tmp_I_Ca_L = outarray[2];
298  Vmath::Sadd(n, -65.0, inarray[0], 1, tmp_I_Ca_L, 1);
299  Vmath::Vmul(n, inarray[10], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
300  Vmath::Vmul(n, inarray[11], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
301  Vmath::Vmul(n, inarray[12], 1, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
302  Vmath::Smul(n, C_m * g_Ca_L, tmp_I_Ca_L, 1, tmp_I_Ca_L, 1);
303  Vmath::Vsub(n, outarray[0], 1, tmp_I_Ca_L, 1, outarray[0], 1);
304 
305  // Na-K Pump Current
306  Array<OneD, NekDouble> &tmp_f_Na_k = outarray[14];
307  Vmath::Smul(n, -F / R / T, inarray[0], 1, tmp_f_Na_k, 1);
308  Vmath::Vexp(n, tmp_f_Na_k, 1, tmp, 1);
309  Vmath::Smul(n, 0.0365 * sigma, tmp, 1, tmp, 1);
310  Vmath::Smul(n, -0.1 * F / R / T, inarray[0], 1, tmp_f_Na_k, 1);
311  Vmath::Vexp(n, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
312  Vmath::Smul(n, 0.1245, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
313  Vmath::Vadd(n, tmp_f_Na_k, 1, tmp, 1, tmp_f_Na_k, 1);
314  Vmath::Sadd(n, 1.0, tmp_f_Na_k, 1, tmp_f_Na_k, 1);
315 
316  Array<OneD, NekDouble> &tmp_I_Na_K = outarray[15];
317  Vmath::Sdiv(n, K_m_Na_i, inarray[16], 1, tmp_I_Na_K, 1);
318  Vmath::Vpow(n, tmp_I_Na_K, 1, 1.5, tmp_I_Na_K, 1);
319  Vmath::Sadd(n, 1.0, tmp_I_Na_K, 1, tmp_I_Na_K, 1);
320  Vmath::Vmul(n, tmp_f_Na_k, 1, tmp_I_Na_K, 1, tmp_I_Na_K, 1);
321  Vmath::Sdiv(n, C_m * I_Na_K_max * K_o / (K_o + K_i), tmp_I_Na_K, 1,
322  tmp_I_Na_K, 1);
323  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na_K, 1, outarray[0], 1);
324  Vmath::Svtvp(n, -3.0, tmp_I_Na_K, 1, outarray[16], 1, outarray[16], 1);
325  Vmath::Svtvp(n, 2.0, tmp_I_Na_K, 1, outarray[18], 1, outarray[18], 1);
326 
327  // Na-Ca exchanger current
328  Array<OneD, NekDouble> &tmp_I_Na_Ca = outarray[3];
329  Vmath::Smul(n, (gamma - 1) * F / R / T, inarray[0], 1, tmp, 1);
330  Vmath::Vexp(n, tmp, 1, tmp, 1);
331  Vmath::Smul(n, K_sat, tmp, 1, tmp_I_Na_Ca, 1);
332  Vmath::Sadd(n, 1.0, tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
333  Vmath::Smul(
334  n, (K_m_Na * K_m_Na * K_m_Na + Na_o * Na_o * Na_o) * (K_m_Ca + Ca_o),
335  tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
336 
337  Vmath::Smul(n, Na_o * Na_o * Na_o, tmp, 1, tmp2, 1);
338  Vmath::Vmul(n, tmp2, 1, inarray[17], 1, tmp2, 1);
339  Vmath::Smul(n, gamma * F / R / T, inarray[0], 1, tmp, 1);
340  Vmath::Vexp(n, tmp, 1, tmp, 1);
341  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
342  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
343  Vmath::Vmul(n, inarray[16], 1, tmp, 1, tmp, 1);
344  Vmath::Svtvm(n, Ca_o, tmp, 1, tmp2, 1, tmp, 1);
345  Vmath::Smul(n, C_m * I_NaCa_max, tmp, 1, tmp, 1);
346  Vmath::Vdiv(n, tmp, 1, tmp_I_Na_Ca, 1, tmp_I_Na_Ca, 1);
347  Vmath::Vsub(n, outarray[0], 1, tmp_I_Na_Ca, 1, outarray[0], 1);
348  Vmath::Svtvp(n, -3.0, tmp_I_Na_Ca, 1, outarray[16], 1, outarray[16], 1);
349 
350  // Calcium Pump current
351  Array<OneD, NekDouble> &tmp_I_p_Ca = outarray[4];
352  Vmath::Sadd(n, 0.0005, inarray[17], 1, tmp_I_p_Ca, 1);
353  Vmath::Vdiv(n, inarray[17], 1, tmp_I_p_Ca, 1, tmp_I_p_Ca, 1);
354  Vmath::Smul(n, C_m * I_p_Ca_max, tmp_I_p_Ca, 1, tmp_I_p_Ca, 1);
355  Vmath::Vsub(n, outarray[0], 1, tmp_I_p_Ca, 1, outarray[0], 1);
356 
357  // Scale currents by capacitance
358  Vmath::Smul(n, 1.0 / C_m, outarray[0], 1, outarray[0], 1);
359 
360  // Scale sodium and potassium by FV_i
361  Vmath::Smul(n, 1.0 / F / V_i, outarray[16], 1, outarray[16], 1);
362  Vmath::Smul(n, 1.0 / F / V_i, outarray[18], 1, outarray[18], 1);
363 
364  // I_tr
365  Array<OneD, NekDouble> &tmp_I_tr = outarray[5];
366  Vmath::Vsub(n, inarray[20], 1, inarray[19], 1, tmp_I_tr, 1);
367  Vmath::Smul(n, 1.0 / tau_tr, tmp_I_tr, 1, tmp_I_tr, 1);
368 
369  // I_up_leak
370  Array<OneD, NekDouble> &tmp_I_up_leak = outarray[6];
371  Vmath::Smul(n, NSR_I_up_max / NSR_I_Ca_max, inarray[20], 1, tmp_I_up_leak,
372  1);
373 
374  // I_up
375  Array<OneD, NekDouble> &tmp_I_up = outarray[7];
376  Vmath::Sdiv(n, NSR_K_up, inarray[17], 1, tmp_I_up, 1);
377  Vmath::Sadd(n, 1.0, tmp_I_up, 1, tmp_I_up, 1);
378  Vmath::Sdiv(n, NSR_I_up_max, tmp_I_up, 1, tmp_I_up, 1);
379 
380  // I_rel
381  Array<OneD, NekDouble> &tmp_I_rel = outarray[8];
382  Vmath::Vsub(n, inarray[19], 1, inarray[17], 1, tmp_I_rel, 1);
383  Vmath::Vmul(n, tmp_I_rel, 1, inarray[13], 1, tmp_I_rel, 1);
384  Vmath::Vmul(n, tmp_I_rel, 1, inarray[13], 1, tmp_I_rel, 1);
385  Vmath::Vmul(n, tmp_I_rel, 1, inarray[14], 1, tmp_I_rel, 1);
386  Vmath::Vmul(n, tmp_I_rel, 1, inarray[15], 1, tmp_I_rel, 1);
387  Vmath::Smul(n, JSR_K_rel, tmp_I_rel, 1, tmp_I_rel, 1);
388 
389  // B1
390  Array<OneD, NekDouble> &tmp_B1 = outarray[9];
391  Vmath::Svtvm(n, 2.0, tmp_I_Na_Ca, 1, tmp_I_p_Ca, 1, tmp_B1, 1);
392  Vmath::Vsub(n, tmp_B1, 1, tmp_I_Ca_L, 1, tmp_B1, 1);
393  Vmath::Vsub(n, tmp_B1, 1, tmp_I_b_Ca, 1, tmp_B1, 1);
394  Vmath::Smul(n, 0.5 / F, tmp_B1, 1, tmp_B1, 1);
395  Vmath::Svtvp(n, JSR_V_up, tmp_I_up_leak, 1, tmp_B1, 1, tmp_B1, 1);
396  Vmath::Svtvp(n, -JSR_V_up, tmp_I_up, 1, tmp_B1, 1, tmp_B1, 1);
397  Vmath::Svtvp(n, JSR_V_rel, tmp_I_rel, 1, tmp_B1, 1, tmp_B1, 1);
398  Vmath::Smul(n, 1.0 / V_i, tmp_B1, 1, tmp_B1, 1);
399 
400  // B2
401  Array<OneD, NekDouble> &tmp_B2 = outarray[10];
402  Vmath::Sadd(n, Km_Cmdn, inarray[17], 1, tmp_B2, 1);
403  Vmath::Vmul(n, tmp_B2, 1, tmp_B2, 1, tmp_B2, 1);
404  Vmath::Sdiv(n, Cmdn_max * Km_Cmdn, tmp_B2, 1, tmp_B2, 1);
405  Vmath::Sadd(n, Km_Trpn, inarray[17], 1, tmp, 1);
406  Vmath::Vmul(n, tmp, 1, tmp, 1, tmp, 1);
407  Vmath::Sdiv(n, Trpn_max * Km_Trpn, tmp, 1, tmp, 1);
408  Vmath::Vadd(n, tmp, 1, tmp_B2, 1, tmp_B2, 1);
409  Vmath::Sadd(n, 1.0, tmp_B2, 1, tmp_B2, 1);
410 
411  // Calcium concentration (18)
412  Vmath::Vdiv(n, tmp_B1, 1, tmp_B2, 1, outarray[17], 1);
413 
414  // Calcium up (21)
415  Vmath::Vsub(n, tmp_I_up, 1, tmp_I_up_leak, 1, outarray[20], 1);
416  Vmath::Svtvp(n, -JSR_V_rel / JSR_V_up, tmp_I_tr, 1, outarray[20], 1,
417  outarray[20], 1);
418 
419  // Calcium rel (20)
420  Vmath::Vsub(n, tmp_I_tr, 1, tmp_I_rel, 1, tmp, 1);
421  Vmath::Sadd(n, Km_Csqn, inarray[19], 1, outarray[19], 1);
422  Vmath::Vmul(n, outarray[19], 1, outarray[19], 1, outarray[19], 1);
423  Vmath::Sdiv(n, Csqn_max * Km_Csqn, outarray[19], 1, outarray[19], 1);
424  Vmath::Sadd(n, 1.0, outarray[19], 1, outarray[19], 1);
425  Vmath::Vdiv(n, tmp, 1, outarray[19], 1, outarray[19], 1);
426 
427  // Process gating variables
428  const NekDouble *v;
429  const NekDouble *x;
430  NekDouble *x_tau;
431  NekDouble *x_new;
432  // m
433  for (i = 0, v = &inarray[0][0], x = &inarray[1][0], x_new = &outarray[1][0],
434  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
438  : (0.32 * (*v + 47.13)) /
439  (1.0 - exp((-0.1) * (*v + 47.13)));
440  beta = 0.08 * exp(-(*v) / 11.0);
441  *x_tau = 1.0 / (alpha + beta);
442  *x_new = alpha * (*x_tau);
443  }
444  // h
445  for (i = 0, v = &inarray[0][0], x = &inarray[2][0], x_new = &outarray[2][0],
446  x_tau = &m_gates_tau[1][0];
447  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
448  {
449  alpha = (*v >= -40.0) ? 0.0 : 0.135 * exp(-((*v) + 80.0) / 6.8);
450  beta = (*v >= -40.0)
451  ? 1.0 / (0.13 * (1.0 + exp(-(*v + 10.66) / 11.1)))
452  : 3.56 * exp(0.079 * (*v)) + 310000.0 * exp(0.35 * (*v));
453  *x_tau = 1.0 / (alpha + beta);
454  *x_new = alpha * (*x_tau);
455  }
456  // j
457  for (i = 0, v = &inarray[0][0], x = &inarray[3][0], x_new = &outarray[3][0],
458  x_tau = &m_gates_tau[2][0];
459  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
460  {
461  alpha =
462  (*v >= -40.0)
463  ? 0.0
464  : (-127140.0 * exp(0.2444 * (*v)) -
465  3.474e-05 * exp(-0.04391 * (*v))) *
466  (((*v) + 37.78) / (1.0 + exp(0.311 * ((*v) + 79.23))));
467  beta = (*v >= -40.0) ? (0.3 * exp(-2.535e-07 * (*v)) /
468  (1.0 + exp(-0.1 * (*v + 32.0))))
469  : 0.1212 * exp(-0.01052 * (*v)) /
470  (1.0 + exp(-0.1378 * (*v + 40.14)));
471  *x_tau = 1.0 / (alpha + beta);
472  *x_new = alpha * (*x_tau);
473  }
474  // oa
475  for (i = 0, v = &inarray[0][0], x = &inarray[4][0], x_new = &outarray[4][0],
476  x_tau = &m_gates_tau[3][0];
477  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
478  {
479  alpha = 0.65 / (exp(-(*v + 10.0) / 8.5) + exp(-(*v - 30.0) / 59.0));
480  beta = 0.65 / (2.5 + exp((*v + 82.0) / 17.0));
481  *x_tau = 1.0 / K_Q10 / (alpha + beta);
482  *x_new = (1.0 / (1.0 + exp(-(*v + 20.47) / 17.54)));
483  }
484  // oi
485  for (i = 0, v = &inarray[0][0], x = &inarray[5][0], x_new = &outarray[5][0],
486  x_tau = &m_gates_tau[4][0];
487  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
488  {
489  alpha = 1.0 / (18.53 + exp((*v + 113.7) / 10.95));
490  beta = 1.0 / (35.56 + exp(-(*v + 1.26) / 7.44));
491  *x_tau = 1.0 / K_Q10 / (alpha + beta);
492  *x_new = (1.0 / (1.0 + exp((*v + 43.1) / 5.3)));
493  }
494  // ua
495  for (i = 0, v = &inarray[0][0], x = &inarray[6][0], x_new = &outarray[6][0],
496  x_tau = &m_gates_tau[5][0];
497  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
498  {
499  alpha = 0.65 / (exp(-(*v + 10.0) / 8.5) + exp(-(*v - 30.0) / 59.0));
500  beta = 0.65 / (2.5 + exp((*v + 82.0) / 17.0));
501  *x_tau = 1.0 / K_Q10 / (alpha + beta);
502  *x_new = 1.0 / (1.0 + exp(-(*v + 30.3) / 9.6));
503  }
504  // ui
505  for (i = 0, v = &inarray[0][0], x = &inarray[7][0], x_new = &outarray[7][0],
506  x_tau = &m_gates_tau[6][0];
507  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
508  {
509  alpha = 1.0 / (21.0 + exp(-(*v - 185.0) / 28.0));
510  beta = exp((*v - 158.0) / 16.0);
511  *x_tau = 1.0 / K_Q10 / (alpha + beta);
512  *x_new = 1.0 / (1.0 + exp((*v - 99.45) / 27.48));
513  }
514  // xr
515  for (i = 0, v = &inarray[0][0], x = &inarray[8][0], x_new = &outarray[8][0],
516  x_tau = &m_gates_tau[7][0];
517  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
518  {
519  alpha = 0.0003 * (*v + 14.1) / (1 - exp(-(*v + 14.1) / 5.0));
520  beta = 7.3898e-5 * (*v - 3.3328) / (exp((*v - 3.3328) / 5.1237) - 1.0);
521  *x_tau = 1.0 / (alpha + beta);
522  *x_new = 1.0 / (1 + exp(-(*v + 14.1) / 6.5));
523  }
524  // xs
525  for (i = 0, v = &inarray[0][0], x = &inarray[9][0], x_new = &outarray[9][0],
526  x_tau = &m_gates_tau[8][0];
527  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
528  {
529  alpha = 4e-5 * (*v - 19.9) / (1.0 - exp(-(*v - 19.9) / 17.0));
530  beta = 3.5e-5 * (*v - 19.9) / (exp((*v - 19.9) / 9.0) - 1.0);
531  *x_tau = 0.5 / (alpha + beta);
532  *x_new = 1.0 / sqrt(1.0 + exp(-(*v - 19.9) / 12.7));
533  }
534  // d
535  for (i = 0, v = &inarray[0][0], x = &inarray[10][0],
536  x_new = &outarray[10][0], x_tau = &m_gates_tau[9][0];
537  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
538  {
539  *x_tau = (1 - exp(-(*v + 10.0) / 6.24)) /
540  (0.035 * (*v + 10.0) * (1 + exp(-(*v + 10.0) / 6.24)));
541  *x_new = 1.0 / (1.0 + exp(-(*v + 10) / 8.0));
542  }
543  // f
544  for (i = 0, v = &inarray[0][0], x = &inarray[11][0],
545  x_new = &outarray[11][0], x_tau = &m_gates_tau[10][0];
546  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
547  {
548  // alpha = 1.0/(1.0 + exp((*v+28.0)/6.9));
549  *x_tau =
550  9.0 /
551  (0.0197 * exp(-0.0337 * 0.0337 * (*v + 10.0) * (*v + 10.0)) + 0.02);
552  *x_new = exp((-(*v + 28.0)) / 6.9) / (1.0 + exp((-(*v + 28.0)) / 6.9));
553  }
554  // f_Ca
555  for (i = 0, v = &inarray[0][0], x = &inarray[12][0],
556  x_new = &outarray[12][0], x_tau = &m_gates_tau[11][0];
557  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
558  {
559  *x_tau = 2.0;
560  *x_new = 1.0 / (1.0 + inarray[17][i] / 0.00035);
561  }
562 
563  Array<OneD, NekDouble> &tmp_Fn = outarray[15];
564  Vmath::Svtsvtp(n, 0.5 * 5e-13 / F, tmp_I_Ca_L, 1, -0.2 * 5e-13 / F,
565  tmp_I_Na_Ca, 1, tmp_Fn, 1);
566  Vmath::Svtvm(n, 1e-12 * JSR_V_rel, tmp_I_rel, 1, tmp_Fn, 1, tmp_Fn, 1);
567 
568  // u
569  for (i = 0, v = &tmp_Fn[0], x = &inarray[13][0], x_new = &outarray[13][0],
570  x_tau = &m_gates_tau[12][0];
571  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
572  {
573  *x_tau = 8.0;
574  *x_new = 1.0 / (1.0 + exp(-(*v - 3.4175e-13) / 1.367e-15));
575  }
576  // v
577  for (i = 0, v = &tmp_Fn[0], x = &inarray[14][0], x_new = &outarray[14][0],
578  x_tau = &m_gates_tau[13][0];
579  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
580  {
581  *x_tau = 1.91 + 2.09 / (1.0 + exp(-(*v - 3.4175e-13) / 13.67e-16));
582  *x_new = 1.0 - 1.0 / (1.0 + exp(-(*v - 6.835e-14) / 13.67e-16));
583  }
584  // w
585  for (i = 0, v = &inarray[0][0], x = &inarray[15][0],
586  x_new = &outarray[15][0], x_tau = &m_gates_tau[14][0];
587  i < n; ++i, ++v, ++x, ++x_new, ++x_tau)
588  {
589  *x_tau = 6.0 * (1.0 - exp(-(*v - 7.9) / 5.0)) /
590  (1.0 + 0.3 * exp(-(*v - 7.9) / 5.0)) / (*v - 7.9);
591  *x_new = 1.0 - 1.0 / (1.0 + exp(-(*v - 40.0) / 17.0));
592  }
593 }
594 
595 /**
596  *
597  */
598 void CourtemancheRamirezNattel98::v_GenerateSummary(SummaryList &s)
599 {
600  SolverUtils::AddSummaryItem(s, "Cell model", "CourtemancheRamirezNattel98");
601  SolverUtils::AddSummaryItem(s, "Cell model var.", lookupIds[model_variant]);
602 }
603 
604 void CourtemancheRamirezNattel98::v_SetInitialConditions()
605 {
606  Vmath::Fill(m_nq, -81.0, m_cellSol[0], 1);
607  Vmath::Fill(m_nq, 2.908e-03, m_cellSol[1], 1);
608  Vmath::Fill(m_nq, 9.649e-01, m_cellSol[2], 1);
609  Vmath::Fill(m_nq, 9.775e-01, m_cellSol[3], 1);
610  Vmath::Fill(m_nq, 3.043e-02, m_cellSol[4], 1);
611  Vmath::Fill(m_nq, 9.992e-01, m_cellSol[5], 1);
612  Vmath::Fill(m_nq, 4.966e-03, m_cellSol[6], 1);
613  Vmath::Fill(m_nq, 9.986e-01, m_cellSol[7], 1);
614  Vmath::Fill(m_nq, 3.296e-05, m_cellSol[8], 1);
615  Vmath::Fill(m_nq, 1.869e-02, m_cellSol[9], 1);
616  Vmath::Fill(m_nq, 1.367e-04, m_cellSol[10], 1);
617  Vmath::Fill(m_nq, 9.996e-01, m_cellSol[11], 1);
618  Vmath::Fill(m_nq, 7.755e-01, m_cellSol[12], 1);
619  Vmath::Fill(m_nq, 2.35e-112, m_cellSol[13], 1);
620  Vmath::Fill(m_nq, 1.0, m_cellSol[14], 1);
621  Vmath::Fill(m_nq, 0.9992, m_cellSol[15], 1);
622  Vmath::Fill(m_nq, 1.117e+01, m_cellSol[16], 1);
623  Vmath::Fill(m_nq, 1.013e-04, m_cellSol[17], 1);
624  Vmath::Fill(m_nq, 1.39e+02, m_cellSol[18], 1);
625  Vmath::Fill(m_nq, 1.488, m_cellSol[19], 1);
626  Vmath::Fill(m_nq, 1.488, m_cellSol[20], 1);
627 }
628 
629 std::string CourtemancheRamirezNattel98::v_GetCellVarName(unsigned int idx)
630 {
631  switch (idx)
632  {
633  case 0:
634  return "u";
635  case 1:
636  return "m";
637  case 2:
638  return "h";
639  case 3:
640  return "j";
641  case 4:
642  return "o_a";
643  case 5:
644  return "o_i";
645  case 6:
646  return "u_a";
647  case 7:
648  return "u_i";
649  case 8:
650  return "x_r";
651  case 9:
652  return "x_s";
653  case 10:
654  return "d";
655  case 11:
656  return "f";
657  case 12:
658  return "f_Ca";
659  case 13:
660  return "U";
661  case 14:
662  return "V";
663  case 15:
664  return "W";
665  case 16:
666  return "Na_i";
667  case 17:
668  return "Ca_i";
669  case 18:
670  return "K_i";
671  case 19:
672  return "Ca_rel";
673  case 20:
674  return "Ca_up";
675  default:
676  return "unknown";
677  }
678 }
679 
680 } // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::CellModel
Cell model base class.
Definition: CellModel.h:66
Nektar::CellModel::m_cellSol
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:126
Nektar::CellModel::m_nq
int m_nq
Number of physical points.
Definition: CellModel.h:117
Nektar::CellModel::m_concentrations
std::vector< int > m_concentrations
Indices of cell model variables which are concentrations.
Definition: CellModel.h:139
Nektar::CellModel::m_gates
std::vector< int > m_gates
Indices of cell model variables which are gates.
Definition: CellModel.h:141
Nektar::CellModel::m_gates_tau
Array< OneD, Array< OneD, NekDouble > > m_gates_tau
Storage for gate tau values.
Definition: CellModel.h:143
Nektar::CellModel::m_nvar
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:119
Nektar::CourtemancheRamirezNattel98::v_GetCellVarName
virtual std::string v_GetCellVarName(unsigned int idx)
Definition: CourtemancheRamirezNattel98.cpp:629
Nektar::CourtemancheRamirezNattel98::v_GenerateSummary
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
Definition: CourtemancheRamirezNattel98.cpp:598
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:160
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
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:48
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:49
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:751
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:209
Vmath::Vlog
void Vlog(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:114
Vmath::Vexp
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:125
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:622
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:664
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:359
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:248
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:324
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:284
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
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:384
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:419
Vmath::Vpow
void Vpow(int n, const T *x, const int incx, const T f, T *y, const int incy)
Definition: Vmath.hpp:136
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Generated on Tue Aug 23 2022 07:27:12 for Nektar++ by doxygen 1.9.1