diff --git a/Assets/Scripts/C2M2/CellData.cs b/Assets/Scripts/C2M2/CellData.cs index ba89a72d..f7923808 100644 --- a/Assets/Scripts/C2M2/CellData.cs +++ b/Assets/Scripts/C2M2/CellData.cs @@ -21,6 +21,10 @@ public class CellData public int refinementLevel; public double timeStep; public double endTime; + // Gating variable data + public GatingVariableData[] gates; + // Ion‐channel data + public ChannelData[] channels; // Clamp data [System.Serializable] @@ -35,14 +39,22 @@ public struct ClampData // Simulation state public double[] U; - public double[] M; - public double[] N; - public double[] H; - public double[] Upre; - public double[] Mpre; - public double[] Npre; - public double[] Hpre; + + // Gating Variables + [System.Serializable] + public class GatingVariableData { + public string name; + public double[] current; + public double[] previous; + } + + // Ion Channels + [System.Serializable] + public class ChannelData { + public string name; + public bool isActive; + } // Graph data [System.Serializable] diff --git a/Assets/Scripts/C2M2/GameManager.cs b/Assets/Scripts/C2M2/GameManager.cs index 8adc038f..b1aad9b3 100644 --- a/Assets/Scripts/C2M2/GameManager.cs +++ b/Assets/Scripts/C2M2/GameManager.cs @@ -24,14 +24,10 @@ public class GameManager : MonoBehaviour // vectors used in SparseSolve; for loading public double[] U; - public double[] M; - public double[] N; - public double[] H; - public double[] Upre; - public double[] Mpre; - public double[] Npre; - public double[] Hpre; + public Dictionary currentStates; + public Dictionary previousStates; + // for loading a file private bool loading = false; diff --git a/Assets/Scripts/C2M2/Menu.cs b/Assets/Scripts/C2M2/Menu.cs index 8abc1f26..0b58ebdf 100644 --- a/Assets/Scripts/C2M2/Menu.cs +++ b/Assets/Scripts/C2M2/Menu.cs @@ -2,6 +2,7 @@ using System.Collections; using System.Collections.Generic; using System.IO; +using System.Linq; using UnityEngine; using TMPro; @@ -139,14 +140,34 @@ public void Save() data = new CellData(); data.U = sim.Get1DValues(); // voltage at every node - data.M = sim.getM(); // M vector - data.N = sim.getN(); // N vector - data.H = sim.getH(); // H vector - data.Upre = sim.getUpre(); // Upre vector - data.Mpre = sim.getMpre(); // Mpre vector - data.Npre = sim.getNpre(); // Npre vector - data.Hpre = sim.getHpre(); // Hpre vector + + // Gating Variable Vectors + var curr = sim.getCurrentStates(); + var prev = sim.getPreviousStates(); + + var gatingVariableList = new List(); + foreach (var kvp in curr) + { + gatingVariableList.Add(new CellData.GatingVariableData { + name = kvp.Key, + current = kvp.Value, + previous = prev[kvp.Key], // lookup by the same key + }); + } + + data.gates = gatingVariableList.ToArray(); + + // Save active ion channels + var channelList = new List(); + foreach (var channel in sim.ionChannels) { + bool active = sim.activeIonChannels.Contains(channel); + channelList.Add(new CellData.ChannelData{ + name = channel.Name, + isActive = active + }); + } + data.channels = channelList.ToArray(); data.simID = sim.simID; data.pos = sim.transform.position; @@ -230,7 +251,7 @@ IEnumerator ShowSaveButton() SaveButtonVisible(true); } - /// +/// /// Load a file. Return true if successful. /// public bool Load(String f) @@ -277,6 +298,7 @@ public bool Load(String f) int ID = 0; // this will be the current ID when placing a new cell after loading + for (; i <= limit; i++) { // retrieve saved data @@ -288,16 +310,23 @@ public bool Load(String f) // restore vectors gm.U = data.U; - gm.M = data.M; - gm.N = data.N; - gm.H = data.H; - gm.Upre = data.Upre; - gm.Mpre = data.Mpre; - gm.Npre = data.Npre; - gm.Hpre = data.Hpre; + + // Build Gating Variables and Ion Channels + + var currStates = new Dictionary(); + var prevStates = new Dictionary(); + + foreach(var g in data.gates) { + currStates[g.name] = g.current; + prevStates[g.name] = g.previous; + } + + gm.currentStates = currStates; + gm.previousStates = prevStates; GameObject go; + try { go = loader.Load(new RaycastHit()); // load the cell @@ -313,9 +342,20 @@ public bool Load(String f) finishedLoading = true; return false; } - + SparseSolverTestv1 sim = go.GetComponent(); + // 1) initialize ion channels (so sim.ionChannels is populated) + sim.InitializeIonChannel(); + // 2) clear active channels then enable/disable each one according to what was saved + sim.activeIonChannels.Clear(); + // Load Ion Channels + foreach (var cd in data.channels) { + // find the matching IonChannel object + var match = sim.ionChannels.FirstOrDefault(ch => ch.Name == cd.name); + if (match != null && cd.isActive) sim.activeIonChannels.Add(match); + } + // restore cell ID sim.simID = data.simID; ID = sim.simID; @@ -352,12 +392,16 @@ public bool Load(String f) for (int j = 0; j < data.graphs.Length; j++) { var graphObj = Instantiate(graphPrefab); + NDLineGraph g = graphObj.GetComponent(); g.ndgraph.FocusVert = data.graphs[j].vertex; g.ndgraph.simulation = sim; graphM.graphs.Add(g.ndgraph); foreach (Vector3 v in data.graphs[j].positions) + { g.positions.Add(v); + + } } } @@ -385,7 +429,7 @@ public bool Load(String f) } syn = Instantiate(GameManager.instance.synapseManagerPrefab.GetComponent().synapsePrefab, ndsim.transform).GetComponentInChildren(); syn.AttachToSimulation(ndsim, synD.syns[j].synVert); - syn.SwitchModel(synD.syns[j].model); + // syn.SwitchModel(synD.syns[j].model); } finishedLoading = true; // this is for ChangeGradient diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs new file mode 100644 index 00000000..8671db30 --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs @@ -0,0 +1,58 @@ +using System; +using System.Collections.Generic; +using Vector = MathNet.Numerics.LinearAlgebra.Vector; +namespace C2M2.NeuronalDynamics.Simulation +{ + public class GatingVariable + { + // Name of the gating variable (i.e: n, m, h) + public string Name { get; } + + // Alpha and Beta functions for the gating variable + public Func Alpha { get; } + public Func Beta { get; } + + public double Exponent { get; } + + // Probability is the initial state probability for the gating variable + public double Probability { get; set; } + + // New flag: if true, this gating variable is computed instantaneously. Default: set to false. + public bool IsInstant { get; set; } + + public GatingVariable(string name, Func alpha, Func beta, double exponent, double probability, int nodeCount, bool isInstant = false) + { + Name = name; + Alpha = alpha; + Beta = beta; + Exponent = exponent; + Probability = probability; + IsInstant = isInstant; + } + } + + public class IonChannel + { + public string Name { get; set; } // Name of the ion channel (e.g., K+ channel) + public double Conductance { get; set; } // g (conductance) + public double ReversalPotential { get; set; } // E (reversal potential) + public List GatingVariables { get; set; } // List of gating variables (can be dynamic) + + public IonChannel(string name, double conductance, double reversalPotential) + { + Name = name; + Conductance = conductance; + ReversalPotential = reversalPotential; + GatingVariables = new List(); + } + public void AddGatingVariable(GatingVariable gatingVariable) + { + GatingVariables.Add(gatingVariable); + } + + public void RemoveGatingVariable(GatingVariable gatingVariable) + { + GatingVariables.Remove(gatingVariable); + } + } +} diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs.meta b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs.meta new file mode 100644 index 00000000..920e86d0 --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannel.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 110e4eedd8e3ead4fadd2d1421eb7005 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs new file mode 100644 index 00000000..f6d3c846 --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs @@ -0,0 +1,445 @@ +using System; +using Vector = MathNet.Numerics.LinearAlgebra.Vector; +using MathNet.Numerics.LinearAlgebra; +using UnityEngine; + +namespace C2M2.NeuronalDynamics.Simulation +{ + /// + /// This file consists of the channels for each Ion Channel. The equations follow the class architecture listed in IonChannels.cs + /// + public static class IonChannelModels + { + + /// + /// Compute steady-state initial value for a gating variable given + /// its alpha and beta rate functions and a reference voltage. + /// The alpha/beta functions accept a Vector and return a Vector; + /// this helper evaluates them at a constant voltage and returns + /// the scalar steady-state value a/(a+b). If alpha or beta is + /// null or a+b is zero/NaN/Inf, a safe fallback is returned. + /// + private static double InitializeGate( + Func alpha, + Func beta, + int nodeCount, + double voltage) + { + // Guard: if alpha is null we cannot initialize; return 0.0 + if (alpha == null) return 0.0; + + // Build a vector of length nodeCount with the requested voltage + var v = Vector.Build.Dense(nodeCount, voltage); + + double a = alpha(v)[0]; + double b = 0.0; + if (beta != null) b = beta(v)[0]; + + return a / (a + b); + } + + /// + /// Potassium Channel + /// gK = 5.0e1, eK = -90e-3 + /// Rate equations taken from Pospischil et al., 2008 + /// + public static IonChannel PotassiumChannel(int nodeCount, double restingV) + { + /// + /// [S/m2] potassium conductance per unit area, this is the Potassium conductance per unit area, it is used in this term + /// ḡKd n^4 (V − ek) + /// where n is the state variable, and ek is the reversal potential. + /// + double gk = 5.0 * 1.0E1; + /// + /// [V] potassium reversal potential + /// + double ek = -90.0 * 1.0E-3; + /// + /// [V] voltage threshold + /// + double vT = 0.0 * 1.0E-3; + /// Declare new ion channel, "potassiumChannel" + IonChannel potassiumChannel = new IonChannel("Potassium Channel", gk, ek); + /// + /// This is \f$\alpha_n\f$ rate function, the rate functions take the form of + /// αn = −0.032(V − VT −15) / exp[−(V − VT − 15)/5] − 1 + /// + /// this is the input voltage + /// alpha_n this function returns the rate at the given voltage + Func alpha_n = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.032) * (15.0 + vT - Vin).PointwiseDivide(((15.0 + vT - Vin) / 5.0).PointwiseExp() - 1.0); + }; + /// + /// This is \f$\beta_n\f$ rate function, the rate functions take the form of + /// βn = 0.5 exp[−(V − VT − 10)/40], + /// + /// this is the input voltage + /// beta_n this function returns the rate at the given voltage + Func beta_n = voltage => + { + var Vin = voltage.Clone(); + // Convert voltage from [V] to [mV] + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.5) * ((10.0 + vT - Vin) / 40.0).PointwiseExp(); + }; + + double n_init = InitializeGate(alpha_n, beta_n, nodeCount, restingV); // n initial probability + + /// Format for returning gating variables + /// new GatingVariable("variable name", alpha_function, beta_function, exponenet, initial probability, nodeCount) + potassiumChannel.AddGatingVariable( + new GatingVariable("n", alpha_n, beta_n, 4, n_init, nodeCount) + ); + return potassiumChannel; + } + + /// + /// Sodium Channel + /// gNa = 50.0e1, eNa = 50.0e-3 + /// Rate equations taken from Pospischil et al., 2008 + /// + public static IonChannel SodiumChannel(int nodeCount, double restingV) + { + /// + /// [S/m2] sodium conductance per unit area, this is the Sodium conductance per unit area, it is used in this term + /// INa = ḡNa m^3h (V − enaa) + /// where \f$m,h\f$ are the state variables, and \f$V_{Na}\f$ is the reversal potential for sodium. + /// + double gna = 50.0 * 1.0E1; + /// + /// [V] sodium reversal potential + /// + double ena = 50.0 * 1.0E-3; + /// + /// [V] voltage threshold + /// + double vT = 0.0 * 1.0E-3; + /// Declaration of new ion channel sodiumChannel. + /// Declared with the struct: (Name, Coductance, Reversal) + IonChannel sodiumChannel = new IonChannel("Sodium Channel", gna, ena); + /// + /// This is \f$\alpha_m\f$ rate function, the rate functions take the form of + /// αm = −0.32(V − VT −13) / exp[−(V − VT − 13)/4] − 1 + /// + /// this is the input voltage + /// alpha_m this function returns the rate at the given voltage + Func alpha_m = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.32) * (13.0 + vT - Vin).PointwiseDivide(((13.0 + vT - Vin) / 4.0).PointwiseExp() - 1.0); + }; + /// + /// This is \f$\beta_m\f$ rate function, the rate functions take the form of + /// βm = 0.28(V − VT −40) / exp[(V − VT − 40)/5] − 1 + /// + /// this is the input voltage + /// beta_m this function returns the rate at the given voltage + Func beta_m = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.28) * (Vin - vT - 40.0).PointwiseDivide(((Vin - vT - 40.0) / 5.0).PointwiseExp() - 1.0); + }; + /// + /// This is \f$\alpha_h\f$ rate function, the rate functions take the form of + /// αh = 0.128 exp[−(V − VT − 17)/18] + /// + /// this is the input voltage + /// alpha_h this function returns the rate at the given voltage + Func alpha_h = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.128) * ((17.0 + vT - Vin) / 18.0).PointwiseExp(); + }; + /// + /// This is \f$\beta_h\f$ rate function, the rate functions take the form of + /// βh=4 / 1 + exp[−(V − VT − 40)/5] + /// + /// this is the input voltage + /// beta_h this function returns the rate at the given voltage + Func beta_h = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * 4.0 / (((40.0 + vT - Vin) / 5.0).PointwiseExp() + 1.0); + }; + + + double h_init = InitializeGate(alpha_h, beta_h, nodeCount, restingV); // h initial probability + double m_init = InitializeGate(alpha_m, beta_m, nodeCount, restingV); // m initial probability + + /// Format for returning gating variables + /// new GatingVariable("variable name", alpha_function, beta_function, exponenet, initial probability, nodeCount) + sodiumChannel.AddGatingVariable( + new GatingVariable("m", alpha_m, beta_m, 3, m_init, nodeCount) + ); + sodiumChannel.AddGatingVariable( + new GatingVariable("h", alpha_h, beta_h, 1, h_init, nodeCount) + ); + + return sodiumChannel; + } + + /// + /// Leakage Channel + /// + + public static IonChannel LeakageChannel(int nodeCount, double restingV) + { + /// + /// [S/m2] leak conductance per unit area, this is the leak conductance per unit area, it is used in this term + /// \f[\bar{g}_{l}(V-V_l)\f] + /// + double gl = 1e-4; // S/m² + /// + /// [V] leak reversal potential + /// + double el = -70.0 * 1.0E-3; + return new IonChannel("Leakage Channel", gl, el); + } + + + /// + /// Calcium Channel + /// gCa = 0.0001e4, eCa = 120.0e-3 + /// Rate equations taken from Pospischil et al., 2008 + /// Conductance taken from Fig. 5, Model of intrinsically bursting cell + /// + public static IonChannel CalciumChannel(int nodeCount, double restingV) + { + double gca = 1.0; // S/m^2 + /// + /// [V] calcium reversal potential + /// + double eca = 120.0 * 1.0E-3; + + IonChannel calciumChannel = new IonChannel("Calcium Channel", gca, eca); + /// + /// This is \f$\alpha_q\f$ rate function, the rate functions take the form of + /// αq = 0.055(−27 − V) / exp[(−27 − V )/3.8] − 1 + /// + /// this is the input voltage + /// alpha_q this function returns the rate at the given voltage + Func alpha_q = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.055) * (-27.0 - Vin) / (((-27.0 - Vin) / 3.8).PointwiseExp() - 1.0); + }; + /// + /// This is \f$\beta_q\f$ rate function, the rate functions take the form of + /// βq = 0.94 exp[(−75 − V) / 17] + /// + /// this is the input voltage + /// beta_q this function returns the rate at the given voltage + Func beta_q = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.94) * (((-75.0 - Vin) / 17.0).PointwiseExp()); + }; + /// + /// This is \f$\alpha_r\f$ rate function, the rate functions take the form of + /// αr = 0.000457 exp[(−13 − V )/50] + /// + /// this is the input voltage + /// alpha_r this function returns the rate at the given voltage + Func alpha_r = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.000457) * (((-13.0 - Vin) / 50.0).PointwiseExp()); + }; + /// + /// This is \f$\beta_r\f$ rate function, the rate functions take the form of + /// βr = 0.0065 / exp[(−15 − V )/28] + 1 + /// + /// this is the input voltage + /// beta_r this function returns the rate at the given voltage + Func beta_r = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return (1.0E3) * (0.0065) / (((-15.0 - Vin) / 28.0).PointwiseExp() + 1.0); + }; + + double q_init = InitializeGate(alpha_q, beta_q, nodeCount, restingV); // q initial probability + double r_init = InitializeGate(alpha_r, beta_r, nodeCount, restingV); // r initial probability + + /// Format for returning gating variables + /// new GatingVariable("variable name", alpha_function, beta_function, exponenet, initial probability, nodeCount) + calciumChannel.AddGatingVariable( + new GatingVariable("q", alpha_q, beta_q, 2, q_init, nodeCount) + ); + calciumChannel.AddGatingVariable( + new GatingVariable("r", alpha_r, beta_r, 1, r_init, nodeCount) + ); + + return calciumChannel; + } + + /// + /// Slow Potassium Channel + /// p∞(V) = 1 / [1 + exp(-(V+35)/10)] + /// τp(V) = τmax / [3.3 * exp((V+35)/20) + exp(-(V+35)/20)] + /// dp/dt = ( p∞(V) - p ) / τp(V) + /// We use alpha_p(V) = p∞(V/)τp(V), beta_p(V) = [1 - p∞(V)]/τp(V). + /// Default gM = 0.004 mS/cm² = 0.004 * 10 = 0.04 S/m², τmax = 4.0 s + /// + public static IonChannel SlowPotassiumChannel(int nodeCount, double restingV) + { + // [S/m²] slow sotassium conductance + double gM = 7.5E-5 * 1.0E4; + /// + /// [V] slow potassium reversal potential + /// + double eK = -90.0 * 1.0E-3; + /// + /// [S] maximum time constant, tmax, for the p-gate. + /// + double tMax = 4; + + IonChannel slowKChannel = new IonChannel("Slow Potassium Channel", gM, eK); + + // alpha_p(V) = p∞(V) / τp(V) + Func alpha_p = voltage => + { + // Convert voltage from [V] to [mV] + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + + // Compute steady-state value p∞(V) + Vector pInf = 1.0 / (1.0 + (-(Vin + 35) / 10).PointwiseExp()); + + // Compute time constant τp(V) + Vector tauP = tMax / (3.3 * ((Vin + 35) / 20).PointwiseExp() + (-(Vin + 35) / 20).PointwiseExp()); + + // Return alpha_p(V) = p∞(V) / τp(V) + return pInf.PointwiseDivide(tauP); + }; + + // Define beta_p(V) = [1 - p∞(V)] / τp(V) + Func beta_p = voltage => + { + // Convert voltage from [V] to [mV] + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + + // Compute steady-state value p∞(V) + Vector pInf = 1.0 / (1.0 + (-(Vin + 35) / 10).PointwiseExp()); + + // Compute time constant τp(V) + Vector tauP = tMax / (3.3 * ((Vin + 35) / 20).PointwiseExp() + (-(Vin + 35) / 20).PointwiseExp()); + + // Return beta_p(V) = [1 - p∞(V)] / τp(V) + return (1 - pInf).PointwiseDivide(tauP); + }; + + // Based on the context of voltage-gated calcium channel kinetics, + // tau_p represents the time constant of inactivation, + // while p_inf represents the steady-state open probability + // (or inactivation) of the channel.  + + double p_init = 1.0 / (1.0 + System.Math.Exp(-(restingV + 35) / 10)); + + // Add the gating variable 'p' with exponent = 1 and initial probability 0.0 + slowKChannel.AddGatingVariable( + // new GatingVariable("p", alpha_p, beta_p, 1, p_init, nodeCount) + new GatingVariable("p", alpha_p, beta_p, 1, p_init, nodeCount) + ); + + return slowKChannel; + } + + /// + /// Low Threshold Calcium Channel (IT): + /// + /// IT = gT * [ s∞(V) ]^2 * u * (V - Eca) + /// + /// s∞(V) = 1 / [1 + exp(-(V + Vx + 57)/6.2)] + /// u∞(V) = 1 / [1 + exp((V + Vx + 81)/4)] + /// τu(V) = [ 30.8 + 211.4 + exp((V + Vx + 113.2)/5) ] + /// / [ 3.7 * (1 + exp((V + Vx + 84)/3.2)) ] + /// + /// The activation s∞(V) is "instantaneous," so we do not create + /// a separate gating variable for s. Instead, we define only + /// the inactivation gating variable "u" with alpha_u/beta_u + /// so that du/dt = [u∞(V) - u]/τu(V). + /// + /// Vx is a uniform shift of the voltage dependence, e.g. +2 mV. + /// + public static IonChannel LowThresholdCalciumChannel(int nodeCount, double restingV) + { + // [S/m²] Low Threshold Calcium Conductance + double gT = 0.0004 * 1.0E4; + // Reversal potential for Ca²⁺ in volts + double eCa = 120.0 * 1.0E-3; + // Voltage shift (in mV) + double Vx = 2.0; + + // Create the IonChannel object with name, conductance, and reversal potential + IonChannel lowTCalciumChannel = new IonChannel("Low Threshold Calcium Channel", gT, eCa); + + Func alpha_u = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); // Convert voltage from [V] to [mV] + + // Compute u∞(V) + Vector uInf = 1.0 / (1.0 + ((Vin + Vx + 81.0) / 4.0).PointwiseExp()); + // Compute τu(V) + Vector tauU = (30.8 + 211.4 + ((Vin + Vx + 113.2) / 5.0).PointwiseExp()) + .PointwiseDivide(3.7 * (1.0 + ((Vin + Vx + 84.0) / 3.2).PointwiseExp())); + // Return alpha_u(V) = u∞(V) / τu(V) + return uInf.PointwiseDivide(tauU); + }; + + Func beta_u = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + + Vector uInf = 1.0 / (1.0 + ((Vin + Vx + 81.0) / 4.0).PointwiseExp()); + Vector tauU = (30.8 + 211.4 + ((Vin + Vx + 113.2) / 5.0).PointwiseExp()) + .PointwiseDivide(3.7 * (1.0 + ((Vin + Vx + 84.0) / 3.2).PointwiseExp())); + // Return beta_u(V) = [1 - u∞(V)] / τu(V) + return (1.0 - uInf).PointwiseDivide(tauU); + }; + + + double u_init = 1.0 / (1.0 + System.Math.Exp((restingV + Vx + 81.0) / 4.0)); // u initial probability + + // Define an instantaneous gating variable "s" that represents s∞(V). + // For example, we assume: + // s∞(V) = 1 / (1 + exp(-(Vin + Vx + 50)/10)) + // This function is computed directly from voltage. + Func sInf = voltage => + { + var Vin = voltage.Clone(); + Vin.Multiply(1.0E3, Vin); + return 1.0 / (1.0 + (-(Vin + Vx + 57.0) / 6.2).PointwiseExp()); + }; + + // Add the gating variable "u" with exponent 1 and initial probability 0.0 + lowTCalciumChannel.AddGatingVariable( + new GatingVariable("u", alpha_u, beta_u, 1, u_init, nodeCount) + ); + + // Create the instantaneous gating variable "s" with exponent 1. + // The instantaneous variable is set via "true" as the last parameter passed + // Initialize instantaneous gating variable 's' to sInf(restingV) + double sInit = sInf(Vector.Build.Dense(nodeCount, restingV))[0]; + lowTCalciumChannel.AddGatingVariable( + new GatingVariable("s", sInf, null, 1, 0.0, nodeCount, true) + ); + + return lowTCalciumChannel; + } + } +} \ No newline at end of file diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs.meta b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs.meta new file mode 100644 index 00000000..9f1559ac --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelModels.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 27703378dfd18224395aaa132631d991 +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs new file mode 100644 index 00000000..031cec2f --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs @@ -0,0 +1,187 @@ +using System.Collections.Generic; +using UnityEngine; +using UnityEngine.UI; +using TMPro; +using C2M2.Interaction; +using C2M2.NeuronalDynamics.Simulation; + +namespace C2M2.NeuronalDynamics.IonChannels +{ + public class IonChannelPanelConfig : MonoBehaviour + { + public Transform channelListParent; + + private static readonly Color ActiveColor = new Color(0.18f, 0.80f, 0.18f); + private static readonly Color InactiveColor = new Color(0.85f, 0.15f, 0.15f); + private static Sprite circleSprite; + + private const float ButtonSize = 18f; + private const float RowHeight = 38f; + private const float RowSpacing = 6f; + + private SparseSolverTestv1 solver; + private readonly List<(IonChannel channel, Image buttonImage)> rows = new List<(IonChannel, Image)>(); + + private void OnEnable() + { + if (GameManager.instance?.activeSims == null || GameManager.instance.activeSims.Count == 0) + { + Debug.LogWarning("IonChannelConfig: no SparseSolverTestv1 found in activeSims."); + return; + } + + solver = GameManager.instance.activeSims[0] as SparseSolverTestv1; + + if (solver == null || solver.ionChannels == null || solver.ionChannels.Count == 0) + { + Debug.LogWarning("IonChannelConfig: ionChannels not yet initialized."); + return; + } + + channelListParent = null; + EnsureListParent(); + PopulateChannels(); + } + + private void EnsureListParent() + { + if (channelListParent != null) return; + + Transform bg = transform.Find("Background"); + if (bg == null) + { + Debug.LogError("IonChannelConfig: channelListParent not assigned and no 'Background' child found."); + return; + } + + Transform existing = bg.Find("ChannelList"); + if (existing != null) + { + channelListParent = existing; + return; + } + + GameObject listGO = new GameObject("ChannelList"); + listGO.transform.SetParent(bg, false); + + RectTransform listRect = listGO.AddComponent(); + listRect.anchorMin = new Vector2(0.5f, 1f); + listRect.anchorMax = new Vector2(0.5f, 1f); + listRect.pivot = new Vector2(0.5f, 1f); + listRect.anchoredPosition = new Vector2(-5f, -45f); + listRect.sizeDelta = new Vector2(270f, 280f); + + VerticalLayoutGroup vlg = listGO.AddComponent(); + vlg.spacing = RowSpacing; + vlg.childAlignment = TextAnchor.UpperLeft; + vlg.childForceExpandWidth = true; + vlg.childForceExpandHeight = false; + vlg.padding = new RectOffset(8, 8, 4, 4); + + channelListParent = listGO.transform; + } + + private void Update() + { + if (solver == null) + { + if (GameManager.instance?.activeSims == null || GameManager.instance.activeSims.Count == 0) return; + solver = GameManager.instance.activeSims[0] as SparseSolverTestv1; + if (solver == null || solver.ionChannels == null || solver.ionChannels.Count == 0) return; + channelListParent = null; + EnsureListParent(); + PopulateChannels(); + return; + } + foreach (var (channel, btnImage) in rows) + btnImage.color = solver.activeIonChannels.Contains(channel) ? ActiveColor : InactiveColor; + } + + private void PopulateChannels() + { + foreach (Transform child in channelListParent) + Destroy(child.gameObject); + rows.Clear(); + + if (circleSprite == null) circleSprite = MakeCircleSprite(64); + + foreach (IonChannel channel in solver.ionChannels) + CreateRow(channel); + } + + private void CreateRow(IonChannel channel) + { + GameObject rowGO = new GameObject(channel.Name + "_Row"); + rowGO.transform.SetParent(channelListParent, false); + rowGO.AddComponent().sizeDelta = new Vector2(0f, RowHeight); + + HorizontalLayoutGroup hlg = rowGO.AddComponent(); + hlg.spacing = 6f; + hlg.childAlignment = TextAnchor.MiddleLeft; + hlg.childForceExpandWidth = false; + hlg.childForceExpandHeight = false; + hlg.padding = new RectOffset(2, 2, 2, 2); + + GameObject btnGO = new GameObject("ToggleBtn"); + btnGO.transform.SetParent(rowGO.transform, false); + btnGO.AddComponent().sizeDelta = new Vector2(ButtonSize, ButtonSize); + + LayoutElement btnLayout = btnGO.AddComponent(); + btnLayout.minWidth = ButtonSize; + btnLayout.minHeight = ButtonSize; + btnLayout.preferredWidth = ButtonSize; + btnLayout.preferredHeight = ButtonSize; + + Image btnImage = btnGO.AddComponent(); + btnImage.sprite = circleSprite; + btnImage.color = solver.activeIonChannels.Contains(channel) ? ActiveColor : InactiveColor; + + btnGO.layer = LayerMask.NameToLayer("Raycast"); + + BoxCollider col = btnGO.AddComponent(); + col.size = new Vector3(ButtonSize, ButtonSize, 1f); + + RaycastPressEvents pressEvents = btnGO.AddComponent(); + btnGO.AddComponent().LRTrigger = pressEvents; + pressEvents.OnPress.AddListener(_ => OnToggleClicked(channel, btnImage)); + + GameObject labelGO = new GameObject("Label"); + labelGO.transform.SetParent(rowGO.transform, false); + + TextMeshProUGUI label = labelGO.AddComponent(); + label.text = channel.Name; + label.fontSize = 20f; + label.color = Color.white; + label.alignment = TextAlignmentOptions.MidlineLeft; + + rows.Add((channel, btnImage)); + } + + private static Sprite MakeCircleSprite(int diameter) + { + Texture2D tex = new Texture2D(diameter, diameter, TextureFormat.ARGB32, false); + Color[] pixels = new Color[diameter * diameter]; + float r = diameter / 2f; + Vector2 center = new Vector2(r - 0.5f, r - 0.5f); + for (int i = 0; i < pixels.Length; i++) + { + float x = i % diameter; + float y = i / diameter; + pixels[i] = Vector2.Distance(new Vector2(x, y), center) <= r ? Color.white : Color.clear; + } + tex.SetPixels(pixels); + tex.Apply(); + return Sprite.Create(tex, new Rect(0, 0, diameter, diameter), new Vector2(0.5f, 0.5f)); + } + + private void OnToggleClicked(IonChannel channel, Image buttonImage) + { + if (solver == null) return; + bool newState = !solver.activeIonChannels.Contains(channel); + solver.ToggleChannel(channel, newState); + if (newState) SparseSolverTestv1.SavedActiveChannelNames.Add(channel.Name); + else SparseSolverTestv1.SavedActiveChannelNames.Remove(channel.Name); + buttonImage.color = newState ? ActiveColor : InactiveColor; + } + } +} diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs.meta b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs.meta new file mode 100644 index 00000000..d0e795ae --- /dev/null +++ b/Assets/Scripts/C2M2/NeuronalDynamics/IonChannels/IonChannelPanelConfig.cs.meta @@ -0,0 +1,11 @@ +fileFormatVersion: 2 +guid: 9695a182ed19dbf4ba2f5d563d80507d +MonoImporter: + externalObjects: {} + serializedVersion: 2 + defaultReferences: [] + executionOrder: 0 + icon: {instanceID: 0} + userData: + assetBundleName: + assetBundleVariant: diff --git a/Assets/Scripts/C2M2/NeuronalDynamics/Simulation/SparseSolverTestv1.cs b/Assets/Scripts/C2M2/NeuronalDynamics/Simulation/SparseSolverTestv1.cs index ac609608..c1d3c204 100644 --- a/Assets/Scripts/C2M2/NeuronalDynamics/Simulation/SparseSolverTestv1.cs +++ b/Assets/Scripts/C2M2/NeuronalDynamics/Simulation/SparseSolverTestv1.cs @@ -51,8 +51,13 @@ public class SparseSolverTestv1 : NDSimulation /// /// This is the voltage for the voltage clamp, this is primarily used for when we do the convergence analysis of the code using a /// soma clamp at 50 [mV], the units for voltage in the solver is [V] that is why vstart is set to 0.05 + /// I REALLY THINK THIS IS NOT THE VOLTAGE FOR THE CLAMP + public double vstart = 0.050; + /// + /// This is the starting voltage of the cells. All indices of U (Voltage) are intialized to the startingVoltage quantity. + /// -0.05 [V] equates to -50 mV. /// - public double vstart = 0.050; + public double v0 = -0.07; /// /// [ohm.m] resistance.length, this is the axial resistence of the neuron, increasing this value has the effect of making the AP waves more localized and slower conduction speed /// decreasing this value has the effect of make the AP waves larger and have a faster conduction speed @@ -63,69 +68,43 @@ public class SparseSolverTestv1 : NDSimulation /// private double cap = 1.0 * 1.0E-2; /// - /// [S/m2] potassium conductance per unit area, this is the Potassium conductance per unit area, it is used in this term - /// \f[\bar{g}_{K}n^4(V-V_k)\f] - /// where \f$n\f$ is the state variable, and \f$V_k\f$ is the reversal potential. - /// - private double gk = 6.0 * 1.0E1; - // private double gk = 36; - /// - /// [S/m2] sodium conductance per unit area, this is the Sodium conductance per unit area, it is used in this term - /// \f[\bar{g}_{Na}m^3h(V-V_{Na})\f] - /// where \f$m,h\f$ are the state variables, and \f$V_{Na}\f$ is the reversal potential for sodium. - /// - private double gna = 56.0 * 1.0E1; - - /// - /// [S/m2] leak conductance per unit area, this is the leak conductance per unit area, it is used in this term - /// \f[\bar{g}_{l}(V-V_l)\f] - /// \f$V_l\f$ is the leak reversal potential. + /// leakConductance is used when setting the target time step. This value is updated in InitializeIonChannel /// - private double gl = 0.1 * 1.0E1; + private double leakConductance = 0.0; /// - /// [V] potassium reversal potential - /// - private double ek = -90.0 * 1.0E-3; - - /// - /// [V] sodium reversal potential + /// These are the solution vectors for the voltage U /// - private double ena = 50.0 * 1.0E-3; - + private Vector U; /// - /// [V] leak reversal potential + /// This is the U that gets modified during the step before U is set to it. /// - private double el = -70.0 * 1.0E-3; + private Vector U_Active; /// - /// [V] threshold voltage + /// this is for storing gating variables previous states /// - private static double Vt = -50.0 * 1.0E-3; /// - /// [] potassium channel state probability, unitless + /// this is for storing previous states /// - private double ni = 0.0009648121738618698; + private Vector Upre; + private Dictionary currentStates; /// - /// [] sodium channel state probability, unitless + /// this is for storing gating variable previous states /// - private double mi = 0.00016423459289788037; + private Dictionary previousStates; /// - /// [] sodium channel state probability, unitless + /// Declaration for the list of IonChannels /// - private double hi = 0.999975419740687; + public static HashSet SavedActiveChannelNames; + public List ionChannels; /// - /// These are the solution vectors for the voltage U - /// the state M, state N, and state H + /// Declaration for the list of active IonChannels /// - private Vector U; + public List activeIonChannels; /// - /// This is the U that gets modified during the step before U is set to it. + /// Total conductance of all channels + /// Used in SetTargetTimestep + public double totalConductance = 0; /// - private Vector U_Active; - /// - /// These are the solution vectors for the voltage U - /// the state M, state N, and state H - /// - private Vector M, N, H; /// /// This is for the synaptic current. It contains: /// [0]: The current at the active time step. @@ -137,10 +116,6 @@ public class SparseSolverTestv1 : NDSimulation /// private Vector surfaceArea; /// - /// this is for storing previous states - /// - private Vector Upre, Npre, Mpre, Hpre; - /// /// This is a vector the Reaction terms /// private Vector R; @@ -152,8 +127,6 @@ public class SparseSolverTestv1 : NDSimulation /// Temporary state vector /// private Vector tempState; - - List reactConst; //This is for passing the reaction function constants List> sparse_stencils; CompressedColumnStorage r_csc; //This is for the rhs sparse matrix CompressedColumnStorage l_csc; //This is for the lhs sparse matrix @@ -316,19 +289,16 @@ protected override void PreSolve() GameManager g = GameManager.instance; // if loading, the values from file will be set in BuildVectors and Set1DValues if (!g.Loading) InitializeNeuronCell(); - else BuildVectors(g.U, g.M, g.N, g.H, g.Upre, g.Mpre, g.Npre, g.Hpre); + else BuildVectors(g.U, g.Upre, g.currentStates, g.previousStates); ///R this is the reaction vector for the reaction solve R = Vector.Build.Dense(Neuron.nodes.Count); tempState = Vector.Build.Dense(Neuron.nodes.Count, 0); - ///reactConst this is a small list for collecting the conductances and reversal potential which is sent to the reaction solve routine - reactConst = new List { gk, gna, gl, ek, ena, el }; /// this sets the target time step size - timeStep = SetTargetTimeStep(cap, 2 * Neuron.MaxRadius,2*Neuron.MinRadius, Neuron.TargetEdgeLength, gna, gk,gl, res, 1.0); + timeStep = SetTargetTimeStep(cap, 2 * Neuron.MaxRadius,2*Neuron.MinRadius, Neuron.TargetEdgeLength, totalConductance, leakConductance, res, 1.0); ///UnityEngine.Debug.Log("Target Time Step = " + timeStep); - ///List> sparse_stencils = makeSparseStencils(Neuron, res, cap, k); Construct sparse RHS and LHS in coordinate storage format, no zeros are stored \n /// sparse_stencils this is a list which contains only two matrices the LHS and RHS matrices for the Crank-Nicolson solve sparse_stencils = makeSparseStencils(Neuron, res, cap, timeStep); @@ -349,7 +319,7 @@ protected override void PreSolve() protected override void SolveStep(int t) { R = U_Active.Clone(); - explicitUpdate(R, Upre, reactF(reactConst, U_Active, N, M, H, cap), reactF(reactConst, Upre, Npre, Mpre, Hpre, cap), timeStep, 1); + explicitUpdate(R, Upre, reactF(activeIonChannels, U_Active, currentStates, cap), reactF(activeIonChannels, Upre, previousStates, cap), timeStep, 1); var Rsyn = Vector.Build.Dense(Neuron.nodes.Count, 0.0); explicitUpdate(Rsyn, Rsyn, Isyn[0], Isyn[1], timeStep, surfaceArea); @@ -361,22 +331,29 @@ protected override void SolveStep(int t) R.Add(Rsyn, R); SBDF_implicit_decomp.Solve(R.ToArray(), b); - tempState = N.Clone(); - explicitUpdate(N, Npre, fS(N, an(U_Active), bn(U_Active), timeStep), fS(Npre, an(Upre), bn(Upre), timeStep), timeStep, 1); - Npre = tempState.Clone(); - - tempState = M.Clone(); - explicitUpdate(M, Mpre, fS(M, am(U_Active), bm(U_Active), timeStep), fS(Mpre, am(Upre), bm(Upre), timeStep), timeStep, 1); - Mpre = tempState.Clone(); - - tempState = H.Clone(); - explicitUpdate(H, Hpre, fS(H, ah(U_Active), bh(U_Active), timeStep), fS(Hpre, ah(Upre), bh(Upre), timeStep), timeStep, 1); - Hpre = tempState.Clone(); - - N.Map(x => System.Math.Max(0.0, System.Math.Min(1.0, x)), N); - M.Map(x => System.Math.Max(0.0, System.Math.Min(1.0, x)), M); - H.Map(x => System.Math.Max(0.0, System.Math.Min(1.0, x)), H); - + foreach (var channel in activeIonChannels) + { + foreach (var gatingVariable in channel.GatingVariables) + { + if (!gatingVariable.IsInstant){ + tempState = currentStates[gatingVariable.Name].Clone(); + + // Use stateexplicitSBDF2 to update the gating variable + explicitUpdate( + currentStates[gatingVariable.Name], + previousStates[gatingVariable.Name], + fS(currentStates[gatingVariable.Name], gatingVariable.Alpha(U_Active), gatingVariable.Beta(U_Active), timeStep), + fS(previousStates[gatingVariable.Name], gatingVariable.Alpha(Upre), gatingVariable.Beta(Upre), timeStep), + timeStep, 1 + ); + + // Update previous state for the next time step + previousStates[gatingVariable.Name] = tempState.Clone(); + + currentStates[gatingVariable.Name].Map(x => System.Math.Max(0.0, System.Math.Min(1.0, x)),currentStates[gatingVariable.Name]); + } + } + } Upre = U_Active.Clone(); U_Active.SetSubVector(0, Neuron.nodes.Count, Vector.Build.DenseOfArray(b)); @@ -418,7 +395,7 @@ internal override void SetOutputValues() /// this is the axial resistance /// this is membrane resistance scale factor, since this is only a fraction of theoretical maximum /// - public static double SetTargetTimeStep(double cap, double maxDiameter, double minDiameter,double edgeLength ,double gna, double gk, double gl, double res, double cfl) + public static double SetTargetTimeStep(double cap, double maxDiameter, double minDiameter, double edgeLength, double totalConductance, double gl, double res, double cfl) { /// here we set the minimum time step size and maximum time step size /// the dtmin is based on prior numerical experiments that revealed that for each refinement level the @@ -442,7 +419,6 @@ public static double SetTargetTimeStep(double cap, double maxDiameter, double mi //GameManager.instance.DebugLogSafe("lower_bound = " + lower_bound.ToString()); return dt; } - /// /// This function initializes the voltage vector U and the state vectors /// M, N, and H \n @@ -456,15 +432,11 @@ private void InitializeNeuronCell() { lock (visualizationValuesLock) { - U = Vector.Build.Dense(Neuron.nodes.Count, -70.0 * 1.0E-3); + U = Vector.Build.Dense(Neuron.nodes.Count, v0); U_Active = U.Clone(); } Upre = U_Active.Clone(); - M = Vector.Build.Dense(Neuron.nodes.Count, mi); - N = Vector.Build.Dense(Neuron.nodes.Count, ni); - H = Vector.Build.Dense(Neuron.nodes.Count, hi); - Isyn = new List(); for (int i = 0; i < 2; i++) { @@ -474,13 +446,84 @@ private void InitializeNeuronCell() } surfaceArea = Vector.Build.Dense(Neuron.nodes.Count, 0.0); - Mpre = M.Clone(); Npre = N.Clone(); Hpre = H.Clone(); + // Initialize ion channels + InitializeIonChannel(); + totalConductance = activeIonChannels.Sum(ch => ch.Conductance); + + // Dictionary's for Gating Variable States + currentStates = new Dictionary(); + previousStates = new Dictionary(); + + foreach (var channel in activeIonChannels) + { + foreach (var gatingVariable in channel.GatingVariables) + { + currentStates[gatingVariable.Name] = Vector.Build.Dense(Neuron.nodes.Count, gatingVariable.Probability); + previousStates[gatingVariable.Name] = currentStates[gatingVariable.Name].Clone(); + } + } // Set color bar range directly on ColorLUT to bypass serialized field defaults - ColorLUT.GlobalMin = -0.07f; // -100 mV + ColorLUT.GlobalMin = 0.0f; // 0 mV ColorLUT.GlobalMax = 0.05f; // 50 mV } /// + /// Initializes activeionchannel list and all ionchannel list. + /// This spawns each Ion Channel from IonChannels folder + /// + public void InitializeIonChannel() + { + int n = Neuron.nodes.Count; + + ionChannels = new List + { + IonChannelModels.PotassiumChannel(n, v0), + IonChannelModels.SodiumChannel(n, v0), + IonChannelModels.LeakageChannel(n, v0), + IonChannelModels.CalciumChannel(n, v0), + IonChannelModels.SlowPotassiumChannel(n, v0), + IonChannelModels.LowThresholdCalciumChannel(n, v0), + }; + + if (SavedActiveChannelNames == null) + SavedActiveChannelNames = new HashSet { "Potassium Channel", "Sodium Channel", "Leakage Channel" }; + + activeIonChannels = new List(); + foreach (var channel in ionChannels) + { + if (SavedActiveChannelNames.Contains(channel.Name)) activeIonChannels.Add(channel); + } + + leakConductance = ionChannels.Find(ch => ch.Name == "Leakage Channel").Conductance; + } + /// + /// Immediately toggles a channel on or off. + /// Turning off removes it from activeIonChannels; state vectors are left intact. + /// Turning on adds it back and reinitializes its gating variable states from Probability. + /// + public void ToggleChannel(IonChannel channel, bool active) + { + lock (visualizationValuesLock) + { + if (active) + { + if (!activeIonChannels.Contains(channel)) + { + activeIonChannels.Add(channel); + foreach (var gv in channel.GatingVariables) + { + currentStates[gv.Name] = Vector.Build.Dense(Neuron.nodes.Count, gv.Probability); + previousStates[gv.Name] = currentStates[gv.Name].Clone(); + } + } + } + else + { + activeIonChannels.Remove(channel); + } + } + } + /// /// This is for constructing the lhs and rhs of system matrix \n /// This will construct a HINES matrix (symmetric), it should be tridiagonal with some off /// diagonal entries corresponding to a branch location in the neuron graph \n @@ -575,39 +618,43 @@ public static List> makeSparseStencils(Neuron myCell, /// /// these are the conductances and reversal potentials defined by List reactConst = new List { gk, gna, gl, ek, ena, el }; /// this is the voltage vector - /// this is the state vector n - /// this is the state vector m - /// this is the state vector h + /// this is the list of ion channels active in the simulation + /// this is the list of state vectors for all gating variables /// this is the capacitance /// - private static Vector reactF(List reactConst, Vector V, Vector NN, Vector MM, Vector HH, double cap) + + private static Vector reactF(List activeIonChannels, Vector V, Dictionary gatingStates, double cap) { - /// initialize the output vector and prod vector these will be used to assemble the different parts of the reaction calculation \n - /// Vector output = Vector.Build.Dense(V.Count, 0.0); initializes the output vector of length equation to number of entries in voltage vector, initialized to 0 \n - /// Vector prod = Vector.Build.Dense(V.Count, 0.0); initializes the product vector of length equation to number of entries in voltage vector, initialized to 0 \n - /// double ek, ena, el, gk, gna, gl; these are the conductances and reversal potentials that we need to assign using reactConst parameter that is sent \n Vector output = Vector.Build.Dense(V.Count, 0.0); - Vector prod = Vector.Build.Dense(V.Count, 0.0); - double ek, ena, el, gk, gna, gl; - /// this sets the constants for the conductances \n - /// gk = reactConst[0]; gna = reactConst[1]; gl = reactConst[2]; - gk = reactConst[0]; gna = reactConst[1]; gl = reactConst[2]; - /// this sets constants for reversal potentials \n - /// ek = reactConst[3]; ena = reactConst[4]; el = reactConst[5]; - ek = reactConst[3]; ena = reactConst[4]; el = reactConst[5]; - /// output.Add(prod.Multiply(gk), output); this adds current due to potassium - prod.SetSubVector(0, V.Count, NN.PointwisePower(4.0)); - prod.SetSubVector(0, V.Count, (V.Subtract(ek)).PointwiseMultiply(prod)); - output.Add(prod.Multiply(gk), output); - /// output.Add(prod.Multiply(gna), output); this adds current due to sodium - prod.SetSubVector(0, V.Count, MM.PointwisePower(3.0)); - prod.SetSubVector(0, V.Count, HH.PointwiseMultiply(prod)); prod.SetSubVector(0, V.Count, (V.Subtract(ena)).PointwiseMultiply(prod)); - output.Add(prod.Multiply(gna), output); - /// output.Add((V.Subtract(el)).Multiply(gl), output); this adds leak current - output.Add((V.Subtract(el)).Multiply(gl), output); - /// Return the negative of the total - output.Multiply(-1.0 / cap, output); + + foreach (var channel in activeIonChannels) + { + Vector prod = Vector.Build.Dense(V.Count, 1.0); + // Only apply leak current directly + if (channel.Name.Contains("Leak")) + { + output.Add(V.Subtract(channel.ReversalPotential).Multiply(channel.Conductance), output); + } + else + { + foreach (var gatingVariable in channel.GatingVariables) + { + Vector state = gatingStates[gatingVariable.Name]; + // Adding a state to Gating Variables specifically for Low Threshold Calcium From Pospischil_Minimal_HH_2008 (s-instantaneous gating variable). + // A better solution may be available. For my current use case this is sufficient. + if (gatingVariable.IsInstant) { + state = gatingVariable.Alpha(V); + } + // Calculate contribution for this channel and gating variable + prod.SetSubVector(0, V.Count, state.PointwisePower(gatingVariable.Exponent).PointwiseMultiply(prod)); + } + prod.SetSubVector(0, V.Count, V.Subtract(channel.ReversalPotential).PointwiseMultiply(prod)); + output.Add(prod.Multiply(channel.Conductance), output); + } + + } + output.Multiply(-1.0 / cap, output); return output; } @@ -626,7 +673,6 @@ private void explicitUpdate(Vector S, Vector Spre, Vector F, Vector Fpre, double S.Add(Spre.Multiply(-1.0 / 3.0), S); S.Add(Fpre.Multiply(-2.0 * dt * scale / 3.0), S); } - /// /// This is the function for the right hand side of the ODE on state S, which is given by: /// \f[\frac{dS}{dt}=\alpha_S(V)(1-S)-\beta_S(V)S\f] @@ -650,200 +696,41 @@ private static double ScalingFactor(double v, double dt) return System.Math.Max(dt / v, 1.0); } - /// - /// This is \f$\alpha_n\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// an this function returns the rate at the given voltage - private static Vector an(Vector V) - { - var V_arr = Vector.Build.DenseOfVector(V); - V_arr.Multiply(1.0E3, V_arr); - double Vt_scaled = 1.0E3 * Vt; - - var d_alpha = V_arr - Vt_scaled - 15.0; - var alpha = Vector.Build.Dense(V.Count); - - for (int i = 0; i < V.Count; i++) - { - double da = d_alpha[i]; - - /// L'hospital's rule for limit form when da is close to 0, - /// this is to avoid numerical instability and NaN values - if (System.Math.Abs(da) < 1.0E-6) - { - // limit form - alpha[i] = 1.0E3 * -0.032 / - ((-1.0 / 5.0) * System.Math.Exp(-da / 5.0)); - } - else - { - alpha[i] = 1.0E3 * -0.032 * da / - (System.Math.Exp(-da / 5.0) - 1.0); - } - } - - return alpha; - } - /// - /// This is \f$\beta_n\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// bn this function returns the rate at the given voltage - private static Vector bn(Vector V) - { - Vector Vin = Vector.Build.DenseOfVector(V); - Vin.Multiply(1.0E3, Vin); - return (1.0E3) * (0.5) * ((10.0 - Vin + 1.0E3 * Vt) / 40.0).PointwiseExp(); - } - /// - /// This is \f$\alpha_m\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// am this function returns the rate at the given voltage - private static Vector am(Vector V) - { - var V_arr = Vector.Build.DenseOfVector(V); - V_arr.Multiply(1.0E3, V_arr); - double Vt_scaled = 1.0E3 * Vt; - - var d_alpha = V_arr - Vt_scaled - 13.0; - - var alpha = Vector.Build.Dense(V.Count); - - for (int i = 0; i < V.Count; i++) - { - double da = d_alpha[i]; - - /// L'hospital's rule for limit form when da is close to 0, - /// this is to avoid numerical instability and NaN values - if (System.Math.Abs(da) < 1.0E-6) - { - // limit form - alpha[i] = -1.0E3 * 0.32 / - ((-1.0 / 4.0) * System.Math.Exp(-da / 4.0)); - } - else - { - alpha[i] = -1e3 * 0.32 * da / - (System.Math.Exp(-da / 4.0) - 1.0); - } - } + // used by save/load functions in Menu.cs - return alpha; + // Returns a map from each gating-variable name to its current values [V.Count] + // returns kvp (key-value pair) + public Dictionary getCurrentStates() { + return currentStates.ToDictionary( + kvp => kvp.Key, kvp => kvp.Value.AsArray() + ); } - /// - /// This is \f$\beta_m\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// bm this function returns the rate at the given voltage - private static Vector bm(Vector V) - { - var V_arr = Vector.Build.DenseOfVector(V); - V_arr.Multiply(1.0E3, V_arr); - double Vt_scaled = 1.0E3 * Vt; - - var d_beta = V_arr - Vt_scaled - 40.0; - - var beta = Vector.Build.Dense(V.Count); - - for (int i = 0; i < V.Count; i++) - { - double db = d_beta[i]; - - /// L'hospital's rule for limit form when da is close to 0, - /// this is to avoid numerical instability and NaN values - if (System.Math.Abs(db) < 1.0E-6) - { - // limit form - beta[i] = 1.0E3 * 0.28 / - ((1.0 / 5.0) * System.Math.Exp(db / 5.0)); - } - else - { - beta[i] = 1.0E3 * 0.28 * db / - (System.Math.Exp(db / 5.0) - 1.0); - } - } - return beta; - } - /// - /// This is \f$\alpha_h\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// ah this function returns the rate at the given voltage - private static Vector ah(Vector V) - { - Vector Vin = Vector.Build.DenseOfVector(V); - Vin.Multiply(1.0E3, Vin); - return (1.0E3) * (0.128) * ((17.0 - Vin + 1.0E3 * Vt) / 18.0).PointwiseExp(); + // Returns a map from each gating-variable name to its current values [V.Count] + // returns kvp (key-value pair) + public Dictionary getPreviousStates() { + return previousStates.ToDictionary( + kvp => kvp.Key, kvp => kvp.Value.AsArray() + ); } - /// - /// This is \f$\beta_h\f$ rate function, the rate functions take the form of - /// \f[ - /// \frac{a_p(V-B_p)}{\exp(\frac{V-B_p}{C_p})-D_p} - /// \f] - /// the constants \f$a_p,B_p,C_p,D_p\f$ are manually coded in for this version of the simulation\n - /// TODO: come up with an implementation where the user can enter in their own parameters (Yale Neuron has this capability) - /// - /// this is the input voltage - /// bh this function returns the rate at the given voltage - private static Vector bh(Vector V) - { - Vector Vin = Vector.Build.DenseOfVector(V); - Vin.Multiply(1.0E3, Vin); - return (1.0E3) * 4.0 / (((40.0 - Vin + 1.0E3 * Vt) / 5.0).PointwiseExp() + 1.0); - } - - // used by save/load functions in Menu.cs - public double[] getM() { return M.AsArray(); } - public double[] getN() { return N.AsArray(); } - public double[] getH() { return H.AsArray(); } public double[] getUpre() { return Upre.AsArray(); } - public double[] getMpre() { return Mpre.AsArray(); } - public double[] getNpre() { return Npre.AsArray(); } - public double[] getHpre() { return Hpre.AsArray(); } - public void BuildVectors(double[] u, double[] m, double[] n, double[] h, - double[] upre, double[] mpre, double[] npre, double[] hpre) + public void BuildVectors(double[] u, double[] upre, + Dictionary currStates, Dictionary prevStates) { lock (visualizationValuesLock) U = Vector.Build.DenseOfArray(u); lock (visualizationValuesLock) U_Active = U.Clone(); Upre = Vector.Build.DenseOfArray(upre); - M = Vector.Build.DenseOfArray(m); - N = Vector.Build.DenseOfArray(n); - H = Vector.Build.DenseOfArray(h); + // convert every double[] -> Vector + + currentStates = currStates.ToDictionary( + kvp => kvp.Key, + kvp => Vector.Build.DenseOfArray(kvp.Value)); - Mpre = Vector.Build.DenseOfArray(mpre); - Npre = Vector.Build.DenseOfArray(npre); - Hpre = Vector.Build.DenseOfArray(hpre); + previousStates = prevStates.ToDictionary( + kvp => kvp.Key, + kvp => Vector.Build.DenseOfArray(kvp.Value)); Isyn = new List(); for (int i = 0; i < 2; i++)