with A83_014_Vector; use A83_014_Vector; with ConsAKD ; use ConsAKD; with CursAKD ; use CursAKD; with DOS_011_Plot ; use DOS_011_Plot; -- with FileAK ; use FileAK ; with InteAKD; use InteAKD; with FloaAKD; use FloaAKD; with MathAK ; -- for Interpolate with MathAKD; use MathAKD; with SounAKD; -- for Speaker_Thump with StriAK ; use StriAK; -- for Image with TextAKD; use TextAKD; package body HHNM is ---------------------------------------------------------------------- ---------------------------------------------------------------------- -- Copyright (C) 1994 David Wallace Croft. All rights reserved. ---------------------------------------------------------------------- ---------------------------------------------------------------------- -- procedure Bandpass_Filter ( -- Freq_Min : in float := Freq_Min_Default; -- Freq_Max : in float := Freq_Max_Default; -- Freq_Step : in float := Freq_Step_Default; -- AC_Ampli : in float := AC_Ampli_Default; -- DC_Value : in float := DC_Value_Default; -- Threshold : in float := Threshold_Default; -- Time_Max : in float := Time_Max_Default; -- Time_Delta : in float := Time_Delta_Default; -- Print_Step : in boolean := false; -- TPS : in float := TPS_Default; -- CM : in float := CM_Default; -- V_K : in float := V_K_Default; -- V_Na : in float := V_Na_Default; -- V_Rest : in float := V_Rest_Default; -- G_K_Mul : in float := G_K_Mul_Default; -- G_Na_Mul : in float := G_Na_Mul_Default; -- J_L_Mul : in float := J_L_Mul_Default; -- J_L_Sub : in float := J_L_Sub_Default; -- Log_File : in string := ""; -- Interpolate_On : in boolean := false ) is -- ---------------------------------------------------------------------- -- T_Print : float := 0.0; -- Time : float := 0.0; -- Freq : float; -- stimulus frequency in Hertz (Hz) = cycles/sec. -- Soma : Soma_Type; -- J_Stim : float; -- current of stimulus -- V_Old : float; -- Bins : long_integer := 1 -- + long_integer ( ( Freq_Max - Freq_Min ) / Freq_Step ); -- Spikes : natural; -- Spikes_DC : natural := 0; -- File : File_Type; -- Rec : MathAK.Interpolate_Type; -- Completed : boolean; -- End_Points_Done : boolean := false; -- Temp_File : constant string := "Bandpass.Tmp"; -- Bin : long_integer; -- begin -- if ( Log_File /= "" ) then -- Create ( File, Out_File, Log_File ); -- if not Interpolate_On then -- Put ( File, "Log_File : " ); Put_Line ( File, Log_File ); -- Put ( File, "Freq_Min : " ); Put_Line ( File, Freq_Min ); -- Put ( File, "Freq_Max : " ); Put_Line ( File, Freq_Max ); -- Put ( File, "Freq_Step : " ); Put_Line ( File, Freq_Step ); -- Put ( File, "AC_Ampli : " ); Put_Line ( File, AC_Ampli ); -- Put ( File, "DC_Value : " ); Put_Line ( File, DC_Value ); -- Put ( File, "Threshold : " ); Put_Line ( File, Threshold ); -- Put ( File, "Time_Max : " ); Put_Line ( File, Time_Max ); -- Put ( File, "Time_Delta: " ); Put_Line ( File, Time_Delta ); -- Put_Line ( File, "" ); -- end if; -- end if; -- Freq := Freq_Min; -- Bin := 1; -- loop -- if Interpolate_On and End_Points_Done then -- MathAK.Interpolate ( Completed, Freq, Rec, -- " " & Image ( Spikes ), Log_File, Temp_File ); -- exit when Completed; -- end if; -- if Print_Step then -- Print_Header; -- end if; -- Initialize ( Soma, V_Rest ); -- Time := 0.0; -- T_Print := 0.0; -- Spikes := 0; -- loop -- J_Stim := DC_Value -- + AC_Ampli * Sin ( 2.0 * MathAK.Pi * Freq * Time ); -- V_Old := Soma.V_Memb; -- Stimulate ( -- Soma, Time, T_Print, -- in out -- J_Stim, Print_Step, -- TPS, Time_Delta, CM, V_K, V_Na, V_Rest, -- in -- G_K_Mul, G_Na_Mul, J_L_Mul, J_L_Sub ); -- in -- if ( V_Old < Threshold ) and ( Soma.V_Memb >= Threshold ) then -- Spikes := Spikes + 1; -- SounAKD.Speaker_Thump; -- if Freq = 0.0 then -- Spikes_DC := Spikes_DC + 1; -- end if; -- end if; -- exit when Time > Time_Max; -- end loop; -- Put ( Freq, 2, 2, 4 ); -- Put ( " " ); -- Put ( Spikes ); -- Put ( " " ); -- if Freq /= 0.0 then -- Put ( float ( Spikes ) -- / ( Freq * Time_Max ), 2, 2, 4 ); -- Put ( " " ); -- Put ( float ( Spikes - Spikes_DC ) -- / ( Freq * Time_Max ), 2, 2, 4 ); -- end if; -- Put_Line ( "" ); -- if ( Log_File /= "" ) and not Interpolate_On then -- Put ( File, Freq, 4, 5, 4 ); -- Put ( File, " " ); -- Put ( File, Spikes ); -- Put ( File, "" ); -- if Freq /= 0.0 then -- Put ( File, float ( Spikes ) -- / ( Freq * Time_Max ), 2, 2, 4 ); -- Put ( File, " " ); -- Put ( File, float ( Spikes - Spikes_DC ) -- / ( Freq * Time_Max ), 2, 2, 4 ); -- end if; -- Put_Line ( File, "" ); -- end if; -- if Interpolate_On then -- if not End_Points_Done then -- Put ( File, Freq ); -- Put_Line ( File, " " & Image ( Spikes ) ); -- if Freq = Freq_Min then -- Freq := Freq_Max; -- else -- Close ( File ); -- End_Points_Done := true; -- end if; -- end if; -- else -- Freq := Freq + Freq_Step; -- Bin := Bin + 1; -- exit when Bin > Bins; -- end if; -- end loop; -- if ( Log_File /= "" ) and not Interpolate_On then -- Close ( File ); -- end if; -- exception -- when others => -- if Log_File /= "" then -- Close ( File ); -- end if; -- raise; -- end Bandpass_Filter; -- -- procedure Bandpass_Filter_Demo is -- ---------------------------------------------------------------------- -- AC_Ampli : float; -- DC_Value : float; -- Freq_Min : float; -- Freq_Max : float; -- Freq_Step : float; -- Threshold : float; -- Time_Max : float; -- Time_Delta : float; -- Print_Step : boolean; -- Log_File : string ( 1..12 ) := ( others => ASCII.Nul ); -- Bad_Delta : boolean; -- Interpolate_On : boolean; -- begin -- Freq_Min := Ask ( "Freq_Min " , Freq_Min_Default, -- -Freq_Max_Default, +Freq_Max_Default, 1, 1, 0 ); -- Freq_Max := Ask ( "Freq_Max " , Freq_Max_Default, -- Freq_Min, Freq_Max_Default, 1, 1, 0 ); -- Freq_Step := Ask ( "Freq_Step ", Freq_Step_Default, -- float'small, Freq_Max - Freq_Min, 1, 1, 0 ); -- AC_Ampli := Ask ( "AC_Ampli ", AC_Ampli_Default, -- Fore => 0, Aft => 0, Exp => 3 ); -- DC_Value := Ask ( "DC_Value ", DC_Value_Default, -- Fore => 0, Aft => 0, Exp => 3 ); -- Threshold := Ask ( "Threshold ", Threshold_Default, -- Fore => 0, Aft => 0, Exp => 3 ); -- Time_Max := Ask ( "Time_Max (seconds) ", Time_Max_Default, -- 0.0, float'large, 0, 0, 3 ); -- Time_Delta := Ask ( "Time_Delta (seconds) ", -- Time_Delta_Default, 0.0, Time_Max, 0, 0, 3 ); -- Print_Step := Ask_Bool ( "Show steps? ", true ); -- Ask ( Log_File, "Log_File: ", "" ); -- Interpolate_On := Ask_Bool ( "Interpolate_On? ", false ); -- loop -- Bad_Delta := false; -- begin -- Bandpass_Filter ( Freq_Min, Freq_Max, Freq_Step, AC_Ampli, -- DC_Value, Threshold, Time_Max, Time_Delta, Print_Step, -- Log_File => Log_File, Interpolate_On => Interpolate_On ); -- exception -- when constraint_error => -- Bad_Delta := true; -- Time_Delta := Time_Delta / 2.0; -- Put ( "Changing Time_Delta to" ); -- Put ( Time_Delta ); -- Put_Line ( "." ); -- delay 5.0; -- end; -- exit when not Bad_Delta; -- end loop; -- end Bandpass_Filter_Demo; procedure Calc_AB ( AN : out float; AM : out float; AH : out float; BN : out float; BM : out float; BH : out float; V_Memb : in float; V_Rest : in float ) is ---------------------------------------------------------------------- -- These equations are written with V_Memb and V_Rest in volts and with -- AN, AM, AH, BN, BM, BH in units of per second. Original HH -- equations were in milliVolts and in per millisecond. ---------------------------------------------------------------------- begin AN := 10.0 * ( -Tightness * ( V_Memb - V_Rest ) + 10.0 ) / ( Exp ( ( -Tightness * ( V_Memb - V_Rest ) + 10.0 ) / 10.0 ) - 1.0 ); AM := 100.0 * ( -Tightness * ( V_Memb - V_Rest ) + 25.0 ) / ( Exp ( ( -Tightness * ( V_Memb - V_Rest ) + 25.0 ) / 10.0 ) - 1.0 ); AH := 70.0 * Exp ( -Tightness * ( V_Memb - V_Rest ) / 20.0 ); BN := 125.0 * Exp ( -Tightness * ( V_Memb - V_Rest ) / 80.0 ); BM := 4000.0 * Exp ( -Tightness * ( V_Memb - V_Rest ) / 18.0 ); BH := Tightness / ( Exp ( ( -Tightness * ( V_Memb - V_Rest ) + 30.0 ) / 10.0 ) + 1.0 ); end Calc_AB; procedure Comparison ( Soma_Count : in positive := 1; Current_Mul : in float := J_Stim_Default ) is ---------------------------------------------------------------------- Time_Delta : constant float := Time_Delta_Default; V_Rest : constant float := V_Rest_Default; File_Name : constant string := "Compare.Dat"; Print_Time : constant float := Time_Delta * 100.0; Time_Max : constant float := 0.2; File : File_Type; Somas : array ( 1..Soma_Count ) of Soma_Type; Time : float := 0.0; J_Stim : float; Print_Elapse : float := Print_Time + Time_Delta; Print_On : boolean; T_Print : float; -- dummy begin Create ( File, Out_File, File_Name ); Put_Line ( "Sending output to file " & File_Name & "." ); for Soma in Somas'range loop Initialize ( Somas ( Soma ), V_Rest ); end loop; loop if Print_Elapse >= Print_Time then Print_On := true; Print_Elapse := 0.0; else Print_On := false; end if; if Print_On then Put ( Time ); Put ( " " ); Put ( File, Time ); Put ( File, " " ); end if; for Soma in Somas'range loop J_Stim := float ( Soma ) * Current_Mul; Stimulate ( Somas ( Soma ), Time, T_Print, -- in out J_Stim => J_Stim, Print_On => false, Time_Delta => Time_Delta ); if Print_On then Put ( Somas ( Soma ).V_Memb ); Put ( " " ); Put ( File, Somas ( Soma ).V_Memb ); Put ( File, " " ); end if; end loop; if Print_On then New_Line; New_Line ( File ); end if; exit when Time > Time_Max; Print_Elapse := Print_Elapse + Time_Delta; end loop; Close ( File ); end Comparison; procedure Demo is ---------------------------------------------------------------------- Default : natural := 1; Min : constant natural := 0; Max : constant natural := 10; Option : natural; begin loop Put_Line ( Copyright ); Put_Line ( "" ); Put_Line ( "Hodgkin-Huxley model for propagating nerve impulses." ); Put_Line ( "" ); Put_Line ( " 0 = Quit" ); Put_Line ( " 1 = Bandpass_Filter" ); Put_Line ( " 2 = Calc_AB" ); Put_Line ( " 3 = Fire" ); Put_Line ( " 4 = Soma_HH" ); Put_Line ( " 5 = Initialize" ); Put_Line ( " 6 = Print" ); Put_Line ( " 7 = Print_Header" ); Put_Line ( " 8 = Stimulate" ); Put_Line ( " 9 = Firing vs. Current (F-I) Curve" ); Put_Line ( "10 = Voltage vs. Current Injection Comparisons" ); New_Line; Option := Ask_Nat ( "Option ", Default, Min, Max ); Default := Option + 1; if Default > Max then Default := Min; end if; case Option is when 0 => exit; -- when 1 => Bandpass_Filter_Demo; -- when 2 => Calc_AB; when 3 => declare Time_Delta : float; Save_File : string ( 1..12 ); Time_Max : float := Time_Max_Default; begin Ask ( Save_File, "Save_File ", "HH_Fire.Dat " ); Time_Max := Ask ( "Time_Max (seconds) ", Time_Max, 0.0, Time_Max_Default * 10.0, 0, 0, 3 ); Time_Delta := Ask ( "Time_Delta (seconds) ", Time_Delta_Default, 0.0, Time_Max, 0, 0, 3 ); Fire ( Save_File, Time_Max, Time_Delta => Time_Delta ); end; -- when 4 => Soma_HH; -- when 5 => Initialize; -- when 6 => Print; when 7 => Print_Header; -- when 8 => Stimulate; when 9 => F_I_Curve; when 10 => Comparison ( Ask_Nat ( "Soma Count", 1 ), Ask ( "Current Multiplier", J_Stim_Default ) ); when others => Put_Line ( "The demonstration of that option is currently " & "not available." ); end case; Pause; New_Line; end loop; end Demo; procedure F_I_Curve ( -- TPS : in float := TPS_Default; Time_Delta : in float := Time_Delta_Default; CM : in float := CM_Default; V_K : in float := V_K_Default; V_Na : in float := V_Na_Default; V_Rest : in float := V_Rest_Default; G_K_Mul : in float := G_K_Mul_Default; G_Na_Mul : in float := G_Na_Mul_Default; J_L_Mul : in float := J_L_Mul_Default; J_L_Sub : in float := J_L_Sub_Default ) is ---------------------------------------------------------------------- Time_Max : constant float := 1.0; -- seconds T_Print : float := 0.0; Time : float; -- seconds Soma : Soma_Type; J_Stim : float; -- current of stimulus Spikes : natural; Firing : boolean; TPS : float := 0.0; File : File_Type; File_Name : constant string := "HH_FI.Dat"; begin Create ( File, Out_File, File_Name ); Put_Line ( "Sending output to file " & File_Name & "." ); for J_Stim_Mul in 1..100 loop Spikes := 0; Firing := false; Time := 0.0; T_Print := 0.0; Initialize ( Soma, V_Rest ); loop J_Stim := float ( J_Stim_Mul ) * 1.0e-7; Stimulate ( Soma, Time, T_Print, -- in out J_Stim, false, -- in TPS, Time_Delta, CM, V_K, V_Na, V_Rest, -- in G_K_Mul, G_Na_Mul, J_L_Mul, J_L_Sub ); -- in exit when Time > Time_Max; if Soma.V_Memb >= 0.0 then if not Firing then Spikes := Spikes + 1; SounAKD.Speaker_Thump; Firing := true; end if; else Firing := false; end if; end loop; Put ( J_Stim ); Put ( " " ); Put ( File, J_Stim ); Put ( File, " " ); Put ( Spikes ); New_Line; Put ( File, Spikes ); New_Line ( File ); end loop; Close ( File ); end F_I_Curve; procedure Fire ( Save_File : in string := ""; Time_Max : in float := Time_Max_Default; TPS : in float := TPS_Default; Time_Delta : in float := Time_Delta_Default; CM : in float := CM_Default; V_K : in float := V_K_Default; V_Na : in float := V_Na_Default; V_Rest : in float := V_Rest_Default; G_K_Mul : in float := G_K_Mul_Default; G_Na_Mul : in float := G_Na_Mul_Default; J_L_Mul : in float := J_L_Mul_Default; J_L_Sub : in float := J_L_Sub_Default ) is ---------------------------------------------------------------------- T_Print : float := 0.0; Time : float := 0.0; Soma : Soma_Type; J_Stim : float; -- current of stimulus X_Data, Y_Data : Vector_Type ( 0..79 ) := ( others => 0.0 ); Plot_Index : natural range X_Data'range := 0; Plot_Last : float := Time_Max / float ( X_Data'last ); begin Print_Header; Initialize ( Soma, V_Rest ); loop J_Stim := 0.0; if ( Time >= 0.5E-3 ) and ( Time < 0.6E-3 ) then J_Stim := J_Stim_Default; end if; Stimulate ( Soma, Time, T_Print, -- in out J_Stim, true, -- in TPS, Time_Delta, CM, V_K, V_Na, V_Rest, -- in G_K_Mul, G_Na_Mul, J_L_Mul, J_L_Sub ); -- in if Plot_Last >= Time_Max / float ( X_Data'last ) then X_Data ( Plot_Index ) := Time; Y_Data ( Plot_Index ) := Soma.V_Memb; exit when Plot_Index = X_Data'last; Plot_Index := Plot_Index + 1; Plot_Last := 0.0; end if; Plot_Last := Plot_Last + Time_Delta; end loop; Pause; Clear_Screen; Plot_XY ( X_Data, Y_Data ); if Save_File /= "" then Save ( Y_Data, Save_File ); end if; Move ( 24, 0 ); end Fire; procedure Initialize ( Soma : out Soma_Type; V_Rest : in float := V_Rest_Default ) is ---------------------------------------------------------------------- AN, AM, AH, BN, BM, BH : float; begin Soma.V_Memb := V_Rest; -- Calculate initial values of N, M, H -- N := AN / (AN + BN), etc. Calc_AB ( AN, AM, AH, BN, BM, BH, V_Rest, V_Rest ); Soma.N := AN / ( AN + BN ); Soma.M := AM / ( AM + BM ); Soma.H := AH / ( AH + BH ); end Initialize; procedure Print ( Time : in float; V_Memb : in float; Status : in Status_Type ) is ---------------------------------------------------------------------- begin Put ( Time , 2, 2, 0 ); Put ( " " ); Put ( V_Memb , 4, 1, 0 ); Put ( " " ); Put ( Status.J_Memb, 2, 2, 3 ); Put ( " " ); Put ( Status.G_Na , 2, 2, 3 ); Put ( " " ); Put ( Status.J_Na , 2, 2, 3 ); Put ( " " ); Put ( Status.G_K , 2, 2, 3 ); Put ( " " ); Put ( Status.J_K , 2, 2, 3 ); Put ( " " ); Put ( Status.J_L , 2, 2, 3 ); Put_Line ( "" ); end Print; procedure Print_Header is ---------------------------------------------------------------------- begin Put_Line ( " Time V_Memb J_Memb G_Na J_Na " & " G_K J_K J_L" ); Put_Line ( " (ms) (mV) (A cm-2) (Ohm-1 cm-2)" ); Put_Line ( "" ); end Print_Header; procedure Soma_HH ( Status : out Status_Type; Soma : in out Soma_Type; J_Stim : in float; CM : in float := CM_Default; Time_Delta : in float := Time_Delta_Default; V_K : in float := V_K_Default; V_Na : in float := V_Na_Default; V_Rest : in float := V_Rest_Default; G_K_Mul : in float := G_K_Mul_Default; G_Na_Mul : in float := G_Na_Mul_Default; J_L_Mul : in float := J_L_Mul_Default; J_L_Sub : in float := J_L_Sub_Default ) is ---------------------------------------------------------------------- dVdT : float; -- time-rate-change of Voltage dNdT, dMdT, dHdT : float; AN, AM, AH, BN, BM, BH : float; G_K_Temp : float; -- Potassium conductance G_Na_Temp : float; -- Sodium conductance J_K_Temp : float; -- current of Potassium J_Na_Temp : float; -- current of Sodium J_L_Temp : float; -- current of leakage J_Memb_Temp : float; -- current of membrane begin Calc_AB ( AN, AM, AH, BN, BM, BH, Soma.V_Memb, V_Rest ); -- Calculate derivatives using Eqs. 6.63, 6.67, and 6.68 dNdT := AN * ( 1.0 - Soma.N ) - BN * Soma.N; dMdT := AM * ( 1.0 - Soma.M ) - BM * Soma.M; dHdT := AH * ( 1.0 - Soma.H ) - BH * Soma.H; Soma.N := Soma.N + dNdT * Time_Delta; Soma.M := Soma.M + dMdT * Time_Delta; Soma.H := Soma.H + dHdT * Time_Delta; -- Conductances in Mho per sq. cm G_K_Temp := ( Soma.N ** 4 ) * G_K_Mul; G_Na_Temp := ( Soma.M ** 3 ) * G_Na_Mul * Soma.H; J_K_Temp := G_K_Temp * ( Soma.V_Memb - V_K ); J_Na_Temp := G_Na_Temp * ( Soma.V_Memb - V_Na ); J_L_Temp := J_L_Mul * ( Soma.V_Memb - V_Rest - J_L_Sub ); J_Memb_Temp := J_K_Temp + J_Na_Temp + J_L_Temp; -- Eq. 6.53 with no X dependence and with J_Stim current dVdT := ( -J_Memb_Temp + J_Stim ) / CM; Soma.V_Memb := Soma.V_Memb + dVdT * Time_Delta; Status.G_K := G_K_Temp; Status.G_Na := G_Na_Temp; Status.J_K := J_K_Temp; Status.J_Na := J_Na_Temp; Status.J_L := J_L_Temp; Status.J_Memb := J_Memb_Temp; end Soma_HH; procedure Stimulate ( Soma : in out Soma_Type; Time : in out float; T_Print : in out float; J_Stim : in float; Print_On : in boolean := true; TPS : in float := TPS_Default; Time_Delta : in float := Time_Delta_Default; CM : in float := CM_Default; V_K : in float := V_K_Default; V_Na : in float := V_Na_Default; V_Rest : in float := V_Rest_Default; G_K_Mul : in float := G_K_Mul_Default; G_Na_Mul : in float := G_Na_Mul_Default; J_L_Mul : in float := J_L_Mul_Default; J_L_Sub : in float := J_L_Sub_Default ) is ---------------------------------------------------------------------- Status : Status_Type; begin Soma_HH ( Status, -- out Soma, -- in out J_Stim, CM, Time_Delta, V_K, V_Na, V_Rest, -- in G_K_Mul, G_Na_Mul, J_L_Mul, J_L_Sub ); -- in if Print_On and ( Time >= T_Print ) then Print ( Time * 1000.0, Soma.V_Memb * 1000.0, Status ); T_Print := T_Print + TPS; end if; Time := Time + Time_Delta; end Stimulate; ---------------------------------------------------------------------- ---------------------------------------------------------------------- end HHNM;