Work done to IonChannel generic function

Author: Carlos Eckert

Assets/Scripts/C2M2/NeuronalDynamics/Simulation/IonChannel.cs

@@ -3,17 +3,25 @@
 
 public class GatingVariable
 {
+    // Alpha function for the gating variable
     public Func<double, double> Alpha { get; set; }
+
+    // Beta function for the gating variable
     public Func<double, double> Beta { get; set; }
+
+    // Name of the gating variable (i.e: n, m, h)
     public string Name { get; set; }
-    public List<double> Coefficients { get; set; }
+    // 
+    public List<double> ChannelConstants { get; set; }
+    public double InitialValue { get; set; }
 
     public GatingVariable(string name, Func<double, double> alpha, Func<double, double> beta)
     {
         Name = name;
         Alpha = alpha;
         Beta = beta;
-        Coefficients = new List<double>();
+        ChannelConstants = new List<double>();
+        InitialValue = InitialValue;
     }
 }
 
@@ -31,7 +39,6 @@ public IonChannel(string name, double conductance, double reversalPotential)
         ReversalPotential = reversalPotential;
         GatingVariables = new List<GatingVariable>();
     }
-
     public void AddGatingVariable(GatingVariable gatingVariable)
     {
         GatingVariables.Add(gatingVariable);

Assets/Scripts/C2M2/NeuronalDynamics/Simulation/IonChannel/obj/Debug/net8.0/IonChannel.AssemblyInfo.cs

@@ -13,7 +13,7 @@
 [assembly: System.Reflection.AssemblyCompanyAttribute("IonChannel")]
 [assembly: System.Reflection.AssemblyConfigurationAttribute("Debug")]
 [assembly: System.Reflection.AssemblyFileVersionAttribute("1.0.0.0")]
-[assembly: System.Reflection.AssemblyInformationalVersionAttribute("1.0.0+bed9a50b52f2757c8e2930f61c4b8fd3ac41c29d")]
+[assembly: System.Reflection.AssemblyInformationalVersionAttribute("1.0.0+361b14b31f5b2aa50307dc72e44a360504076790")]
 [assembly: System.Reflection.AssemblyProductAttribute("IonChannel")]
 [assembly: System.Reflection.AssemblyTitleAttribute("IonChannel")]
 [assembly: System.Reflection.AssemblyVersionAttribute("1.0.0.0")]

Assets/Scripts/C2M2/NeuronalDynamics/Simulation/IonChannel/obj/Debug/net8.0/IonChannel.AssemblyInfoInputs.cache

@@ -1 +1 @@
-0f2750cddfbea5680889400aec062eef631265d0509eb50adc500f3b464c949f
+3510c39156aab660920e3d3e16bb304c6d20e39712337b1fef7b4171c4045928

Assets/Scripts/C2M2/NeuronalDynamics/Simulation/SparseSolverTestv1.cs

@@ -64,36 +64,6 @@ public class SparseSolverTestv1 : NDSimulation
         /// </summary>
         private double cap = 1.0 * 1.0E-2;
         /// <summary>
-        /// [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.
-        /// </summary>
-        private double gk = 5.0 * 1.0E1;
-        /// <summary>
-        /// [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.
-        /// </summary>
-        private double gna = 50.0 * 1.0E1;
-        /// <summary>
-        /// [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.
-        /// </summary>
-        private double gl = 0.0 * 1.0E1;
-        /// <summary>
-        /// [V] potassium reversal potential
-        /// </summary>
-        private double ek = -90.0 * 1.0E-3;
-        /// <summary>
-        /// [V] sodium reversal potential
-        /// </summary>
-        private double ena = 50.0 * 1.0E-3;
-        /// <summary>
-        /// [V] leak reversal potential
-        /// </summary>
-        private double el = -70.0 * 1.0E-3;
-        /// <summary>
         /// [] potassium channel state probability, unitless
         /// </summary>
         private double ni = 0.0376969;
@@ -380,28 +350,51 @@ protected override void PreSolve()
             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.M, g.N, g.H, g.Upre, g.Mpre, g.Npre, g.Hpre);
+
+            ///<c>R</c> this is the reaction vector for the reaction solve
+            R = Vector.Build.Dense(Neuron.nodes.Count);            
+
+            tempState = Vector.Build.Dense(Neuron.nodes.Count, 0);
+            ///<c>reactConst</c> this is a small list for collecting the conductances and reversal potential which is sent to the reaction solve routine
+            // reactConst = new List<double> { 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);
+            ///UnityEngine.Debug.Log("Target Time Step = " + timeStep);
+
+            ///<c>List<CoordinateStorage<double>> sparse_stencils = makeSparseStencils(Neuron, res, cap, k);</c> Construct sparse RHS and LHS in coordinate storage format, no zeros are stored \n
+            /// <c>sparse_stencils</c> 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);
+            ///<c>CompressedColumnStorage</c> call Compresses the sparse matrices which are stored in <c>sparse_stencils[0]</c> and <c>sparse_stencils[1]</c>
+            r_csc = CompressedColumnStorage<double>.OfIndexed(sparse_stencils[0]); //null;
+            l_csc = CompressedColumnStorage<double>.OfIndexed(sparse_stencils[1]); //null;
+            ///<c>double [] b</c> we define storage for the diffusion solve part
+            b = new double[Neuron.nodes.Count];
+            ///<c>var lu = SparseLU.Create(l_csc, ColumnOrdering.MinimumDegreeAtA, 0.1);</c> this creates the LU decomposition of the HINES matrix which is defined by <c>l_csc</c>
+            lu = SparseLU.Create(l_csc, ColumnOrdering.MinimumDegreeAtA, 0.1);
 
-            // Define alpha and beta functions for the gating variable (example for potassium channel)
-            Func<double, double> alpha_n = (voltage) => 0.01 * (voltage + 55) / (1 - Math.Exp(-(voltage + 55) / 10));
-            Func<double, double> beta_n = (voltage) => 0.125 * Math.Exp(-(voltage + 65) / 80);
+            // // Define alpha and beta functions for the gating variable (example for potassium channel)
+            // Func<double, double> alpha_n = (voltage) => 0.01 * (voltage + 55) / (1 - Math.Exp(-(voltage + 55) / 10));
+            // Func<double, double> beta_n = (voltage) => 0.125 * Math.Exp(-(voltage + 65) / 80);
 
-            // Create a gating variable for potassium channel (n gate)
-            GatingVariable potassiumGate = new GatingVariable("n", alpha_n, beta_n)
-            {
-                Coefficients = new List<double> { 1.0 }
-            };
+            // // Create a gating variable for potassium channel (n gate)
+            // GatingVariable potassiumGate = new GatingVariable("n", alpha_n, beta_n)
+            // {
+            //     Coefficients = new List<double> { 1.0 }
+            // };
 
-            // Create the potassium ion channel
-            IonChannel potassiumChannel = new IonChannel("Potassium", , ); 
+            // // Create the potassium ion channel
+            // IonChannel potassiumChannel = new IonChannel("Potassium", , ); 
 
-            // Add the gating variable to the potassium channel
-            potassiumChannel.AddGatingVariable(potassiumGate);
+            // // Add the gating variable to the potassium channel
+            // potassiumChannel.AddGatingVariable(potassiumGate);
 
-            // Calculate the current at a certain voltage
-            double voltage = -65;
-            double current = potassiumChannel.CalculateCurrent(voltage);
+            // // Calculate the current at a certain voltage
+            // double voltage = -65;
+            // double current = potassiumChannel.CalculateCurrent(voltage);
 
-            Console.WriteLine($"Potassium channel current at {voltage} mV: {current} A");
+            // Console.WriteLine($"Potassium channel current at {voltage} mV: {current} A");
 
 
 
@@ -428,6 +421,61 @@ protected override void PreSolve()
             lu = SparseLU.Create(l_csc, ColumnOrdering.MinimumDegreeAtA, 0.1);
         }
 
+        public void InitializeIonChannels()
+        {
+            /// <summary>
+            /// [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.
+            /// </summary>
+            private double gk = 5.0 * 1.0E1;
+            /// <summary>
+            /// [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.
+            /// </summary>
+            private double gna = 50.0 * 1.0E1;
+            /// <summary>
+            /// [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.
+            /// </summary>
+            private double gl = 0.0 * 1.0E1;
+            /// <summary>
+            /// [V] potassium reversal potential
+            /// </summary>
+            private double ek = -90.0 * 1.0E-3;
+            /// <summary>
+            /// [V] sodium reversal potential
+            /// </summary>
+            private double ena = 50.0 * 1.0E-3;
+            /// <summary>
+            /// [V] leak reversal potential
+            /// </summary>
+            private double el = -70.0 * 1.0E-3;
+
+            // 
+            private double n_alphaFunction;
+            private double m_alphaFunction;
+            private double h_alphaFunction;
+            private double n_betaFunction;
+            private double m_betaFunction;
+            private double h_betaFunction;
+
+
+
+            IonChannel potassiumChannel = new IonChannel("Potassium Channel", gk, ek);
+            potassiumChannel.AddGatingVariable(new GatingVariable("n", alpha: n_alphaFunction, beta: n_betaFunction));
+
+            IonChannel sodiumChannel = new IonChannel("Sodium Channel", gna, ena);
+            sodiumChannel.AddGatingVariable(new GatingVariable("m", alpha: m_alphaFunction, beta: m_betaFunction));
+            sodiumChannel.AddGatingVariable(new GatingVariable("h", alpha: h_alphaFunction, beta: h_betaFunction));
+
+            IonChannel leakageChannel = new IonChannel("Leakage Channel", gl, el);
+            leadChannel.AddGatingVariable(new GatingVariable("n"));
+            
+        }
+
         /// <summary>
         /// This is the main solver, it is running on it own thread.
         /// The solver using SBDF2 for time steping, the implicit part is used for the diffusion and the explicit