/
ChR2H134R.mod
198 lines (151 loc) · 7.4 KB
/
ChR2H134R.mod
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
COMMENT
-----------------------------------------------------------------------------
ChR2H134R.mod
Model of Channelrhodopsin-2 (mutant H134R)
==========================================
from:
Williams JC, Xu J, Lu Z, Klimas A, Chen X, et al. (2013) Computational Optogenetics: Empirically-Derived Voltage- and Light-Sensitive Channelrhodopsin- 2 Model.
PLoS Comput Biol 9(9): e1003220. doi:10.1371/journal.pcbi.1003220
Original (MATLAB) code by: John C. Williams
Implemented by: Michele Giugliano, SISSA Trieste, 8/8/2018, mgiugliano@gmail.com
This mechanism includes the voltage, light, and temperature dependences of CHR2-H134R.
Set the desired temperature by the hoc assignment statement ==> e.g. celsius = 37
-----------------------------------------------------------------------------
ENDCOMMENT
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...
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 (mV) : Be careful with implementing eqs. 1 and 12!
B = -14.6408 (mV) :
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 : scaling factor for losses of photons due to scattering or absorption
sigma_retinal = 12.E-20 : retinal cross-sectional area, m^2
tauChR2 = 1.3 (ms) : time constant for ChR2 activation
temp = 22 (degC) : original temperature
Q10_Gd1 = 1.97 : Q10 value for the temperature sensitivity
Q10_Gd2 = 1.77 : Q10 value for the temperature sensitivity
Q10_Gr = 2.56 : Q10 value for the temperature sensitivity
Q10_e12dark = 1.1 : Q10 value for the temperature sensitivity
Q10_e21dark = 1.95 : Q10 value for the temperature sensitivity
Q10_epsilon1 = 1.46 : Q10 value for the temperature sensitivity
Q10_epsilon2 = 2.77 : Q10 value for the temperature sensitivity
}
ASSIGNED {
: variables calculated by the present mechanism or by NEURON
v (mV) : Membrane potential
celsius (degC) : Temperature
i (mA/cm2) : Membrane current
Gd1 (1./ms) : rate constant for O1->C1 transition
Gd2 (1./ms) : rate constant for O2->C2 transition
Gr (1./ms) : rate constant for C2->C1 transition
e12 (1./ms) : rate constant for O1->O2 transition
e21 (1./ms) : rate constant for O2->O1 transition
epsilon1 : quantum efficiency for photon absorption from C1
epsilon2 : quantum efficiency for photon absorption from C2
F (1./ms) : photon flux: number of photons per molecule per second
S0 : time- & irradiance-dep. activation func. (post-isomerization)
}
STATE {
: Let's declare the state variables
: Note: in any kinetic scheme with N state, you have N-1 *independent* state variables and
: this is the reason why C2 is not defined below, but rather expressed as (1 - O1 - O2 - C1).
:
: O1, O2, C1 and C2 are the fractions of channels in open and closed states, while p is
: an additional state variable, capturing the kinetics of ChR2 activation
O1 O2 C1 p
}
BREAKPOINT {
SOLVE states METHOD cnexp
Irradiance = 0. : typically with values in the range 0 - 10 mW/mm^2
: The control of the Irradiance waveform within Neuron is very rudimentary and
: it is made explicit below. In other words, the time course of Irradiance is
: upfront defined - in this example - as a single-pulse photoactivation protocol.
: Modify it to fit your needs!
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 = (0.001) * (gmax * (A + B * exp(v /C))/v * (O1 + gamma * O2) * (v - EChR2))
: note: the above expression includes the voltage-dep. rectification of ChR2
}
INITIAL {
: Let's set the state variables to their initial values...
rates(v)
C1 = 1.
O1 = 0.
O2 = 0.
p = 0.
}
DERIVATIVE states {
: The state variables are computed here...
rates(v)
: Note: these are the full equations:
: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
: However, only 3 of them are independent. Let's then define C2 as (1. - C1 - O1 - O2)
O1' = -(Gd1+e12) * O1 + e21 * O2 + epsilon1*F*p * C1
O2' = -(Gd2+e21) * O2 + e12 * O1 + epsilon2*F*p * (1. - C1 - O1 - O2)
C1' = -epsilon1*F*p * C1 + Gd1 * O1 + Gr * (1. - C1 - O1 - 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
: These values are 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