Metabolic module: add Metabolic HH to neuronmodels

.gitignore

@@ -5,4 +5,4 @@ Manifest.toml
 .vscode
 neuroblox.code-workspace
 *.mat
-.DS_Store
+*.DS_Store

src/Neuroblox.jl

@@ -222,7 +222,7 @@ export JansenRitSPM12, next_generation, qif_neuron, if_neuron, hh_neuron_excitat
     hh_neuron_inhibitory, van_der_pol, Generic2dOscillator
 export HHNeuronExciBlox, HHNeuronInhibBlox, IFNeuron, LIFNeuron, QIFNeuron, IzhikevichNeuron, LIFExciNeuron, LIFInhNeuron,
     CanonicalMicroCircuitBlox, WinnerTakeAllBlox, CorticalBlox, SuperCortical, HHNeuronInhib_MSN_Adam_Blox, HHNeuronInhib_FSI_Adam_Blox, HHNeuronExci_STN_Adam_Blox,
-    HHNeuronInhib_GPe_Adam_Blox, Striatum_MSN_Adam, Striatum_FSI_Adam, GPe_Adam, STN_Adam, LIFExciCircuitBlox, LIFInhCircuitBlox
+    HHNeuronInhib_GPe_Adam_Blox, Striatum_MSN_Adam, Striatum_FSI_Adam, GPe_Adam, STN_Adam, LIFExciCircuitBlox, LIFInhCircuitBlox, MetabolicHHNeuron
 export LinearNeuralMass, HarmonicOscillator, JansenRit, WilsonCowan, LarterBreakspear, NextGenerationBlox, NextGenerationResolvedBlox, NextGenerationEIBlox, KuramotoOscillator
 export Matrisome, Striosome, Striatum, GPi, GPe, Thalamus, STN, TAN, SNc
 export HebbianPlasticity, HebbianModulationPlasticity

src/blox/connections.jl

@@ -960,6 +960,40 @@ function (bc::BloxConnector)(
     end
 end
 
+function (bc::BloxConnector)(
+    HH_out::Union{MetabolicHHNeuron},
+    HH_in::Union{MetabolicHHNeuron};
+    kwargs...
+)
+    """
+    Connection rule specific to the metabolic HH neurons based on Dutta et al:
+
+		Dutta, Shrey, et al. "Mechanisms underlying pathological cortical bursts
+		during metabolic depletion." Nature Communications 14.1 (2023): 4792.
+
+		https://www.nature.com/articles/s41467-023-40437-0
+		https://zenodo.org/records/8013692
+
+    """
+    sys_out = get_namespaced_sys(HH_out)
+    sys_in = get_namespaced_sys(HH_in)
+
+    w = generate_weight_param(HH_out, HH_in; kwargs...)
+    push!(bc.weights, w)
+
+    eq = if HH_out.neurontype == :excitatory
+        sys_in.I_syn ~ -w * sys_in.G_exc * (sys_in.V - sys_in.E_syn_exc) *
+            sys_out.S * exp(-sys_out.χ/5)
+    elseif HH_out.neurontype == :inhibitory
+        sys_in.I_syn ~ -w * sys_in.G_inh * (sys_in.V - sys_in.E_syn_inh) *
+            sys_out.S * exp(-sys_out.χ/5)
+    else
+        warning("Unknown neuron type. Assuming excitatory neuron type.")
+    end
+
+    accumulate_equation!(bc, eq)
+end
+
 function (bc::BloxConnector)(
     bloxout::PYR_Izh, 
     bloxin::PYR_Izh; 
@@ -993,4 +1027,4 @@ function (bc::BloxConnector)(
     eq = sys_in.jcn ~ w*x
 
     accumulate_equation!(bc, eq)
-end
+end

src/blox/neuron_models.jl

@@ -310,8 +310,6 @@ struct HHNeuronInhib_FSI_Adam_Blox <: AbstractInhNeuronBlox
 			a = a
 			b = b
 			T = T
-            τ = τ
-            τₛ = τₛ
 		end
 
         @brownian χ
@@ -335,6 +333,7 @@ struct HHNeuronInhib_FSI_Adam_Blox <: AbstractInhNeuronBlox
 			   D(hD)~(hD_inf(V)-hD)/τₕD(V),
 			   D(G)~(-1/τ)*G + G_asymp(V,a,b)*(1-G),
 			   D(Gₛ)~(-1/τₛ)*Gₛ + G_asymp(V,a,b)*(1-Gₛ)
+
 		]
 
 		sys = System(
@@ -863,3 +862,165 @@ struct IzhikevichNeuron <: AbstractNeuronBlox
 		new(p, sts[2], sts[5], sts[1], sys, namespace)
 	end
 end
+
+struct MetabolicHHNeuron <: AbstractNeuronBlox
+
+	"""
+
+	A Hodgkin-Huxley model expanded with:
+
+	- dynamic ion concentrations
+	- ATPase kinetic rate
+	- dynamic oxygen concentration
+	- astrocytic potassium buffering
+
+	Based on Dutta et al:
+
+		Dutta, Shrey, et al. "Mechanisms underlying pathological cortical bursts
+		during metabolic depletion." Nature Communications 14.1 (2023): 4792.
+
+		https://www.nature.com/articles/s41467-023-40437-0
+		https://zenodo.org/records/8013692
+
+	"""
+
+	odesystem
+    output
+    namespace
+	neurontype::Symbol
+
+	function MetabolicHHNeuron(
+		;name,
+		namespace=nothing,
+		neurontype=:excitatory,
+		Naᵢᵧ = 18.0,  # Intracellular Naconcentration, in mM
+		ρₘₐₓ = 1.25,  # Maximum pump rate, in mM/s
+		α = 5.3,  # Conversion factor from pump current to O2 consumption rate, in g/mol
+		λ = 1.,  # *Relative cell density
+		ϵ₀ = 0.17,  # O2 diffusion rate, in s^-1
+		O₂ᵦ = 32.,  # O2 buffer concentration, in mg/L #TODO: potentially unrealistic value (found values are ~0.5)
+		γ = 0.0445,  # conversion factor from current to concentration, in (mM/s)/(uA/cm2)
+		β = 7.,  # Ratio of intracellular vs extracellular volume
+		ϵₖ = 0.33,  # K+ diffusion rate, in 1/s
+		Kₒᵦ = 3.5,  # K+ buffer concentration, in mM
+		Gᵧ = 8.0,  # Glia uptake strength of K+, in mM/s
+		Clᵢ = 6.0, # Intracellular Cl- concentration, in mM
+		Clₒ = 130.0, # Extracellular Cl- concentration, in mM
+		R = 8.314,  # Ideal gas constant, in J/(mol*K)
+		T = 310.0,  # Temperature, in K
+		F = 96485.0,  # Faraday's constant, in C/mol
+		Gₙₐ = 30.,  # Na+ maximum conductance, in mS/cm^2
+		Gₖ = 25.,  # K+ maximum conductance, in mS/cm^2
+		Gₙₐ_L = 0.0175,  # Na+ leak conductance, in mS/cm^2
+		Gₖ_L = 0.05,  # K+ leak conductance, in mS/cm^2
+		G_cl_L = 0.05,  # Cl- leak conductance, in mS/cm^2
+		C_m = 1.,  # Membrane capacitance, in uF/cm^2
+		I_in = 0.,  # External current input, in uA/cm^2
+		G_exc = 0.022,  # Conductance of excitatory synapses, in mS/cm^2
+		G_inh = 0.374,  # Conductance of inhibitory synapses, in mS/cm^2
+		E_syn_exc = 0., # Excitatory synaptic reversal potential, in mV
+		E_syn_inh = -80.,  # Inhibitory synaptic reversal potential, in mV
+		τ = 4.,  # *Time constant for synapse, in ms
+		#*: significantly varies between excitatory and inhibitory neurons
+	)
+
+		# Parameters
+		ps = @parameters begin
+			Naᵢᵧ=Naᵢᵧ
+			ρₘₐₓ=ρₘₐₓ
+			α=α
+			λ=λ
+			ϵ₀=ϵ₀
+			O₂ᵦ=O₂ᵦ
+			γ=γ
+			β=β
+			ϵₖ=ϵₖ
+			Kₒᵦ=Kₒᵦ
+			Gᵧ=Gᵧ
+			R=R
+			T=T
+			F=F
+			Gₙₐ=Gₙₐ
+			Gₖ=Gₖ
+			Gₙₐ_L=Gₙₐ_L
+			Gₖ_L=Gₖ_L
+			G_cl_L=G_cl_L
+			C_m=C_m
+			I_in=I_in
+			G_exc=G_exc
+			G_inh=G_inh
+			E_syn_exc=E_syn_exc
+			E_syn_inh=E_syn_inh
+			τ=τ
+		end
+
+		# State variables
+		sts = @variables begin
+			V(t)=-60.0
+			O₂(t)=25.0
+			Kₒ(t)=3.0
+			Naᵢ(t)=15.0
+			m(t)=0.0
+			h(t)=0.0
+			n(t)=0.0
+			I_syn(t)=0.0 
+			[input=true] 
+			S(t)=0.1
+			[output=true] 
+			χ(t)=0.0
+			[output=true] 
+		end
+
+		# Pump currents
+		ρ = ρₘₐₓ / (1.0 + exp((20.0 - O₂)/3.0))
+		I_pump = ρ / (1.0 + exp((25.0 - Naᵢ)/3.0)*(1.0 + exp(5.5 - Kₒ)))
+		I_gliapump = ρ / (3.0*(1.0 + exp((25.0 - Naᵢᵧ)/3.0))*(1.0 + exp(5.5 - Kₒ)))
+
+		# Glia current
+		I_glia = Gᵧ / (1.0 + exp((18.0 - Kₒ)/2.5))
+
+		# Ion concentrations
+		Kᵢ = 140.0 + (18.0 - Naᵢ)
+		Naₒ = 144.0 - β*(Naᵢ - 18.0)
+
+		# Ion reversal potentials
+		Eₙₐ = R*T/F * log(Naₒ/Naᵢ) * 1000.0
+		Eₖ = R*T/F * log(Kₒ/Kᵢ) * 1000.0
+		E_cl = R*T/F * log(Clᵢ/Clₒ) * 1000.0
+
+		# Ion currents
+		Iₙₐ = Gₙₐ*m^3.0*h*(V - Eₙₐ) + Gₙₐ_L*(V - Eₙₐ)
+		Iₖ = Gₖ*n^4.0*(V - Eₖ) + Gₖ_L*(V - Eₖ)
+		I_cl = G_cl_L*(V - E_cl)
+
+		# Ion channel gating rate equations
+		aₘ = 0.32*(V + 54.0)/(1.0 - exp(-0.25*(V + 54.0)))
+		bₘ = 0.28*(V + 27.0)/(exp(0.2*(V + 27.0)) - 1.0)
+		aₕ = 0.128*exp(-(V + 50.0)/18.0)
+		bₕ = 4.0/(1.0 + exp(-0.2*(V + 27.0)))
+		aₙ = 0.032*(V + 52.0)/(1.0 - exp(-0.2*(V + 52.0)))
+		bₙ = 0.5*exp(-(V + 57.0)/40.0)
+
+		# Depolarization factor, as continuous variable
+		η = 0.4/(1.0 + exp(-10.0*(V + 30.0)))/(1.0 + exp(10.0*(V + 10.0)))
+
+		# Differential equations
+		eqs = [
+			D(O₂) ~ -α*λ*(I_pump + I_gliapump) + ϵ₀*(O₂ᵦ - O₂),
+			D(Kₒ) ~ γ*β*Iₖ - 2.0*β*I_pump - I_glia - 2.0*I_gliapump - ϵₖ*(Kₒ - Kₒᵦ),
+			D(Naᵢ) ~ -γ*Iₙₐ - 3.0*I_pump,
+			D(m) ~ aₘ * (1.0 - m) - bₘ*m,
+			D(h) ~ aₕ * (1.0 - h) - bₕ*h,
+			D(n) ~ aₙ * (1.0 - n) - bₙ*n,
+			D(V) ~ (-Iₙₐ - Iₖ - I_cl - I_syn - I_in)/C_m,
+			D(S) ~ (20.0/(1.0 + exp(-(V + 20.0)/3.0)) * (1.0 - S) - S)/τ,
+			D(χ) ~ η*(V + 50.0) - 0.4*χ
+			]
+
+		# Define the ODE system
+		sys = ODESystem(eqs, t, sts, ps; name=name)
+
+		# Construct the neuron
+		new(sys, sts[1], namespace, neurontype)
+	end
+end