dsThevEquiv

 Calculates the Th�venin equivalent voltage and conductance for a
 given set of M specified ionic channels.
 Inputs:
   data - DynaSim data structure (see dsCheckData)
   fields_currents - 1xM cell array of field names that
           contain the ionic currents (M entries, one for each ionic channel).
   field_voltage - 1x1 string specifying membrane voltage
   reversals_list - 1xM array containing a list of all reversal
           potential
   output_field_name - string containing the desired base output field
           name (ETH and gTH will be appended for Th�venin reversal potential
           and conductance).
 Outputs:
   data: data structure containing summed data
 Example:
   data2 = dsThevEquiv(data,{'IB_NG_IBaIBdbiSYNseed_ISYN','IB_NG_iGABABAustin_IGABAB','IB_FS_IBaIBdbiSYNseed_ISYN'},'IB_V',[-95,-95,-95]);

 Algorithm:
 Ionic channel's conductance, j, will be estimated by g_j(t) = I_j(t) ./ (V(t) - E_j)
 Th�venin equivalent conductance and reversal will be estimated as follows
        gTh(t) = \sum_j{g_j(t)}
        ETH = \sum_j( g_j(t) * E_j(t) ) / \sum_j( g_j(t) )
 v1.0 David Stanley, Boston University 2016, stanleyd@bu.edu

0001 function data = dsThevEquiv(data, fields_currents, field_voltage, reversals_list, output_field_name,varargin)
0002 % Calculates the Th�venin equivalent voltage and conductance for a
0003 % given set of M specified ionic channels.
0004 % Inputs:
0005 %   data - DynaSim data structure (see dsCheckData)
0006 %   fields_currents - 1xM cell array of field names that
0007 %           contain the ionic currents (M entries, one for each ionic channel).
0008 %   field_voltage - 1x1 string specifying membrane voltage
0009 %   reversals_list - 1xM array containing a list of all reversal
0010 %           potential
0011 %   output_field_name - string containing the desired base output field
0012 %           name (ETH and gTH will be appended for Th�venin reversal potential
0013 %           and conductance).
0014 % Outputs:
0015 %   data: data structure containing summed data
0016 % Example:
0017 %   data2 = dsThevEquiv(data,{'IB_NG_IBaIBdbiSYNseed_ISYN','IB_NG_iGABABAustin_IGABAB','IB_FS_IBaIBdbiSYNseed_ISYN'},'IB_V',[-95,-95,-95]);
0018 %
0019 % Algorithm:
0020 % Ionic channel's conductance, j, will be estimated by g_j(t) = I_j(t) ./ (V(t) - E_j)
0021 % Th�venin equivalent conductance and reversal will be estimated as follows
0022 %        gTh(t) = \sum_j{g_j(t)}
0023 %        ETH = \sum_j( g_j(t) * E_j(t) ) / \sum_j( g_j(t) )
0024 % v1.0 David Stanley, Boston University 2016, stanleyd@bu.edu
0025 
0026     
0027     if nargin < 5
0028         % Default field name
0029         output_field_name = 'thev_equiv';
0030         
0031         % Add default prefix
0032         ind = strfind(fields_currents{1},'_');
0033         prefix = fields_currents{1}(1:ind);
0034         output_field_name = strcat(prefix,output_field_name);
0035     end
0036     
0037     % Add prefix if necessary
0038     if isempty(strfind(output_field_name,'_'))
0039         warning('No population prefix specified (see documentation). Adding a default prefix.');
0040         ind = strfind(fields_currents{1},'_');
0041         prefix = fields_currents{1}(1:ind);
0042         output_field_name = strcat(prefix,output_field_name);
0043     end
0044     
0045     data = dsCheckData(data, varargin{:});
0046     % note: calling dsCheckData() at beginning enables analysis function to
0047     % accept data matrix [time x cells] in addition to DynaSim data structure.
0048 
0049     for i = 1:length(data)
0050         dat = data(i);
0051         gTH = zeros(size(dat.(fields_currents{1})));
0052         sum_gj_times_Ej = zeros(size(dat.(fields_currents{1})));
0053         for j = 1:length(fields_currents)
0054             g_j = dat.(fields_currents{j}) ./ ( dat.(field_voltage) - reversals_list(j) ); % g_j = I ./ (V - E)
0055             
0056             gTH = gTH + g_j;
0057             sum_gj_times_Ej = sum_gj_times_Ej + (g_j * reversals_list(j));
0058         end
0059         ETH = sum_gj_times_Ej ./ gTH;
0060         
0061         
0062         on1 = strcat(output_field_name,'_ETH');
0063         on2 = strcat(output_field_name,'_gTH');
0064         data(i).(on1) = ETH;
0065         data(i).(on2) = gTH;
0066         if all(strcmp(data(i).labels,on1) == false); data(i).labels(end+1) = {on1}; end
0067         if all(strcmp(data(i).labels,on2) == false); data(i).labels(end+1) = {on2}; end
0068 
0069     end
0070 
0071 end