electrophysiology_cell_routines.f90

Go to the documentation of this file.

module electrophysiology_cell_routines
  USE base_routines
  use field_routines
  use kinds
  use strings
  use types

#include "macros.h"  

  implicit none
  private

  !interfaces

  public  bueno_orovio_initialize, bueno_orovio_integrate
  public  tentusscher06_initialize, tentusscher06_integrate

  interface pow
    module procedure powr
    module procedure powi
  end  interface pow

contains
  ! for auto generated code

  real(dp) function powr(a,b)
    real(dp), intent(in) :: a,b
    powr = a**b
  end function powr
  real(dp) function powi(a,b)
    real(dp), intent(in) :: a
    integer(intg), intent(in) :: b
    powi = a**b
  end function powi

  ! Bueno-Orovio (2008) minimal cell model. Membrane potential rescaled, partial evaluation applied.

  ! initialize field
  subroutine bueno_orovio_initialize(field,err,error,*)
    type(field_type), intent(inout), pointer :: field
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    
    integer(intg) :: i
    REAL(DP), PARAMETER, DIMENSION(1:4) :: y0 = (/ -8.09242e+01,  9.99993e-01, 2.16366e-02, 9.84320e-01 /) ! paced initial condition

    enters('bueno_orovio_init',err,error,*999)
    DO i=1,4
      CALL field_component_values_initialise(field,field_u_variable_type,field_values_set_type,i,y0(i),err,error,*999)
    END DO

    exits('bueno_orovio_init')
    return
999 errorsexits('bueno_orovio_init',err,error)
    return 1

  end subroutine bueno_orovio_initialize

  subroutine bueno_orovio_dydt(t,y,dydt,activ)
    real(dp), dimension(:), intent(out) :: dydt(:)
    real(dp), dimension(:), intent(in) :: y(:)
    real(dp), intent(in) :: t
    real(dp), intent(in) :: activ

    real(dp),parameter :: u_m=3.0e-01, u_p=1.3e-01, u_q=6.0e-03, u_r=6.0e-03, vscale=81.0, period = 1000.0
    real(dp) :: J_fi, J_si, J_so, Vm, bv, m, p, q, r, tau_o, tau_s, tau_so, tau_v_minus, tau_w_minus, w_inf, J_stim
  
    if (mod(t,period)>=0 .and. mod(t,period)<=0.999999) then ! apply stimulus in first ms
      j_stim = -0.65 * activ ! stimulus: ~52 mV for 1 ms
    else
     j_stim  = 0
    end if

    bv = (1+(y(1)/vscale)); 
    if (bv < u_m) then
      m = 0
    else
      m = 1
    end if
    if (bv < u_q) then
      q = 0
    else
      q = 1
    end if
    if (bv < u_p) then
      p = 0
    else
      p = 1
    end if

    if (bv < u_r) then
      r = 0
    else
      r = 1
    end if

    j_fi = (9.09090909090909e+00*(-m)*y(2)*(-3.0e-01+bv)*(1.55e+00-bv))
    tau_v_minus = ((1150*q)+(60*(1-q)))
    dydt(2) = (((1-m)*(0-y(2))/tau_v_minus)-(6.89655172413793e-01*m*y(2)))
    tau_o = ((400*(1-r))+(6*r))
    tau_so = (3.002e+01+(-1.4512e+01*(1+tanh((2.046e+00*(-6.5e-01+bv))))))
    j_so = ((bv*(1-p)/tau_o)+(p/tau_so))
    j_si = (5.29801324503311e-01*(-p)*y(4)*y(3))
    dydt(1) = ((-j_fi-j_so-j_si-j_stim)*vscale)
    vm = (-83+(8.57e+01*bv))
    tau_s = ((2.7342e+00*(1-p))+(16*p))
    dydt(3) = (((5.0e-01*(1+tanh((2.0994e+00*(-9.087e-01+bv)))))-y(3))/tau_s)
    w_inf = (((1-r)*(1-(1.42857142857143e+01*bv)))+(9.4e-01*r))
    tau_w_minus = (60+(-2.25e+01*(1+tanh((65*(-3.0e-02+bv))))))
    dydt(4) = (((1-r)*(w_inf-y(4))/tau_w_minus)-(5.0e-03*r*y(4)))
  end subroutine bueno_orovio_dydt   

  ! integrates all cell models
  subroutine bueno_orovio_integrate(cells,materials,t0,t1,err,error,*)
    type(field_type), intent(inout), pointer :: cells, materials
    real(dp), intent(in)    :: t0, t1
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    integer, parameter :: celldim = 4
    real(dp), dimension(1:celldim) :: y, dydt
    real(dp) :: t, dt, activ
    real(dp), dimension(:), pointer :: celldata, activdata

    integer(intg) :: ncells, i, d
    type(domain_ptr_type), pointer :: domain
    TYPE(field_variable_type), POINTER :: CELLS_VARIABLE, ACTIV_VARIABLE

    enters('bueno_orovio_integrate',err,error,*999)

    domain=>cells%decomposition%domain(1)
    ncells = domain%ptr%topology%nodes%number_of_nodes ! local

    cells_variable=>cells%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
    activ_variable=>materials%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR

    CALL field_parameter_set_data_get(cells,field_u_variable_type,field_values_set_type,celldata,err,error,*999)
    CALL field_parameter_set_data_get(materials,field_u_variable_type,field_values_set_type,activdata,err,error,*999)

    do i=1,ncells
      !   field ->   y
      do d=1,celldim
        !Default to version 1 of each node derivative
        y(d) = celldata(cells_variable%COMPONENTS(d)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1))
      end do
      !Default to version 1 of each node derivative
      activ = activdata(activ_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1))
      
      t = t0
      do while (t < t1 - 1e-6)
        ! integrate one cell, one time step. 
        ! just use adaptive forward euler, tested in matlab
        call bueno_orovio_dydt(t,y,dydt,activ)
        dt = min(t1-t,0.5)  ! default max
        dt = min(dt,1 / abs(dydt(1))) ! maximum increase in Vm
!        IF(DT < 1e-6)  WRITE(*,*) 'CELL ',i, 'DVDT = ',dydt(1),'dt = ',dt,' activ=',activ,'y=',y
        y = y + dt * dydt ! FWE
        t = t + dt
      end do
      !   y -> field  
      do d=1,celldim
 !       call field_parameter_set_update_local_node(cells,field_u_variable_type,field_values_set_type,1,i,d,y(d), err,error,*999)
        !Default to version 1 of each node derivative
        celldata(cells_variable%COMPONENTS(d)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1)) = y(d)
      end do
    end do
    CALL field_parameter_set_data_restore(cells,field_u_variable_type,field_values_set_type,celldata,err,error,*999)
    CALL field_parameter_set_data_restore(materials,field_u_variable_type,field_values_set_type,activdata,err,error,*999)


    exits('bueno_orovio_integrate')
    return
999 errorsexits('bueno_orovio_integrate',err,error)
    return 1
  end subroutine bueno_orovio_integrate
  


  ! Ten Tusscher and Panfilov (2006) human ventricular cell model

  ! initialize field
  subroutine tentusscher06_initialize(field,err,error,*)
    type(field_type), intent(inout), pointer :: field
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    
    integer(intg) :: i
    REAL(DP), PARAMETER, DIMENSION(1:19) :: y0 = (/ -85.23, 0.9755, 0.9953, 0.7888, 3.64, 0.000126, 0.00036, 0.9073, 0.7444,&
        & 0.7045, 0.00172,  3.373e-5, 136.89, 0.00621, 0.4712, 0.0095, 8.604, 2.42e-8, 0.999998 /) ! paced initial condition

    enters('tentusscher06_init',err,error,*999)
    DO i=1,19
      CALL field_component_values_initialise(field,field_u_variable_type,field_values_set_type,i,y0(i),err,error,*999)
    END DO

    exits('tentusscher06_init')
    return
999 errorsexits('tentusscher06_init',err,error)
    return 1

  end subroutine tentusscher06_initialize

  subroutine tentusscher06_dydt(t,y,dydt,activ)
    real(dp), dimension(:), intent(out) :: dydt(:)
    real(dp), dimension(:), intent(in) :: y(:)
    real(dp), intent(in) :: t
    real(dp), intent(in) :: activ

    real(dp),parameter :: period=1000
    real(dp) :: Ca_i_bufc, Ca_sr_bufsr, Ca_ss_bufss, E_Ca, E_K, E_Ks, E_Na, O, alpha_K1, alpha_d, alpha_m, alpha_xr1, alpha_xr2,&
    & alpha_xs, beta_K1, beta_d, beta_m, beta_xr1, beta_xr2, beta_xs, d_inf, f2_inf, fCass_inf, f_inf, gamma_d, h_inf, i_CaL,&
    & i_K1, i_Kr, i_Ks, i_Na, i_NaCa, i_NaK, i_Stim, i_b_Ca, i_b_Na, i_leak, i_p_Ca, i_p_K, i_rel, i_to, i_up, i_xfer, j_inf, k1,&
    & k2, kcasr, m_inf, r_inf, s_inf, tau_d, tau_f, tau_f2, tau_fCass, tau_h, tau_j, tau_m, tau_r, tau_s, tau_xr1, tau_xr2,tau_xs,&
    & xK1_inf, xr1_inf, xr2_inf, xs_inf

    if (mod(t,period)>=0 .and. mod(t,period)<=1.999999) then ! apply stimulus in 2 ms
      i_stim = -3.57142857142857e+01 * activ
    else
      i_stim = 0
    end if

    i_cal = (5.75002020961245e-01*y(12)*y(4)*y(2)*y(3)*(-15+y(1))*(-2+(2.5e-01*y(7)*exp((7.48677816454909e-02*(-15+y(1))))))/&
          & (-1+exp((7.48677816454909e-02*(-15+y(1))))))
    
    d_inf = (1/(1+exp((1.33333333333333e-01*(-8-y(1))))))
    alpha_d = (2.5e-01+(1.4e+00/(1+exp((7.69230769230769e-02*(-35-y(1)))))))
    beta_d = (1.4e+00/(1+exp((2.0e-01*(5+y(1))))))
    gamma_d = (1/(1+exp((5.0e-02*(50-y(1))))))
    tau_d = ((1*alpha_d*beta_d)+gamma_d)
    f2_inf = (3.3e-01+(6.7e-01/(1+exp((1.42857142857143e-01*(35+y(1)))))))
    tau_f2 = ((562*exp((4.16666666666667e-03*(-pow((27+y(1)),2)))))+(31/(1+exp((1.0e-01*(25-y(1))))))+&
           & (80/(1+exp((1.0e-01*(30+y(1)))))))
    dydt(2) = ((f2_inf-y(2))/tau_f2)
    fcass_inf = (4.0e-01+(6.0e-01/(1+pow((20*y(7)),2))))
    tau_fcass = (2+(80/(1+pow((20*y(7)),2))))
    dydt(3) = ((fcass_inf-y(3))/tau_fcass)
    f_inf = (1/(1+exp((1.42857142857143e-01*(20+y(1))))))
    tau_f = (20+(1.1025e+03*exp((4.44444444444444e-03*(-pow((27+y(1)),2)))))+(200/(1+exp((1.0e-01*(13-y(1))))))+&
          &(180/(1+exp((1.0e-01*(30+y(1)))))))
    dydt(4) = ((f_inf-y(4))/tau_f)
    e_ca = (1.33568803298478e+01*log((2/y(6))))
    i_b_ca = (5.92e-04*(y(1)-e_ca))
    kcasr = (2.5e+00-(1.5e+00/(1+pow((1.5e+00/y(5)),2))))
    k1 = (1.5e-01/kcasr)
    o = (k1*pow(y(7),2)*y(8)/(6.0e-02+(k1*pow(y(7),2))))
    i_rel = (1.02e-01*o*(y(5)-y(7)))
    i_up = (6.375e-03/(1+(6.25e-08/pow(y(6),2))))
    i_leak = (3.6e-04*(y(5)-y(6)))
    i_xfer = (3.8e-03*(y(7)-y(6)))
    k2 = (4.5e-02*kcasr)
    dydt(8) = (((-k2)*y(7)*y(8))+(5.0e-03*(1-y(8))))
    ca_i_bufc = (1/(1+(2.0e-04/pow((1.0e-03+y(6)),2))))
    ca_sr_bufsr = (1/(1+(3/pow((3.0e-01+y(5)),2))))
    ca_ss_bufss = (1/(1+(1.0e-04/pow((2.5e-04+y(7)),2))))
    i_p_ca = (1.238e-01*y(6)/(5.0e-04+y(6)))
    i_naca = (8.66622022994244e-05*((2*exp((1.31018617879609e-02*y(1)))*pow(y(17),3))-(6860000*&
           &exp((-2.43320290347846e-02*y(1)))*y(6)))/(1+(1.0e-01*exp((-2.43320290347846e-02*y(1))))))
    dydt(6) = (ca_i_bufc*((6.66910509631797e-02*(i_leak-i_up))+i_xfer-(5.84427487220097e-05*(i_b_ca+i_p_ca-(2*i_naca)))))
    dydt(5) = (ca_sr_bufsr*(i_up-i_rel-i_leak))
    dydt(7) = (ca_ss_bufss*((-1.75328246166029e-02*i_cal)+(2.00073152889539e+01*i_rel)-(300*i_xfer)))
    e_na = (2.67137606596956e+01*log((140/y(17))))
    i_na = (1.4838e+01*pow(y(11),3)*y(9)*y(10)*(y(1)-e_na))
    h_inf = (1/pow((1+exp((1.34589502018843e-01*(7.155e+01+y(1))))),2))
    if (y(1) < -40.0) then
      tau_h = (1/((5.7e-02*exp((1.47058823529412e-01*(-80-y(1)))))+(2.7e+00*exp((7.9e-02*y(1))))+(310000*exp((3.485e-01*y(1))))))
    else
      tau_h = (1.68831168831169e-01*(1+exp((-9.00900900900901e-02*(1.066e+01+y(1))))))
    end if
    dydt(9) = ((h_inf-y(9))/tau_h)
    j_inf = (1/pow((1+exp((1.34589502018843e-01*(7.155e+01+y(1))))),2))
    if (y(1) < -40.0) then
     tau_j = (1/((1*((-25428*exp((2.444e-01*y(1))))-(6.948e-06*exp((-4.391e-02*y(1)))))*(3.778e+01+y(1))/&
          &(1+exp((3.11e-01*(7.923e+01+y(1))))))+(2.424e-02*exp((-1.052e-02*y(1)))/(1+exp((-1.378e-01*(4.014e+01+y(1))))))))
    else
     tau_j = (1.66666666666667e+00/exp((5.7e-02*y(1)))*(1+exp((-1.0e-01*(32+y(1))))))
    end if
    dydt(10) = ((j_inf-y(10))/tau_j)
    m_inf = (1/pow((1+exp((1.10741971207087e-01*(-5.686e+01-y(1))))),2))
    alpha_m = (1/(1+exp((2.0e-01*(-60-y(1))))))
    beta_m = ((1.0e-01/(1+exp((2.0e-01*(35+y(1))))))+(1.0e-01/(1+exp((5.0e-03*(-50+y(1)))))))
    tau_m = (1*alpha_m*beta_m)
    dydt(11) = ((m_inf-y(11))/tau_m)
    e_k = (2.67137606596956e+01*log((5.4e+00/y(13))))
    alpha_k1 = (1.0e-01/(1+exp((6.0e-02*(-200+y(1)-e_k)))))
    beta_k1 = (((3*exp((2.0e-04*(100+y(1)-e_k))))+exp((1.0e-01*(-10+y(1)-e_k))))/(1+exp((-5.0e-01*(y(1)-e_k)))))
    xk1_inf = (alpha_k1/(alpha_k1+beta_k1))
    i_k1 = (5.405e+00*xk1_inf*(y(1)-e_k))
    i_to = (2.94e-01*y(18)*y(19)*(y(1)-e_k))
    i_kr = (1.53e-01*y(14)*y(15)*(y(1)-e_k))
    e_ks = (2.67137606596956e+01*log((9.6e+00/(y(13)+(3.0e-02*y(17))))))
    i_ks = (3.92e-01*pow(y(16),2)*(y(1)-e_ks))
    i_nak = (2.298375e+00*y(17)/(40+y(17))/(1+(1.245e-01*exp((-3.74338908227455e-03*y(1))))+&
          &(3.53e-02*exp((3.74338908227455e-02*(-y(1)))))))
    i_b_na = (2.9e-04*(y(1)-e_na))
    i_p_k = (1.46e-02*(y(1)-e_k)/(1+exp((1.67224080267559e-01*(25-y(1))))))
    dydt(12) = ((d_inf-y(12))/tau_d)
    dydt(1) = (-1*(i_k1+i_to+i_kr+i_ks+i_cal+i_nak+i_na+i_b_na+i_naca+i_b_ca+i_p_k+i_p_ca+i_stim))
    dydt(13) = (-1.16885497444019e-04*(i_k1+i_to+i_kr+i_ks+i_p_k+i_stim-(2*i_nak)))
    xr1_inf = (1/(1+exp((1.42857142857143e-01*(-26-y(1))))))
    alpha_xr1 = (450/(1+exp((1.0e-01*(-45-y(1))))))
    beta_xr1 = (6/(1+exp((8.69565217391304e-02*(30+y(1))))))
    tau_xr1 = (1*alpha_xr1*beta_xr1)
    dydt(14) = ((xr1_inf-y(14))/tau_xr1)
    xr2_inf = (1/(1+exp((4.16666666666667e-02*(88+y(1))))))
    alpha_xr2 = (3/(1+exp((5.0e-02*(-60-y(1))))))
    beta_xr2 = (1.12e+00/(1+exp((5.0e-02*(-60+y(1))))))
    tau_xr2 = (1*alpha_xr2*beta_xr2)
    dydt(15) = ((xr2_inf-y(15))/tau_xr2)
    xs_inf = (1/(1+exp((7.14285714285714e-02*(-5-y(1))))))
    alpha_xs = (1400/sqrt((1+exp((1.66666666666667e-01*(5-y(1)))))))
    beta_xs = (1/(1+exp((6.66666666666667e-02*(-35+y(1))))))
    tau_xs = (80+(1*alpha_xs*beta_xs))
    dydt(16) = ((xs_inf-y(16))/tau_xs)
    dydt(17) = (-1.16885497444019e-04*(i_na+i_b_na+(3*i_nak)+(3*i_naca)))
    r_inf = (1/(1+exp((1.66666666666667e-01*(20-y(1))))))
    tau_r = (8.0e-01+(9.5e+00*exp((5.55555555555556e-04*(-pow((40+y(1)),2))))))
    dydt(18) = ((r_inf-y(18))/tau_r)
    s_inf = (1/(1+exp((2.0e-01*(20+y(1))))))
    tau_s = (3+(85*exp((3.125e-03*(-pow((45+y(1)),2)))))+(5/(1+exp((2.0e-01*(-20+y(1)))))))
    dydt(19) = ((s_inf-y(19))/tau_s)
  end subroutine tentusscher06_dydt   

  ! integrates all cell models
  subroutine tentusscher06_integrate(cells,materials,t0,t1,err,error,*)
    type(field_type), intent(inout), pointer :: cells, materials
    real(dp), intent(in)    :: t0, t1
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    integer, parameter :: celldim = 19
    real(dp), dimension(1:celldim) :: dydt
    real(dp), dimension(:), pointer :: y
    real(dp) :: t, dt, activ, m_inf, d_inf, m_inf0, d_inf0, prev_v
    real(dp), dimension(:), pointer :: celldata, activdata

    integer(intg) :: ncells, i, d
    type(domain_ptr_type), pointer :: domain
    TYPE(field_variable_type), POINTER :: CELLS_VARIABLE, ACTIV_VARIABLE

    enters('tentusscher06_integrate',err,error,*999)

    domain=>cells%decomposition%domain(1)
    ncells = domain%ptr%topology%nodes%number_of_nodes ! local

    cells_variable=>cells%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR
    activ_variable=>materials%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR

    CALL field_parameter_set_data_get(cells,field_u_variable_type,field_values_set_type,celldata,err,error,*999)
    CALL field_parameter_set_data_get(materials,field_u_variable_type,field_values_set_type,activdata,err,error,*999)

    do i=1,ncells
      !Default to version 1 of each node derivative
      d = cells_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1)
      y => celldata(d:d+celldim-1)
      prev_v = y(1)
      !Default to version 1 of each node derivative
      activ = activdata(activ_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1))
      
  
      t = t0
      do while (t < t1 - 1e-6)
        ! integrate one cell, one time step. 
        ! just use adaptive forward euler, tested in matlab
        call tentusscher06_dydt(t,y,dydt,activ)
        dt = min(t1-t,1.0/3)  ! default max
        dt = min(dt,1 / abs(dydt(1))) ! maximum increase in Vm

        m_inf0 = (1/pow((1+exp((1.10741971207087e-01*(-5.686e+01-y(1))))),2))
        d_inf0 = (1/(1+exp((1.33333333333333e-01*(-8-y(1))))))
        
        y = y + dt * dydt ! FWE
        
        m_inf = (1/pow((1+exp((1.10741971207087e-01*(-5.686e+01-y(1))))),2))
        d_inf = (1/(1+exp((1.33333333333333e-01*(-8-y(1))))))
        
        ! FWEoo on m,d gates : prevent gates from crossing midpoint steady state
        
        m_inf = (m_inf0 + m_inf)/2
        d_inf = (d_inf0 + d_inf)/2
        if(dydt(11) * (y(11) - m_inf) > 0) then
          y(11) = m_inf
        end if
        if(dydt(12) * (y(12) - d_inf) > 0) then
          y(12) = d_inf
        end if        
        
        t = t + dt
      end do
      if(prev_v < 0 .and. y(1) > 0) then ! store activation times, where else?
        !Default to version 1 of each node derivative
        activdata(activ_variable%COMPONENTS(7)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(i)%DERIVATIVES(1)%VERSIONS(1))=t1
      end if
    end do ! for cells
    CALL field_parameter_set_data_restore(cells,field_u_variable_type,field_values_set_type,celldata,err,error,*999)
    CALL field_parameter_set_data_restore(materials,field_u_variable_type,field_values_set_type,activdata,err,error,*999)

    exits('tentusscher06_integrate')
    return
999 errorsexits('tentusscher06_integrate',err,error)
    return 1
  end subroutine tentusscher06_integrate
end module electrophysiology_cell_routines