Merge pull request #65 from openworm/development

Option to show extra plots of driving force and m3h etc.

Author: pgleeson

Tutorial/Source/HodgkinHuxley.py

@@ -15,7 +15,8 @@ def __init__(self, C_m=1, gmax_Na=120, gmax_K=36, gmax_L=0.3, E_Na=50,
                  I_inj_amplitude=0, I_inj_duration=0, I_inj_delay=0,
                  vc_delay=10, vc_duration=30, vc_condVoltage=-65,
                  vc_testVoltage=10, vc_returnVoltage=-65, runMode='iclamp',
-                 injected_current_plot=True, gating_plot=True, cond_dens_plot=True,
+                 injected_current_plot=True, gating_plot=True, cond_scaling_plot=False,
+                 cond_dens_plot=True, driving_force_plot=False,
                  current_plot=True, memb_pot_plot=True):
 
         self.C_m  = C_m
@@ -78,12 +79,15 @@ def __init__(self, C_m=1, gmax_Na=120, gmax_K=36, gmax_L=0.3, E_Na=50,
         # plotting conditionals
         self.injected_current_plot = injected_current_plot
         self.gating_plot = gating_plot
+        self.cond_scaling_plot = cond_scaling_plot
         self.cond_dens_plot = cond_dens_plot
+        self.driving_force_plot = driving_force_plot
         self.current_plot = current_plot
         self.memb_pot_plot = memb_pot_plot
 
         self.num_plots = (int(self.injected_current_plot) +
-                          int(self.gating_plot) + int(self.cond_dens_plot) +
+                          int(self.gating_plot)+ int(self.cond_scaling_plot) +
+                          int(self.cond_dens_plot) + int(self.driving_force_plot) +
                           int(self.current_plot) + int(self.memb_pot_plot))
 
         self.plot_count = 0
@@ -310,10 +314,24 @@ def Main(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]):
                 except NameError:
                     ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)
                     plt.title('Simulation of Hodgkin Huxley model neuron')
-                plt.plot(self.t, m, 'r', label='m')
-                plt.plot(self.t, h, 'g', label='h')
-                plt.plot(self.t, n, 'b', label='n')
-                plt.ylabel('Gating Variable')
+                plt.plot(self.t, m, 'r', label='$m$')
+                plt.plot(self.t, h, 'g', label='$h$')
+                plt.plot(self.t, n, 'b', label='$n$')
+                plt.ylabel('Gating variable')
+                plt.legend()
+                self.plot_count += 1
+
+            if self.cond_scaling_plot:
+                try:
+                    plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)
+                except NameError:
+                    ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)
+                    plt.title('Simulation of Hodgkin Huxley model neuron')
+                scale_na = m*m*m*h
+                scale_k = n*n*n*n
+                plt.plot(self.t, scale_na, 'c', label='$m^{3}h$')
+                plt.plot(self.t, scale_k, 'y', label='$n^{4}$')
+                plt.ylabel('Cond scaling')
                 plt.legend()
                 self.plot_count += 1
 
@@ -325,8 +343,34 @@ def Main(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]):
                     plt.title('Simulation of Hodgkin Huxley model neuron')
                 plt.plot(self.t, gna, 'c', label='$g_{Na}$')
                 plt.plot(self.t, gk, 'y', label='$g_{K}$')
-                plt.ylabel('Cond. dens ($mS/cm^2$)')
+                plt.ylabel('Cond dens ($mS/cm^2$)')
+                plt.legend()
+                self.plot_count += 1
+
+
+            if self.driving_force_plot:
+                try:
+                    ax_here = plt.subplot(self.num_plots,1,self.plot_count+1, sharex = ax1)
+                except NameError:
+                    ax1 = plt.subplot(self.num_plots,1,self.plot_count + 1)
+                    plt.title('Simulation of Hodgkin Huxley model neuron')
+                    ax_here = ax1
+
+                dna = V - self.E_Na
+                dk = V - self.E_K
+                zero = [0 for v in V]
+
+                #plt.plot(self.t, dna, 'c', label='$V - E_{Na}$')
+                ax_here.fill_between(self.t, dna, color='c', alpha=0.5)
+                ax_here.fill_between(self.t, dk, color='y', alpha=0.5)
+
+                plt.plot(self.t, dna, 'c', label='$V_{m} - E_{Na}$', linewidth=0.8)
+                plt.plot(self.t, dk, 'y', label='$V_{m} - E_{K}$', linewidth=0.8)
+                plt.plot(self.t, zero, 'k', linestyle='dashed', linewidth=0.5)
+                plt.ylabel('Driving force (mV)')
                 plt.legend()
+                #if not self.is_vclamp(): plt.ylim(-85,60)
+                #plt.ylim(-1, 40)
                 self.plot_count += 1
 
             if self.current_plot:
@@ -338,7 +382,7 @@ def Main(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]):
                 plt.plot(self.t, ina, 'c', label='$I_{Na}$')
                 plt.plot(self.t, ik, 'y', label='$I_{K}$')
                 plt.plot(self.t, il, 'm', label='$I_{L}$')
-                plt.ylabel('Current ($\\mu{A}/cm^2$)')
+                plt.ylabel('Curr dens ($\\mu{A}/cm^2$)')
                 plt.legend()
                 self.plot_count += 1
 
@@ -350,7 +394,7 @@ def Main(self, init_values=[-64.99584, 0.05296, 0.59590, 0.31773]):
                     plt.title('Simulation of Hodgkin Huxley model neuron')
                 plt.plot(self.t, V, 'k')
                 plt.ylabel('$V_{m}$ (mV)')
-                plt.xlabel('t (ms)')
+                plt.xlabel('Time (ms)')
                 if not self.is_vclamp(): plt.ylim(-85,60)
                 #plt.ylim(-1, 40)
                 self.plot_count += 1

Tutorial/Source/hhcell.cell.nml

@@ -2,15 +2,15 @@
 
 <neuroml xmlns="http://www.neuroml.org/schema/neuroml2"
          xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
-         xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2  https://raw.githubusercontent.com/NeuroML/NeuroML2/master/Schemas/NeuroML2/NeuroML_v2beta3.xsd"   
+         xsi:schemaLocation="http://www.neuroml.org/schema/neuroml2  https://raw.githubusercontent.com/NeuroML/NeuroML2/master/Schemas/NeuroML2/NeuroML_v2.2.xsd"
          id="hhcell">
 
     <include href="passiveChan.channel.nml"/> <!-- Include the channel definitions -->
     <include href="naChan.channel.nml"/>
     <include href="kChan.channel.nml"/>
 
     <cell id="hhcell">
-        
+
         <notes>Conductance based cell model NeuroML2 format: standard Hodgkin Huxley model cell with Na, K and passive conductances</notes>
 
         <morphology id="morphology">
@@ -28,7 +28,7 @@
         <biophysicalProperties id="bioPhys1">
 
             <membraneProperties>
-                        
+
                 <channelDensity id="leak" ionChannel="passiveChan" condDensity="0.3 mS_per_cm2" erev="-54.387mV" ion="non_specific"/>
                 <channelDensity id="naChans" ionChannel="naChan" condDensity="120.0 mS_per_cm2" erev="50.0 mV" ion="na"/>
                 <channelDensity id="kChans" ionChannel="kChan" condDensity="36 mS_per_cm2" erev="-77mV" ion="k"/>
@@ -49,4 +49,3 @@
 
 
 </neuroml>
-

notebooks/Python_HH_version/Python_Notebook_HH.ipynb

@@ -37,7 +37,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "012cf4db99d64e63b2aba86949411ced",
+       "model_id": "c131ca61e49a441e95dff9c74a0b29f9",
        "version_major": 2,
        "version_minor": 0
       },
@@ -51,7 +51,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "76a0b6d59696460dbbf15f755c87eee9",
+       "model_id": "3080e0df7a9645299d1e6f2f7992ef53",
        "version_major": 2,
        "version_minor": 0
       },
@@ -88,11 +88,15 @@
     "\n",
     "2) The second plot shows the **activation/inactivation variables** of the ion channels in the neuron: activation variable **m** and inactivation variable **h** for *Na*, and activation variable **n** for *K*. Activation variables increase as the cell becomes depolarised (**V<sub>m</sub>** increases) and the inactivation variable decreases. **m** responds to these changes quicker than **h** or **n**.\n",
     "\n",
-    "3) The third plot shows the **conductance densities** for Sodium (**g<sub>Na</sub> = gmax<sub>Na</sub> m<sup>3</sup> h**) and Potassium (**g<sub>K</sub> = gmax<sub>K</sub> n<sup>4</sup>**). Note, the leak conductance density value, **g<sub>L</sub>** is not voltage dependent, and so does not vary with time. The maximum conductance densities for the 3 conductances can be changed or set to zero.\n",
+    "3) The third (optional) plot shows the **conductance scaling factors** for Sodium (**m<sup>3</sup> * h**) and Potassium (**n<sup>4</sup>**). Note: these only ever have a value between 0 and 1. \n",
     "\n",
-    "4) The fourth plot shows the influx (negative y-axis) and outflux (positive y-axis) of **ionic currents (I<sub>Na</sub> = g<sub>Na</sub> (V<sub>m</sub> - E<sub>Na</sub>), I<sub>K</sub> = g<sub>K</sub> (V<sub>m</sub> - E<sub>K</sub>) and I<sub>L</sub> = g<sub>L</sub> (V<sub>m</sub> - E<sub>L</sub>**) passing through each type of ion channel being modeled. The reversal potentials of each of these (**E<sub>X</sub>**) can be altered above. \n",
+    "4) The fourth plot shows the **conductance densities** for Sodium (**g<sub>Na</sub> = gmax<sub>Na</sub> m<sup>3</sup> h**) and Potassium (**g<sub>K</sub> = gmax<sub>K</sub> n<sup>4</sup>**). Note, the leak conductance density value, **g<sub>L</sub>** is not voltage dependent, and so does not vary with time. The maximum conductance densities for the 3 conductances can be changed or set to zero.\n",
     "\n",
-    "5) The bottom plot shows the **neural membrane potential activity, V<sub>m</sub>**. The spikes here are called \"action potentials\". This variable is governed by the expression: **dV<sub>m</sub>/dt = (I<sub>inj</sub> - I<sub>Na</sub> - I<sub>K</sub> - I<sub>L</sub>)/C<sub>m</sub>**\n",
+    "5) The fifth (optional) plot shows the **driving force** for the Sodium (**V<sub>m</sub> - E<sub>Na</sub>**) and Potassium (**V<sub>m</sub> - E<sub>K</sub>**) currents. Note, the important value here is how different these are from the zero line, i.e. the magnitude of the force driving ions in or out of the cell. \n",
+    "\n",
+    "6) The sixth plot shows the influx (negative y-axis) and outflux (positive y-axis) of **ionic currents** (**I<sub>Na</sub> = g<sub>Na</sub> (V<sub>m</sub> - E<sub>Na</sub>), I<sub>K</sub> = g<sub>K</sub> (V<sub>m</sub> - E<sub>K</sub>)** and **I<sub>L</sub> = g<sub>L</sub> (V<sub>m</sub> - E<sub>L</sub>)**) passing through each type of ion channel being modeled. The reversal potentials of each of these (**E<sub>X</sub>**) can be altered above. \n",
+    "\n",
+    "7) The bottom plot shows the **neural membrane potential activity, V<sub>m</sub>**. The spikes here are called \"action potentials\". This variable is governed by the expression: **dV<sub>m</sub>/dt = (I<sub>inj</sub> - I<sub>Na</sub> - I<sub>K</sub> - I<sub>L</sub>)/C<sub>m</sub>**\n",
     "\n",
     "\n",
     "![hh](../../Tutorial/_media/equivalentCircuit.png)\n"
@@ -104,7 +108,9 @@
    "source": [
     "<font size=\"5\">Notes</font>\n",
     "\n",
-    ">**Note: Sign of channel current $I_{Na}$, $I_{K}$ and $I_{L}$ are opposite to NeuroML tutorial<br />**\n",
+    ">**Note 1: Python source of the HH model<br />** The source of the HH model in Python as used in this notebook can be found [here](https://github.com/openworm/hodgkin_huxley_tutorial/blob/master/Tutorial/Source/HodgkinHuxley.py).\n",
+    "\n",
+    ">**Note 2: Sign of channel current $I_{Na}$, $I_{K}$ and $I_{L}$ are opposite to NeuroML tutorial<br />**\n",
     "The Python code is interested in the current coming through the channel (so positive if flowing \"out\" like K, negative if flowing \"in\", like Na), whereas for the NeuroML channels (more specifically the channel density element) is what's flowing into the cell (so positive for Na)."
    ]
   },

notebooks/Python_HH_version/ui_widget.py

@@ -171,7 +171,8 @@ def highlight_slider():
 #define toggle button for default values and connect to fucntion call
 showValue_togglebtn = ipywidgets.ToggleButton(value=False,description='Default Values',disabled=False,button_style='info',tooltip='Show/Hide default value below') # 'success', 'info', 'warning', 'danger' or ''
 showValue_togglebtn.observe(showDefault)
-defalultValues = ipywidgets.HTMLMath(value=r"\(C = 1.0\)<br>\(G_{Na} = 120, G_{K} =  36, G_{L} =  0.3\)<br>\(V_{Na} =  50, V_{K} = -77, G_{L} = -54.387\)")
+defalultValues = ipywidgets.HTMLMath(value=r"\(C = %s\)<br>\(G_{Na} = %s, G_{K} = %s, G_{L} =  %s\)<br>\(E_{Na} =  %s, E_{K} = %s, E_{L} = %s\)" % \
+                     (default_capacitance, default_cond_Na, default_cond_K, default_cond_L, default_E_Na, default_E_K, default_E_L))
 defalultValues.layout.display = 'none'
 
 #define toggle buttons for iclamp/vclamp run mode
@@ -205,13 +206,15 @@ def highlight_slider():
 button_row=ipywidgets.HBox([reset_button,showValue_togglebtn])
 
 #plot selectors
-header_plotting = ipywidgets.HTMLMath(value=r"<b> Select plots</b>")
-injected_current_plot_value = ipywidgets.Checkbox(value=True, description="Current injection", disabled=False,)
-gating_plot_value = ipywidgets.Checkbox(value=True, description="Gating variables", disabled=False,)
-cond_dens_plot_value = ipywidgets.Checkbox(value=True, description="Conductance densities", disabled=False,)
-current_plot_value = ipywidgets.Checkbox(value=True, description="Current densities", disabled=False,)
-memb_pot_plot_value = ipywidgets.Checkbox(value=True, description="Membrane potential", disabled=False,)
-plot_selection_row=ipywidgets.VBox([header_plotting, ipywidgets.HBox([injected_current_plot_value, gating_plot_value, cond_dens_plot_value]), ipywidgets.HBox([current_plot_value, memb_pot_plot_value])])
+header_plotting = ipywidgets.HTMLMath(value=r"<b> Select plots to show</b>")
+injected_current_plot_value = ipywidgets.Checkbox(value=True, description="1) Current injection", disabled=False,)
+gating_plot_value = ipywidgets.Checkbox(value=True, description="2) Gating variables", disabled=False,)
+cond_scaling_value = ipywidgets.Checkbox(value=False, description="3) Conductance scaling", disabled=False,)
+cond_dens_plot_value = ipywidgets.Checkbox(value=True, description="4) Conductance densities", disabled=False,)
+driving_force_value = ipywidgets.Checkbox(value=False, description="5) Driving force", disabled=False,)
+current_plot_value = ipywidgets.Checkbox(value=True, description="6) Current densities", disabled=False,)
+memb_pot_plot_value = ipywidgets.Checkbox(value=True, description="7) Membrane potential", disabled=False,)
+plot_selection_row=ipywidgets.VBox([header_plotting, ipywidgets.HBox([injected_current_plot_value, gating_plot_value, cond_scaling_value, cond_dens_plot_value]), ipywidgets.HBox([driving_force_value, current_plot_value, memb_pot_plot_value])])
 
 #layout vertically all the widgets defined above
 modelInputs=ipywidgets.VBox([h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,runMode_iclamp,runMode_vclamp,button_row,defalultValues,plot_selection_row])
@@ -229,7 +232,7 @@ def launch_interactive_widget():
     def runHH(C_m, g_Na, g_K, g_L, E_Na, E_K, E_L, t_0, t_n, delta_t,
               I_inj_max, I_inj_width, I_inj_trans, vc_delay, vc_duration,
               vc_condVoltage, vc_testVoltage, vc_returnVoltage, runMode,
-              injected_current_plot, gating_plot, cond_dens_plot, current_plot, memb_pot_plot):
+              injected_current_plot, gating_plot, cond_scaling_plot, cond_dens_plot, driving_force_plot, current_plot, memb_pot_plot):
 
         highlight_slider()
         runner = HHmodel.HodgkinHuxley(C_m, g_Na, g_K, g_L, E_Na, E_K, E_L,
@@ -240,7 +243,9 @@ def runHH(C_m, g_Na, g_K, g_L, E_Na, E_K, E_L, t_0, t_n, delta_t,
                                        runMode,
                                        injected_current_plot=injected_current_plot,
                                        gating_plot=gating_plot,
+                                       cond_scaling_plot=cond_scaling_plot,
                                        cond_dens_plot=cond_dens_plot,
+                                       driving_force_plot=driving_force_plot,
                                        current_plot=current_plot,
                                        memb_pot_plot=memb_pot_plot)
         # init_values are the steady state values for v,m,h,n at zero current injection
@@ -256,7 +261,9 @@ def runHH(C_m, g_Na, g_K, g_L, E_Na, E_K, E_L, t_0, t_n, delta_t,
                                             'vc_testVoltage':textBox_testVoltage,'vc_returnVoltage':textBox_returnVoltage,
                                             'runMode':runMode_togglebtns, 'injected_current_plot': injected_current_plot_value,
                                                       'gating_plot':gating_plot_value,
+                                                      'cond_scaling_plot':cond_scaling_value,
                                                       'cond_dens_plot':cond_dens_plot_value,
+                                                      'driving_force_plot':driving_force_value,
                                                       'current_plot':current_plot_value,
                                                       'memb_pot_plot':memb_pot_plot_value})