001 package com.croftsoft.core.ai.neuro.imp;
002
003 import com.croftsoft.core.ai.neuro.Channel;
004 import com.croftsoft.core.ai.neuro.HhNeuronMut;
005 import com.croftsoft.core.lang.NullException;
006 import com.croftsoft.core.math.MathLib;
007 import com.croftsoft.core.sim.DeltaClock;
008 import com.croftsoft.core.sim.Sim;
009 import com.croftsoft.core.util.seq.Seq;
010
011 /***********************************************************************
012 * Hodgkin-Huxley (HH) neuron.
013 *
014 * Voltages are in units of volts instead of millivolts.
015 * Area-based units are in meters squared instead of centimeters squared.
016 *
017 * For a good description of how the HH equations are derived, see
018 * Nelson, "Chapter 17: Electrophysiological Models", p285, in
019 * Koslow and Subramaniam, "Databasing the Brain: From Data to
020 * Knowledge (Neuroinformatics)", 2005.
021 *
022 * Other references used in implementing the HH equations include:
023 *
024 * Doi and Kumagai, "Nonlinear Dynamics of Small-Scale Biophysical Neural
025 * Networks", Chapter 10, p263, in Poznanski, "Biophysical Neural
026 * Networks: Foundations of Integrative Neuroscience", 2001.
027 *
028 * Feng, "Computational Neuroscience: A Comprehensive Approach", 2004,
029 * p4 and p30.
030 *
031 * @version
032 * $Id: HhNeuronImp.java,v 1.14 2008/09/07 01:48:31 croft Exp $
033 * @since
034 * 2008-08-08
035 * @author
036 * <a href="https://www.croftsoft.com/">David Wallace Croft</a>
037 ***********************************************************************/
038
039 public final class HhNeuronImp
040 implements HhNeuronMut, Sim
041 ////////////////////////////////////////////////////////////////////////
042 ////////////////////////////////////////////////////////////////////////
043 {
044
045 // Area-based units have been converted from per centimeters squared to
046 // per meters squared. Voltages converted from millivolts to volts.
047
048 public static final double
049 H_INIT = 0.596111046, // steady-state resting
050 LEAK_CONDUCTIVITY = 0.3e+1, // 0.3 mS/cm^2
051 LEAK_REVERSAL_POTENTIAL = -54.4e-3, // -54.4 mV
052 LHOPITAL_THRESHOLD = 1e-9, // L'Hopital's Rule
053 M_INIT = 0.052934218, // steady-state resting
054 MEMBRANE_CAPACITIVITY = 1e-2, // 1 microF/cm^2
055 MEMBRANE_VOLTAGE_INIT = -0.064999722, // steady-state resting
056 N_INIT = 0.317681168, // steady-state resting
057 POTASSIUM_CONDUCTIVITY = 36e+1, // 36 mS/cm^2
058 POTASSIUM_REVERSAL_POTENTIAL = -77e-3, // -77 mV
059 SODIUM_CONDUCTIVITY = 120e+1, // 120 mS/cm^2
060 SODIUM_REVERSAL_POTENTIAL = 50e-3, // 50 mV
061 THRESHOLD = 0; // crossing for spike
062
063 /** 10 microns in diameter in meters. */
064 private static final double DIAMETER = 10e-6;
065
066 /** Area of sphere of 10 microns in diameter in meters squared */
067 public static final double
068 MEMBRANE_AREA = Math.PI * DIAMETER * DIAMETER;
069
070 private final Seq<Channel> channelSeq;
071
072 private final DeltaClock deltaClock;
073
074 //
075
076 private double
077 deltaH,
078 deltaMembraneVoltage,
079 deltaM,
080 deltaN,
081 h,
082 leakConductivity,
083 leakReversalPotential,
084 m,
085 membraneArea,
086 membraneCapacitivity,
087 membraneVoltage,
088 n,
089 potassiumConductivity,
090 potassiumReversalPotential,
091 sodiumConductivity,
092 sodiumReversalPotential,
093 threshold;
094
095 private boolean
096 spiking,
097 willSpike;
098
099 ////////////////////////////////////////////////////////////////////////
100 ////////////////////////////////////////////////////////////////////////
101
102 public HhNeuronImp (
103 final Seq<Channel> channelSeq,
104 final DeltaClock deltaClock,
105 final double h,
106 final double leakConductivity,
107 final double leakReversalPotential,
108 final double m,
109 final double membraneArea,
110 final double membraneCapacitivity,
111 final double membraneVoltage,
112 final double n,
113 final double potassiumConductivity,
114 final double potassiumReversalPotential,
115 final double sodiumConductivity,
116 final double sodiumReversalPotential,
117 final double threshold )
118 ////////////////////////////////////////////////////////////////////////
119 {
120 NullException.check (
121 this.channelSeq = channelSeq,
122 this.deltaClock = deltaClock );
123
124 this.h = h;
125
126 this.leakConductivity = leakConductivity;
127
128 this.leakReversalPotential = leakReversalPotential;
129
130 this.m = m;
131
132 this.membraneArea = membraneArea;
133
134 this.membraneCapacitivity = membraneCapacitivity;
135
136 this.membraneVoltage = membraneVoltage;
137
138 this.n = n;
139
140 this.potassiumConductivity = potassiumConductivity;
141
142 this.potassiumReversalPotential = potassiumReversalPotential;
143
144 this.sodiumConductivity = sodiumConductivity;
145
146 this.sodiumReversalPotential = sodiumReversalPotential;
147
148 this.threshold = threshold;
149 }
150
151 public HhNeuronImp (
152 final Seq<Channel> channelSeq,
153 final DeltaClock deltaClock )
154 ////////////////////////////////////////////////////////////////////////
155 {
156 this (
157 channelSeq,
158 deltaClock,
159 H_INIT,
160 LEAK_CONDUCTIVITY,
161 LEAK_REVERSAL_POTENTIAL,
162 M_INIT,
163 MEMBRANE_AREA,
164 MEMBRANE_CAPACITIVITY,
165 MEMBRANE_VOLTAGE_INIT,
166 N_INIT,
167 POTASSIUM_CONDUCTIVITY,
168 POTASSIUM_REVERSAL_POTENTIAL,
169 SODIUM_CONDUCTIVITY,
170 SODIUM_REVERSAL_POTENTIAL,
171 THRESHOLD );
172 }
173
174 ////////////////////////////////////////////////////////////////////////
175 // accessor methods
176 ////////////////////////////////////////////////////////////////////////
177
178 public double getH ( )
179 ////////////////////////////////////////////////////////////////////////
180 {
181 return h;
182 }
183
184 public double getLeakConductivity ( )
185 ////////////////////////////////////////////////////////////////////////
186 {
187 return leakConductivity;
188 }
189
190 public double getM ( )
191 ////////////////////////////////////////////////////////////////////////
192 {
193 return m;
194 }
195
196 public double getMembraneArea ( )
197 ////////////////////////////////////////////////////////////////////////
198 {
199 return membraneArea;
200 }
201
202 public double getMembraneCapacitivity ( )
203 ////////////////////////////////////////////////////////////////////////
204 {
205 return membraneCapacitivity;
206 }
207
208 public double getMembraneVoltage ( )
209 ////////////////////////////////////////////////////////////////////////
210 {
211 return membraneVoltage;
212 }
213
214 public double getN ( )
215 ////////////////////////////////////////////////////////////////////////
216 {
217 return n;
218 }
219
220 public double getPotassiumConductivity ( )
221 ////////////////////////////////////////////////////////////////////////
222 {
223 return potassiumConductivity;
224 }
225
226 public double getSodiumConductivity ( )
227 ////////////////////////////////////////////////////////////////////////
228 {
229 return sodiumConductivity;
230 }
231
232 public double getThreshold ( )
233 ////////////////////////////////////////////////////////////////////////
234 {
235 return threshold;
236 }
237
238 public boolean isSpiking ( )
239 ////////////////////////////////////////////////////////////////////////
240 {
241 return spiking;
242 }
243
244 ////////////////////////////////////////////////////////////////////////
245 // mutator methods
246 ////////////////////////////////////////////////////////////////////////
247
248 public void setLeakConductivity ( final double leakConductivity )
249 ////////////////////////////////////////////////////////////////////////
250 {
251 this.leakConductivity = leakConductivity;
252 }
253
254 public void setMembraneArea ( final double membraneArea )
255 ////////////////////////////////////////////////////////////////////////
256 {
257 this.membraneArea = membraneArea;
258 }
259
260 public void setMembraneCapacitivity (
261 final double membraneCapacitivity )
262 ////////////////////////////////////////////////////////////////////////
263 {
264 this.membraneCapacitivity = membraneCapacitivity;
265 }
266
267 public void setMembraneVoltage ( final double membraneVoltage )
268 ////////////////////////////////////////////////////////////////////////
269 {
270 this.membraneVoltage = membraneVoltage;
271 }
272
273 public void setPotassiumConductivity (
274 final double potassiumConductivity )
275 ////////////////////////////////////////////////////////////////////////
276 {
277 this.potassiumConductivity = potassiumConductivity;
278 }
279
280 public void setSodiumConductivity ( final double sodiumConductivity )
281 ////////////////////////////////////////////////////////////////////////
282 {
283 this.sodiumConductivity = sodiumConductivity;
284 }
285
286 public void setThreshold ( final double threshold )
287 ////////////////////////////////////////////////////////////////////////
288 {
289 this.threshold = threshold;
290 }
291
292 ////////////////////////////////////////////////////////////////////////
293 // lifecycle methods
294 ////////////////////////////////////////////////////////////////////////
295
296 public void access ( )
297 ////////////////////////////////////////////////////////////////////////
298 {
299 final double deltaTime = deltaClock.getDeltaTime ( );
300
301 final int size = channelSeq.size ( );
302
303 double capacitiveCurrent = 0;
304
305 for ( int i = 0; i < size; i++ )
306 {
307 final Channel channel = channelSeq.get ( i );
308
309 if ( channel.isOpen ( ) )
310 {
311 final double channelCurrent = channel.getConductance ( )
312 * ( membraneVoltage - channel.getReversalPotential ( ) );
313
314 capacitiveCurrent += channelCurrent;
315 }
316 }
317
318 // Multiplying by 1000 to adjust for HH values given in millivolts.
319
320 final double v = ( membraneVoltage - -65e-3 ) * 1000;
321
322
323 // L'Hopital's Rule is applied when the dividend is nearly zero.
324
325 final double alphaM
326 = Math.abs ( 25.0 - v ) < LHOPITAL_THRESHOLD ? 1.0 // -0.1 / -0.1
327 : 0.1 * ( 25.0 - v )
328 / ( Math.exp ( ( 25.0 - v ) / 10.0 ) - 1.0 );
329
330 final double alphaN
331 = Math.abs ( 10.0 - v ) < LHOPITAL_THRESHOLD ? 0.1 // -0.01 / -0.1
332 : 0.01 * ( 10.0 - v )
333 / ( Math.exp ( ( 10.0 - v ) / 10.0 ) - 1.0 );
334
335 final double alphaH
336 = 0.07 * Math.exp ( -v / 20.0 );
337
338 final double betaM
339 = 4.0 * Math.exp ( -v / 18.0 );
340
341 final double betaN
342 = 0.125 * Math.exp ( -v / 80.0 );
343
344 final double betaH
345 = 1.0 / ( Math.exp ( ( 30.0 - v ) / 10.0 ) + 1.0 );
346
347 // Alternate equations for deltaM, deltaN, and deltaH.
348 //
349 // final double
350 // hInf = alphaH / ( alphaH + betaH ),
351 // mInf = alphaM / ( alphaM + betaM ),
352 // nInf = alphaN / ( alphaN + betaN );
353 //
354 // final double
355 // tauH = 1 / ( alphaH + betaH ),
356 // tauM = 1 / ( alphaM + betaM ),
357 // tauN = 1 / ( alphaN + betaN );
358 //
359 // deltaM = deltaTime * ( mInf - m ) / tauM;
360 //
361 // deltaN = deltaTime * ( nInf - n ) / tauN;
362 //
363 // deltaH = deltaTime * ( hInf - h ) / tauH;
364
365 // Multiplied by 1000 to adjust for HH rates in per millisecond units.
366
367 deltaM = deltaTime * 1000 * ( alphaM * ( 1.0 - m ) - betaM * m );
368
369 deltaN = deltaTime * 1000 * ( alphaN * ( 1.0 - n ) - betaN * n );
370
371 deltaH = deltaTime * 1000 * ( alphaH * ( 1.0 - h ) - betaH * h );
372
373 final double
374 leakConductance = membraneArea * leakConductivity,
375 membraneCapacitance = membraneArea * membraneCapacitivity,
376 potassiumConductance = membraneArea * potassiumConductivity,
377 sodiumConductance = membraneArea * sodiumConductivity;
378
379 final double hhCurrent
380 = sodiumConductance * m * m * m * h
381 * ( membraneVoltage - sodiumReversalPotential )
382 + potassiumConductance * n * n * n * n
383 * ( membraneVoltage - potassiumReversalPotential )
384 + leakConductance
385 * ( membraneVoltage - leakReversalPotential );
386
387 capacitiveCurrent += hhCurrent;
388
389 final double deltaCharge = deltaTime * -capacitiveCurrent;
390
391 deltaMembraneVoltage = deltaCharge / membraneCapacitance;
392
393 willSpike = ( membraneVoltage < threshold )
394 && ( membraneVoltage + deltaMembraneVoltage >= threshold );
395 }
396
397 public void mutate ( )
398 ////////////////////////////////////////////////////////////////////////
399 {
400 h = MathLib.clip ( h + deltaH, 0, 1 );
401
402 m = MathLib.clip ( m + deltaM, 0, 1 );
403
404 n = MathLib.clip ( n + deltaN, 0, 1 );
405
406 membraneVoltage += deltaMembraneVoltage;
407
408 spiking = willSpike;
409 }
410
411 ////////////////////////////////////////////////////////////////////////
412 ////////////////////////////////////////////////////////////////////////
413 }