subcellular.f95

!        -*- f95 -*-
! (c) 2016 - author information redacted by request 2026-08-07
! INRIA Rhone-Alpes
! STDP model : This module defines models of mGluR-dependent signaling, CICR, calcium buffering, production of endocannabinoids,
! phenomenological model of endocannabinoid-dependent plasticity
module subcellular

    use pars_mod
    use general_math

    implicit none

    logical, save :: subcellular_compute_once = .true.

    contains

    real*8 function tauCa(C, pars)
        implicit none
        real*8 :: C
        type(pars_type) :: pars
        tauCa = (1+pars%CaBuff%BT/pars%CaBuff%KdB/(1+C/pars%CaBuff%KdB)**2)
    end function tauCa

    real*8 function JIP3R_CICR_func(IP3, Ca_cyt, Ca_ER, h, pars)
        implicit none
        real*8 :: IP3, Ca_cyt, Ca_ER, h
        type(pars_type) :: pars
        real*8 :: minf, ninf
        minf = IP3/(IP3+pars%CICR%d1)
        ninf = Ca_cyt/(Ca_cyt+pars%CICR%d5)
        JIP3R_CICR_func = pars%CICR%rc * (minf*ninf*h)**3 * (Ca_ER-Ca_cyt)
    end function JIP3R_CICR_func

    real*8 function Jserca_CICR_func(Ca_cyt, pars)
        implicit none
        real*8 :: Ca_cyt
        type(pars_type) :: pars
        Jserca_CICR_func = pars%CICR%ver * Hill(Ca_cyt,pars%CICR%ker,2)
    end function Jserca_CICR_func

    real*8 function Jleak_CICR_func(Ca_cyt, Ca_ER, pars)
        implicit none
        real*8 :: Ca_cyt, Ca_ER
        type(pars_type) :: pars
        Jleak_CICR_func = pars%CICR%rl * (Ca_ER-Ca_cyt)
    end function Jleak_CICR_func

    real*8 function dh_CICR(Ca_cyt, IP3, h, pars)
        implicit none
        real*8 :: Ca_cyt, IP3, h
        type(pars_type) :: pars
        type(CICR_type) :: p
        p = pars%CICR
        dh_CICR = (p%a2 * p%d2 * (IP3+p%d1)/(IP3+p%d3))*(1-h) - p%a2*Ca_cyt*h
    end function dh_CICR

    real*8 function dCa_cyt_func(JCa, Ca_cyt, pars)
        implicit none
        real*8 :: JCa, Ca_cyt
        type(pars_type) :: pars
        dCa_cyt_func = ( JCa - ( Ca_cyt - pars%CaBuff%Cab )/pars%CaBuff%tauCab )/tauCa(Ca_cyt, pars)
    end function dCa_cyt_func

    real*8 function dCa_ER_func(JCaER, Ca_ER, pars)
        implicit none
        real*8 :: JCaER, Ca_ER
        ! JCaER = JIP3R-Jserca+Jleak
        type(pars_type) :: pars
        dCa_ER_func = -JCaER*pars%CICR%rhoER/tauCa(Ca_ER, pars)
    end function dCa_ER_func

    real*8 function vglu_func(Glu,Ca_cyt,pars)
        implicit none
        real*8 :: Glu,Ca_cyt
        type(pars_type) :: pars
        real*8 :: Hill1
        Hill1=Ca_cyt/(Ca_cyt+pars%IP3%kpi)
        vglu_func = pars%IP3%vbeta * Glu / (Glu+( pars%IP3%kr * (1+pars%IP3%kp/pars%IP3%kr*Hill1) ))
    end function vglu_func

    real*8 function vplcg_func(IP3,Ca_cyt,pars)
        implicit none
        real*8 :: IP3,Ca_cyt
        type(pars_type) :: pars
        vplcg_func = pars%IP3%vdelta/(1+IP3/pars%IP3%kappad)*Hill(Ca_cyt,pars%IP3%kdelta,2)
    end function vplcg_func

    real*8 function v3k_func(IP3,CaMKIIact,pars)
        implicit none
        real*8 :: IP3,CaMKIIact
        type(pars_type) :: pars
        v3k_func = pars%IP3%v3k*CaMKIIact * Hill(IP3,pars%IP3%k3,int(pars%IP3%n3))
    end function v3k_func

    real*8 function dIP3_func(IP3, vip3prod, pars,    v3k)
        implicit none
        real*8 :: IP3, vip3prod,    v3k
        type(pars_type) :: pars
        dIP3_func = vip3prod - v3k - pars%IP3%r5p*IP3
    end function dIP3_func

    real*8 function dDAG_func(DAG, DAGLP, vip3prod, pars)
        implicit none
        real*8 :: DAG, DAGLP, vip3prod
        type(pars_type) :: pars
        type(DGLandDAG_type) :: p
        p = pars%DGLandDAG
        dDAG_func = vip3prod-p%rDGL*DAGLP*DAG/(DAG+p%KDGL)-p%kDAGK*DAG
    end function dDAG_func

    real*8 function dDAGLP_simple_func(DAGLP, Ca_cyt, pars)
        implicit none
        real*8 :: DAGLP, Ca_cyt
        type(pars_type) :: pars
        type(KandP_on_DAGLP_type) :: p
        p = pars%KandP_on_DAGLP
        dDAGLP_simple_func = p%rK*Ca_cyt**int(p%nK)*(1-DAGLP) - p%rP*DAGLP
    end function dDAGLP_simple_func

    subroutine dtwoAG_dAEA_ECb(twoAG, AEA, DAG, DAGLP, Ca_cyt, pars,    dtwoAG, dAEA)
        implicit none
        real*8, intent(in) :: twoAG, AEA, DAG, DAGLP, Ca_cyt
        type(pars_type), intent(in) :: pars
        real*8, intent(out) :: dtwoAG, dAEA
        type(ECb_type) :: pECb
        type(DGLandDAG_type) :: parDnD
        pECb = pars%ECb
        parDnD = pars%DGLandDAG
        !d2-AG/dt
        dtwoAG = parDnD%rDGL * DAGLP * DAG / (DAG + parDnD%KDGL) - parDnD%kMAGL * twoAG
        !dAEA/dt
        dAEA = pECb%vATAEA * Ca_cyt - pECb%vFAAH * AEA / (pECb%KFAAH + AEA)
    end subroutine dtwoAG_dAEA_ECb

    subroutine ctrl1_ctrl2_ECb(x, pars,    ctrl1, ctrl2)
        implicit none
        real*8, intent(in) :: x
        real*8, intent(out) :: ctrl1, ctrl2
        type(pars_type), intent(in) :: pars
        ctrl1 = x + pars%DA%gamma1DA*pars%DA%DA
        ctrl2 = x + pars%DA%gamma2DA*pars%DA%DA
    end subroutine ctrl1_ctrl2_ECb

    subroutine Smooth_Om_ECb(x, fpre, pars, Omega)
        implicit none
        real*8, intent(in) :: x, fpre
        type(pars_type), intent(in) :: pars
        real*8, intent(out) :: Omega
        ! {to tabulate LTPwin}
        integer, parameter :: n_x=100
        real*8, save :: xst, xstep, LTPwin_tab(n_x)
        real*8 :: v(n_x), ka
        integer :: i
        ! end {to tabulate LTPwin}
        real*8 :: x0, LTDw, w, LTDwin, LTPwin
        x0 = 0.5*(pars%ECb%LTDstart + pars%ECb%LTDstop)
        LTDw = x0-pars%ECb%LTDstart
        w = pars%ECb_smooth%kw*LTDw
        ka = pars%ECb_smooth%K + pars%ECb_smooth%kadd * (1-abs(x-x0)/LTDw)**int(pars%ECb_smooth%kn) * Rect(x, x0-LTDw, x0+LTDw)
        LTDwin = (Hill(x-x0, ka, int(pars%ECb_smooth%n))-1)*Rect(x,x0-w,x0+w)

        ! tabulate for speed
        if (subcellular_compute_once) then
            subcellular_compute_once = .false.
            xst = pars%ECb%LTPstart-6*pars%ECb_smooth%tau
            xstep = (2*6*pars%ECb_smooth%tau)/n_x
            forall (i=1:n_x)
                v(i) = xst+xstep*(i-1)
                LTPwin_tab(i) = 1/(1+exp(-(v(i)-pars%ECb%LTPstart)/pars%ECb_smooth%tau))
            end forall
        end if
        LTPwin = lin_interpTable_brds(LTPwin_tab, n_x, (x-xst)/xstep)

        Omega = 1.0 + pars%ECb%LTDMax*LTDwin + pars%ECb%LTPMax*LTPwin
    end subroutine Smooth_Om_ECb

    subroutine Sharp_Om_ECb(x, fpre, pars, Omega)
        implicit none
        real*8, intent(in) :: x, fpre
        type(pars_type), intent(in) :: pars
        real*8, intent(out) :: Omega
        Omega = 1 - pars%ECb%LTDMax*(heaviside(x-pars%ECb%LTDstart)-heaviside(x-pars%ECb%LTDstop)) + heaviside(x-pars%ECb%LTPstart)*pars%ECb%LTPMax
    end subroutine Sharp_Om_ECb

    subroutine dfpre_ECb(Om_func, ctrl1, ctrl2, fpre, pars,    dfpre)
        implicit none
        real*8, intent(in) :: ctrl1, ctrl2, fpre
        type(pars_type), intent(in) :: pars
        real*8, intent(out) :: dfpre
        type(ECb_type) :: pECb
        real*8 :: Omega, taufpre
        interface
            subroutine Om_func(x, fpre, pars, Omega)
                use pars_mod
                implicit none
                real*8, intent(in) :: x, fpre
                type(pars_type), intent(in) :: pars
                real*8, intent(out) :: Omega
            end subroutine Om_func
        end interface
        pECb = pars%ECb
        call Om_func(ctrl1, fpre, pars, Omega)
        !make it static
        taufpre = pECb%P1 / ( pECb%P2**pECb%P3 + ctrl2**pECb%P3 ) + pECb%P4
        !dfpre/dt
        dfpre = (Omega-fpre)/taufpre
    end subroutine dfpre_ECb

end module subcellular