Monitor ion channel current with built-in models

[PikaPei]

Hello,

I used NEURON simulator before. In NEURON, I can use something like record(ina) to record sodium current. Is there a way to monitor the current with built-in models instead of constructing models from scratch?

One thing I can think of is this example (https://brainpy.readthedocs.io/en/latest/tutorial_building/build_conductance_neurons_v2.html#build-a-hh-model-by-composing-existing-ion-channels), it will be great if I can monitor ['INa', 'IK', 'IL'] separately. And this should be useful for plotting the results of voltage clamp simulation.

Another small question: When I search the documentation of ion channels, I am confused about some types of channels. For example, INa_HH1952 and INa_HH1952v2. As far as I know, one master_type is HHTypedNeuron, and the other is Sodium. I don't understand what is the difference between them.

Thank you!

[ztqakita]

Thanks for the question! All the built-in models contain the current() function that can return the current of the ion channel, but users cannot access the return value by using the monitor.

 def current(self, V):
    return self.g_max * self.p ** 3 * self.q * (self.E - V)

Instead, you can use the following methods to get the current values:

1. brainpy.math.for_loop(): Use for loop to iterate the model update and return the current:

   
   def step_run(i, x):
   model.update(i, x)
   INa_current = model.INa.current(self.model.V)
   IK_current = model.IK.current(self.model.V)
   IL_current = model.IL.current(self.model.V)
   return INa_current, IK_current, IL_current

INa_current, IK_current, IL_current = bm.for_loop(step_run, [indices, inputs], progress_bar=True, jit=True)

2. Customize the ion channel model to save the current as `brainpy.math.Variable`. Then you can access the current value by using monitor in `DSRunner`.

The difference between `INa_HH1952` and `INa_HH1952v2` is the form of the `update()` and `reset_state()` functions. For `INa_HH1952` and its `master_type` `HHTypedNeuron`, these functions only receive `V` as parameters:

class IK(bp.dyn.IonChannel): master_type = bp.dyn.HHTypedNeuron

def update(self, V, *args, **kwargs): pass

def reset_state(self, V, *args, **kwargs): pass

This version is compatible with the old BrainPy ion channels.

For `INa_HH1952v2 ` and `Sodium`, these functions receive `V, C, E` as parameters:

class INa(bp.dyn.IonChannel): master_type = bp.dyn.Sodium

def update(self, V, C, E, *args, **kwargs): pass

def reset_state(self, V, C, E, *args, **kwargs): pass

[PikaPei]

Thanks for providing the workaround and details of INa_HH1952 and INa_HH1952v2!

As a practice, I tried to implement with customized DSRunner. My goal is to simulate the voltage clamp. If I separate the runner.run into two times, for example, (1) V=-70mv with 100ms, and (2) use hh.Vclamp[:] = 0.0 to set V=0mV with 10ms, the results are correct. However, if I want to send another array as values for the voltage clamp to DSRunner, the results are totally wrong.

Here is my code:

# Use a new class to create variable for current
class IonChannelCurrent:
    def __init__(self, size, **kwargs):
        super().__init__(size, **kwargs)
        self.I = bm.Variable(bm.zeros(self.num))

    def update(self, V):
        super().update(V)
        self.I.value = self.current(V)

class INa_HH1952_Current(IonChannelCurrent, bp.dyn.INa_HH1952):
    pass

class IK_HH1952_Current(IonChannelCurrent, bp.dyn.IK_HH1952):
    pass

class IL_Current(IonChannelCurrent, bp.dyn.IL):
    pass

# HH model
class HH(bp.dyn.CondNeuGroup):
    def __init__(self, size):
        super().__init__(size=size)

        self.INa = INa_HH1952_Current(size)
        self.IK = IK_HH1952_Current(size)
        self.IL = IL_Current(size, E=-54.387, g_max=0.03)
        self.Vclamp = bm.Variable(bm.ones(self.num) * -70)

    def update(self):
        super().update()
        self.V.value = self.Vclamp

# Run the model (method 1, correct results)
hh = HH(1)

runner = bp.DSRunner(
    hh,
    monitors={"V": hh.V, "INa": hh.INa.I},
)

runner.run(100)
hh.Vclamp[:] = 0.0
runner.run(10)

# Run the model (method 2, wrong results)
hh = HH(1)

v_clamp, duration = bp.inputs.section_input(
    values=[-70, 0],
    durations=[100, 10],
    return_length=True,
)

runner = bp.DSRunner(
    hh,
    monitors={"V": hh.V, "INa": hh.INa.I},
    inputs=("Vclamp", v_clamp, "iter")
)

runner.run(duration)

I still think this should work, but I don't know how to keep going. Could you help me check if I did something wrong? Thank you!

[chaoming0625]

Let me first figure out what you want.

[chaoming0625]

I understand your problem.

Actually, this is very easy to fix. Instead of using inputs as tuple, I suggest to use function:

def give_input():
  i = bp.share['i']
  hh.V[:] = v_clamp[i]

runner = bp.DSRunner(hh, monitors={"V": hh.V, "INa": hh.INa.I}, inputs=give_input)

However, if you would like to use tuple, the correct inputs should be:

runner = bp.DSRunner(
    hh,
    monitors={"V": hh.V, "INa": hh.INa.I},
    inputs=("Vclamp", v_clamp, "iter", '=')  # the defaul is "+"
)

[chaoming0625]

Overall, the code will be:

import brainpy as bp
import brainpy.math as bm
import matplotlib.pyplot as plt

# Use a new class to create variable for current
class IonChannelCurrent:
  def __init__(self, size, **kwargs):
    super().__init__(size, **kwargs)
    self.I = bm.Variable(bm.zeros(self.num))

  def update(self, V):
    super().update(V)
    self.I.value = self.current(V)

class INa_HH1952_Current(IonChannelCurrent, bp.dyn.INa_HH1952):
  pass

class IK_HH1952_Current(IonChannelCurrent, bp.dyn.IK_HH1952):
  pass

class IL_Current(IonChannelCurrent, bp.dyn.IL):
  pass

# HH model
class HH(bp.dyn.CondNeuGroup):
  def __init__(self, size):
    super().__init__(size=size)

    self.INa = INa_HH1952_Current(size)
    self.IK = IK_HH1952_Current(size)
    self.IL = IL_Current(size, E=-54.387, g_max=0.03)
    self.Vclamp = bm.Variable(bm.ones(self.num) * -70)

  def update(self):
    super().update()
    self.V.value = self.Vclamp

# Run the model (method 1, correct results)
hh = HH(1)

runner = bp.DSRunner(
    hh,
    monitors={"V": hh.V, "INa": hh.INa.I},
)

runner.run(100)
hh.Vclamp[:] = 0.0
runner.run(10)

# Run the model (method 2, wrong results)
hh = HH(1)
v_clamp, duration = bp.inputs.section_input(
    values=[-70, 0],
    durations=[100, 10],
    return_length=True,
)

def give_input():
  i = bp.share['i']
  hh.V[:] = v_clamp[i]

runner = bp.DSRunner(hh, monitors={"V": hh.V, "INa": hh.INa.I}, inputs=give_input)

runner.run(duration)
plt.plot(runner.mon.ts, runner.mon["V"], label="V")
plt.show()

[chaoming0625]

For the second question:

When I search the documentation of ion channels, I am confused about some types of channels. For example, INa_HH1952 and INa_HH1952v2. As far as I know, one master_type is HHTypedNeuron, and the other is Sodium. I don't understand what is the difference between them.

This is a very great question!

Actually, INa_HH1952v2 is more flexible to be used in other cases.

For example, if your neuron has a sodium ion, just like calcium which has a dynamically changed reversal potential (E_Na), in this case, you need INa_HH1952v2 since it relies on Sodium to compute its current.

[PikaPei]

Thank you so much for your help! Sorry for the unclear question and bothering you. I've just found there's actually documentation on the operator of DSRunner.

And I also understand the idea of INa_HH1952v2. BrainPy is a great tool. Keep going!