Reproducing Neuron Channels in matlab

[aluchko]

I'm working on a script to export QUB (http://www.qub.buffalo.edu) models to NEURON and am having trouble getting the models to match. As an intermediate step I've also been trying to reproduce them in matlab, I can get the same traces in matlab as in QUB, but when I get to NEURON (particularly when I have more than 2 states) the traces don't match.

The instance I'm looking at now is a 3 state model with an SEClamp holding at -20 mV till 100ms, than jumping to 60 mV for 400 ms. The channel states probabilities end up quite divergent, the NEURON open channel hovers around .06 for the opening 100 ms (apparently I'm not setting the STEADYSTATE properly as it hasn't settled) and drops to .004 by the end of the simulation.

In matlab it opens at .2724 and drops to only .1583, I'm not sure why the NEURON trace should be so different.

thanks for any assistance.


This is the channel I create in NEURON

Code: Select all

TITLE channel model imported from QUB

NEURON {
  SUFFIX QUB_K_CHANNEL
  USEION k READ ek WRITE ik
  RANGE gbar, gk, ik, k0_C1_O1, k1_C1_O1, k0_O1_C1, k1_O1_C1, k0_C2_C1, k1_C2_C1, k0_C1_C2, k1_C1_C2
}

UNITS {
  (mA) = (milliamp)
  (mV) = (millivolt)
  (pS) = (picosiemens)
  (um) = (micron)
}

PARAMETER {
  gbar = 0.001  (pS/um2)
  k0_C1_O1 = 0.0005 (/ms)
  k1_C1_O1 = 0.08 (/mV)
  k0_O1_C1 = 0.0006 (/ms)
  k1_O1_C1 = 0.04 (/mV)
  k0_C2_C1 = 0.00081  (/ms)
  k1_C2_C1 = 0.57 (/mV)
  k0_C1_C2 = 0.002  (/ms)
  k1_C1_C2 = 3  (/mV)
}

STATE {
  O1 C2 C1 }

ASSIGNED {
  v (mV)
  ek (mV)
  ik (mA/cm2)
  gk (pS/um2)
  C2C1  (/ms)
  C1C2  (/ms)
  C1O1  (/ms)
  O1C1  (/ms)
}

INITIAL {
  SOLVE states STEADYSTATE sparse
}

BREAKPOINT {
  SOLVE states METHOD sparse
  gk = gbar*O1
  ik = (1e-4)*gk*(v - ek)
}

KINETIC states {
  rates(v)
  CONSERVE O1 + C2 + C1 = 1
  ~ C2 <-> C1 (C2C1, C1C2)
  ~ C1 <-> O1 (C1O1, O1C1)
}

PROCEDURE rates(v(millivolt)) {
  C2C1 = k0_C2_C1 * exp(k1_C2_C1*v)
  C1C2 = k0_C1_C2 * exp(k1_C1_C2*v)
  C1O1 = k0_C1_O1 * exp(k1_C1_O1*v)
  O1C1 = k0_O1_C1 * exp(k1_O1_C1*v)
}

And here is the matlab code

Code: Select all

time_end = 500; % 500 ms

v_rev = 0;
V = zeros(time_end,1);
V(1:100) = -20;
V(101:time_end) = 60;
C = [0, .001, 0]; % conductance vector in pS/um2
K = zeros(3,3,2);

K(:,:,1) = [0, 0.0006, 0.00081; 0.0005, 0, 0; .002, 0, 0]; % the k0
K(:,:,2) = [0, 0.04, 0.57; 0.08, 0, 0; 3, 0, 0]; % k1


Q = zeros(size(K,2),size(K,2));
p = zeros(size(K,2),time_end+1);
I = zeros(time_end);

% calculate the inital Q matrix
for x=1:size(K,2)
  for y=1:size(K,2)
    Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(1)); % q_ij = k0*exp(k1*v)
  end
end
for x=1:size(K,2) % set the diagonals
  Q(x,x) = 0;
  Q(x,x) = - sum(Q(:,x));
end

% solve for steady state
p(:,1) = [Q; ones(1, length(Q))] \ [zeros(length(Q), 1); 1];

for t=1:time_end
    % generate the new transition matrix
    for x=1:size(K,2)
      for y=1:size(K,2)
        Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(t)); % q_ij = k0*exp(k1*v)
      end
    end
    for x=1:size(K,2) % set the diagonals
      Q(x,x) = 0;
      Q(x,x) = - sum(Q(:,x));
    end
    A = expm(Q);
    p(:,t+1) = A*p(:,t);
    p(:,t+1) = p(:,t+1)/sum(p(:,t+1)); % re-normalize
    
    I(t) =  1e-4 * C * p(:,t+1) * (V(t)-v_rev); % current in mA/cm2
end
figure(1)
plot(V)
figure(2)
plot(I)
figure(3)
plot(p(2,:))

[aluchko]

So I'm having 2 issues.

First, the initial steadystate probabilities seem to be off. The neuron channel starts at a much higher percentage of open channels and starts to decline despite the voltage being held constant.

When I use a 2 state matlab model, and set the simulation to have the same starting probabilities, than I can reproduce the neuron trace, but if I calculate my own steady state activation (which seems to be the right one) it doesn't match neuron's.


The second issue is that when I increase to 3 or more channels even if I match the initial probabilities the traces end up diverging, however at the moment I just want to understand why the initial probabilities are wrong.


This is the full 2 state channel model I've been using

Code: Select all

TITLE channel model imported from QUB

NEURON {
	SUFFIX QUB_K_CHANNEL2
	USEION k READ ek WRITE ik
	RANGE gbar, gk, ik, k0_C1_O1, k1_C1_O1, k0_O1_C1, k1_O1_C1
}

UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
	(pS) = (picosiemens)
	(um) = (micron)
}

PARAMETER {
	gbar = 0.001	(pS/um2)
	k0_C1_O1 = 0.0006	(/ms)
	k1_C1_O1 = 0.04	(/mV)
	k0_O1_C1 = 0.0005	(/ms)
	k1_O1_C1 = 0.08	(/mV)
}

STATE {
	O1 C1 }

ASSIGNED {
	v (mV)
	ek (mV)
	ik (mA/cm2)
	gk (pS/um2)
	C1O1	(/ms)
	O1C1	(/ms)
}

INITIAL {
	SOLVE states STEADYSTATE sparse
}

BREAKPOINT {
	SOLVE states METHOD sparse
	gk = gbar*O1
	ik = (1e-4)*gk*(v - ek)
}

KINETIC states {
	rates(v)
	CONSERVE O1 + C1 = 1
	~ C1 <-> O1	(C1O1, O1C1)
}

PROCEDURE rates(v(millivolt)) {
	C1O1 = k0_C1_O1 * exp(k1_C1_O1*v)
	O1C1 = k0_O1_C1 * exp(k1_O1_C1*v)
}

And the session file I've been using to create a patch clamp and some graphs to compare with matlab

Code: Select all

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[6]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}

//Begin SingleCompartment
{
load_file("single.hoc")
}
ocbox_ = new SingleCompartment(0)
ocbox_.inserter = new Inserter(0)
{object_push(ocbox_.inserter)}
{
mt.select("QUB_K_CHANNEL2") i = mt.selected()
ms[i] = new MechanismStandard("QUB_K_CHANNEL2")
ms[i].set("gbar_QUB_K_CHANNEL2", 0.001, 0)
ms[i].set("k0_C1_O1_QUB_K_CHANNEL2", 0.0006, 0)
ms[i].set("k1_C1_O1_QUB_K_CHANNEL2", 0.04, 0)
ms[i].set("k0_O1_C1_QUB_K_CHANNEL2", 0.0005, 0)
ms[i].set("k1_O1_C1_QUB_K_CHANNEL2", 0.08, 0)
mstate[i]= 1
maction(i)
}
{object_pop() doNotify()}
{object_push(ocbox_)}
{inserter.v1.map()}
{endbox()}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.vbox
ocbox_.map("SingleCompartment", 71, 124, 113.28, 115.2)
}
objref ocbox_
//End SingleCompartment

//Begin PointProcessManager
{
load_file("pointman.hoc")
}
{
soma ocbox_ = new PointProcessManager(0)
}
{object_push(ocbox_)}
{
mt.select("VClamp") i = mt.selected()
ms[i] = new MechanismStandard("VClamp")
ms[i].set("dur", 100, 0)
ms[i].set("dur", 400, 1)
ms[i].set("dur", 0, 2)
ms[i].set("amp", -20, 0)
ms[i].set("amp", 60, 1)
ms[i].set("amp", 0, 2)
ms[i].set("gain", 100000, 0)
ms[i].set("rstim", 1, 0)
ms[i].set("tau1", 0.001, 0)
ms[i].set("tau2", 0, 0)
ms[i].set("e0", -0.000599988, 0)
ms[i].set("vo0", 59.9988, 0)
ms[i].set("vi0", 59.9988, 0)
ms[i].set("fac", 0, 0)
mt.select("VClamp") i = mt.selected() maction(i)
hoc_ac_ = 0.5
sec.sec move() d1.flip_to(0)
}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.v1
ocbox_.map("PointProcessManager", 75, 264, 208.32, 326.4)
}
objref ocbox_
//End PointProcessManager
{
xpanel("RunControl", 0)
v_init = -65
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 500
xvalue("t","t", 2 )
tstop = 500
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 0.6
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(296,217)
}
{
save_window_ = new Graph(0)
save_window_.size(0,500,-70,60)
scene_vector_[3] = save_window_
{save_window_.view(0, -70, 500, 130, 576, 97, 300.48, 200.32)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
}
{
save_window_ = new Graph(0)
save_window_.size(-10,500,0,1.3e-05)
scene_vector_[4] = save_window_
{save_window_.view(-10, 0, 510, 1.3e-05, 576, 358, 539.52, 201.28)}
graphList[1].append(save_window_)
save_window_.save_name("graphList[1].")
save_window_.addvar("soma.ik( 0.5 )", 1, 1, 0.8, 0.9, 2)
}
{
save_window_ = new Graph(0)
save_window_.size(0,500,0,1)
scene_vector_[5] = save_window_
{save_window_.view(0, 0, 500, 1, 897, 97, 300.48, 200.32)}
graphList[2].append(save_window_)
save_window_.save_name("graphList[2].")
save_window_.addvar("soma.O1_QUB_K_CHANNEL2( 0.5 )", 1, 1, 0.8, 0.9, 2)
}
objectvar scene_vector_[1]
{doNotify()}

And this is the matlab function that attempts to do the same experiment with the same model as neuron. Calling simplePatch_2state(true) uses the same starting probabilities that neuron uses, and should reproduce the same trace. Setting it to false and the script will use its own steadystate propabilities (which seem to be correct) and thus get a different trace than neuron.

Code: Select all

function simplePatch_2state(use_neuron_prob)

time_end = 500; % 500 ms

v_rev = -77;
V = zeros(time_end,1);
V(1:100) = -20;
V(101:time_end) = 60;
C = [.001, 0]; % conductance vector in pS/um2
K = zeros(2,2,2);

gbar = 0.001;
k0_c1_o1 = 0.0006;
k1_c1_o1 = 0.04;
k0_o1_c1 = 0.0005;
k1_o1_c1 = 0.08;

% put in matrix form so we can calculate the steady state
K(:,:,1) = [0, k0_c1_o1; k0_o1_c1, 0;]; % the k0
K(:,:,2) = [0, k1_c1_o1; k1_o1_c1, 0;]; % k1

Q = zeros(size(K,2),size(K,2));
p = zeros(size(K,2),time_end+1);
I = zeros(time_end);

% calculate the inital Q matrix
for x=1:size(K,2)
  for y=1:size(K,2)
    Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(1)); % q_ij = k0*exp(k1*v)
  end
end
for x=1:size(K,2) % set the diagonals
  Q(x,x) = 0;
  Q(x,x) = - sum(Q(:,x));
end

% solve for steady state
p(:,1) = [Q; ones(1, length(Q))] \ [zeros(length(Q), 1); 1];

if (use_neuron_prob) % use the neuron starting probabilities instead
  p(1,1) = 0.941693;
  p(2,1) = 1 - 0.941693;
end

for t=1:time_end
    % generate the new transition matrix
    c1o1 = k0_c1_o1 * exp(k1_c1_o1*V(t));
    o1c1 = k0_o1_c1 * exp(k1_o1_c1*V(t));
    p(1,t+1) = p(1,t) - p(1,t) * o1c1 + p(2,t) * c1o1;
    p(2,t+1) = p(2,t) + p(1,t) * o1c1 - p(2,t) * c1o1;
    I(t) =  1e-4 * C * p(:,t+1) * (V(t)-v_rev); % current in mV/cm2
end
figure(1)
plot(V)
figure(2)
plot(I)
figure(3)
plot(p(1,:))

I'm sorry for the big code dump but it seemed to be the simplest way to make everything clear. I'd really appreciate if anyone could shed some light on what's going on as I don't understand where neuron is getting those steadystate probabilities from.

thanks,
Aaron

[ted]

First, the initial steadystate probabilities seem to be off. The neuron channel starts at a much higher percentage of open channels and starts to decline despite the voltage being held constant.

The first thing to do is examine the voltage dependence of the forward and reverse rates and of the open fraction itself. Resorting to a compact notation

c: closed state fraction
o: open state fraction
Kco and Koc: forward and reverse rates for c<->o

and taking parameters and equations from the NMODL code, a plot of Kco and Koc and mathematical analysis of the undelying equations reveals that Kco = Koc at "v1/2" = 25 ln(6/5) which is about 4.56 mV, and Kco > Koc for more hyperpolarized v. So steady state open fraction must be 0.5 at about 4.56 mV, approaches 1 as v hyperpolarizes, and falls toward 0 as v depolarizes. And sure enough, that's what the NMODL code does. Is that what your Matlab implementation predicts? Sorry, I'm unfamiliar with Matlab.

[aluchko]

So maybe we're thinking of different things with steadystate.

What I mean is hold at -20 mV for a long period, than what probability vector does the channel settle on?
oc>koc=.0005*exp(.08*-20)
oc>kco=.0006*exp(.04*-20)
oc>koc
0.00010094826
oc>kco
0.00026959738

Put these in a matrix so if you multiply by p it will give the change in p, dp, in the next timestep.
[-koc koc;
kco -kco]

Then set
A=[-koc koc;
kco -kco;
1 1]
y = [0;0;1]

and solve Ax=y, where x is p, p doesn't change when multiplied by the transition matrix, and it sums to 1.

This is what I'm doing in the matlab code and what I thought the steadystate would do in the neuron code, though I may have misunderstood its purpose.

[ted]

What I mean is hold at -20 mV for a long period, than what probability vector does the channel settle on?
oc>koc=.0005*exp(.08*-20)
oc>kco=.0006*exp(.04*-20)
oc>koc
0.00010094826
oc>kco
0.00026959738

Exactly what the code in the mod file computes.

Put these in a matrix so if you multiply by p it will give the change in p, dp, in the next timestep.

In the next timestep the open fraction remains about 0.728 and p doesn't change if the membrane is still being held at -20 mV. If the membrane is not clamped at -20, and no other currents are present, the open fraction of channels gradually increases and membrane potential gradually approaches ek--takes a bit less than 100 seconds to flatten out. Is that what your Matlab code predicts?

[aluchko]

My matlab code predicts the same, that is hits steady state at .728, but when neuron starts the simulation it has an open fraction of 0.941693.

So I guess my question is two parts, first where does 0.941693 come from, and second, does neuron have a capability to calculate the steady state probability, .728?

[ted]

My matlab code predicts the same, that is hits steady state at .728, but when neuron starts the simulation it has an open fraction of 0.941693.

Ah, I see the source of confusion: the RunControl panel specifies that membrane potential will be initialized to -65 mV--see the value assigned to v_init in the ses file, and look at the numeric field next to the RunControl's Init button. So at time t=0, voltage-dependent state variables will have steady state values appropriate for -65 mV. You can see for yourself that v starts out at -65 mV by changing the color of v(.5) to red, or even just by clicking on the v(.5) trace when that Graph is in "Crosshair" mode and dragging the cursor down and to the left . . . sure enough, at t=0 v is -65, at 0.025 it's -21.73..., and only at -0.05 does it finally reach -20.066... Change the value next to the Init button to -20 and you'll see that your mod file works fine. Which should be good news because it suggests that your script to export QUB models to NEURON is working.

I first tried reproducing the problem you reported by setting up my own "experimental rig" but of course I knew to change the value in the field next to the Init button to -20 mv, and so had absolutely no idea what the problem was. It wasn't until I used your ses file that it dawned on me.

To avoid this kind of problem in the future, just make sure that v_init and the voltage clamp's first command potential are identical.

Further comments--

A suggestion: Use SEClamp instead of VClamp viewtopic.php?f=28&t=505

Where are you finding QUB models? Is there an archive of them somewhere?

[aluchko]

Ohh, I never even thought to look at the v_init :(

With the v_init set properly the steady state probabilities are working as they should, and the 2-state matlab model is behaving the same as the neuron model.

However, when I go up to 3-states my matlab model and neuron model start behaving very differently. In matlab I get an initial opening fraction of 0.05397, while in neuron I get 0.461538. In both cases this initial open fraction is maintained, it's just that the trace in matlab is very different from the trace in NEURON.

Note that the matlab model very closely resembles what the same model does in QUB, so it's apparently doing something correct, just a different correct than NEURON. Am I misunderstanding how these NEURON channels are running?

Here's KChannel.mod

Code: Select all

TITLE channel model imported from QUB

NEURON {
        SUFFIX QUB_K_CHANNEL
        USEION k READ ek WRITE ik
        RANGE gbar, gk, ik, k0_C1_O1, k1_C1_O1, k0_O1_C1, k1_O1_C1, k0_C2_C1, k1_C2_C1, k0_C1_C2, k1_C1_C2
}

UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
        (pS) = (picosiemens)
        (um) = (micron)
}

PARAMETER {
        gbar = 0.001    (pS/um2)
        k0_C1_O1 = 0.0005       (/ms)
        k1_C1_O1 = 0.08 (/mV)
        k0_O1_C1 = 0.0006       (/ms)
        k1_O1_C1 = 0.04 (/mV)
        k0_C2_C1 = 0.0008       (/ms)
        k1_C2_C1 = 0.07 (/mV)
        k0_C1_C2 = 0.002        (/ms)
        k1_C1_C2 = .03  (/mV)
}

STATE {
        O1 C2 C1 }

ASSIGNED {
        v (mV)
        ek (mV)
        ik (mA/cm2)
        gk (pS/um2)
        C2C1    (/ms)
        C1C2    (/ms)
        C1O1    (/ms)
        O1C1    (/ms)
}

INITIAL {
        SOLVE states STEADYSTATE sparse
}

BREAKPOINT {
        SOLVE states METHOD sparse
        gk = gbar*O1
        ik = (1e-4)*gk*(v - ek)
}

KINETIC states {
        rates(v)
        CONSERVE O1 + C2 + C1 = 1
        ~ C1 <-> C2     (C2C1, C1C2)
        ~ O1 <-> C1     (C1O1, O1C1)
}

PROCEDURE rates(v(millivolt)) {
        C2C1 = k0_C2_C1 * exp(k1_C2_C1*v)
        C1C2 = k0_C1_C2 * exp(k1_C1_C2*v)
        C1O1 = k0_C1_O1 * exp(k1_C1_O1*v)
        O1C1 = k0_O1_C1 * exp(k1_O1_C1*v)
}

init.hoc

Code: Select all

{load_file("nrngui.hoc")}
objectvar save_window_, rvp_
objectvar scene_vector_[6]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{pwman_place(0,0,0)}

//Begin SingleCompartment
{
load_file("single.hoc")
}
ocbox_ = new SingleCompartment(0)
ocbox_.inserter = new Inserter(0)
{object_push(ocbox_.inserter)}
{
mt.select("QUB_K_CHANNEL") i = mt.selected()
ms[i] = new MechanismStandard("QUB_K_CHANNEL")
mstate[i]= 1
maction(i)
}
{object_pop() doNotify()}
{object_push(ocbox_)}
{inserter.v1.map()}
{endbox()}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.vbox
ocbox_.map("SingleCompartment", 80, 193, 113.28, 115.2)
}
objref ocbox_
//End SingleCompartment

//Begin PointProcessManager
{
load_file("pointman.hoc")
}
{
soma ocbox_ = new PointProcessManager(0)
}
{object_push(ocbox_)}
{
mt.select("SEClamp") i = mt.selected()
ms[i] = new MechanismStandard("SEClamp")
ms[i].set("rs", 1, 0)
ms[i].set("dur1", 100, 0)
ms[i].set("amp1", -20, 0)
ms[i].set("dur2", 400, 0)
ms[i].set("amp2", 60, 0)
ms[i].set("dur3", 0, 0)
ms[i].set("amp3", 0, 0)
mt.select("SEClamp") i = mt.selected() maction(i)
hoc_ac_ = 0.5
sec.sec move() d1.flip_to(0)
}
{object_pop() doNotify()}
{
ocbox_ = ocbox_.v1
ocbox_.map("PointProcessManager", 84, 333, 208.32, 326.4)
}
objref ocbox_
//End PointProcessManager

{
xpanel("RunControl", 0)
v_init = -20
xvalue("Init","v_init", 1,"stdinit()", 1, 1 )
xbutton("Init & Run","run()")
xbutton("Stop","stoprun=1")
runStopAt = 5
xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 )
runStopIn = 1
xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 )
xbutton("Single Step","steprun()")
t = 500
xvalue("t","t", 2 )
tstop = 500
xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 )
dt = 0.025
xvalue("dt","dt", 1,"setdt()", 0, 1 )
steps_per_ms = 40
xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 )
screen_update_invl = 0.05
xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 )
realtime = 0.33
xvalue("Real Time","realtime", 0,"", 0, 1 )
xpanel(305,286)
}
{
save_window_ = new Graph(0)
save_window_.size(0,500,-70,60)
scene_vector_[3] = save_window_
{save_window_.view(0, -70, 500, 130, 585, 166, 300.48, 200.32)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
}
{
save_window_ = new Graph(0)
save_window_.size(-10,500,0,1.3e-05)
scene_vector_[4] = save_window_
{save_window_.view(-10, 0, 510, 1.3e-05, 585, 427, 539.52, 201.28)}
graphList[1].append(save_window_)
save_window_.save_name("graphList[1].")
save_window_.addvar("soma.ik( 0.5 )", 1, 1, 0.8, 0.9, 2)
}
{
save_window_ = new Graph(0)
save_window_.size(0,500,0,1)
scene_vector_[5] = save_window_
{save_window_.view(0, 0, 500, 1, 906, 166, 300.48, 200.32)}
graphList[2].append(save_window_)
save_window_.save_name("graphList[2].")
save_window_.addexpr("soma.O1_QUB_K_CHANNEL( 0.5 )", 1, 1, 0.8, 0.9, 2)
}
objectvar scene_vector_[1]
{doNotify()}

And
matchannel.m

Code: Select all

function matchannel(use_neuron_prob)

time_end = 500; % 500 ms
dt=1;
v_rev = -77;
V = zeros(time_end,1);
V(1:100) = -20;
V(101:time_end) = 60;
C = [.001, 0, 0]; % conductance vector in pS/um2
K = zeros(3,3,2);

gbar = 0.001;
k0_c1_o1 = 0.0005;
k1_c1_o1 = 0.08;
k0_o1_c1 = 0.0006;
k1_o1_c1 = 0.04;
k0_c1_c2 = 0.002;
k1_c1_c2 = .03;
k0_c2_c1 = 0.0008;
k1_c2_c1 = .07;


% put in matrix form so we can calculate the steady state
K(:,:,1) = [0, k0_c1_o1, 0; k0_o1_c1, 0, k0_c2_c1; 0, k0_c1_c2, 0]; % the k0
K(:,:,2) = [0, k1_c1_o1, 0; k1_o1_c1, 0, k1_c2_c1; 0, k1_c1_c2, 0]; % the k1

Q = zeros(size(K,2),size(K,2));
p = zeros(size(K,2),time_end+1);
I = zeros(time_end);

% calculate the inital Q matrix
for x=1:size(K,2)
  for y=1:size(K,2)
    Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(1)); % q_ij = k0*exp(k1*v)
  end
end
for x=1:size(K,2) % set the diagonals
  Q(x,x) = 0;
  Q(x,x) = - sum(Q(:,x));
end

% solve for steady state
p(:,1) = [Q; ones(1, length(Q))] \ [zeros(length(Q), 1); 1];

if (use_neuron_prob) % use the neuron starting probabilities instead
  p(1,1) = 0.6932;
  p(2,1) = 0.259565;
  p(3,1) = p(3,1) - p(1,1) - p(2,1);
end

for t=1:time_end
  % generate the new transition matrix
  c2c1 = k0_c2_c1 * exp(k1_c2_c1*V(t))*dt;
  c1c2 = k0_c1_c2 * exp(k1_c1_c2*V(t))*dt;
  c1o1 = k0_c1_o1 * exp(k1_c1_o1*V(t))*dt;
  o1c1 = k0_o1_c1 * exp(k1_o1_c1*V(t))*dt;
    

  p(1,t+1) = p(1,t) - p(1,t) * o1c1 + p(2,t) * c1o1;
  p(2,t+1) = p(2,t) + p(1,t) * o1c1 + p(3,t) * c2c1 - (p(2,t) * c1o1 + p(2,t) * c1c2);
  p(3,t+1) = p(3,t) + p(2,t) * c1c2 - p(3,t) * c2c1;
  I(t) =  1e-4 * C * p(:,t+1) * (V(t)-v_rev); % current in mV/cm2
end
figure(1)
plot(V)
figure(2)
plot(I)
figure(3)
plot(p(1,:))

Is there some way to put print statements into the nmodl file to help me debug things? I could tweak and recompile the generated C, but I thought there might be a simpler way to do it.

thanks a bunch for the help so far,
Aaron

ps. I don't know of any repository of QUB files but they're pretty easy to generate by hand.

[ted]

I never even thought to look at the v_init

It's _the_ parameter that governs NEURON's default initialization.

when I go up to 3-states my matlab model and neuron model start behaving very differently. In matlab I get an initial opening fraction of 0.05397, while in neuron I get 0.461538.

The NMODL syntax for state transitions is

Code: Select all

~ A <-> B (kAB, kBA)

where kAB and kBA are the rates for the forward and backward transitions. In your two state model, transitions between the two states were specified by

Code: Select all

       ~ C1 <-> O1     (C1O1, O1C1)

Assuming that you interpreted C1O1 and O1C1 as the rates of the forward and backward transitions, there would be a close match between your conceptual model and the equations that the computer would be solving.

However, in your three state model the transitions are specified by

Code: Select all

        ~ C1 <-> C2     (C2C1, C1C2)
        ~ O1 <-> C1     (C1O1, O1C1)

which would be incorrect if C1C2 and O1C1 are your forward rates and the other two are your backward rates. Fix that and your NMODL code should work as you expect.

Is there some way to put print statements into the nmodl file to help me debug things?

NMODL allows insertion of
printf("format string", variable_names_separated_by_commas)
statements that use the customary C syntax for formatting.

[aluchko]

It's _the_ parameter that governs NEURON's default initialization.

I'd been spending a lot of the time with the multiple run fitter which doesn't show the run control panel and I forgot about v_init and ended up thinking that the SEClamp was setting the initial voltage instead (as happens in QUB).

which would be incorrect if C1C2 and O1C1 are your forward rates and the other two are your backward rates. Fix that and your NMODL code should work as you expect.

Sorry about that, it seems that when I was trying to make a small test case, and I went from my script generated NMODL to hand generated, I must have introduced some typos at the same time.

So it looks like my reproductions in matlab were working the whole time, as well as the QUB export script, it was only the v_init that was causing problems and making the rest of the simulations go awry.

thanks for all the help.

[ted]

I'd been spending a lot of the time with the multiple run fitter which doesn't show the run control panel

It is possible, even handy, to use a RunControl panel to specify simulation parameters such as v_init, tstop, dt. Of course, these parameters can also be specified with assignment statements in hoc (or Python).

ended up thinking that the SEClamp was setting the initial voltage instead (as happens in QUB).

The convenience that QUB users may gain from that design decision comes at a price. Works great if there is to be only one voltage clamp, but what if there is more than one, and what happens if the user wants to initialize the model to a different potential than the clamp's first command potential? On the other hand, NEURON users have to remember to ensure that the clamp's amp1 has the same value as v_init, if that's what they want. Then there's the more general issue of whether it makes sense for instrumentation (a voltage clamp) to have the side effect of altering one of the standard run system's major initialization parameters.

it seems that when I was trying to make a small test case, and I went from my script generated NMODL to hand generated, I must have introduced some typos at the same time.

Easy to do. I would probably have chosen quite different names for the NMODL code's parameters and variables--shorter (means less typing, so fewer chances to mistype), and less likely to result in hiccups or reduplication. For example, I might have called this mechanism
k3st : three state k channel
or even
kqub3st
and defined
RANGE gbar, g, i, . . .
STATE { c0 c1 o }
so the BREAKPOINT block statements become
g = gbar*o
i = (1e-4)*g*(v - ek)
ik = i
Then the hoc name for thie conductance at any moment is g_kqub3st, and the current it generates is i_kqub3st

So it looks like my reproductions in matlab were working the whole time, as well as the QUB export script, it was only the v_init that was causing problems and making the rest of the simulations go awry.

Yep. That's the kind of problem I like to have: code that works fine if it weren't for something minor that is easy to fix.

[aluchko]

It is possible, even handy, to use a RunControl panel to specify simulation parameters such as v_init, tstop, dt. Of course, these parameters can also be specified with assignment statements in hoc (or Python).

Yeah, I've written up a control system around the MRF to import data and set up the appropriate experimental conditions (SEClamp settings, and v_init), I'm guessing this is part of the reason I forgot about some of the things I shoved behind the hood.

The convenience that QUB users may gain from that design decision comes at a price. Works great if there is to be only one voltage clamp, but what if there is more than one, and what happens if the user wants to initialize the model to a different potential than the clamp's first command potential? On the other hand, NEURON users have to remember to ensure that the clamp's amp1 has the same value as v_init, if that's what they want. Then there's the more general issue of whether it makes sense for instrumentation (a voltage clamp) to have the side effect of altering one of the standard run system's major initialization parameters.

QUB (or at least QUB-express which I've been using), is designed largely for fitting ion channel models while NEURON is a much more general tool so I think the design decisions fit those purposes.

Easy to do. I would probably have chosen quite different names for the NMODL code's parameters and variables--shorter (means less typing, so fewer chances to mistype), and less likely to result in hiccups or reduplication. For example, I might have called this mechanism

We'll probably be able to make the mechanism name user configurable, as for the state names they're already taken from whatever label the user gives them in QUB.

[ted]

We'll probably be able to make the mechanism name user configurable, as for the state names they're already taken from whatever label the user gives them in QUB.

That makes sense for an automated tool. After verifying that the converted code works properly, users are free to do whatever they want with names.