tnnp_mod.f90

Go to the documentation of this file.

module tnnp_mod

  use real_type, only: rp

  implicit none
  private

  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: tnnp_ni, tnnp_ny, tnnp_nw
  public :: tnnp_y0, tnnp_ab_0, tnnp_ab_w, tnnp_ilist


  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!                              !!
  !!      IList   DESCRIPTION     !!
  !!                              !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  integer, parameter :: tnnp_ni=  15    ! Number of currents
  integer, parameter :: jna    =  1
  integer, parameter :: jnaca  =  2  
  integer, parameter :: jnak   =  3  
  integer, parameter :: jbna   =  4  
  integer, parameter :: jcal   =  5  
  integer, parameter :: jpca   =  6  
  integer, parameter :: jbca   =  7  
  integer, parameter :: jto    =  8  
  integer, parameter :: jks    =  9  
  integer, parameter :: jkr    =  10 
  integer, parameter :: jk1    =  11 
  integer, parameter :: jpk    =  12
  integer, parameter :: jup    =  13 
  integer, parameter :: jrel   =  14 
  integer, parameter :: jleak  =  15


  !!!!!!!!!!!!!!!!!!!!!!
  !!                   
  !!      Y DESCRIPTION
  !!                   
  !!!!!!!!!!!!!!!!!!!!!!!
  ! Y(1:NW)    =   W      = gating variables
  ! Y(NW+1:NY) = ( C, V ) = ( concentrations, potentiel)
  integer, parameter :: tnnp_ny  =  17     ! Number of vaiables 
  integer, parameter :: tnnp_nw  =  12

  integer, parameter :: ym      =  1
  integer, parameter :: yh      =  2  
  integer, parameter :: yj      =  3  
  integer, parameter :: yd      =  4  
  integer, parameter :: yf      =  5  
  integer, parameter :: yfca    =  6  
  integer, parameter :: yr      =  7  
  integer, parameter :: ys      =  8  
  integer, parameter :: yxs     =  9  
  integer, parameter :: yxr1    =  10 
  integer, parameter :: yxr2    =  11 
  integer, parameter :: yg      =  12 
  integer, parameter :: ynai    =  13 
  integer, parameter :: ycai    =  14 
  integer, parameter :: yki     =  15 
  integer, parameter :: ycasr   =  16 
  integer, parameter :: yv      =  tnnp_ny


  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!                        !!
  !!      CONSTANTES TNNP   !!
  !!                        !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 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 :: f     = 96485.3415_rp
  real(RP), parameter :: rtsf  = 26.713760659695648_rp              
  real(RP), parameter :: fsrt  = 1._rp/(rtsf)         

  real(RP), parameter :: vrest    =  -86.424_rp    ! [mV] cellML = -86.2
  real(RP), parameter :: cm       = 0.185_rp       ! [muF]  

  ! ordre de grandeur taille cellule = 7 000 - 45 000 [micro m^3]
  ! "A healthy adult cardiomyocyte has a cylindrical shape that is approximately 100μm 
  ! long and 10-25μm  in diameter" (wikipedia) 
  real(RP), parameter :: vc       = 0.016404_rp    ! [mu L] = [mm^3]  cellML
  real(RP), parameter :: vsr      = 0.001094_rp    ! [mu L] = [mm^3]  cellML
!!$  ! cellML = 1000 fois trop grand pour les volumes  
!!$  real(RP), parameter :: Vc       = 16404E-9_RP   ! [mu L] = [mm^3]  
!!$  real(RP), parameter :: Vsr      =  1094E-9_RP   ! [mu L] = [mm^3]  

  real(RP), parameter :: nao      = 140.0_rp       ! [mM]
  real(RP), parameter :: cao      = 2.0_rp         ! [mM]
  real(RP), parameter :: ko       = 5.4_rp         ! [mM]
  real(RP), parameter :: gna      = 14.838_rp      ! [S/mF] 
  real(RP), parameter :: gk1      = 5.405_rp       ! [S/mF]
  real(RP), parameter :: gto      = 0.294_rp       ! [S/mF] 
  real(RP), parameter :: gkr      = 0.096_rp       ! [S/mF]
  real(RP), parameter :: gks      = 0.245_rp       ! [S/mF]
  real(RP), parameter :: pkna     = 0.03_rp        ! sans dimension
  real(RP), parameter :: gcal     = 1.75e-4_rp     ! [cm^3.muF^{-1}.s^{-1}]
  real(RP), parameter :: knaca    = 1000._rp       ! [A/F]
  real(RP), parameter :: gamma    = 0.35_rp        ! sans dimension 
  real(RP), parameter :: kmca     = 1.38_rp        ! [mM]    
  real(RP), parameter :: kmnai    = 87.5_rp        ! [mM]
  real(RP), parameter :: ksat     = 0.1_rp         ! sans dimension
  real(RP), parameter :: alpha    = 2.5_rp         ! sans dimension
  real(RP), parameter :: pnak     = 1.362_rp       ! [A/F]
  real(RP), parameter :: kmk      = 1._rp          ! [mM]
  real(RP), parameter :: kmna     = 40._rp         ! [mM]
  real(RP), parameter :: gpk      = 0.0146_rp      ! [S/mF]
  real(RP), parameter :: gpca     = 0.825_rp       ! [nS/pF] ??? unite = [A/F]
  real(RP), parameter :: kpca     = 0.0005_rp      ! [mM]
  real(RP), parameter :: gbna     = 2.9e-4_rp      ! [S/mF]
  real(RP), parameter :: gbca     = 5.92e-4_rp     ! [S/mF]
  real(RP), parameter :: vmaxup   = 4.25e-4_rp     ! [mM/ms]
  real(RP), parameter :: kup      = 2.5e-4_rp      ! [mM]
  real(RP), parameter :: arel     = 0.016464_rp    ! [mM]
  real(RP), parameter :: brel     = 0.25_rp        ! [mM]
  real(RP), parameter :: crel     = 0.008232_rp    ! [mM/s]
  real(RP), parameter :: vleak    = 8.e-5_rp       ! [ms{-1}]
  real(RP), parameter :: bufc     = 0.15_rp        ! [mM]
  real(RP), parameter :: kbufc    = 1.e-3_rp       ! [mM]
  real(RP), parameter :: bufsr    = 10._rp         ! [mM]
  real(RP), parameter :: kbufsr   = 0.3_rp         ! [mM]


contains

  subroutine tnnp_y0(Y) 
    real(RP), dimension(TNNP_NY), intent(out) :: Y

    ! potentiel de repos
    y(yv) =  vrest

    ! concentrations en [mM]  (M=mol/L)
    y(ynai)  =  11.212_rp      ! cellML = 11.6
    y(yki)   = 137.614_rp      ! cellML = 138.3
    y(ycai)  =   0.445e-4_rp   ! cellML = 2e-4
    y(ycasr) =   0.513_rp      ! cellML = 0.2      

    ! variables de porte sans dimensions
    y(ym)   = 0.133e-2_rp   ! cellML =0.00  
    y(yh)   = 0.776_rp       ! cellML = 0.75  
    y(yj)   = 0.776_rp       ! cellML =  0.75  
    y(yd)   = 0.193e-4_rp    ! cellML =  0.00  
    y(yf)   = 1.00_rp 
    y(yfca) = 1.00_rp  
    y(yr)   = 0.00_rp  
    y(ys)   = 1.00_rp  
    y(yxs)  = 0.297e-2_rp    ! cellML =  0.00  
    y(yxr1) = 0.00_rp  
    y(yxr2) = 0.485_rp       ! cellML =  1.00
    y(yg)   = 1.00_rp  

  end subroutine tnnp_y0


  subroutine tnnp_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(TNNP_NW) :: WInf, tau

    call tnnp_winf_tau(winf, tau, y) 
    b(1:tnnp_nw)    = (winf - y(1:tnnp_nw)) * (1._rp / tau)

    call  tnnp_g(b(tnnp_nw+1:tnnp_ny), y, i)
    
  end subroutine tnnp_ab_0


  subroutine tnnp_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 

    real(RP), dimension(TNNP_NW) :: WInf, tau

    call tnnp_winf_tau(winf, tau, y) 

    tau     = 1._rp/tau
    a(1:tnnp_nw) = - tau

    b(1:tnnp_nw) = winf * tau

    call  tnnp_g(b(tnnp_nw+1:tnnp_ny), y, i)
    
  end subroutine tnnp_ab_w


  !!   IList
  !!                        
  subroutine tnnp_ilist(IList, Y) 
    real(RP), dimension(TNNP_nI), intent(out) :: IList
    real(RP), dimension(TNNP_nY), intent(in)  :: Y

    real(RP) :: a, b, ENa, ECa, EK, EKs, U, V, VmE
    real(RP) :: Nai, Cai, Ki, Cas
    real(RP) :: a_k1, b_k1, xK1_inf

    ! concentrations
    nai = y(ynai)
    cai = y(ycai)
    ki  = y(yki)
    cas = y(ycasr)

    ! Potentiels de Nernst
    ena = rtsf * log( nao / nai )
    eca = rtsf * log( cao / cai ) * 0.5_rp
    ek  = rtsf * log( ko / ki )
    eks =       ( ko + pkna * nao ) 
    eks = eks / ( ki + pkna * nai ) 
    eks = rtsf* log( eks ) 

    v = y(yv)
    u = v * fsrt

    !-------------!
    !   SODIUM    !
    !-------------!
    vme = v - ena

    !  fast sodium
    ilist(jna) = gna * y(ym)**3 * y(yh) * y(yj) * vme

    !  Echangeur Na+/Ca2+
    a = exp(gamma * u ) * nai**3 * cao
    a = a - exp( (gamma-1._rp)* u ) * nao**3 * cai * alpha
    b = ( kmnai**3 + nao**3 )*( kmca + cao )
    b = b * ( 1._rp+ ksat * exp( (gamma-1._rp)* u ) )
    ilist(jnaca) = knaca * a / b

    !  Echangeur Na+/K+ 
    a = ( ko + kmk ) * (nai + kmna)
    a = a * ( 1._rp + 0.1245_rp * exp( -0.1_rp* u ) + 0.0353_rp*exp(- u ) )
    ilist(jnak) = pnak * ko * nai / a

    !  background Na
    ilist(jbna) = gbna * vme 


    !-------------!
    !  CALCIUM    !
    !-------------!
    
    !  Ca canal L
    a = gcal * y(yd) * y(yf) * y(yfca)
    a = a * 4._rp * u * f 
    a = a * ( cai*exp( 2._rp* u ) - 0.341_rp* cao )
    ilist(jcal) = a / ( exp( 2._rp* u ) - 1._rp )

    !  pompe  Ca2+
    ilist(jpca)  =  gpca * cai / (kpca + cai)

    !  background Ca
    ilist(jbca)  =  gbca * (v - eca) 


    !-------------!
    !  POTASSIUM  !
    !-------------!
    
    vme = v - ek

    !  Transient outward current
    ilist(jto) = gto * y(yr) * y(ys) * vme
 
    !  Slow delayed rectifier current
    ilist(jks) = gks * y(yxs)**2 * ( v - eks)

    !  Rapid Delayed Rectifier Current
    !!$    a = GKr * sqrt( Ko / 5.4 ) * VmE  !! modifie car Ko = 5.4
    a = gkr *  vme
    ilist(jkr)  = a * y(yxr1) * y(yxr2)

    !  Inward Rectifier K+ Current
    ! calcul xkl_inf 
    a_k1    = 0.1_rp / ( 1._rp + exp( 0.06_rp * ( vme -200._rp ) ) )
    b_k1    = 3._rp * exp( 2e-4_rp * ( vme + 100._rp ) ) + exp( 0.1_rp * ( vme - 10._rp ) )
    b_k1    = b_k1 / ( 1._rp + exp( -0.5_rp * ( vme ) ) )
    xk1_inf = a_k1 / ( a_k1 + b_k1 )
    !!$    IList(jK1) = GK1 * sqrt(Ko/5.4) * xk1_inf * VmE  !! modifie car Ko = 5.4
    ilist(jk1) = gk1 * xk1_inf * vme

    !  pompe  K+  
    ilist(jpk) = gpk * vme / ( 1._rp + exp( ( 25._rp - v ) / 5.98_rp ) )


    !---------------------!
    !  COURANTS INTERNES  !
    !---------------------!

    ! Iup   en [mM/ms]
    ilist(jup) = vmaxup / ( 1._rp + (kup / cai)**2 )

    ! Irel  en [mM/ms] 
    a = arel * cas**2 / ( brel**2 + cas**2 ) + crel
    ilist(jrel) = a * y(yd) * y(yg)

    ! Ileak en [mM/ms]
    ilist(jleak) = vleak * (cas - cai)

  end subroutine tnnp_ilist


  !!   Winf_tau             
  !!                        
  subroutine tnnp_winf_tau(Winf, tau, Y) 
    real(RP), dimension(TNNP_NY), intent(in)  :: Y
    real(RP), dimension(TNNP_NW), intent(out) :: Winf, tau

    real(RP) :: V, Cai, a, b, c

    v   = y(yv)
    cai = y(ycai)

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TAU

    !  tau_m
    a = 1._rp   / ( 1._rp + exp( (-60._rp    - v  ) * 0.2_rp   ) )
    b = 0.1_rp / ( 1._rp + exp( ( v     + 35._rp ) * 0.2_rp   ) ) 
    b = 0.1_rp / ( 1._rp + exp( ( v     - 50._rp ) * 0.005_rp ) ) + b
    tau(ym) = a * b

    !  tau_h
    if ( v >= -40._rp ) then
       a = 0._rp
       b = 0.77_rp / ( 0.13_rp  * ( 1._rp + exp( -( v + 10.66_rp ) / 11.1_rp) ) )
    else
       a = 0.057_rp * exp( -( v + 80._rp ) / 6.8_rp )
       b = 2.7_rp   * exp( 0.079_rp*v ) + 3.1e5_rp * exp( 0.3485_rp*v )
    end if
    tau(yh) = 1._rp / (a + b)

    !  tau_j
    if ( v >= -40._rp ) then
       a = 0._rp
       b = 0.6_rp * exp(0.057_rp*v) / (1._rp + exp( -0.1_rp*(v+32._rp) ) )
    else   
       a = - 25428._rp * exp(0.2444_rp*v) - 6.948e-6_rp*exp(-0.04391_rp*v) 
       a = a * (v+37.78_rp) / ( 1._rp + exp( 0.311_rp*(v+79.23_rp) ) ) 
       b = 0.02424_rp * exp(-0.01052_rp*v) / ( 1._rp + exp( -0.1378_rp * ( v + 40.14_rp ) ) )
    end if
    tau(yj) = 1._rp / (a + b)

    !  tau_d
    a = 1.4_rp / ( 1._rp + exp( (-35._rp - v) / 13._rp  ) ) + 0.25_rp
    b = 1.4_rp / ( 1._rp + exp( ( v + 5._rp ) * 0.20_rp ) )
    c = 1.0_rp / ( 1._rp + exp( (50._rp - v ) * 0.05_rp ) )
    tau(yd) = a * b + c

    !  tau_f
    a = 165._rp / ( 1._rp + exp( (25._rp - v) * 0.1_rp  ) ) 
    tau(yf) = a + 1125._rp * exp( -( v + 27._rp)**2 / 240._rp ) +  80._rp

    !  tau_fCa
    tau(yfca) = 2._rp

    !  tau_r
    tau(yr) = 9.5_rp * exp( -( v + 40._rp )**2 / 1800._rp ) + 0.8_rp

    !  tau_s
    a = 85._rp * exp( - ( v + 45._rp )**2 / 320._rp ) 
    tau(ys) = a + 5._rp / ( 1._rp + exp( ( v - 20._rp ) * 0.2_rp)) + 3._rp

    !  tau_xs
    a = 1100._rp / sqrt( 1._rp+ exp( (-10._rp - v ) / 6._rp ) )
    b  = 1._rp   / ( 1._rp+ exp( ( v - 60._rp ) * 0.05_rp ) )
    tau(yxs)  = a * b

    ! tau_xr1
    a = 450._rp / ( 1._rp + exp( ( -45._rp - v ) * 0.1_rp ) )
    b = 6._rp   / ( 1._rp + exp( ( v   + 30._rp) / 11.5_rp) )
    tau(yxr1) = a * b

    !  tau_xr2
    a = 3._rp    / ( 1._rp + exp( ( -60._rp - v  ) * 0.05_rp ) )
    b = 1.12_rp / ( 1._rp + exp( ( v   - 60._rp ) * 0.05_rp ) )
    tau(yxr2)  = a * b

    ! tau_g
    tau(yg)   = 2._rp


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WInf

    !  m_inf
    winf(ym) = 1._rp  / ( 1._rp + exp( (-56.86_rp - v  ) / 9.03_rp  ) )**2

    !  h_inf
    winf(yh) = 1._rp / ( 1._rp + exp( ( v + 71.55_rp ) / 7.43_rp ) )**2
    
    !  j_inf = h_inf
    winf(yj) = winf(yh) 

    !  d_inf
     winf(yd) = 1._rp   / ( 1._rp + exp( (-5._rp - v ) / 7.5_rp ) )

    !  f_inf
     winf(yf) = 1._rp   / ( 1._rp + exp( ( v + 20._rp) / 7._rp   ) )

    !  fcs_inf
    a = 1._rp  / ( 1._rp +      ( cai / 0.000325_rp )**8       )
    b = 0.1_rp / ( 1._rp + exp( ( cai - 0.0005_rp   ) * 1e4_rp  ) )
    c = 0.2_rp / ( 1._rp + exp( ( cai - 0.00075_rp  ) * 1250._rp ) )
    winf(yfca)   = ( a + b + c + 0.23_rp ) / 1.46_rp

    !  r_inf
    winf(yr) = 1._rp/ ( 1._rp + exp( ( 20._rp- v ) / 6._rp ) )

    !  s_inf
    winf(ys) = 1._rp / ( 1._rp + exp( ( v + 20._rp ) * 0.2_rp ) )

    !  xs_inf
    winf(yxs) = 1._rp    / ( 1._rp + exp( ( -5._rp - v ) / 14._rp ) )

    
    !  xr1_inf
    winf(yxr1) = 1._rp   / ( 1._rp + exp( ( -26._rp - v ) / 7._rp   ) )

    !  xr2_inf
    winf(yxr2)  = 1._rp / ( 1._rp + exp( ( v   + 88._rp ) / 24._rp   ) )

    !  g_inf
    if ( cai < 3.5e-4_rp  ) then
       winf(yg) = 1._rp / ( 1._rp + ( cai/3.5e-4_rp )**6 )
    else
       winf(yg) = 1._rp / ( 1._rp + ( cai/3.5e-4_rp )**16 )
    end if


    !!!!!!!!!!!!!!!!!!!!!! APPENDING TAU TO +INFINITY


    ! tau_fCa
    if ( ( winf(yfca) > y(yfca) ) .and. ( v > -60._rp) )  then
       tau(yfca) = 1e127_rp

    end if

    ! tau_g
    if ( ( winf(yg) > y(yg) ) .and. ( v > -60._rp) )  then
       tau(yg) = 1e127_rp
    end if

  end subroutine tnnp_winf_tau


  subroutine tnnp_g(G, Y, I) 
    real(RP), dimension(TNNP_NY)          , intent(in)  :: Y
    real(RP)                              , intent(in)  :: I
    real(RP), dimension(TNNP_NW+1:TNNP_NY), intent(out) :: G

    real(RP), dimension(TNNP_NI)  :: IList
    real(RP) :: INa, IK, ICa, D1, D2, J

    ! calcul des courants ioniques
    call tnnp_ilist(ilist, y)

    ! courant sodium total
    ina = ilist(jna) + ilist(jbna)  + 3._rp*( ilist(jnak) + ilist(jnaca) )

    ! Nai
    g(ynai) = -ina * cm / (vc*f)

    ! courant potassium total
    ik = sum( ilist(jto:jpk) ) - 2._rp* ilist(jnak) 

    ! Ki
    g(yki) = -( ik + i ) * cm / (vc*f)

    ! courant calcium total 
    ica = sum( ilist(jcal:jbca) ) - 2._rp* ilist(jnaca) 

    ! flux intracellulaire
    j  = ilist(jleak) - ilist(jup) + ilist(jrel)

    ! Cai
    ! calcul de D1 = d Cai-bufC / d Cai 
    d1 = bufc * kbufc / (y(ycai) + kbufc)**2   ! [sans dimension]
    ! calcul de D2 = d Cai_total / dt
    d2 = j - ica / (2._rp*vc*f) * cm
    g(ycai) = d2 / ( 1._rp + d1 )

    ! Casr
    ! calcul de D1 = d CaSR_bufrs / d CaSR 
    d1 = bufsr * kbufsr / (y(ycasr) + kbufsr)**2    ! [sans dimension]
    ! calcul de D2 = d CaSR_total / dt
    d2 = - j * (vc/vsr)
    g(ycasr) = d2 / ( 1._rp + d1 )

    ! vm
    g(yv) =  - (ina + ik + ica) + i

  end subroutine tnnp_g

end module tnnp_mod