Creation of channelrhodopsin2 using Neuron

[ram1712]

Hi,
I wish to create a ChannelRhodopsin2 channel using Neuron. I have all the parameters required but i don't know how to proceed further as in i have to create a .mod file



Can somebody please help me out regarding this as to how to proceed with this??

[ted]

Presumably this would be a re-implementation with NEURON of a previously published model. Several such models have been published--which one do you have in mind? Are you familiar with
Williams JC, Xu J, Lu Z, Klimas A, Chen X, Ambrosi CM, Cohen IS, and Entcheva E
Computational Optogenetics: Empirically-Derived Voltage- and Light-Sensitive Channelrhodopsin-2 Model
PLoS Computational Biology, 2013 Sep;9(9):e1003220.
http://www.ploscompbiol.org/article/inf ... bi.1003220

[ram1712]

Actually i'm working on a model of C-fibre that i have currently and i plan to create a ChR2 channel in it using the parameters specified by :--

Theoretical principles underlying optical stimulation of myelinated axons expressing channelrhodopsin-2
R.L. Arlow, T.J. Foutz, C.C. McIntyre
http://www.sciencedirect.com/science/ar ... 221300537X

[ted]

You would have read the article, then, and would know whether or not the authors used NEURON to implement the model they describe. If they did, you might ask them if they will be making their source code available through ModelDB. That's the easiest way to avoid reinventing the wheel.

[ram1712]

No they haven't created using NEURON hence i tried to create the .mod file by my own using the parameters mentioned

Could you please help me figure out if the code that i have written is correct or not??

There is a 4 state markov chain containing 2 open and 2 closed states

Ka1:- Transisition rate from closed state 1 to open state 1 ; Kd1:- Transisition rate from open state 1 to closed state 1

e12:- Transisition rate from open state 1 to open state 2 ; e21:- Transisition rate from open state 2 to open state 1

Ka2:- Transisition rate from open state 2 to closed state 2 ; Kd2:- Transisition rate from closed state 2 to open state 2

Kr:- Transisition rate from closed state 2 to closed state 1


The code is as follows:--

Code: Select all

TITLE  CHANNELRHODOPSIN 2

NEURON  {
    SUFFIX CHR2
	NONSPECIFIC_CURRENT ichr2
	RANGE ichr2
	}
	
UNITS  {
    (um) = (micron)
    (mA) = (milliamp)
	(mV) = (millivolt)	
	(uA) = (microamp)
    (S) = (siemens)
	}
	
PARAMETER {
    Kd1 = 0.130        (/ms)	
	Kd2 = 0.025        (/ms)
	e12 = 0.053        (/ms)
	e21 = 0.023        (/ms)
	Ka1 = 0.369        (/ms)
	Ka2 = 0.369        (/ms)
	Kr  = 0.0004       (/ms)
	g1 = 0.75           (S)
	g2 = 0.75           (S)
	
	}
	
ASSIGNED {
   ichr2       (mA/cm2)
   v       (mV)

}

BREAKPOINT {
	SOLVE states METHOD sparse
		ichr2=V*(g1*O1 + g2*O2)
}   
	

STATE {
    C1 O1 C2 O2
}   
	
	
INITIAL { SOLVE states STEADYSTATE sparse}

KINETIC states {

        ~C1 <->O1       (Ka1,Kd1)
        ~O1 <->O2       (e12,e21)
        ~O2 <->C2       (Ka2,Kd2)
        ~C2 <-> C1      (Kr,0)
        
	CONSERVE C1+O1+C2+O2=1
	}

[ted]

No they haven't created using NEURON

I wonder why the "Experimental" section of that paper would start with this sentence:

Our PONS model system was created in NEURON (v 7.1) using Python (2.7.1)

However, if you prefer to wing it,
1. Check your equations to make sure that the transition rates are correct (they aren't).
2. "SUFFIX" means this will be a density mechanism, which is a reasonable choice. However, a density mechanism must have a conductance parameter that is in density units e.g. S/cm2, and it must generate a current specified in units of mA/cm2.
3. Use modlunit to check your mod file for consistency of units and syntax errors. If you are using UNIX/Linux/OS X

Code: Select all

REPEAT
  at the command line execute modlunit chr2.mod (or whatever you have called this file)
  correct the errors that it reports
UNTIL there are no more errors

If you are using OS X, you could just drag the mod file and drop it onto the modlunit icon.
If you are using MSWin, double click on the modlunit icon in the NEURON program group, then browse to the directory that contains the mod file and double click on the mod file. For either OS X or MSWin, fix the error that is detected, and check again with modlunit.

[ram1712]

I modified the code and tried to compile using nmodll but i got the following error:---

Illegal or unimplemented SOLVE type :states at line 63 (i.e final line of the code)

the code is:---

Code: Select all

TITLE  CHANNELRHODOPSIN 2

NEURON  {
    SUFFIX CHR2
	NONSPECIFIC_CURRENT ichr2
    RANGE gchr2
	}
	
UNITS  {
    (um) = (micron)
    (mA) = (milliamp)
	(mV) = (millivolt)	
	(uA) = (microamp)
    (S) = (siemens)
	(fS)= (femtosiemens)
	}
	
PARAMETER {

    ghcr2=0.75         (mS/cm^2)
    Kd1 = 0.130        (/ms)	
	Kd2 = 0.025        (/ms)
	e12 = 0.053        (/ms)
	e21 = 0.023        (/ms)
	Ka1 = 0.369        (/ms)
	Ka2 = 0.369        (/ms)
	Kr  = 0.0004       (/ms)
	g1 = 50             (fS)
	g2 = 2.5             (fS)
	
	}
	

STATE {
    C1 O1 C2 O2
}   
	
BREAKPOINT {
	SOLVE states METHOD sparse
		ichr2 = gchr2*(v)
}	
	
INITIAL { SOLVE states METHOD sparse}
		
		

KINETIC states {

        ~C1 <->O1       (Ka1,Kd1)
        ~O1 <->O2       (e12,e21)
        ~O2 <->C2       (Ka2,Kd2)
        ~C2 <-> C1      (Kr,0)
        
	CONSERVE C1+O1+C2+O2=1
	}
	
DERIVATIVE states {
	
	C1' = -Ka1*C1 + Kd1*O1 + Kr*C2
	C2' =  Kd2*O2 - Kd*C2 - Kr*C2
	O1' = Ka1*C1 - Kd1*O1 - e12*O1 + e21*O2
	O2' = Ka2*C2 + e12*O1 - Kd2*O2 - e21*O2
   }

[ted]

Before compiling a mod file, it is important to test the file with modlunit. If you're using Linux or OS X, at the command line type
modlunit yourfilename.mod
Under MSWin, click on the modlunit item in the NEURON program group, use the file browser to navigate to the directory that contains your mod file, click on your mod file, then click on the "Selected directory" button. Next select your mod file, then click on "Check Units."

You'll find that modlunit returns the following informative error message:
Multiple declaration of states at line 71 in file whatever.mod

Which means that a mod file can't contain two code blocks that have the same name. The solution is to get rid of one of the "states" blocks, then check the revised code with modlunit again.

modlunit checks for syntax errors and consistency of units. It stops after the first error it finds, so you may have to
test, then revise
many times before all problems have been detected and repaired.

One minor convenience for MSWin users: after you have used modlunit to check a file in a particular directory, the path to that directory is saved so you can go back to it by just clicking on modlunit's "Recent directories" button and selecting the path.

[ram1712]

Thank you very much for letting me know about modllunit.I fixed a lot of errors but there's this one error which just keeps on appearing repeatedly and i couldn't figure it out

The previous primary expression with units :0.01 coul/m2-sec is missing a conversion factor and should read :
(0.001)*()
at line 46 ie ichr2=gchr2*(v)<<ERROR>>

I think its some problem with the units.Is it ???

i have assigned mS/cm2 for gchr2 and mv for v so ichr2 should be mA/cm2.


I have posted the code again for reference:----

Code: Select all

TITLE  CHANNELRHODOPSIN 2

NEURON  {
    SUFFIX CHR2
	NONSPECIFIC_CURRENT ichr2
    RANGE gchr2
	}
	
UNITS  {
    (um) = (micron)
    (mA) = (milliamp)
	(mV) = (millivolt)	
	(uA) = (microamp)
    (S) = (siemens)
	(fS)= (femtosiemens)
	(mS)= (millisiemens)
	
	}
	
PARAMETER {

    gchr2 = 0.75        (mS/cm2)
    Kd1 = 0.130        (/ms)	
	Kd2 = 0.025        (/ms)
	e12 = 0.053        (/ms)
	e21 = 0.023        (/ms)
	Ka1 = 0.369        (/ms)
	Ka2 = 0.369        (/ms)
	Kr  = 0.0004       (/ms)
	g1 = 50             (fS)
	g2 = 2.5             (fS)
	
	}
	
ASSIGNED {
	v       (mV)
	ichr2    (mA/cm2)
}	

STATE {
    C1 O1 C2 O2
}   
	
BREAKPOINT {
	SOLVE states METHOD sparse
		ichr2 = gchr2*(v)
}	
	
INITIAL { SOLVE states METHOD sparse}
		
		

KINETIC states {

        ~C1 <->O1       (Ka1,Kd1)
        ~O1 <->O2       (e12,e21)
        ~O2 <->C2       (Ka2,Kd2)
        ~C2 <-> C1      (Kr,0)
        
	CONSERVE C1+O1+C2+O2=1
	}

[ted]

Indeed, inconsistency of units is the problem. But before we get to that, here are a couple of suggestions.

First, you may have noticed that NEURON's point process and artificial spiking cell classes have names that are capitalized, but that the names of density mechanisms use lower case. This reflects a common practice of using capitalized names for class names (density mechanisms are not treated as classes). Since your NMODL code defines a density mechanism, you might want to follow this practice and use a lower case name for your mechanism, such as
SUFFIX chr2.

Second, you might want to change the name of the current variable from ichr2 to plain old i, so its name in hoc will be i_chr2. Otherwise its hoc name will be something like ichr2_chr2 or ichr2_CHR2. That that would be would be redundant redundant, no no? Besides, i_chr2 is easier to type. For similar reasons, you might also want to change the name of the channel density parameter from gchr2 to plain old g.

i have assigned mS/cm2 for gchr2 and mv for v so ichr2 should be mA/cm2.

Figuring out how to fix inconsistent units is a recurring headache. A quick and easy trick that often works is to start with something you are absolutely certain about, then make one change at a time until the answer jumps out at you. Here's how to get the right answer for this one:
1 S * 1 volt = 1 ___
1 S * 1 mV = 1 ___
1 mS * 1 mV = 1 ___
So if you use units of ___/cm2 for g, i (the product of g and v) will be in units of mA/cm2, and no scale factor will be necessary.

Another trick is to look at the source code for a mod file that has no inconsistency of units, like hh.mod.
gl (S/cm2) : leak conductance
v (mV)
el (mV)
il (mA/cm2)
il = gl*(v - el)
It won't tell you why units are consistent, but at least it'll show you what units to use.
"But that's cheating!"
Only if you cheat yourself by not taking the time to figure out why these units are consistent.


Final hints:
1. Make sure that your kinetic scheme accurately reproduces the authors' original model description.
2. In the INTIAL block you'll want to use
SOLVE states STEADYSTATE sparse
instead of
SOLVE states METHOD sparse
This is discussed in the NEURON book (see Appendix for entries related to
INITIAL block / SOLVE )
3. In the BREAKPOINT block are you sure that
i = g*(v)
? Isn't channel conductance a function of the fraction of open channels, and isn't current proportional to the electrochemical gradient (not just the electric potential difference across the membrane)?

[ram1712]

Yes
Determination of the ChR2 transmembrane channel conductance depends on the transmembrane voltage (Vm), the reversal potential (Ecat, set to 0 mV), and the channel conductance (gChR2).

The ChR2current during illumination (imax) was determined by Ohm’s law:
imax= (Vm-Ecat)gChR2

ChR2 channel conductance depends on the state of the channel, with zero conductance in states C1 and C2, low conductance (g2) in state O2, and high conductance (g1) in state O1

[mgiugliano]

Presumably this would be a re-implementation with NEURON of a previously published model. Several such models have been published--which one do you have in mind? Are you familiar with
Williams JC, Xu J, Lu Z, Klimas A, Chen X, Ambrosi CM, Cohen IS, and Entcheva E
Computational Optogenetics: Empirically-Derived Voltage- and Light-Sensitive Channelrhodopsin-2 Model
PLoS Computational Biology, 2013 Sep;9(9):e1003220.
http://www.ploscompbiol.org/article/inf ... bi.1003220

Hi Ted, Hi everyone!

I join this thread in the hope to be helpful and to receive help and feedback.

I read the above 2013 PLoS paper, obtained the original (MATLAB) code, and attempted its porting to NMODL.

Compilation by nrnivmodl works but the actual simulation silently fails.
No error message appears, but when I insert my new mechanism into a passive single-compartment soma and run the simulation, the voltage trace is plotted at rest for just few samples and then it appears truncated.
I suspect "NaN" are generated. After spending a day on it, I still cannot locate where or why they arise. I am stuck.

May I ask you or anybody else if you spot any obvious mistake in the (NMODL) syntax/style of my code below? Thanks in advance.

Code: Select all

TITLE Channelrhodopsin-2 (mutant H134R) current density 


UNITS {
: Convenient aliases for the units...
	(mS) = (millisiemens)
	(mV) = (millivolt)
	(mA) = (milliamp)
}


NEURON {
: Public interface of the present mechanism...
    :THREADSAFE
	SUFFIX ChR2
	NONSPECIFIC_CURRENT i

	RANGE i, gmax, Irradiance
    RANGE light_intensity, light_delay, pulse_width

    GLOBAL A, B, C, gamma
	GLOBAL wavelength, hc, wloss, sigma_retinal
	GLOBAL Q10_Gd1, Q10_Gd2, Q10_Gr, Q10_e12dark, Q10_e21dark, Q10_epsilon1, Q10_epsilon2
}



PARAMETER {
: Below are constants, variables that are not changed within
: this mechanism, and other variables changed by the user through the hoc code...
	Irradiance      = 0.
	gmax  			= 0.4  		(mS/cm2)	  : maximal conductance
	EChR2 			= 0.     	(mV)          : reversal potential

	light_delay     = 100.		(ms)		  : initial delay, before switching on the light pulse
	pulse_width     = 100.		(ms)		  : width of the light pulse
	light_intensity = 5.					  : mW/mm^2, intensity of the light pulse

    gamma           = 0.1					  : ratio of conductances in states O2/O1, unit-less
	A 				= 10.6408
	B 				= -14.6408
	C 				= -42.7671  (mV)

	wavelength 		= 470		    		  : wavelength of max absorption for retinal, nm
	hc       		= 1.986446E-25  		  : Planck's constant * speed of light, kg m^3/s^2
	wloss    		= 1.3           		  : unit-less
	sigma_retinal 	= 12.E-20       		  : retinal cross-sectional area, m^2

	tauChR2  		= 1.3		(ms)          :

	temp 			= 22	  (degC)		  : original temp 
	Q10_Gd1      	= 1.97					  : Q10 values, temperature sensitivity
	Q10_Gd2      	= 1.77					  : Q10 values, temperature sensitivity
	Q10_Gr       	= 2.56					  : Q10 values, temperature sensitivity
	Q10_e12dark  	= 1.1					  : Q10 values, temperature sensitivity
	Q10_e21dark  	= 1.95					  : Q10 values, temperature sensitivity
	Q10_epsilon1 	= 1.46					  : Q10 values, temperature sensitivity
	Q10_epsilon2 	= 2.77					  : Q10 values, temperature sensitivity
}


ASSIGNED {
: variables calculated by the present mechanism or by NEURON
	v      		(mV)
	celsius		(degC)

	i    		(mA/cm2)

	Gd1    		(1./ms)
	Gd2    		(1./ms)
	Gr     		(1./ms)
	e12    		(1./ms)
	e21    		(1./ms)

	epsilon1
	epsilon2
	F   		(1./ms)
	S0
}


STATE {
: Let's declare the state variables
	O1 O2 C1 C2 p
}


BREAKPOINT {
	SOLVE states METHOD cnexp

	Irradiance = 0.
:    if (t < light_delay)                      { Irradiance = 0. }
:    else if (t < (light_delay + pulse_width)) { Irradiance = light_intensity }
:    else if (t > (light_delay + pulse_width)) { Irradiance = 0. }

	i = gmax * (A + B * exp(v /C)) * (O1 + gamma * O2) * (v - EChR2)

	: includes the voltage-dep. rectification of ChR2
}


INITIAL {
: Let's set the state variables to their initial values... 
	rates(v)
	C1 = 1.
	C2 = 0.
	O1 = 0.
	O2 = 0.
	p  = 0.
}


DERIVATIVE states {
: The state variables are computed here...
    rates(v)
	O1' = -(Gd1+e12)           * O1 + e21 * O2 + epsilon1*F*p * C1
	O2' = -(Gd2+e21)           * O2 + e12 * O1 + epsilon2*F*p * C2
	C1' = -epsilon1*F*p        * C1 + Gd1 * O1 + Gr  * C2   
	C2' = -(epsilon2*F*p + Gr) * C2 + Gd2 * O2 
	p'  =  (S0 - p) / tauChR2
}


UNITSOFF
PROCEDURE rates(v (mV)) {
    LOCAL e12dark, e21dark, logphi0, Ephoton, flux

	: Values at 22 degrees celsius...
	e12dark  = 0.011                                 : ms^-1
	e21dark  = 0.008                                 : ms^-1
	epsilon1 = 0.8535
	epsilon2 = 0.14
	Gd1 = 0.075 + 0.043 * tanh( -(v+20.) / 20.)    	 : dark-adapted deactivation rate, ms^-1
	Gd2 = 0.05                                  	 : ms^-1
	Gr  = 0.0000434587 * exp(-0.0211539274 * v)    	 : recovery rate ms^-1

	: Values adjusted to the temperature specified by the user...
	e12dark  = e12dark  * Q10_e12dark^((celsius-temp)/10.)    : scale with temp, using Q10
	e21dark  = e21dark  * Q10_e21dark^((celsius-temp)/10.)    : scale with temp, using Q10
	epsilon1 = epsilon1 * Q10_epsilon1^((celsius-temp)/10.)   : scale with temp, using Q10
	epsilon2 = epsilon2 * Q10_epsilon2^((celsius-temp)/10.)   : scale with temp, using Q10
	Gd1 	 = Gd1           * Q10_Gd1^((celsius-temp)/10.)	  : scale with temp, using Q10
	Gd2 	 = Gd2           * Q10_Gd2^((celsius-temp)/10.)	  : scale with temp, using Q10
	Gr  	 = Gr             * Q10_Gr^((celsius-temp)/10.)   : scale with temp, using Q10

	if (Irradiance>0) {
		logphi0  = log(1. + Irradiance / 0.024)      : unit-less
	}
	else {
		logphi0 = 0.	 							 : for consistency
	}
	e12      = e12dark + 0.005 * logphi0             : ms^-1
	e21      = e21dark + 0.004 * logphi0             : ms^-1

	S0       = 0.5 * (1. + tanh(120.*(100. * Irradiance - 0.1))) : unit-less
	Ephoton  = 1.E9 * hc / wavelength          : J, scaling wavelength from nm to m	
	flux     = 1000. * Irradiance / Ephoton    : 1/(s*m^2), scaling irradiance from mW/mm^2 to W/m^2
	F        = flux  * sigma_retinal / (wloss * 1000.) : ms^-1, scaling F from 1/s to 1/ms
}
UNITSON

[mgiugliano]

Posting in a prestigious forum as this one has strange (positive) effects: debugging occurs telepathically ;-))))))

I have now solved my (presumed) numerical instability upon removing one of the o.d.e. describing the kinetic scheme: they are NOT all independent and NEURON probably was annoyed by my mistake.

I have also found that the definition of the current i = .... must be prefixed by a minus sign (and possibly also by a 0.001 coefficient, pointed out by modlunit) to give a depolarising response.

I will further validate the mechanism against the PLoS Comp Bio. paper and will post it as an entry on ModelDB, as soon as I am fully confident.

I hope that my previous mistake might be of use to someone.