dsUnitMake_autogenData

% Make DynaSim Autogen Files

0001 %% Make DynaSim Autogen Files
0002 
0003 warning('off','MATLAB:lang:cannotClearExecutingFunction');
0004 
0005 cwd = pwd;
0006 
0007 % Set where to save outputs
0008 output_directory = fullfile(dsGetConfig('ds_data_path'), 'autoGenSave_temp');
0009 
0010 % Erase output_directory
0011 if isdir(output_directory)
0012   rmdir(output_directory, 's')
0013 end
0014 
0015 % move to root directory where outputs will be saved
0016 mkdirSilent(output_directory);
0017 cd(output_directory);
0018 
0019 autogenOptions = {'auto_gen_test_data_flag',1, 'random_seed', 1, 'verbose_flag',0, 'visible', 'off'};
0020 
0021 %% Lorenz equations
0022 eqns={
0023   's=10; r=27; b=2.666';
0024   'dx/dt=s*(y-x)';
0025   'dy/dt=r*x-y-x*z';
0026   'dz/dt=-b*z+x*y';
0027 };
0028 data=dsSimulate(eqns, 'tspan',[0 100], 'ic',[1 2 .5], 'solver','rk4', 'study_dir','demo_lorenz', autogenOptions{:});
0029 dsPlot(data, autogenOptions{:}); close gcf
0030 
0031 %% Izhikevich neuron with noisy drive
0032 eqns={
0033   'C=100; vr=-60; vt=-40; k=.7; Iapp=70; ton=10; toff=20';
0034   'a=.03; b=-2; c=-50; d=100; vpeak=35';
0035   'dv/dt=(k*(v-vr)*(v-vt)-u+I(t))/C; v(0)=vr';
0036   'du/dt=a*(b*(v-vr)-u); u(0)=0';
0037   'if(v>vpeak)(v=c; u=u+d)';
0038   'I(t)=Iapp*(t>ton&t<toff)*(1+.5*rand)'; % define applied input using reserved variables 't' for time and 'dt' for fixed time step of numerical integration
0039   'monitor I';                            % indicate to store applied input during simulation
0040 };
0041 data=dsSimulate(eqns, 'tspan',[0 100], 'study_dir','demo_izhikevich', autogenOptions{:});
0042 dsPlot(data, autogenOptions{:}); close gcf
0043 
0044 %% RUNNING SETS OF Izhikevich SIMULATIONS
0045 eqns={
0046   'a=.02; b=.2; c=-65; d=6; I=14';
0047   'dv/dt=.04*v^2+5*v+140-u+I; v(0)=-70';
0048   'du/dt=a*(b*v-u); u(0)=-20';
0049   'if(v>=30)(v=c;u=u+d)';
0050   };
0051 P='pop1'; % name of population
0052 vary={
0053   {P,'a',.02; P,'b',.2 ; P,'c',-50; P,'d',2;  P,'I',15} % tonic bursting
0054   {P,'a',.01; P,'b',.2 ; P,'c',-65; P,'d',8;  P,'I',30} % spike frequency adaptation
0055   {P,'a',.02; P,'b',.2 ; P,'c',-65; P,'d',6;  P,'I',7}  % spike latency
0056   {P,'a',.03; P,'b',.25; P,'c',-52; P,'d',0;  P,'I',0}  % rebound burst
0057   {P,'a',1;   P,'b',1.5; P,'c',-60; P,'d',0;  P,'I',-65}% bistability
0058   {P,'a',.02; P,'b',1  ; P,'c',-55; P,'d',4;  P,'I',1}  % accomodation
0059   {P,'a',-.02;P,'b',-1 ; P,'c',-60; P,'d',8;  P,'I',80} % inhibition-induced spiking
0060   {P,'a',-.026;P,'b',-1; P,'c',-45; P,'d',0;  P,'I',70} % inhibition-induced bursting
0061   };
0062 
0063 %% normal
0064 data=dsSimulate(eqns, 'tspan',[0 100], 'vary',vary, 'study_dir','demo_izhikevich_vary', autogenOptions{:});
0065 
0066 % TODO: fix dsPlot bug
0067 % dsPlot(data, autogenOptions{:}); close gcf
0068 
0069 %% parfor
0070 data=dsSimulate(eqns, 'tspan',[0 100], 'vary',vary, 'study_dir','demo_izhikevich_vary_parfor',...
0071   'parallel_flag',1, autogenOptions{:});
0072 
0073 %% parfor compile
0074 data=dsSimulate(eqns, 'tspan',[0 100], 'vary',vary, 'study_dir','demo_izhikevich_vary_parfor_compile',...
0075   'parallel_flag',1, 'compile_flag',1, autogenOptions{:});
0076 
0077 %% compile
0078 data=dsSimulate(eqns, 'tspan',[0 100], 'vary',vary, 'study_dir','demo_izhikevich_vary_compile',...
0079   'compile_flag',1, autogenOptions{:});
0080 
0081 %% Hodgkin-Huxley neuron equations (without predefined mechanisms)
0082 eqns={
0083   'gNa=120; gK=36; Cm=1';
0084   'INa(v,m,h) = gNa.*m.^3.*h.*(v-50)';
0085   'IK(v,n) = gK.*n.^4.*(v+77)';
0086   'dv/dt = (10-INa(v,m,h)-IK(v,n))/Cm; v(0)=-65';
0087   'dm/dt = aM(v).*(1-m)-bM(v).*m; m(0)=.1';
0088   'dh/dt = aH(v).*(1-h)-bH(v).*h; h(0)=.1';
0089   'dn/dt = aN(v).*(1-n)-bN(v).*n; n(0)=0';
0090   'aM(v) = (2.5-.1*(v+65))./(exp(2.5-.1*(v+65))-1)';
0091   'bM(v) = 4*exp(-(v+65)/18)';
0092   'aH(v) = .07*exp(-(v+65)/20)';
0093   'bH(v) = 1./(exp(3-.1*(v+65))+1)';
0094   'aN(v) = (.1-.01*(v+65))./(exp(1-.1*(v+65))-1)';
0095   'bN(v) = .125*exp(-(v+65)/80)';
0096 };
0097 data=dsSimulate(eqns, 'study_dir','demo_hh_1', autogenOptions{:});
0098 
0099 dsPlot(data, autogenOptions{:}); close gcf
0100 
0101 % Equivalent Hodgkin-Huxley neuron with predefined mechanisms
0102 data=dsSimulate('dv/dt=10+@current/Cm; Cm=1; v(0)=-65; {iNa,iK}', 'study_dir','demo_hh_2', autogenOptions{:});
0103 
0104 % Example of a bursting neuron model using three active current mechanisms:
0105 eqns='dv/dt=5+@current; {iNaF,iKDR,iM}; gNaF=100; gKDR=5; gM=1.5; v(0)=-70';
0106 data=dsSimulate(eqns, 'tspan',[0 100], 'study_dir','demo_hh_3', autogenOptions{:});
0107 dsPlot(data, autogenOptions{:}); close gcf
0108 
0109 %% Sparse Pyramidal-Interneuron-Network-Gamma (sPING)
0110 
0111 % define equations of cell model (same for E and I populations)
0112 eqns={
0113   'dv/dt=Iapp+@current+noise*randn(1,N_pop)';
0114   'monitor iGABAa.functions, iAMPA.functions'
0115 };
0116 % Tip: monitor all functions of a mechanism using: monitor MECHANISM.functions
0117 
0118 % create DynaSim specification structure
0119 s=[];
0120 s.populations(1).name='E';
0121 s.populations(1).size=8;
0122 s.populations(1).equations=eqns;
0123 s.populations(1).mechanism_list={'iNa','iK'};
0124 s.populations(1).parameters={'Iapp',5,'gNa',120,'gK',36,'noise',40};
0125 s.populations(2).name='I';
0126 s.populations(2).size=2;
0127 s.populations(2).equations=eqns;
0128 s.populations(2).mechanism_list={'iNa','iK'};
0129 s.populations(2).parameters={'Iapp',0,'gNa',120,'gK',36,'noise',40};
0130 s.connections(1).direction='I->E';
0131 s.connections(1).mechanism_list={'iGABAa'};
0132 s.connections(1).parameters={'tauD',10,'gSYN',.1,'netcon','ones(N_pre,N_post)'};
0133 s.connections(2).direction='E->I';
0134 s.connections(2).mechanism_list={'iAMPA'};
0135 s.connections(2).parameters={'tauD',2,'gSYN',.1,'netcon',ones(8,2)};
0136 
0137 data=dsSimulate(s, 'study_dir','demo_sPING_0', autogenOptions{:});
0138 dsPlot(data, autogenOptions{:}); close all
0139 
0140 %% Covary
0141 % Vary two parameters (run a simulation for all combinations of values)
0142 vary={'(E,I)','(EK1,EK2)',[-80 -60]};
0143 %       This sets modifications:
0144 %           * E_EK1, E_EK2, I_EK1, I_EK2 = -80
0145 %           * E_EK1, E_EK2, I_EK1, I_EK2 = -60
0146 data = dsSimulate(s, 'study_dir','demo_sPING_covary1', 'vary',vary, autogenOptions{:});
0147               
0148 vary={'(E,I)','(EK1,EK2)',[-80 -60; -85 -65]};
0149 %       This sets modifications:
0150 %           * E_EK1, I_EK1 = -80 and E_EK2, I_EK2 = -85
0151 %           * E_EK1, I_EK1 = -60 and E_EK2, I_EK2 = -65
0152 data = dsSimulate(s, 'study_dir','demo_sPING_covary2', 'vary',vary, autogenOptions{:});
0153               
0154 vary={'(E,I)','(EK1,EK2)',cat(3,[-80 -60], [-85 -65])};
0155 %       This sets modifications:
0156 %           * E_EK1, E_EK2 = -80 and I_EK1, I_EK2 = -85
0157 %           * E_EK1, E_EK2 = -60 and I_EK1, I_EK2 = -65
0158 data = dsSimulate(s,'study_dir','demo_sPING_covary3', 'vary',vary, autogenOptions{:});
0159               
0160 vary={'(E,I)','(EK1,EK2)',cat(3, [-75 -55; -80 -60], [-85 -65; -90 -70])};
0161 %       This sets modifications:
0162 %           * E_EK1 = -75, E_EK2 = -80, I_EK1 = -85, I_EK2 = -90.
0163 %           * E_EK1 = -55, E_EK2 = -60, I_EK1 = -65, I_EK2 = -70.
0164 data = dsSimulate(s,'study_dir','demo_sPING_covary4', 'vary',vary, autogenOptions{:});
0165 
0166 
0167 %% SAVING SIMULATED DATA
0168 
0169 %% Save data from a single simulation
0170 % Example using the previous sPING model:
0171 data=dsSimulate(s, 'save_data_flag',1, 'study_dir','demo_sPING_1', autogenOptions{:});
0172 
0173 %% Save data from a set of simulations
0174 
0175 % Specify what to vary
0176 % Tip: use 'vary' Syntax 2 to systematically vary a parameter
0177 vary={'E','Iapp',[0 10 20]}; % vary the amplitude of tonic input to E-cells
0178 data=dsSimulate(s, 'save_data_flag',1, 'study_dir','demo_sPING_2',...
0179                      'vary',vary, autogenOptions{:});
0180 
0181 % load and plot the saved data
0182 % data_from_disk = dsImport('demo_sPING_2', 'auto_gen_test_data_flag',1);
0183 % dsPlot(data_from_disk, autogenOptions{:}); close all
0184 % dsPlot(data_from_disk,'variable','E_v', autogenOptions{:}); close all
0185 
0186 % Vary two parameters (run a simulation for all combinations of values)
0187 vary={
0188   'E'   ,'Iapp',[0 10 20];      % amplitude of tonic input to E-cells
0189   'I->E','tauD',[5 10 15]       % inhibition decay time constant from I to E
0190   };
0191 data=dsSimulate(s, 'save_data_flag',1, 'study_dir','demo_sPING_3',...
0192                 'vary',vary, 'verbose_flag',1, autogenOptions{:});
0193 % data=dsImport('demo_sPING_3', autogenOptions{:});
0194 dsPlot(data, autogenOptions{:}); close all
0195 dsPlot(data,'plot_type','rastergram', autogenOptions{:}); close all
0196 dsPlot(data,'plot_type','power', autogenOptions{:}); close all
0197 dsPlotFR(data, autogenOptions{:}); close all
0198 
0199 %% RUNNING SIMULATIONS ON THE CLUSTER
0200 
0201 % Run three simulations in parallel jobs and save the simulated data
0202 eqns='dv/dt=@current+I; {iNa,iK}';
0203 vary={'','I',[0 10 20]};
0204 dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_loop',...
0205                    'vary',vary,'cluster_flag',1, autogenOptions{:});
0206 
0207 % Array mode
0208 eqns='dv/dt=@current+I; {iNa,iK}';
0209 vary={'','I',[0 10 20]};
0210 dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_array',...
0211                    'vary',vary,'cluster_flag',1, 'qsub_mode', 'array',...
0212                    autogenOptions{:});
0213 
0214 % Repeat but also save plotted data
0215 eqns='dv/dt=@current+I; {iNa,iK}';
0216 vary={'','I',[0 10 20]};
0217 dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_2',...
0218                    'vary',vary,'cluster_flag',1, autogenOptions{:},...
0219                    'plot_functions',@dsPlot);
0220 
0221 % Save multiple plots and pass custom options to each plotting function
0222 eqns='dv/dt=@current+I; {iNa,iK}';
0223 vary={'','I',[0 10 20]};
0224 dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_3',...
0225                    'vary',vary,'cluster_flag',1, autogenOptions{:},...
0226                    'plot_functions',{@dsPlot,@dsPlot},...
0227                    'plot_options',{{},{'plot_type','power'}});
0228 
0229 %% Compilation
0230 
0231 % Compile
0232 data=dsSimulate(s,'compile_flag',1, 'study_dir','demo_sPING_3_compile', autogenOptions{:});
0233 
0234 % Combine compilation and parallelization
0235 vary={'','I',0:2:14};
0236 data=dsSimulate(s, 'compile_flag',1, 'parallel_flag',1, 'vary', vary, 'study_dir','demo_sPING_3_compile_parallel', autogenOptions{:});
0237 dsPlot(data, autogenOptions{:}); close all
0238 
0239 %% Matlab Solvers
0240 % Run three simulations in parallel jobs and save the simulated data
0241 eqns='dv/dt=@current+I; {iNa,iK}';
0242 vary={'','I',[0 10 20]};
0243 
0244 % od45
0245 data = dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_ml_solver_ode45',...
0246                    'vary',vary,'cluster_flag',0,'overwrite_flag',1,...
0247                    'dt',0.01, 'downsample_factor',100,'solver','ode45',...
0248                     'matlab_solver_options', {'InitialStep', 0.01}, 'compile_flag',0, autogenOptions{:});
0249   
0250 % ode45 compiled
0251 data = dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_ml_solver_ode45_compiled',...
0252                    'vary',vary,'cluster_flag',0,'overwrite_flag',1,...
0253                    'dt',0.01, 'downsample_factor',100,'solver','ode45',...
0254                     'matlab_solver_options', {'InitialStep', 0.01}, 'compile_flag',1, autogenOptions{:});
0255                
0256 % ode23
0257 data = dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_ml_solver_ode23s',...
0258                    'vary',vary,'cluster_flag',0,'overwrite_flag',1,...
0259                    'dt',0.01, 'downsample_factor',100,'solver','ode23s',...
0260                     'matlab_solver_options', {'InitialStep', 0.01}, 'compile_flag',0, autogenOptions{:});
0261                   
0262 % ode23 compiled
0263 data = dsSimulate(eqns,'save_data_flag',1,'study_dir','demo_cluster_1_ml_solver_ode23s_compiled',...
0264                    'vary',vary,'cluster_flag',0,'overwrite_flag',1,...
0265                    'dt',0.01, 'downsample_factor',100,'solver','ode23s',...
0266                     'matlab_solver_options', {'InitialStep', 0.01}, 'compile_flag',1, autogenOptions{:});
0267 
0268 %% one file mode
0269 study_dir=fullfile('.', 'study_HH_varyI_one_file');
0270 eqns='dv/dt=@current+I; {iNa,iK}';
0271 vary={'','I',[0:10:30]};
0272 
0273 dsSimulate(eqns,'vary',vary, 'study_dir',study_dir,'save_data_flag',1,...
0274   'cluster_flag',1, 'qsub_mode','array', 'one_solve_file_flag',1,...
0275   'tspan',[0 100], 'dt',0.01, 'solver','euler', 'downsample_factor', 100,...
0276   'plot_functions',{@dsPlot}, 'plot_options',{{'format', 'jpg', 'visible', 'off'}},...
0277   autogenOptions{:});
0278 
0279 %% Remove output_directory
0280 cd(cwd)
0281 
0282 fprintf('Removing temporary output_directory: %s\n', output_directory)
0283 rmdir(output_directory, 's')
0284 
0285 fprintf('\nDone making autogen data.\n\n')