Development (#30)

*General improvements for  code and documentation

Co-authored-by: João Pedro Gomes dos Santos <[email protected]>
Co-authored-by: João Pedro Gomes dos Santos <[email protected]>
Co-authored-by: Kadri Pajo <[email protected]>
Co-authored-by: Kadri Pajo <[email protected]>
Co-authored-by: dkeller9 <[email protected]>

Author: 5 people

.gitignore

@@ -1,8 +1,4 @@
-
 **/NEURON/single_sec_cell.py
-!**/Matlab/Model/Place models here.txt
-**/Matlab/Model/**
 *.pdf
 *.tiff
-*.png
-!**/Docs/Images/**
+*.png

Code for paper figures/Reproduce paper figure 3/spike.m renamed to Code for paper figures/Common scripts/spike.m

@@ -1,4 +1,6 @@
 function thisAmp = spike(timePoint, basalAmp, jumpAmp, risingOffset, k1, k2)
+% Script provided by Anu G. Nair to reproduce the calcium spike present in
+% his 2016 paper (https://doi.org/10.1371/journal.pcbi.1005080)
 
     tForMax = log(k1/k2)/(k1-k2);
     maxAmpForThisRate = (exp(-k1*tForMax) - exp(-k2*tForMax));

Code for paper figures/Reproduce paper figure 6/spiketraindd_Ca.m renamed to Code for paper figures/Common scripts/spiketraindd_Ca.m

@@ -1,6 +1,10 @@
 function thisAmp = spiketraindd_Ca(t, tStart, frequency, maxSpikeNum, ...
-                            basalAmp, maxAmp, k1, k2, spikelengthLowerThresh, spikeAmpMaxThresh, burstFrequency, maxBurstNum, burstMinInter, regorand)
-                       
+                            basalAmp, maxAmp, k1, k2, ...
+                            spikelengthLowerThresh, spikeAmpMaxThresh,...
+                            burstFrequency, maxBurstNum, burstMinInter,...
+                            regorand)
+% Script provided by Anu G. Nair to reproduce the calcium spike present in
+% his 2016 paper (https://doi.org/10.1371/journal.pcbi.1005080)
 
     t = t * 1000;
     cv = 0.1;

Code for paper figures/Reproduce paper figure 3/Analysis_new.mat renamed to Code for paper figures/Figure 3/Analysis_diagnostics_fig3.mat

@@ -0,0 +1,28 @@
+function thisAmp = Ca_input_fig3(times,n)
+persistent CaX
+persistent CaT
+persistent Cah1
+persistent Cah2
+if isempty(CaX)
+    Data = coder.load('../Subcellular_workflow/Code for paper figures/Figure 3/Nair_2016_calcium.mat','calcium','time');
+    CaX = Data.calcium{n};
+    CaT = Data.time{n};
+    Cah1 = 1;
+end
+thisAmp = CaX(end);
+if times ~= 30
+    if times == 0
+        Cah1 = 1;
+        thisAmp = CaX(1);
+    else
+        while times > CaT(Cah1)
+            Cah1 = Cah1 + 1;
+        end
+        while times < CaT(Cah1 -1)
+            Cah1 = Cah1 - 1;
+        end
+        Cah2 = (CaT(Cah1)-times)*1/(CaT(Cah1)-CaT(Cah1-1));
+        thisAmp = (CaX(Cah1-1)*Cah2 + CaX(Cah1)*(1-Cah2));
+    end
+end
+end

Code for paper figures/Figure 3/Ca_input_fig3.m

@@ -0,0 +1,103 @@
+function f_Reproduce_data_fig3_CD(folder)
+%% Running this script will generate figures 1C and 2C from Nair et al 2016
+%% using the updated model with new parameters governing CaMKII autophosphorylation
+warning('off','SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless')
+
+load(folder + "model_Nair_2016_optimized_fig3.mat",'modelobj');
+obj = modelobj;
+
+%% Get steady state values after equilibrating for 50000 s
+
+cnfst = getconfigset(obj);
+cnfst.SolverType = 'ode15s';
+cnfst.StopTime = 50000;
+cnfst.TimeUnits = 'second';
+cnfst.SolverOptions.AbsoluteTolerance = 1e-7;
+cnfst.SolverOptions.RelativeTolerance = 1e-4;
+cnfst.CompileOptions.UnitConversion = 1;
+cnfst.SolverOptions.AbsoluteToleranceScaling = 1;
+cnfst.RunTimeOptions.StatesToLog = 'all';
+
+[t,species,~] = sbiosimulate(obj);
+
+SteadyState=species(length(t),:);
+
+for i=1:length(SteadyState)
+    if SteadyState(i)<0
+        SteadyState(i)=0;
+    end
+end
+
+for i=1:length(SteadyState)
+    set(obj.species(i), 'InitialAmount', SteadyState(i));
+end
+
+%% Add parameters used by the Calcium input
+
+ruleobj=addrule(obj, 'Spine.Ca = Ca_input_fig3(time,Ca_input)');
+set(ruleobj,'RuleType','RepeatedAssignment');
+addparameter(obj, 'Ca_input', 0);
+set(obj.parameters(end), 'ValueUnits', 'dimensionless');
+
+% Assign DA_expression to dopamine
+ruleobj=addrule(obj, 'Spine.DA = Spine.DA_expression');
+set(ruleobj,'RuleType','RepeatedAssignment');
+
+%% Run simulations
+
+DA = -4:0.2:4;
+CaStart = 4;
+cnfst.StopTime = 30;
+cnfst.SolverOptions.MaxStep = 0.01;
+cnfst.SolverOptions.OutputTimes = 0:0.01:30;
+cnfst.RuntimeOptions.StatesToLog = {'pSubstrate'};
+set(obj.rules(2), 'Active', 1);
+set(obj.rules(3), 'Active', 0);
+
+obj.parameters(end).Value = 1;%Ca_input
+[t_noDA,x_noDA,~] = sbiosimulate(obj);
+activationArea = sum(x_noDA) - x_noDA(1) * length(x_noDA);
+
+set(obj.rules(3), 'Active', 1);
+set(obj.parameters(228), 'ValueUnits', 'second'); % par 228 = DA_start
+obj.parameters(228).Value = CaStart + 1; % par 228 = DA_start
+obj.parameters(end).Value = 2;%Ca_input
+
+[t_DA,x_DA,~] = sbiosimulate(obj);
+
+activationAreaWithMultipleDA = zeros(1,length(DA));
+
+for n = 1:length(DA)
+objm{n} = copyobj(obj);
+objm{n}.parameters(end).Value = n+2;%Ca_input
+objm{n}.parameters(228).Value = CaStart + DA(n); % par 228 = DA_start
+objm_cnfst{n} = getconfigset(objm{n});
+objm_cnfst{n}.StopTime = 30;
+objm_cnfst{n}.SolverOptions.MaxStep = 0.01;
+objm_cnfst{n}.SolverOptions.OutputTimes = 0:0.01:30;
+objm_cnfst{n}.RuntimeOptions.StatesToLog = {'pSubstrate'};
+objm_cnfst{n}.SolverType = 'ode15s';
+objm_cnfst{n}.TimeUnits = 'second';
+objm_cnfst{n}.SolverOptions.AbsoluteTolerance = 1e-7;
+objm_cnfst{n}.SolverOptions.RelativeTolerance = 1e-4;
+objm_cnfst{n}.CompileOptions.UnitConversion = 1;
+objm_cnfst{n}.SolverOptions.AbsoluteToleranceScaling = 1;
+end
+
+parfor n = 1:length(DA)
+    warning('off','SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless')
+    [~,x,~] = sbiosimulate(objm{n},objm_cnfst{n});
+    activationAreaWithMultipleDA(n) = sum(x) - x(1) * length(x);
+%     ME = MException(errID,msgtext)
+end
+
+Nair_2016_optimized_Matlab_data_fig3{1}(:,1) = t_noDA;
+Nair_2016_optimized_Matlab_data_fig3{1}(:,2) = x_noDA;
+Nair_2016_optimized_Matlab_data_fig3{2}(:,1) = t_DA;
+Nair_2016_optimized_Matlab_data_fig3{2}(:,2) = x_DA;
+Nair_2016_optimized_Matlab_data_fig3{3}(:,1) = DA;
+Nair_2016_optimized_Matlab_data_fig3{3}(:,2) = activationAreaWithMultipleDA/activationArea;
+
+save(folder + "Nair_2016_optimized_Matlab_data_fig3.mat",...
+    'Nair_2016_optimized_Matlab_data_fig3')
+end