translate KINETIC block to DERIVATIVE block

[ilariac]

Hi,

I am trying to compare two mod files but one is written using the KINETIC block and the other the DERIVATIVE block. Since I would like to solve the ODEs using the same method and not (SOLVE kin METHOD sparse) or (SOLVE deriv METHOD cnexp), I am thinking to rewrite the KINETIC block as the DERIVATIVE block, but I don't know how to write the conservation law.

In the KINETIC block it is simply CONSERVE A + B + C = 1; is there any way to impose the conservation law when the DERIVATIVE block is used instead of the KINETIC one? Can I do like this?

For example if I have:

KINETIC kin {
rates(v)
A <-> B (v1,v2)
B <-> C (v3,v4)

CONSERVE A + B + C = 1
}

can I 'translate' this block in the following way:

DERIVATIVE deriv {
rates(v)
A' = v2*B - v1*A
B' = v1*A - v4*C - (v2+v3)*B
C = 1 - (A + B)
}

So instead of writing C'= v3*B - v4*C, in order to include the conservation law, can I write C=1-(A+B)?

I also wonder how can I set the initial conditions. In the KINETIC block it is simply:

INITIAL {
SOLVE kin
STEADYSTATE sparse
},

what is the equivalent in the DERIVATIVE block?

Thank you for your help,

Ilaria

[ted]

I am trying to compare two mod files but one is written using the KINETIC block and the other the DERIVATIVE block. Since I would like to solve the ODEs using the same method and not (SOLVE kin METHOD sparse) or (SOLVE deriv METHOD cnexp, I am thinking to rewrite the KINETIC block as the DERIVATIVE block

By rewriting the code are you trying to avoid confounding the comparison by minor numerical differences caused by using different numerical methods? If yes, have you done the simplest thing first, i.e. use the files as is, and compare the results? Or do you have a more substantive reason to do the rewrite?

Code: Select all

is there any way to impose the conservation law when the DERIVATIVE block is used instead of the KINETIC one?
 . . .
For example if I have:

KINETIC kin {
             rates(v)
             A <-> B (v1,v2)
             B <-> C (v3,v4)

             CONSERVE A + B + C = 1
}

can I 'translate' this block in the following way:

DERIVATIVE deriv {
             rates(v)
             A' = v2*B - v1*A
             B' = v1*A - v4*C - (v2+v3)*B
             C  = 1 - (A + B)
}

You're almost there. In the right hand side of B' = . . .
substitute the algebraic formula for C, then gather terms and rearrange. You'll have two ODEs whose right hand sides are entirely in terms of A and B, and an algebraic equation for C. Move the algebraic equation to the BREAKPOINT block right after the SOLVE statement, and you're done.

how can I set the initial conditions.

Set the left hand sides of the two ODEs to 0 and treat the resulting equations as a pair of algebraic equations in A and B. Solve them for A and B. Change your initial block to

Code: Select all

INITIAL {
  rates(v)
  A = . . .
  B = . . .
  C = 1 - A - B
}

[ilariac]

Thank you for the reply. I have been thinking about what you wrote for the "initial conditions' issue". Just to make sure I understood, did you mean to solve the linear system by hand and write in the INITIAL block the analytic form of A and B?

Which if I am not wrong is something like:

A= v2v4 / (v2v4 + v1v3 + v1v4)

B= v1v4 / (v2v4 + v1v3 + v1v4)

My problem is that in reality I have 5 differential equations, and the corresponding linear system is messy. I found the analytical solution using Mathematica and it's 3 pages long.

So what are my options? Is it possible to solve such linear system in the .mod?

I solved the linear system in python, but then I don't know how to pass the values to the .mod. My idea was to find iA, iB (initial value for A and initial value for B) and then write

Code: Select all

INITIAL {
    A = iA
    B = iB
    C = 1 - iA - iB :actually I could also set directly C=iC
}

But how to pass these values to the .mod from python (the same python script I use to insert the mechanism in the soma for example) is beyond me. Could you also please help me with that?

I have tried unsuccessfully to do this:

Code: Select all

NEURON {
	SUFFIX na15od
	USEION na READ ena WRITE ina
	RANGE gbar, ina, g, iA, iB, iC
}

.
.
.

PARAMETER {
	v (mV)
	ena (mV)
	celsius
	gbar  = 0.1	 (mho/cm2)

    iA
    iB
    iC

    .
    .
    .
}

.
.
.

INITIAL {
    A = iA
    B = iB
    C = 1 - iA - iB :actually I could also set directly C=iC
}

but it doesn't read the values iA, iB and iC. They are set to the default value 0.


I know that another solution is to use the VERBATIM, I don't have experience in C, but with the help of one of my colleagues we wrote a .c to solve the linear system, but I haven't figured out yet how to include this function we made in the VERBATIM block.

I would prefer to use python though.

Thank you again.

ps: the reason why I want to rewrite the .mod is because I want to compare also the computational cost of these two mod files. KINETIC block should slow down the running time, so I want to have both mod files written using the DERIVATIVE block to properly compare them.

[ilariac]

I solved the problem by writing in the python file:

Code: Select all

for seg in cell:
    seg.mech.iA=initial_values[0]
    seg.mech.iB=initial_values[1]
    seg.mech.iC=initial_values[2]

Thank you for your help,

Ilaria

[ted]

did you mean to solve the linear system by hand and write in the INITIAL block the analytic form of A and B?

Exactly.

Thank you for describing how you dealt with your problem.

the reason why I want to rewrite the .mod is because I want to compare also the computational cost of these two mod files. KINETIC block should slow down the running time, so I want to have both mod files written using the DERIVATIVE block to properly compare them.

Please be sure to say how your test turns out. Unless this is the only mechanism used in your model, AND it is used in every compartment of your model, my guess is that run time will be about the same regardless of whether the mechanism is expressed by a kinetic scheme or as a family of ODEs.

To some readers of this post, 5 ODEs may sound like a lot of equations for a single ion channel mechanism, but it's actually pretty small compared to the equation sets required to describe the kinds of kinetic schemes that channel biophysicists deal with. Check out these examples in ModelDB and imagine reimplementing them with ODEs:
rsg.mod https://senselab.med.yale.edu/modeldb/s ... %2frsg.mod
Narsg.mod https://senselab.med.yale.edu/modeldb/s ... fNarsg.mod
Na.mod https://senselab.med.yale.edu/modeldb/s ... l%2fNa.mod

There's a point at which programmer time becomes more important than computer time.