: ============================================================================
: Cavitational Capacitive Drive (CCD) Model Mechanism (mW/cm2 units)
: Author: Mithun Padmakumar
: Date: July 2026
: 
: Citation:
:   Padmakumar, M., Rajan, D., & Steephen, J. E. (2026). Cavitational capacitive 
:   drive: A computationally efficient model for ultrasonic neuromodulation. 
:   Journal of Neural Engineering.
: 
: Description:
:   This mechanism models the time-varying intramembrane capacitance changes 
:   induced by Focused Ultrasound Stimulation (FUSS) using the Cavitational 
:   Capacitive Drive (CCD) formulation with US intensity recorded in mW/cm2.
: 
: Mechanism Dynamics:
:   The relation q = c(t) * v leads to total capacitive current:
:     i_cap = c(t) * dv/dt + (dc/dt) * v
:   - c(t) is assigned to the compartment's 'cm' via a POINTER (c) in 
:     BEFORE BREAKPOINT.
:   - (dc/dt) * v current contribution is computed in BREAKPOINT as nonspecific current.
: 
: Parameters:
:   - usi    : Ultrasound intensity (mW/cm2), range 10 - 2000 mW/cm2
:   - usf    : Ultrasound frequency (kHz or /ms), range 100 - 1000 kHz
:   - tbegin : Stimulation start time (ms)
:   - tdur   : Stimulation duration (ms)
:   - c1     : Baseline membrane capacitance (default 1.0 uF/cm2)
: 
: Usage Instructions:
:   1. Load parameter lookup tables first: load_file("ccd_tables.hoc")
:   2. Insert mechanism: insert ccd
:   3. Link pointer to segment capacitance: for (x, 0) setpointer c_ccd(x), cm(x)
:   4. Adjust integration time step: dt = 0.025 / usf
: ============================================================================

UNITS {
  (mA) = (milliamp)
  (mV) = (millivolt)
  (uF) = (microfarad)
  PI = (pi) (1)
}

NEURON {
  SUFFIX ccd
  THREADSAFE
  RANGE c1, usi, usf, tbegin, tdur, dc, PWMperiod, PWMdc, k, se, sc
  POINTER c
  NONSPECIFIC_CURRENT i
}
 
PARAMETER {
  c1 = 1 (uF/cm2)
  usf = 200 (/ms)
  usi = 10 (mW/cm2)
  tbegin = 1e9 (ms) 
  tdur = 0 (ms)
  dc (uF/cm2-ms)
  PWMperiod = 10 (ms)
  PWMdc = 1 (1)
}
 
ASSIGNED {
  c (uF/cm2)
  i (mA/cm2)
  v (mV)
  dt (ms)
  k (1)
  se (uF/cm2)
  sc (uF/cm2)
}

FUNCTION_TABLE k100_table(x (mW/cm2)) (1)
FUNCTION_TABLE k200_table(x (mW/cm2)) (1)
FUNCTION_TABLE k500_table(x (mW/cm2)) (1)
FUNCTION_TABLE k800_table(x (mW/cm2)) (1)
FUNCTION_TABLE k1000_table(x (mW/cm2)) (1)

FUNCTION_TABLE se100_table(x (mW/cm2)) (uF/cm2)
FUNCTION_TABLE se200_table(x (mW/cm2)) (uF/cm2)
FUNCTION_TABLE se500_table(x (mW/cm2)) (uF/cm2)
FUNCTION_TABLE se800_table(x (mW/cm2)) (uF/cm2)
FUNCTION_TABLE se1000_table(x (mW/cm2)) (uF/cm2)

FUNCTION_TABLE sc100_table(x (mW/cm2)) (uF/cm2) 
FUNCTION_TABLE sc200_table(x (mW/cm2)) (uF/cm2) 
FUNCTION_TABLE sc500_table(x (mW/cm2)) (uF/cm2) 
FUNCTION_TABLE sc800_table(x (mW/cm2)) (uF/cm2) 
FUNCTION_TABLE sc1000_table(x (mW/cm2)) (uF/cm2) 

FUNCTION cm(t(ms)) (uF/cm2) {
  if ( (usi >= 10) && (usi <= 2000) && (t >= tbegin) && (t <= (tbegin + tdur)) && ( fmod(t,PWMperiod) <= (PWMperiod*PWMdc) )) {
    	if (sin(2*PI*usf*(t-tbegin)) > 0 ) {
		cm = c1 + sc*tanh( k * sin(2*PI*usf*(t - tbegin)) ) / tanh( k )
	} else {
		cm = c1 + se*tanh( k * sin(2*PI*usf*(t - tbegin)) ) / tanh( k )
	}
  }else{
    cm = c1
  }
}
        
FUNCTION dcmdt(t(ms))(uF/cm2-ms) {
  if ( (usi >= 10) && (usi <= 2000) && (t >= tbegin) && (t <= (tbegin + tdur)) && ( fmod(t,PWMperiod) <= (PWMperiod*PWMdc) )) {
    dcmdt = (cm(t+dt) - cm(t) ) / dt 
  }else{
    dcmdt = 0
  }
}

PROCEDURE update_params() {
    LOCAL alpha

    if (usf < 100) {
        k  = k100_table(usi)
        se = se100_table(usi)
        sc = sc100_table(usi)
    
    } else if (usf < 200) {
        
	alpha = (200 - usf)/(200 - 100)
        
        k  = alpha*k100_table(usi)  + (1-alpha)*k200_table(usi)
        se = alpha*se100_table(usi) + (1-alpha)*se200_table(usi)
        sc = alpha*sc100_table(usi) + (1-alpha)*sc200_table(usi)

    } else if (usf < 500) {

        alpha = (500 - usf)/(500 - 200)

        k  = alpha*k200_table(usi)  + (1-alpha)*k500_table(usi)
        se = alpha*se200_table(usi) + (1-alpha)*se500_table(usi)
        sc = alpha*sc200_table(usi) + (1-alpha)*sc500_table(usi)

    } else if (usf < 800) {

        alpha = (800 - usf)/(800 - 500)

        k  = alpha*k500_table(usi)  + (1-alpha)*k800_table(usi)
        se = alpha*se500_table(usi) + (1-alpha)*se800_table(usi)
        sc = alpha*sc500_table(usi) + (1-alpha)*sc800_table(usi)

    } else if (usf < 1000) {

        alpha = (1000 - usf)/(1000 - 800)

        k  = alpha*k800_table(usi)  + (1-alpha)*k1000_table(usi)
        se = alpha*se800_table(usi) + (1-alpha)*se1000_table(usi)
        sc = alpha*sc800_table(usi) + (1-alpha)*sc1000_table(usi)

    } else {

        k  = k1000_table(usi)
        se = se1000_table(usi)
        sc = sc1000_table(usi)
    }
}

INITIAL {
  dc = 0
  update_params()
}

BEFORE BREAKPOINT {
  c = cm(t)
  dc = dcmdt(t)
}

BREAKPOINT {
  at_time(tbegin)
  at_time(tbegin + tdur)
  
  i = dc*v*(0.001)
}