Ex15: Calcium Dynamics

Intracellular calcium concentration pool and calcium-dependent K channels

Model: Calcium Dynamics

A cerebellar granule cell (Maex & De Schutter, 1998) with 7 ion channel types, a decaying calcium concentration pool, and a calcium-dependent K channel:

\[\frac{d[\text{Ca}^{2+}]}{dt} = -\frac{I_{Ca}}{2 F \cdot \text{shell}} - \frac{[\text{Ca}^{2+}] - [\text{Ca}^{2+}]_{\text{rest}}}{\tau}\]

The full composition tree uses neuroml:cell with detailed morphology, channelDensity for each channel, species for Ca²⁺, and a mix of standard NeuroML rate types and custom ComponentType definitions.

---

1. Define Granule Cell in TVBO

from tvbo import SimulationExperiment

exp = SimulationExperiment.from_string("""
label: "NeuroML Ex15: Ca Dynamics"
dynamics:
  name: Granule_98
  iri: neuroml:cell
  parameters:
    diameter:            { value: 10.0 }
    length:              { value: 10.0 }
    specificCapacitance: { value: 1.0, unit: uF_per_cm2 }
    initMembPotential:   { value: -65.0, unit: mV }
    spikeThresh:         { value: 0, unit: mV }
    resistivity:         { value: 0.1, unit: kohm_cm }
    network_temperature: { value: 32, unit: degC }
  components:
    # ── Passive leak ──
    GranPassiveCond:
      name: GranPassiveCond
      iri: neuroml:ionChannelPassive
      parameters:
        conductance: { value: 10, unit: pS }
        condDensity: { value: 0.0330033, unit: mS_per_cm2 }
        erev:        { value: -65.0, unit: mV }
        ion:         { description: non_specific }

    # ── Na fast inactivating channel (gateHHratesTau with custom timeCourse) ──
    Gran_NaF_98:
      name: Gran_NaF_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: na }
        condDensity: { value: 55.7227, unit: mS_per_cm2 }
        erev:        { value: 55.0, unit: mV }
        ion:         { description: na }
      components:
        m:
          name: m
          iri: neuroml:gateHHratesTau
          parameters:
            instances: { value: 3 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 1500, unit: per_s }
                scale:    { value: 0.012345679, unit: V }
                midpoint: { value: -0.028999999999999998, unit: V }
            reverseRate:
              name: reverseRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 1500, unit: per_s }
                scale:    { value: -0.0151515, unit: V }
                midpoint: { value: -0.028999999999999998, unit: V }
            timeCourse:
              name: timeCourse
              iri: neuroml:Gran_NaF_98_m_tau_tau
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V:     { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                ALPHA: { equation: { rhs: "alpha * TIME_SCALE" } }
                BETA:  { equation: { rhs: "beta * TIME_SCALE" } }
                t:     { equation: { rhs: "Piecewise(( ( 0 ) * TIME_SCALE , (ALPHA + BETA) == 0), ( ( 0.00005 ) * TIME_SCALE , 1/(ALPHA + BETA) < ( 0.00005 ) ), ( ( 1/(ALPHA + BETA)) * TIME_SCALE , True))" } }
        h:
          name: h
          iri: neuroml:gateHHratesTau
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 120, unit: per_s }
                scale:    { value: -0.01123596, unit: V }
                midpoint: { value: -0.04, unit: V }
            reverseRate:
              name: reverseRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 120, unit: per_s }
                scale:    { value: 0.01123596, unit: V }
                midpoint: { value: -0.04, unit: V }
            timeCourse:
              name: timeCourse
              iri: neuroml:Gran_NaF_98_h_tau_tau
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V:     { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                ALPHA: { equation: { rhs: "alpha * TIME_SCALE" } }
                BETA:  { equation: { rhs: "beta * TIME_SCALE" } }
                t:     { equation: { rhs: "Piecewise(( ( 0 ) * TIME_SCALE , (ALPHA + BETA) == 0), ( ( 0.000225 ) * TIME_SCALE , 1/(ALPHA + BETA) < ( 0.000225 ) ), ( ( 1/(ALPHA + BETA)) * TIME_SCALE , True))" } }

    # ── Delayed rectifier K channel (custom rates) ──
    Gran_KDr_98:
      name: Gran_KDr_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: k }
        condDensity: { value: 8.89691, unit: mS_per_cm2 }
        erev:        { value: -90.0, unit: mV }
        ion:         { description: k }
      components:
        m:
          name: m
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 4 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:Gran_KDr_98_m_alpha_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "(170 * (exp (73 *(V - (-0.038))))) / TIME_SCALE" } }
            reverseRate:
              name: reverseRate
              iri: neuroml:Gran_KDr_98_m_beta_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "(170 * (exp ((-18) *(V - (-0.038))))) / TIME_SCALE" } }
        h:
          name: h
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:Gran_KDr_98_h_alpha_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "Piecewise(( ( 0.76 ) / TIME_SCALE , V > ( -0.046 ) ), ( ( 0.7 + 0.065*(exp (-80*(V - (-0.046))))) / TIME_SCALE , True))" } }
            reverseRate:
              name: reverseRate
              iri: neuroml:Gran_KDr_98_h_beta_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "(1.1/(1 + (exp (-80.7 * (V - (-0.044)))))) / TIME_SCALE" } }

    # ── Anomalous inward rectifier H channel (standard rates) ──
    Gran_H_98:
      name: Gran_H_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: h }
        condDensity: { value: 0.03090506, unit: mS_per_cm2 }
        erev:        { value: -42.0, unit: mV }
        ion:         { description: h }
      components:
        n:
          name: n
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 0.8, unit: per_s }
                scale:    { value: -0.01100110011, unit: V }
                midpoint: { value: -0.065, unit: V }
            reverseRate:
              name: reverseRate
              iri: neuroml:HHExpRate
              parameters:
                rate:     { value: 0.8, unit: per_s }
                scale:    { value: 0.01100110011, unit: V }
                midpoint: { value: -0.065, unit: V }

    # ── High-voltage activated Ca channel (custom h-gate rates) ──
    Gran_CaHVA_98:
      name: Gran_CaHVA_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: ca }
        condDensity: { value: 0.9084216, unit: mS_per_cm2 }
        erev:        { value: 80.0, unit: mV }
        ion:         { description: ca }
      components:
        m:
          name: m
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 2 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:HHSigmoidRate
              parameters:
                rate:     { value: 1600, unit: per_s }
                scale:    { value: 0.01388888889, unit: V }
                midpoint: { value: 0.015, unit: V }
            reverseRate:
              name: reverseRate
              iri: neuroml:HHExpLinearRate
              parameters:
                rate:     { value: 100, unit: per_s }
                scale:    { value: -0.005, unit: V }
                midpoint: { value: 0.0011000000000000003, unit: V }
        h:
          name: h
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:Gran_CaHVA_98_h_alpha_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "Piecewise(( ( 5.0 ) / TIME_SCALE , V < ( -0.060 ) ), ( ( 5 * (exp (-50 * (V - (-0.060))))) / TIME_SCALE , V >= ( -0.060 ) ))" } }
            reverseRate:
              name: reverseRate
              iri: neuroml:Gran_CaHVA_98_h_beta_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                r: { equation: { rhs: "Piecewise(( ( 0 ) / TIME_SCALE , V < ( -0.060 ) ), ( ( 5 - (5 * (exp (-50 * (V - (-0.060)))))) / TIME_SCALE , V >= ( -0.060 ) ))" } }

    # ── A-type K channel (gateHHtauInf with custom timeCourse) ──
    Gran_KA_98:
      name: Gran_KA_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: k }
        condDensity: { value: 1.14567, unit: mS_per_cm2 }
        erev:        { value: -90.0, unit: mV }
        ion:         { description: k }
      components:
        m:
          name: m
          iri: neuroml:gateHHtauInf
          parameters:
            instances: { value: 3 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 1 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            timeCourse:
              name: timeCourse
              iri: neuroml:Gran_KA_98_m_tau_tau
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                t: { equation: { rhs: "(0.410e-3 * ((exp (( ((V) + 0.0435) / (-0.0428))))) + 0.167e-3) * TIME_SCALE" } }
            steadyState:
              name: steadyState
              iri: neuroml:HHSigmoidVariable
              parameters:
                rate:     { value: 1 }
                scale:    { value: 0.0198, unit: V }
                midpoint: { value: -0.036699999999999997, unit: V }
        h:
          name: h
          iri: neuroml:gateHHtauInf
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 1 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            timeCourse:
              name: timeCourse
              iri: neuroml:Gran_KA_98_h_tau_tau
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V: { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                t: { equation: { rhs: "(0.001) * (10.8 + (0.03 * V * 1000) + (1 / (57.9 * (exp (V * 0.127 * 1000)) + (134e-6 * (exp (V * (-0.059 * 1000))))))) * TIME_SCALE" } }
            steadyState:
              name: steadyState
              iri: neuroml:HHSigmoidVariable
              parameters:
                rate:     { value: 1 }
                scale:    { value: -0.0084, unit: V }
                midpoint: { value: -0.0688, unit: V }

    # ── Ca-dependent K channel (custom rates depending on caConc) ──
    Gran_KCa_98:
      name: Gran_KCa_98
      iri: neuroml:ionChannelHH
      parameters:
        conductance: { value: 10, unit: pS }
        species:     { description: k }
        condDensity: { value: 17.9811, unit: mS_per_cm2 }
        erev:        { value: -90.0, unit: mV }
        ion:         { description: k }
      components:
        m:
          name: m
          iri: neuroml:gateHHrates
          parameters:
            instances: { value: 1 }
          components:
            q10Settings:
              name: q10Settings
              iri: neuroml:q10ExpTemp
              parameters:
                q10Factor:        { value: 3 }
                experimentalTemp: { value: 17.350264793, unit: degC }
            forwardRate:
              name: forwardRate
              iri: neuroml:Gran_KCa_98_m_alpha_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                CONC_SCALE: { value: 1, unit: mol_per_m3 }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V:      { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                CACONC: { equation: { rhs: "caConc / CONC_SCALE" } }
                r:      { equation: { rhs: "(2500/(1 + ( (1.5e-3 *(exp (-85*V))) / CACONC))) / TIME_SCALE" } }
            reverseRate:
              name: reverseRate
              iri: neuroml:Gran_KCa_98_m_beta_rate
              parameters:
                TIME_SCALE: { value: 1, unit: s }
                VOLT_SCALE: { value: 1, unit: V }
                CONC_SCALE: { value: 1, unit: mol_per_m3 }
                offset:     { value: 0.010, unit: V }
              derived_variables:
                V:      { equation: { rhs: "(v - offset) / VOLT_SCALE" } }
                CACONC: { equation: { rhs: "caConc / CONC_SCALE" } }
                r:      { equation: { rhs: "(1500/(1 + (CACONC / (1.5e-4 * (exp (-77*V)))))) / TIME_SCALE" } }

    # ── Calcium concentration pool ──
    Gran_CaPool:
      name: Gran_CaPool
      iri: neuroml:decayingPoolConcentrationModel
      parameters:
        ion:                    { description: ca }
        restingConc:            { value: 7.55e-11, unit: mol_per_cm3 }
        decayConstant:          { value: 10.0, unit: ms }
        shellThickness:         { value: 8.4e-6, unit: cm }
        initialConcentration:   { value: 7.55e-11, unit: mol_per_cm3 }
        initialExtConcentration: { value: 2.4e-6, unit: mol_per_cm3 }

    # ── Pulse generator input ──
    pulseGen1:
      name: pulseGen1
      iri: neuroml:pulseGenerator
      parameters:
        delay:     { value: 100, unit: ms }
        duration:  { value: 500, unit: ms }
        amplitude: { value: 0.01, unit: nA }
network:
  number_of_nodes: 1
integration:
  method: euler
  step_size: 0.001
  duration: 700.0
  time_scale: ms
""")
print(f"Model: {exp.dynamics.name if exp.dynamics else 'network'}")

Model: Granule_98

2. Render LEMS XML

xml = exp.render("lems")
print(xml[:1500])

<Lems>
  <Target component="sim_NeuroML_Ex15__Ca_Dynamics"/>

  <Include file="Cells.xml"/>
  <Include file="Networks.xml"/>
  <Include file="Inputs.xml"/>
  <Include file="Simulation.xml"/>

    <ComponentType name="Gran_CaHVA_98_h_alpha_rate" extends="baseVoltageDepRate">
        <Constant name="TIME_SCALE" dimension="time" value="1 s"/>
        <Constant name="VOLT_SCALE" dimension="voltage" value="1 V"/>
        <Constant name="offset" dimension="voltage" value="0.01 V"/>
        <Dynamics>
            <DerivedVariable name="V" dimension="none" value="(v - offset) / VOLT_SCALE"/>
            <ConditionalDerivedVariable name="r" exposure="r" dimension="per_time">
                <Case condition="V .lt. ( -0.060 )" value="( 5.0 ) / TIME_SCALE"/>
                <Case condition="V .geq. ( -0.060 )" value="( 5 * (exp (-50 * (V - (-0.060))))) / TIME_SCALE"/>
            </ConditionalDerivedVariable>
        </Dynamics>
    </ComponentType>

    <ComponentType name="Gran_CaHVA_98_h_beta_rate" extends="baseVoltageDepRate">
        <Constant name="TIME_SCALE" dimension="time" value="1 s"/>
        <Constant name="VOLT_SCALE" dimension="voltage" value="1 V"/>
        <Constant name="offset" dimension="voltage" value="0.01 V"/>
        <Dynamics>
            <DerivedVariable name="V" dimension="none" value="(v - offset) / VOLT_SCALE"/>
            <ConditionalDerivedVariable name="r" exposure="r" dimension="per_time">
                <Case condition="V .lt. ( -0.060 )" value="( 0 

3. Run Reference

from tvbo.adapters.neuroml import run_lems_example

ref_outputs = run_lems_example("LEMS_NML2_Ex15_CaDynamics.xml")
for name, arr in ref_outputs.items():
    print(f"  {name}: shape={arr.shape}")

  ex15_cond_dens.dat: shape=(700001, 6)
  ex15_curr_dens.dat: shape=(700001, 6)
  ex15_v.dat: shape=(700001, 2)

4. Run TVBO

result = exp.run("neuroml")
da = result.integration.data
print(f"TVBO: {da.dims}, shape={da.shape}")

TVBO: ('time', 'variable'), shape=(700001, 1)

5. Compare & Plot

from tvbo.adapters.neuroml import plot_lems_comparison
plot_lems_comparison("LEMS_NML2_Ex15_CaDynamics.xml", ref_outputs, result.integration.data, title_prefix="Ex15")