Ex10: Q10 Temperature Dependence

Temperature scaling of ion channel rates via Q10 coefficient

Model: Q10 Temperature Dependence

Demonstrates how ion channel rate constants scale with temperature using the Q10 coefficient: \(k(T) = k(T_{\text{ref}}) \cdot Q_{10}^{(T - T_{\text{ref}})/10}\).

At network temperature \(T = 22°C\) with \(T_{\text{ref}} = 32°C\) and \(Q_{10} = 3\):

\[\text{factor} = 3^{(22 - 32)/10} = 3^{-1} \approx 0.333\]

All gating rates are slowed by a factor of 3 compared to reference temperature. The cell is otherwise standard HH with the same conductances as Ex1/Ex5.

1. Define in TVBO

from tvbo import SimulationExperiment

# Q10 factor = 3^((22-32)/10) = 1/3
# All standard HH rates multiplied by Q10_factor = 0.3333
# Same conductances as Ex1/Ex5: g_Na=1200nS, g_K=360nS, g_L=3nS, C=10pF

exp = SimulationExperiment.from_string("""
label: "NeuroML Ex10: Q10 Temperature Dependence"
dynamics:
  name: HodgkinHuxley_Q10
  parameters:
    C:     { value: 10.0 }
    g_Na:  { value: 1200.0 }
    g_K:   { value: 360.0 }
    g_L:   { value: 3.0 }
    E_Na:  { value: 50.0 }
    E_K:   { value: -77.0 }
    E_L:   { value: -54.3 }
    I_amp: { value: 0.08 }
    pulse_delay:    { value: 100.0, unit: ms }
    pulse_duration: { value: 100.0, unit: ms }
    Q10:   { value: 0.333333, description: "Q10 factor = 3^((22-32)/10)" }
  derived_variables:
    I_ext:
      equation:
        rhs: "Piecewise((I_amp, (t >= pulse_delay) & (t < pulse_delay + pulse_duration)), (0.0, True))"
    alpha_m:
      equation:
        rhs: "Q10 * Piecewise((1.0, Eq(v, -40.0)), (0.1*(v + 40.0)/(1.0 - exp(-(v + 40.0)/10.0)), True))"
    beta_m:
      equation:
        rhs: "Q10 * 4.0*exp(-(v + 65.0)/18.0)"
    alpha_h:
      equation:
        rhs: "Q10 * 0.07*exp(-(v + 65.0)/20.0)"
    beta_h:
      equation:
        rhs: "Q10 * 1.0/(1.0 + exp(-(v + 35.0)/10.0))"
    alpha_n:
      equation:
        rhs: "Q10 * Piecewise((0.1, Eq(v, -55.0)), (0.01*(v + 55.0)/(1.0 - exp(-(v + 55.0)/10.0)), True))"
    beta_n:
      equation:
        rhs: "Q10 * 0.125*exp(-(v + 65.0)/80.0)"
  state_variables:
    v:
      equation:
        rhs: "(-g_Na*m**3*h*(v - E_Na) - g_K*n**4*(v - E_K) - g_L*(v - E_L) + I_ext*1000) / C"
      initial_value: -65.0
      variable_of_interest: true
    m:
      equation: { rhs: "alpha_m*(1 - m) - beta_m*m" }
      initial_value: 0.052932
    h:
      equation: { rhs: "alpha_h*(1 - h) - beta_h*h" }
      initial_value: 0.596121
    n:
      equation: { rhs: "alpha_n*(1 - n) - beta_n*n" }
      initial_value: 0.317677
network:
  number_of_nodes: 1
integration:
  method: euler
  step_size: 0.01
  duration: 300.0
  time_scale: ms
""")
print(f"Model: {exp.dynamics.name}")
print(f"Q10 scaling: rates × {1/3:.4f}")
Model: HodgkinHuxley_Q10
Q10 scaling: rates × 0.3333

2. Render LEMS XML

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

<Lems>

  <!-- Tell jLEMS/jNeuroML which component is the simulation entry point. -->
  <Target component="sim_NeuroML_Ex10__Q10_Temperature_Dependence"/>

  <!-- ════════════════════════════════════════════════════════════════
       Dimensions & Units (inline — no external includes needed)
       ════════════════════════════════════════════════════════════════ -->

  <!-- Dimensions -->
  <Dimension name="none"/>
  <Dimension name="time" t="1"/>
  <Dimension name="voltage" m="1" l="2" t="-3" i="-1"/>
  <Dimension name="per_time" t="-1"/>
  <Dimension name="conductance" m="-1" l="-2" t="3" i="2"/>
  <Dimension name="capacitance" m="-1" l="-2" t="4" i="2"/>
  <Dimension name="current" i="1"/>
  <Dimension name="resistance" m="1" l="2" t="-3" i="-2"/>
  <Dimension name="concentration" l="-3" n="1"/>
  <Dimension name="substance" n="1"/>
  <Dimension name="charge" t="1" i="1"/>
  <Dimension name="temperature" k="1"/>

  <!-- Units -->
  <Unit symbol="s" dimension="time" power="0"/>
  <Unit symbol="ms" dimension="time" power="-3"/>
  <Unit symbol="us" dimension="time" power="-6"/>
  <Unit symbol="V" dimension="voltage" power="0"/>
  <Unit symbol="mV" dimension="voltage" power="-3"/>
  <Unit symbol="A" dimension="current" power="0"/>
  <Unit symbol="mA" dimension="current" power="-3"/>
  <Unit symbol="nA" dimension="current" power="-9"/>
  <Unit symbol="pA" dimension="current" power="-12"/>
  <Unit symbol="S" dimension="conductance" power="0"/>
  <Unit symbol="mS" dimension="con

3. Run Reference

import sys, os
sys.path.insert(0, os.path.dirname(os.path.abspath(".")))
from _nml_helpers import run_lems_example

ref_outputs = run_lems_example("LEMS_NML2_Ex10_Q10.xml")
for name, arr in ref_outputs.items():
    print(f"  {name}: shape={arr.shape}")
  hhq10_v.dat: shape=(30001, 2)

4. Run TVBO Version

import numpy as np

result = exp.run("neuroml")
da = result.integration.data
tvbo_arr = np.column_stack([da.coords['time'].values, da.values])
print(f"TVBO: shape={tvbo_arr.shape}")
TVBO: shape=(30001, 5)

5. Numerical Comparison

from _nml_helpers import compare_traces
import numpy as np

ref_arr = list(ref_outputs.values())[0]
ref_v = ref_arr[:, [0, 1]]
tvbo_v = tvbo_arr[:, [0, 1]]

compare_traces(ref_v, tvbo_v, ref_cols=['time', 'v'], tvbo_cols=['time', 'v'])
  v: RMSE=0.579011  max_err=0.661121  corr=-0.379677  ⚠️
{'v': {'rmse': np.float64(0.5790106114873341),
  'max_err': np.float64(0.6611210000000001),
  'corr': np.float64(-0.37967686883902035),
  'close': False}}

6. Plot

from _nml_helpers import plot_comparison

plot_comparison(
    ref_v, tvbo_v,
    ref_cols=['time', 'v'], tvbo_cols=['time', 'v'],
    title="Ex10: Q10 HH (22°C) — NeuroML vs TVBO",
    time_scale=1.0, time_unit="s",
)