data = calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants, current_suffix, varargin) v1.0 David Stanley, Boston University 2017, stanleyd@bu.edu
0001 function data = dsCalcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants, current_suffix, varargin) 0002 %data = calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants, current_suffix, varargin) 0003 % 0004 % v1.0 David Stanley, Boston University 2017, stanleyd@bu.edu 0005 0006 % 0007 % Purpose: Calculates synaptic or ionic currents post-simulation, given state 0008 % variables. This allows the user to avoid recording synaptic currents 0009 % using monitor functions, which can be slow in large simulations. 0010 % 0011 % Usage: 0012 % calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants, current_suffix, varargin) 0013 % calcCurrentPosthoc(data, mechanism_prefix, current_equation); 0014 % calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields); 0015 % calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants); 0016 % calcCurrentPosthoc(data, mechanism_prefix, current_equation, additional_fields, additional_constants, current_suffix); 0017 % 0018 % Inputs: 0019 % data: DynaSim data structure 0020 % mechanism_prefix: Prefix of the mechanism to use to calculate the current. Most times this will be the same as the 0021 % mechanism filename (e.g. iNa) 0022 % current_equation: Char array containing the formula used to calculate the synaptic or ionic current. 0023 % (E.g. 'm.^3*h.^1*(iNa_V - ENa)'). Generally this can be copied directly from the 0024 % function line in the original mechanism. 0025 % 0026 % Inputs (Optional): 0027 % additional_fields: Additional fields to pull out from data that current_equation 0028 % might need access too. Useful for referencing 0029 % state variables outside of the current 0030 % mechanism, such as the membrane voltage of the 0031 % parent population, or variables in another 0032 % mechanism (e.g. calcium pool). This is provided 0033 % in the form of a cell array of chars, listing 0034 % the desired field names. The full field names 0035 % will then be available in current_equation (see 0036 % examples). 0037 % 0038 % additional_constants: Any additional constants you might be using in 0039 % current_equation. This would likely be things 0040 % like the reversal potential. This is provided 0041 % in the form of a structure: 0042 % 0043 % additional_constants.(constant_name) = constant_value 0044 % 0045 % Constants "constant_name" will then be accessible 0046 % within current_equation (see examples). 0047 % 0048 % current_suffix: Suffix to append to the newly created field 0049 % (default: 'posthoc') 0050 % 0051 % Outputs: 0052 % data: DynaSim data structure 0053 % 0054 % Examples: 0055 % % % % % % % % % % % % % % % % % % Example 1 % % % % % % % % % % % % % % % % % 0056 % tic; data1=dsSimulate('dv[5]/dt=10+@current; {iNa,iK}; monitor iNa.functions','tspan',[0 200]); toc % Simulate with monitoring 0057 % tic; data2=dsSimulate('dv[5]/dt=10+@current; {iNa,iK};','tspan',[0 200]); toc; % Simulate without montioring 0058 % % Add the iNa currents to data2 0059 % mechanism_prefix = 'pop1_iNa'; 0060 % additional_constants = struct; 0061 % additional_constants.ENa = 50; % From mech file, iNa.mech 0062 % additional_constants.gNa = 120; % From mech file, iNa.mech 0063 % current_string = 'gNa.*m.^3.*h.*(pop1_v-ENa)'; % Adapted from mech file, iNa.mech 0064 % % m and h are state variables with prefix matching mechanism_prefix (pop1_iNa) 0065 % % ENa and gNa are supplied as additional constants and hence will be recognized 0066 % % pop1_v is specified below in "additional fields" 0067 % additional_fields = {'pop1_v'}; 0068 % data2b = dsCalcCurrentPosthoc(data2,mechanism_prefix, current_string, additional_fields, additional_constants, 'INa'); 0069 % 0070 % figure; plot(data1.pop1_iNa_INa,data1.pop1_v); xlabel('INa'); ylabel('Vm'); title('Vm vs INa with monitor functions'); 0071 % figure; plot(data2b.pop1_iNa_INa,data2b.pop1_v); xlabel('INa'); ylabel('Vm'); title('Vm vs INa with posthoc calculation'); 0072 % 0073 % % % % % % % % % % % % % % % % % Example 2 - PING network % % % % % % % % % % % % % % % % % 0074 % % Simulate PING network with and without monitor functions 0075 % eqns1={ 0076 % 'dv/dt=Iapp+@current+noise*randn(1,N_pop)'; 0077 % 'monitor functions' 0078 % }; 0079 % 0080 % eqns2={ 0081 % 'dv/dt=Iapp+@current+noise*randn(1,N_pop)'; 0082 % }; 0083 % 0084 % scale_factor = 1; 0085 % N_E = 80*scale_factor; 0086 % N_I = 20*scale_factor; 0087 % % create DynaSim specification structure 0088 % s=[]; 0089 % s.populations(1).name='E'; 0090 % s.populations(1).size=N_E; 0091 % s.populations(1).mechanism_list={'iNa','iK'}; 0092 % s.populations(1).parameters={'Iapp',5,'gNa',120,'gK',36,'noise',0}; 0093 % s.populations(2).name='I'; 0094 % s.populations(2).size=N_I; 0095 % s.populations(2).mechanism_list={'iNa','iK'}; 0096 % s.populations(2).parameters={'Iapp',0,'gNa',120,'gK',36,'noise',0}; 0097 % s.connections(1).direction='I->E'; 0098 % s.connections(1).mechanism_list={'iGABAa'}; 0099 % s.connections(1).parameters={'tauD',10,'gSYN',.1,'netcon','ones(N_pre,N_post)'}; 0100 % s.connections(2).direction='E->I'; 0101 % s.connections(2).mechanism_list={'iAMPA'}; 0102 % s.connections(2).parameters={'tauD',2,'gSYN',.1,'netcon',ones(N_E,N_I)}; 0103 % 0104 % % Simulate Sparse Pyramidal-Interneuron-Network-Gamma (sPING) 0105 % s.populations(1).equations=eqns1; 0106 % s.populations(2).equations=eqns1; 0107 % tic; data1=dsSimulate(s,'tspan',[0 200]); toc; 0108 % 0109 % % Re-simulate without monitor functions 0110 % s.populations(1).equations=eqns2; 0111 % s.populations(2).equations=eqns2; 0112 % tic; data2=dsSimulate(s,'tspan',[0 100]); toc; 0113 % 0114 % % Add synaptic currents to data2 0115 % mechanism_prefix = 'I_E_iAMPA'; 0116 % additional_constants = struct; 0117 % additional_constants.ESYN = 0; % From mech file, iAMPA.mech 0118 % additional_constants.gSYN = 0.1; % From mech file, iAMPA.mech 0119 % additional_constants.netcon = ones(N_E,N_I); % From mech file, iAMPA.mech 0120 % current_string = '(gSYN.*(s*netcon).*(I_v-ESYN))';% Adapted from mech file, iAMPA.mech 0121 % % s is state variables with prefix matching mechanism_prefix (I_E_iAMPA) 0122 % % I_v is the presynaptic voltage. calCurrentPostdoc knows to look for it in data 0123 % % because we specified it as "additional_fields" (see below) 0124 % % gSYN and ESYN are supplied as additional constants and hence will be recognized 0125 % 0126 % additional_fields = {'I_v'}; 0127 % data2b = dsCalcCurrentPosthoc(data2,mechanism_prefix, current_string, additional_fields, additional_constants, 'ISYN'); 0128 % 0129 % figure; plot(data1.I_E_iAMPA_ISYN,data1.I_v); xlabel('iAMPA'); ylabel('Vm'); title('Vm vs iAMPA with monitor functions'); 0130 % figure; plot(data2b.I_E_iAMPA_ISYN,data2b.I_v); xlabel('iAMPA'); ylabel('Vm'); title('Vm vs iAMPA with posthoc calculation'); 0131 % 0132 % 0133 % Note: This could ultimately be built into the core dynasim solver, 0134 % as a way of circumventing monitor functions running on the fly; instead, 0135 % all monitors could be calculated posthoc. 0136 % 0137 % Author: Dave Stanley, Boston University, 2017 0138 % 0139 % See also: thevEquiv 0140 0141 %%Code%% 0142 0143 0144 %% localfn output 0145 if ~nargin 0146 fout1 = localfunctions; 0147 return 0148 end 0149 0150 %% auto_gen_test_data_flag argin 0151 options = dsCheckOptions(varargin,{'auto_gen_test_data_flag',0,{0,1}},false); 0152 if options.auto_gen_test_data_flag 0153 % specific to this function 0154 varargs = varargin; 0155 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0156 argin = [{fin1},{fin2}, varargs]; 0157 end 0158 0159 %% main fn body here 0160 if nargin <= 3 0161 error('Need to supply at least the first 3 arguments. Remainders are optional'); 0162 end 0163 0164 0165 % Pull supplied constants into workspace, if any 0166 if nargin >= 4 0167 if isempty(additional_constants); additional_constants = struct; end 0168 vars_pull(additional_constants); 0169 end 0170 0171 ac_bool = false; 0172 if nargin >= 5 0173 ac_bool = true; % boolean for instructions to pull in additional constants 0174 if ischar(additional_fields); additional_fields = {additional_fields}; end 0175 if ~iscellstr(additional_fields); error('additional_fields must be a char array or cell array of chars'); end 0176 end 0177 0178 if nargin < 6 0179 current_suffix = 'posthoc'; 0180 end 0181 0182 if ~isempty(strfind(current_suffix,'_')) 0183 warning('Suffix should not contain underscores (_). Removing'); 0184 current_suffix = strrep(current_suffix,'_',''); 0185 end 0186 0187 if strcmp(mechanism_prefix(end),'_') 0188 warning('Mechnaism_prefix should not contain a trailing underscore (e.g. RS_). Removing'); 0189 mechanism_prefix = mechanism_prefix(1:end-1); 0190 end 0191 0192 new_fieldname = [mechanism_prefix '_' current_suffix]; 0193 0194 0195 N = length(data); 0196 for i = 1:N 0197 % Pull state variables into workspace, if any 0198 myfield = fieldnames(data(1)); 0199 ind = strfind(myfield,mechanism_prefix); % Get set of state variables matching the current mechanism prefix 0200 ind2 = ~cellfun(@isempty,ind); 0201 chosen_fields = myfield(ind2); 0202 vars_pull(data,chosen_fields); % Pull in the variables 0203 0204 % For each state variable, mechanism_prefix portion of the variable name 0205 for j = 1:length(chosen_fields) 0206 curr_varname = chosen_fields{j}; 0207 curr_varname_shortened = strrep(curr_varname,[mechanism_prefix '_'],''); 0208 eval([curr_varname_shortened ' = ' curr_varname ';']); 0209 end 0210 0211 % Pull additional variables into the workspace 0212 if ac_bool 0213 vars_pull(data,additional_fields) 0214 end 0215 0216 % Evaluate the expression for calculating the current 0217 eval(['temp = ' current_equation ';']); 0218 0219 % Store this in the data structure 0220 data(i).(new_fieldname) = temp; 0221 if isempty(strcmp(data(1).labels,new_fieldname)) 0222 data(i).labels(end+1) = {new_fieldname}; 0223 end 0224 end 0225 0226 %% auto_gen_test_data_flag argout 0227 if options.auto_gen_test_data_flag 0228 % specific to this function 0229 argout = {output, fout2}; 0230 0231 dsSaveAutoGenTestData(argin, argout); 0232 end 0233 0234 end