: Ca local fluxe. 2018, Savtchenko et al., 

NEURON {
    POINT_PROCESS GapCaSt
    NONSPECIFIC_CURRENT icagap
    USEION ca READ cai WRITE ica
    RANGE r, icagap, ica, tau1, tau2, ECa, period, amplit
}

UNITS {
    (molar) =	(1/liter)
    (mM) =	(millimolar)
    (um) =	(micron)
    (mA) =	(milliamp)
    FARADAY =	(faraday)	(10000 coulomb)
    PI = (pi)	(1)
    (uS) = (microsiemens)
}

PARAMETER {
    r = 1e16  (ms)
    ECa = 0.00006 (mM)
    CaBasic = 1 (mM/um)
    tau1=100 (ms) <1e-9,1e9>
    tau2 = 1000 (ms) <1e-9,1e9>
    period = 0.100
    amplit=0.0000001
}

ASSIGNED {
    cai (mM)
    diam (um)
    icagap (nanoamp)
    ica  (nanoamp) : so you can plot the calcium current generated by this mechanism
    factor
}

STATE {
    ACa
    BCa
}

INITIAL {
LOCAL tp
    if (tau1/tau2 > .9999) {
        tau1 = .9999*tau2
    }
    ACa = 0
    BCa = 0
    tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
    factor = -exp(-tp/tau1) + exp(-tp/tau2)
    factor = 1/factor
}

BREAKPOINT {
    SOLVE state METHOD cnexp
    : ica = (1e16)*(BCa-ACa)*(cai - ECa)*(2*FARADAY)/((r*CaBasic)*(PI*diam))
    ica = -(1e16)*(BCa-ACa)*(ECa)*(2*FARADAY)/((r*CaBasic)*(PI*diam))
}

DERIVATIVE state {
    ACa' = -ACa/tau1
    BCa' = -BCa/tau2
    : ACa' = amplit*cos(2*3.14*period*t)
    : BCa' = 0
}

NET_RECEIVE(weight) {
    state_discontinuity(BCa, BCa + weight*factor)
    state_discontinuity(ACa, ACa + weight*factor)
}