br_mod.f90

Go to the documentation of this file.

module br_mod

  use real_type, only: rp

  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: br_ny, br_nw, br_ni

  public :: br_y0, br_ilist
  public :: br_ab_0, br_ab_w, br_ab_wv

  public :: br_sp_y0, br_sp_ilist
  public :: br_sp_ab_w

  !! Dimensions
  integer, parameter :: br_ny = 8        
  integer, parameter :: br_nw = 6
  integer, parameter :: br_ni = 4

  !! IList description
  integer, parameter :: jna = 1
  integer, parameter :: js  = 2
  integer, parameter :: jx1 = 3
  integer, parameter :: jk1 = 4

  !! Y description
  integer, parameter :: ym   = 1
  integer, parameter :: yh   = 2 
  integer, parameter :: yj   = 3
  integer, parameter :: yd   = 4
  integer, parameter :: yf   = 5
  integer, parameter :: yx1  = 6
  integer, parameter :: ycai = 7
  integer, parameter :: yv   = 8

  ! Model constants
  ! R = 8314.472    [1E-3 J.K^{-1}.mol^{-1}]
  ! F =  96485.3415 [A.s.mol^{-1}]
  ! T = 310.0       [K]
  ! RTsF = RT / F   [mV]
  real(RP), parameter :: rtsf  = 26.713760659695648_rp 
  real(RP), parameter :: fsrt  = 1._rp/(rtsf)   
  real(RP), parameter :: cm    = 1.0_rp     ! [muF] 
  real(RP), parameter :: g_na  = 4.0_rp
  real(RP), parameter :: g_nac = 0.003_rp
  real(RP), parameter :: g_s   = 0.09_rp    
  real(RP), parameter :: ena   = 50.0_rp
  real(RP), parameter :: c_ix1 = exp(0.04_rp*(77._rp-35._rp))

contains

  subroutine br_y0(Y) 
    real(RP), dimension(BR_NY), intent(out) :: Y

    y(ym)   = 1.09819687232595410246167255529228788e-0002_rp
    y(yh)   = 0.987721175487584100912586864929795661_rp   
    y(yj)   = 0.974838138981649082196488362602052421_rp
    y(yd)   = 2.97072466320784053758100507598202674e-0003_rp
    y(yf)   = 0.999981333893687946493220122499403334_rp
    y(yx1)  = 5.62865057052210353702214502232095934e-0003_rp
    y(ycai) = 1.78200721561980600333356651988516442e-0007_rp
    y(yv)  = -84.5737561222608724445622259686443887_rp

  end subroutine br_y0


  !       %----------------------------------------%
  !       |                                        |
  !       |          REACTION TERMS DEF.           |
  !       |                                        |
  !       %----------------------------------------%

  subroutine br_gating_variables(a, b, V) 
    real(RP), dimension(BR_NW), intent(out) :: a, b
    real(RP)                  , intent(in)  :: V

    b(ym) = -(v + 47._rp)/(exp(-0.1_rp * (v + 47._rp)) - 1._rp)
    a(ym)  = 40._rp * exp(-0.056_rp*(v + 72._rp))
    
    b(yh) = 0.126_rp * exp(-0.25_rp*(v + 77._rp))
    a(yh)  = 1.7_rp / (exp(-0.082_rp*(v+22.5_rp))+1._rp)
    
    b(yj) = 0.055_rp*exp(-0.25_rp*(v + 78._rp))/&
         &      (exp(-0.2_rp*(v+78._rp))+1._rp)
    a(yj)  = 0.3_rp / (exp( -0.1_rp*(v+32._rp))+1._rp)
    
    b(yd) = 0.095_rp*exp(-0.01_rp*(v-5._rp))/&
         &      (1._rp+exp(-0.072_rp*(v-5._rp)))
    a(yd)  = 0.07_rp*exp(-0.017_rp*(v+44._rp))/&
         &      (1._rp+exp(0.05_rp*(v+44._rp)))
    
    b(yf) = 0.012_rp*exp(-0.008_rp*(v+28._rp))/&
         &      (1._rp+exp(0.15_rp*(v+28._rp)))
    a(yf)  = 0.0065_rp*exp(-0.02_rp*(v+30._rp))/&
         &      (1._rp+exp(-0.2_rp*(v+30._rp)))

    b(yx1) = 0.0005_rp*exp(0.083_rp*(v+50._rp))/&
         &       (1._rp+exp(0.057_rp*(v+50._rp)))
    a(yx1)  = 0.0013_rp*exp(-0.06_rp*(v+20._rp))/&
         &       (1._rp+exp(-0.04_rp*(v+20._rp)))

    a = -( a + b )

  end subroutine br_gating_variables


  subroutine br_ab_0(a, b, I, Y, N, Na) 
    real(RP), dimension(Na), intent(out) :: a
    real(RP), dimension(N) , intent(out) :: b
    real(RP)               , intent(in)  :: I
    real(RP), dimension(N) , intent(in)  :: Y
    integer                , intent(in)  :: N, Na 

    real(RP), dimension(BR_NW) :: alpha

    ! ionic currents temporarily stored in b(1:NI)
    call br_ilist(b(1:br_ni), y(1:br_ny))

    ! Cai
    b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))

    ! vm
    b(yv) = (- sum( b(1:br_ni) ) + i )/cm

    ! gating variables
    call br_gating_variables(alpha, b(1:br_nw), y(yv)) 
    b(1:br_nw) = b(1:br_nw) + alpha*y(1:br_nw)
    
  end subroutine br_ab_0


  subroutine br_ab_w(a, b, I, Y, N, Na) 
    real(RP), dimension(Na), intent(out) :: a
    real(RP), dimension(N) , intent(out) :: b
    real(RP)               , intent(in)  :: I
    real(RP), dimension(N) , intent(in)  :: Y
    integer                , intent(in)  :: N, Na 

    ! ionic currents temporarily stored in b(1:NI)
    call br_ilist(b(1:br_ni), y(1:br_ny))

    ! Cai
    b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))

    ! vm
    b(yv) = (- sum( b(1:br_ni) ) + i ) / cm

    ! gating variables
    call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv)) 
    
  end subroutine br_ab_w


  subroutine br_ab_wv(a, b, I, Y, N, Na) 
    real(RP), dimension(Na), intent(out) :: a
    real(RP), dimension(N) , intent(out) :: b
    real(RP)               , intent(in)  :: I
    real(RP), dimension(N) , intent(in)  :: Y
    integer                , intent(in)  :: N, Na 

    ! Cai
    ! I_s temporarily stored in b(1)
    b(1)    = i_s( y )
    b(ycai) = -1.e-7_rp*b(1) + 0.07_rp*(1.e-7_rp - y(ycai))
    a(ycai) = 0._rp

    ! vm
    call itot_dv(a(yv), b(yv), y)
    a(yv) = - a(yv) / cm
    b(yv) = (- b(yv) + i ) / cm

    ! gating variables
    call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv)) 

  end subroutine br_ab_wv



  !       %----------------------------------------%
  !       |                                        |
  !       |          IONIC CURRENTS                |
  !       |                                        |
  !       %----------------------------------------%
  subroutine br_ilist(IList, Y) 
    real(RP), dimension(BR_NI), intent(out) :: IList
    real(RP), dimension(BR_NY), intent(in)  :: Y

    ilist(jna) = i_na(y)
    ilist(js ) = i_s(y)
    ilist(jk1) = i_k1(y)
    ilist(jx1) = i_x1(y)

  end subroutine br_ilist

  !! Total current
  !! 
  !!   Itot = a*V + b
  subroutine itot_dv(a, b, Y)
    real(RP), intent(in), dimension(BR_NY) :: Y
    real(RP) :: a, b

    ! first def : stabilization is on I_Na only
    call i_na_dv(a, b, y)
    b =  b + i_k1(y) + i_s(y) + i_x1(y) 

  end subroutine itot_dv

  !! current I_K1
  !!
  function i_k1(Y)
    real(RP), dimension(BR_NY), intent(in) :: Y
    real(RP) :: I_K1

    real(RP) :: f

    ! f1 = ff(Y(yv), 0.04_RP, 85._RP)
    ! f2 = ff(Y(yv), 0.04_RP, 53._RP)
    ! f3 = ff(Y(yv),-0.04_RP, 23._RP)

    f    = exp(0.04_rp * (y(yv) + 85._rp) )
    i_k1 = 4._rp*( f  - 1._rp )

    f    = exp(0.04_rp * (y(yv) + 53._rp) )
    i_k1 = i_k1 / ( f + f**2 ) 

    f    = exp(-0.04_rp * (y(yv) + 23._rp) )
    i_k1 = i_k1 + 0.2_rp * ( y(yv) + 23._rp) / ( 1._rp - f )

    i_k1 = i_k1 * 0.35_rp

  end function i_k1

  !! current I_S
  !!
  function i_s(Y)
    real(RP), dimension(BR_NY), intent(in) :: Y
    real(RP) :: I_s

    i_s  = -82.3_rp - 13.0287_rp*log( y(ycai) )
    i_s  = g_s * y(yd) * y(yf) * (y(yv) - i_s)

  end function i_s


  !! current I_x1
  !!
  function i_x1(Y)
    real(RP), dimension(BR_NY), intent(in) :: Y
    real(RP) :: I_x1

    ! I_X1 = C_Ix1 - ff( Y(yv), -0.04_RP, 35._RP)

    i_x1 = c_ix1 - exp( -0.04_rp * (y(yv) + 35._rp) )
    i_x1 = i_x1 * y(yx1) * 0.8_rp

  end function i_x1


  !! current i_na
  !!
  function i_na(Y)
    real(RP), dimension(BR_NY), intent(in) :: Y
    real(RP) :: I_Na

    i_na = g_na * y(ym)**3 * y(yh) * y(yj) + g_nac
    i_na = i_na * ( y(yv) - ena )

  end function i_na

  !!
  !! I_Na = a*V + b, a = d(INa)/dV
  !!
  subroutine i_na_dv(a, b, Y)
    real(RP), intent(in), dimension(BR_NY) :: Y
    real(RP) :: a, b

    a = g_na * y(ym)**3 * y(yh) * y(yj) + g_nac
    b =-a * ena

  end subroutine i_na_dv

  function ff(x, a, b)
    real(RP), intent(in) :: x, a, b
    real(RP) :: ff

    ff = exp(a*(x+b))

  end function ff



  !       %----------------------------------------%
  !       |                                        |
  !       |   MODIFIED BR FOR SPIRAL WAVES         |
  !       |                                        |
  !       %----------------------------------------%
  !! current I_S, modified value for 'g_s'
  !!
  function i_s_sp(Y) result(I_s)
    real(RP), dimension(BR_NY), intent(in) :: Y
    real(RP) :: I_s

    real(RP), parameter :: g_S_modified   = 0.01_rp    

    i_s  = -82.3_rp - 13.0287_rp*log( y(ycai) )
    i_s  = g_s_modified * y(yd) * y(yf) * (y(yv) - i_s)

  end function i_s_sp

  subroutine br_sp_ilist(IList, Y) 
    real(RP), dimension(BR_NI), intent(out) :: IList
    real(RP), dimension(BR_NY), intent(in)  :: Y

    ilist(jna) = i_na(y)
    ilist(js ) = i_s_sp(y)
    ilist(jk1) = i_k1(y)
    ilist(jx1) = i_x1(y)

  end subroutine br_sp_ilist


  subroutine br_sp_y0(Y) 
    real(RP), dimension(BR_NY) :: Y

    y(ym)   = 1.06278410016131944945954637794827772e-0002_rp 
    y(yh)   = 0.988691698520646794656617541431759221_rp
    y(yj)   = 0.975989197046760783619339482293370610_rp  
    y(yd)   = 2.90835646511938874988339682604803243e-0003_rp
    y(yf)   = 0.999982195505595856678047191499808037_rp
    y(yx1)  = 5.49126107877021382608497656435722875e-0003_rp
    y(ycai) = 1.17485016892666208390185442837442507e-0007_rp
    y(yv)   =-84.8250996515381622468412766225817842_rp  

  end subroutine br_sp_y0

  subroutine br_sp_ab_w(a, b, I, Y, N, Na) 
    real(RP), dimension(Na), intent(out) :: a
    real(RP), dimension(N) , intent(out) :: b
    real(RP)               , intent(in)  :: I
    real(RP), dimension(N) , intent(in)  :: Y
    integer                , intent(in)  :: N, Na 

    ! ionic currents temporarily stored in b(1:NI)
    call br_sp_ilist(b(1:br_ni), y)

    ! Cai
    b(ycai) = -1.e-7_rp*b(js) + 0.07_rp*(1.e-7_rp - y(ycai))

    ! vm
    b(yv) = ( - sum(b(1:br_ni)) + i ) / cm 

    ! gating variables
    call br_gating_variables( a(1:br_nw), b(1:br_nw), y(yv)) 
    
  end subroutine br_sp_ab_w


end module br_mod