late night

Author: Carlos Eckert

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

@@ -61,7 +61,7 @@ public class SparseSolverTestv1 : NDSimulation
         /// <summary>
         /// [F/m2] capacitance per unit area, this is the plasma membrane capacitance, this a standard value for the capacitance
         /// </summary>
-                private double cap = 1.0 * 1.0E-2;
+        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]
@@ -157,7 +157,7 @@ public class SparseSolverTestv1 : NDSimulation
         /// </summary>
         private Vector tempState;
 
-        List<double> reactConst;                            //This is for passing the reaction function constants
+        // List<double> reactConst;                            //This is for passing the reaction function constants
         List<CoordinateStorage<double>> sparse_stencils;
         CompressedColumnStorage<double> r_csc;              //This is for the rhs sparse matrix
         CompressedColumnStorage<double> l_csc;              //This is for the lhs sparse matrix
@@ -396,7 +396,6 @@ protected override void PreSolve()
                 Debug.Log("Here");
                 InitializeIonChannel();
 
-
                 if (!g.Loading)
                 {
                     InitializeNeuronCell();
@@ -427,6 +426,7 @@ protected override void PreSolve()
                 // Other PreSolve logic (time step, matrices setup, etc.)
                 R = Vector.Build.Dense(Neuron.nodes.Count);
                 tempState = Vector.Build.Dense(Neuron.nodes.Count, 0);
+                timeStep = SetTargetTimeStep(cap, 2 * Neuron.MaxRadius, 2 * Neuron.MinRadius, Neuron.TargetEdgeLength, ionChannels, res, 1.0);
                 sparse_stencils = makeSparseStencils(Neuron, res, cap, timeStep);
                 r_csc = CompressedColumnStorage<double>.OfIndexed(sparse_stencils[0]);
                 l_csc = CompressedColumnStorage<double>.OfIndexed(sparse_stencils[1]);
@@ -519,7 +519,7 @@ protected override void SolveStep(int t)
                         gatingVariable.CurrentState,
                         gatingVariable.PreviousState,
                         fS(gatingVariable.CurrentState, gatingVariable.Alpha(U_Active), gatingVariable.Beta(U_Active)),
-                        fS(gatingVariable.PreviousState, gatingVariable.Alpha(Upre), gatingVariable.Beta(Upre)),
+                        fS(gatingVariable.PreviousState, gatingVariable.Alpha(Upre), gatingVariable.Beta(Upre)), 
                         timeStep
                     );
 
@@ -563,12 +563,12 @@ internal override void SetOutputValues()
         /// <param name="cap"></param> this is the capacitance
         /// <param name="maxDiameter"></param> this is the maximum diameter
         /// <param name="edgeLength"></param> this is the target edge length of the graph geometry
-        /// <param name="gna"></param> this is the sodium conductance
-        /// <param name="gk"></param> this is the potassium conductance
+        /// <param name="ionChannels"></param> this is the list of ionChannels. Used to grather conductance of each channel
         /// <param name="res"></param> this is the axial resistance
         /// <param name="Rmemscf"></param> this is membrane resistance scale factor, since this is only a fraction of theoretical maximum
+        /// <param name="cfl"></param> this isn't defined or used in the function but appears to be intended as a safety variable for the time step. Example usage: dt = Math.Min(upper_bound, dtmax) * cfl;
         /// <returns></returns>
-        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, List<IonChannel> ionChannels, 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 
@@ -577,18 +577,29 @@ public static double SetTargetTimeStep(double cap, double maxDiameter, double mi
             //double dtmin = 2e-6;  
             double dtmax = 50e-6;
             double dt;
+            double scf = 1E-6; // to convert to micrometer of edgelengths and radii don't forget this!!!!
 
-            double gll = gl; double scf = 1E-6; // to convert to micrometer of edgelengths and radii don't forget this!!!!
-
-            // what happens if the leak conductance is 0
-            if (gll == 0.0) { gll = 1.0; }
+            // Calculate effective conductance based on all ion channels
+            double effectiveConductance = 0.0; 
 
-            double upper_bound = cap * edgeLength*scf * System.Math.Sqrt(res / (gll*minDiameter*scf));
+            foreach (IonChannel channel in ionChannels)
+            {
+                effectiveConductance += channel.Conductance;
+            }
+
+            // if effective conductance is 0
+            if (effectiveConductance == 0.0)
+            {
+                effectiveConductance = 1.0; // Set a minimal conductance
+            }
+
+
+            double upper_bound = cap * edgeLength * scf * System.Math.Sqrt(res / (effectiveConductance * minDiameter * scf));
             //double lower_bound = cap * edgeLength*scf * System.Math.Sqrt(res / (gna + gk + gl) * maxDiameter*scf);
             //GameManager.instance.DebugLogSafe("upper_bound = " + upper_bound.ToString());
 
             // some cells may have an upper bound that is too large for the solver, so choose the smaller of the two dtmax or upper_bound
-            dt = System.Math.Min(upper_bound,dtmax);
+            dt = System.Math.Min(upper_bound, dtmax);
             //GameManager.instance.DebugLogSafe("lower_bound = " + lower_bound.ToString());
             return dt;       
         }
@@ -932,3 +943,4 @@ private void LoadSavedStates(GameManager g)
         }
     }
 }
+