Choral
br_mod.f90
Go to the documentation of this file.
1 
28 module br_mod
29 
30  use real_type, only: rp
31 
32  implicit none
33  private
34 
35 
36  ! %----------------------------------------%
37  ! | |
38  ! | PUBLIC DATA |
39  ! | |
40  ! %----------------------------------------%
41  public :: br_ny, br_nw, br_ni
42 
43  public :: br_y0, br_ilist
44  public :: br_ab_0, br_ab_w, br_ab_wv
45 
46  public :: br_sp_y0, br_sp_ilist
47  public :: br_sp_ab_w
48 
49  !! Dimensions
50  integer, parameter :: br_ny = 8
51  integer, parameter :: br_nw = 6
52  integer, parameter :: br_ni = 4
53 
54  !! IList description
55  integer, parameter :: jna = 1
56  integer, parameter :: js = 2
57  integer, parameter :: jx1 = 3
58  integer, parameter :: jk1 = 4
59 
60  !! Y description
61  integer, parameter :: ym = 1
62  integer, parameter :: yh = 2
63  integer, parameter :: yj = 3
64  integer, parameter :: yd = 4
65  integer, parameter :: yf = 5
66  integer, parameter :: yx1 = 6
67  integer, parameter :: ycai = 7
68  integer, parameter :: yv = 8
69 
70  ! Model constants
71  ! R = 8314.472 [1E-3 J.K^{-1}.mol^{-1}]
72  ! F = 96485.3415 [A.s.mol^{-1}]
73  ! T = 310.0 [K]
74  ! RTsF = RT / F [mV]
75  real(RP), parameter :: rtsf = 26.713760659695648_rp
76  real(RP), parameter :: fsrt = 1._rp/(rtsf)
77  real(RP), parameter :: cm = 1.0_rp ! [muF]
78  real(RP), parameter :: g_na = 4.0_rp
79  real(RP), parameter :: g_nac = 0.003_rp
80  real(RP), parameter :: g_s = 0.09_rp
81  real(RP), parameter :: ena = 50.0_rp
82  real(RP), parameter :: c_ix1 = exp(0.04_rp*(77._rp-35._rp))
83 
84 contains
85 
88  subroutine br_y0(Y)
89  real(RP), dimension(BR_NY), intent(out) :: Y
90 
91  y(ym) = 1.09819687232595410246167255529228788e-0002_rp
92  y(yh) = 0.987721175487584100912586864929795661_rp
93  y(yj) = 0.974838138981649082196488362602052421_rp
94  y(yd) = 2.97072466320784053758100507598202674e-0003_rp
95  y(yf) = 0.999981333893687946493220122499403334_rp
96  y(yx1) = 5.62865057052210353702214502232095934e-0003_rp
97  y(ycai) = 1.78200721561980600333356651988516442e-0007_rp
98  y(yv) = -84.5737561222608724445622259686443887_rp
99 
100  end subroutine br_y0
101 
102 
103  ! %----------------------------------------%
104  ! | |
105  ! | REACTION TERMS DEF. |
106  ! | |
107  ! %----------------------------------------%
108 
118  subroutine br_gating_variables(a, b, V)
119  real(RP), dimension(BR_NW), intent(out) :: a, b
120  real(RP) , intent(in) :: V
121 
122  b(ym) = -(v + 47._rp)/(exp(-0.1_rp * (v + 47._rp)) - 1._rp)
123  a(ym) = 40._rp * exp(-0.056_rp*(v + 72._rp))
124 
125  b(yh) = 0.126_rp * exp(-0.25_rp*(v + 77._rp))
126  a(yh) = 1.7_rp / (exp(-0.082_rp*(v+22.5_rp))+1._rp)
127 
128  b(yj) = 0.055_rp*exp(-0.25_rp*(v + 78._rp))/&
129  & (exp(-0.2_rp*(v+78._rp))+1._rp)
130  a(yj) = 0.3_rp / (exp( -0.1_rp*(v+32._rp))+1._rp)
131 
132  b(yd) = 0.095_rp*exp(-0.01_rp*(v-5._rp))/&
133  & (1._rp+exp(-0.072_rp*(v-5._rp)))
134  a(yd) = 0.07_rp*exp(-0.017_rp*(v+44._rp))/&
135  & (1._rp+exp(0.05_rp*(v+44._rp)))
136 
137  b(yf) = 0.012_rp*exp(-0.008_rp*(v+28._rp))/&
138  & (1._rp+exp(0.15_rp*(v+28._rp)))
139  a(yf) = 0.0065_rp*exp(-0.02_rp*(v+30._rp))/&
140  & (1._rp+exp(-0.2_rp*(v+30._rp)))
141 
142  b(yx1) = 0.0005_rp*exp(0.083_rp*(v+50._rp))/&
143  & (1._rp+exp(0.057_rp*(v+50._rp)))
144  a(yx1) = 0.0013_rp*exp(-0.06_rp*(v+20._rp))/&
145  & (1._rp+exp(-0.04_rp*(v+20._rp)))
146 
147  a = -( a + b )
148 
149  end subroutine br_gating_variables
150 
151 
158  subroutine br_ab_0(a, b, I, Y, N, Na)
159  real(RP), dimension(Na), intent(out) :: a
160  real(RP), dimension(N) , intent(out) :: b
161  real(RP) , intent(in) :: I
162  real(RP), dimension(N) , intent(in) :: Y
163  integer , intent(in) :: N, Na
164 
165  real(RP), dimension(BR_NW) :: alpha
166 
167  ! ionic currents temporarily stored in b(1:NI)
168  call br_ilist(b(1:br_ni), y(1:br_ny))
169 
170  ! Cai
171  b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))
172 
173  ! vm
174  b(yv) = (- sum( b(1:br_ni) ) + i )/cm
175 
176  ! gating variables
177  call br_gating_variables(alpha, b(1:br_nw), y(yv))
178  b(1:br_nw) = b(1:br_nw) + alpha*y(1:br_nw)
179 
180  end subroutine br_ab_0
181 
182 
188  subroutine br_ab_w(a, b, I, Y, N, Na)
189  real(RP), dimension(Na), intent(out) :: a
190  real(RP), dimension(N) , intent(out) :: b
191  real(RP) , intent(in) :: I
192  real(RP), dimension(N) , intent(in) :: Y
193  integer , intent(in) :: N, Na
194 
195  ! ionic currents temporarily stored in b(1:NI)
196  call br_ilist(b(1:br_ni), y(1:br_ny))
197 
198  ! Cai
199  b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))
200 
201  ! vm
202  b(yv) = (- sum( b(1:br_ni) ) + i ) / cm
203 
204  ! gating variables
205  call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv))
206 
207  end subroutine br_ab_w
208 
209 
216  subroutine br_ab_wv(a, b, I, Y, N, Na)
217  real(RP), dimension(Na), intent(out) :: a
218  real(RP), dimension(N) , intent(out) :: b
219  real(RP) , intent(in) :: I
220  real(RP), dimension(N) , intent(in) :: Y
221  integer , intent(in) :: N, Na
222 
223  ! Cai
224  ! I_s temporarily stored in b(1)
225  b(1) = i_s( y )
226  b(ycai) = -1.e-7_rp*b(1) + 0.07_rp*(1.e-7_rp - y(ycai))
227  a(ycai) = 0._rp
228 
229  ! vm
230  call itot_dv(a(yv), b(yv), y)
231  a(yv) = - a(yv) / cm
232  b(yv) = (- b(yv) + i ) / cm
233 
234  ! gating variables
235  call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv))
236 
237  end subroutine br_ab_wv
238 
239 
240 
241  ! %----------------------------------------%
242  ! | |
243  ! | IONIC CURRENTS |
244  ! | |
245  ! %----------------------------------------%
246  subroutine br_ilist(IList, Y)
247  real(RP), dimension(BR_NI), intent(out) :: IList
248  real(RP), dimension(BR_NY), intent(in) :: Y
249 
250  ilist(jna) = i_na(y)
251  ilist(js ) = i_s(y)
252  ilist(jk1) = i_k1(y)
253  ilist(jx1) = i_x1(y)
254 
255  end subroutine br_ilist
256 
257  !! Total current
258  !!
259  !! Itot = a*V + b
260  subroutine itot_dv(a, b, Y)
261  real(RP), intent(in), dimension(BR_NY) :: Y
262  real(RP) :: a, b
263 
264  ! first def : stabilization is on I_Na only
265  call i_na_dv(a, b, y)
266  b = b + i_k1(y) + i_s(y) + i_x1(y)
267 
268  end subroutine itot_dv
269 
270  !! current I_K1
271  !!
272  function i_k1(Y)
273  real(RP), dimension(BR_NY), intent(in) :: Y
274  real(RP) :: I_K1
275 
276  real(RP) :: f
277 
278  ! f1 = ff(Y(yv), 0.04_RP, 85._RP)
279  ! f2 = ff(Y(yv), 0.04_RP, 53._RP)
280  ! f3 = ff(Y(yv),-0.04_RP, 23._RP)
281 
282  f = exp(0.04_rp * (y(yv) + 85._rp) )
283  i_k1 = 4._rp*( f - 1._rp )
284 
285  f = exp(0.04_rp * (y(yv) + 53._rp) )
286  i_k1 = i_k1 / ( f + f**2 )
287 
288  f = exp(-0.04_rp * (y(yv) + 23._rp) )
289  i_k1 = i_k1 + 0.2_rp * ( y(yv) + 23._rp) / ( 1._rp - f )
290 
291  i_k1 = i_k1 * 0.35_rp
292 
293  end function i_k1
294 
295  !! current I_S
296  !!
297  function i_s(Y)
298  real(RP), dimension(BR_NY), intent(in) :: Y
299  real(RP) :: I_s
300 
301  i_s = -82.3_rp - 13.0287_rp*log( y(ycai) )
302  i_s = g_s * y(yd) * y(yf) * (y(yv) - i_s)
303 
304  end function i_s
305 
306 
307  !! current I_x1
308  !!
309  function i_x1(Y)
310  real(RP), dimension(BR_NY), intent(in) :: Y
311  real(RP) :: I_x1
312 
313  ! I_X1 = C_Ix1 - ff( Y(yv), -0.04_RP, 35._RP)
314 
315  i_x1 = c_ix1 - exp( -0.04_rp * (y(yv) + 35._rp) )
316  i_x1 = i_x1 * y(yx1) * 0.8_rp
317 
318  end function i_x1
319 
320 
321  !! current i_na
322  !!
323  function i_na(Y)
324  real(RP), dimension(BR_NY), intent(in) :: Y
325  real(RP) :: I_Na
326 
327  i_na = g_na * y(ym)**3 * y(yh) * y(yj) + g_nac
328  i_na = i_na * ( y(yv) - ena )
329 
330  end function i_na
331 
332  !!
333  !! I_Na = a*V + b, a = d(INa)/dV
334  !!
335  subroutine i_na_dv(a, b, Y)
336  real(RP), intent(in), dimension(BR_NY) :: Y
337  real(RP) :: a, b
338 
339  a = g_na * y(ym)**3 * y(yh) * y(yj) + g_nac
340  b =-a * ena
341 
342  end subroutine i_na_dv
343 
344  function ff(x, a, b)
345  real(RP), intent(in) :: x, a, b
346  real(RP) :: ff
347 
348  ff = exp(a*(x+b))
349 
350  end function ff
351 
352 
353 
354  ! %----------------------------------------%
355  ! | |
356  ! | MODIFIED BR FOR SPIRAL WAVES |
357  ! | |
358  ! %----------------------------------------%
359  !! current I_S, modified value for 'g_s'
360  !!
361  function i_s_sp(Y) result(I_s)
362  real(RP), dimension(BR_NY), intent(in) :: Y
363  real(RP) :: I_s
364 
365  real(RP), parameter :: g_S_modified = 0.01_rp
366 
367  i_s = -82.3_rp - 13.0287_rp*log( y(ycai) )
368  i_s = g_s_modified * y(yd) * y(yf) * (y(yv) - i_s)
369 
370  end function i_s_sp
371 
374  subroutine br_sp_ilist(IList, Y)
375  real(RP), dimension(BR_NI), intent(out) :: IList
376  real(RP), dimension(BR_NY), intent(in) :: Y
377 
378  ilist(jna) = i_na(y)
379  ilist(js ) = i_s_sp(y)
380  ilist(jk1) = i_k1(y)
381  ilist(jx1) = i_x1(y)
382 
383  end subroutine br_sp_ilist
384 
385 
388  subroutine br_sp_y0(Y)
389  real(RP), dimension(BR_NY) :: Y
390 
391  y(ym) = 1.06278410016131944945954637794827772e-0002_rp
392  y(yh) = 0.988691698520646794656617541431759221_rp
393  y(yj) = 0.975989197046760783619339482293370610_rp
394  y(yd) = 2.90835646511938874988339682604803243e-0003_rp
395  y(yf) = 0.999982195505595856678047191499808037_rp
396  y(yx1) = 5.49126107877021382608497656435722875e-0003_rp
397  y(ycai) = 1.17485016892666208390185442837442507e-0007_rp
398  y(yv) =-84.8250996515381622468412766225817842_rp
399 
400  end subroutine br_sp_y0
401 
407  subroutine br_sp_ab_w(a, b, I, Y, N, Na)
408  real(RP), dimension(Na), intent(out) :: a
409  real(RP), dimension(N) , intent(out) :: b
410  real(RP) , intent(in) :: I
411  real(RP), dimension(N) , intent(in) :: Y
412  integer , intent(in) :: N, Na
413 
414  ! ionic currents temporarily stored in b(1:NI)
415  call br_sp_ilist(b(1:br_ni), y)
416 
417  ! Cai
418  b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))
419 
420  ! vm
421  b(yv) = ( - sum(b(1:br_ni)) + i ) / cm
422 
423  ! gating variables
424  call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv))
425 
426  end subroutine br_sp_ab_w
427 
428 
429 end module br_mod
subroutine, public br_sp_y0(Y)
Model IONIC_BR_SP: rest state
Definition: br_mod.f90:389
integer, parameter yj
Definition: br_mod.f90:63
integer, parameter jna
Definition: br_mod.f90:55
integer, parameter yv
Definition: br_mod.f90:68
integer, parameter yd
Definition: br_mod.f90:64
subroutine, public br_y0(Y)
re-computed rest state (solving F(Y0)=0)
Definition: br_mod.f90:89
real(rp), parameter c_ix1
Definition: br_mod.f90:82
integer, parameter, public rp
real(kind=RP) = real precision in the code REAL_TOL = epsilon to test real equality ...
Definition: real_type.F90:90
real(rp) function i_na(Y)
Definition: br_mod.f90:324
integer, parameter ycai
Definition: br_mod.f90:67
integer, parameter, public br_nw
Definition: br_mod.f90:51
real(rp), parameter rtsf
Definition: br_mod.f90:75
integer, parameter ym
Definition: br_mod.f90:61
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
real(rp), parameter cm
Definition: br_mod.f90:77
integer, parameter, public br_ni
Definition: br_mod.f90:52
Beeler Reuter ionic model
Definition: br_mod.f90:28
integer, parameter jx1
Definition: br_mod.f90:57
subroutine, public br_ab_w(a, b, I, Y, N, Na)
Model IONIC_BR: reaction term
Definition: br_mod.f90:189
real(rp) function i_x1(Y)
Definition: br_mod.f90:310
integer, parameter yx1
Definition: br_mod.f90:66
integer, parameter js
Definition: br_mod.f90:56
real(rp) function e(x, v1, v2)
subroutine, public br_sp_ab_w(a, b, I, Y, N, Na)
Model IONIC_BR_SP: reaction term
Definition: br_mod.f90:408
integer, parameter jk1
Definition: br_mod.f90:58
subroutine, public br_sp_ilist(IList, Y)
Model IONIC_BR_SP: ionic current list
Definition: br_mod.f90:375
real(rp) function i_s_sp(Y)
Definition: br_mod.f90:362
real(rp) function i_s(Y)
Definition: br_mod.f90:298
real(rp), parameter fsrt
Definition: br_mod.f90:76
subroutine, public br_ilist(IList, Y)
Definition: br_mod.f90:247
integer, parameter, public br_ny
Definition: br_mod.f90:50
subroutine i_na_dv(a, b, Y)
Definition: br_mod.f90:336
integer, parameter yf
Definition: br_mod.f90:65
subroutine, public br_ab_0(a, b, I, Y, N, Na)
Model IONIC_BR: reaction term
Definition: br_mod.f90:159
real(rp), parameter g_nac
Definition: br_mod.f90:79
real(rp), parameter g_s
Definition: br_mod.f90:80
subroutine itot_dv(a, b, Y)
Definition: br_mod.f90:261
real(rp) function ff(x, a, b)
Definition: br_mod.f90:345
real(rp) function i_k1(Y)
Definition: br_mod.f90:273
subroutine br_gating_variables(a, b, V)
Model IONIC_BR: Gating variables dw_i/dt = (w_i-w_{i,~})/{i} =1,,{ NW} 1/ -w_{i,~})/
Definition: br_mod.f90:119
integer, parameter yh
Definition: br_mod.f90:62
real(rp), parameter g_na
Definition: br_mod.f90:78
real(rp), parameter ena
Definition: br_mod.f90:81
subroutine, public br_ab_wv(a, b, I, Y, N, Na)
Model IONIC_BR: reaction term
Definition: br_mod.f90:217