AxonThreshPairs.m

function Is = AxonThreshPairs(depth,diam);

% calculate axon thresholds based on pairwise combinations of depth and diameter

if size(depth) ~= size(diam)
    error('The two input arguments must be of equal size');
end    

L = 400;                                % length of axon field (cm)
dVmt = 10;                              % depolarization threshold (mV)
dAC = 2.5;                              % distance between anode and cathode (cm)
dz = 0.01;                              % axial spatial resolution (cm)

zA = L/2 + dAC/2;                       % anode position (cm)
zC = zA - dAC;                          % cathode position (cm)
z = 0:dz:L;                             % vector of axial positions (cm)

Rm = 5000;                    % specific membrane resistance (ohm.cm^2)
Ra = 70;                      % specific axial resistance (ohm.cm)
re = 350;                     % extracellular resistance (ohm.cm)

Is = zeros(size(depth));

for i = 1: length(depth)
    
    r = diam(i) / 1e4 / 2;                  % axon radius (cm)
    rm = Rm./(2*pi.*r);                     % membrane resistance (ohm.cm)
    ra = Ra./(pi.*r.^2);                    % axial resistance (ohm/cm)
    k = sqrt(ra/rm);                        % spatial frequency of intracellular potential
                                            % 1/um
    d = depth(i) / 1e4;                     % convert depth to cm
    
    dA = sqrt( (d-r)^2 + (zA-z).^2 );       % absolute distances from electrodes
    dC = sqrt( (d-r)^2 + (zC-z).^2 );       % to axon surface coordinates (cm)

    VA = re./(4*pi*dA);                     % potential outside membrane due to 1mA anode (mV)
    VC = re./(4*pi*dC);                     % potential outside membrane due to 1mA cathode (mV)
    Ve = VA - VC;                           % potential applied to membrane (mV)

    dVi = Visolve(Ve,rm,ra,dz);
    dVm1 = max(dVi - Ve);

    Is(i) = dVmt ./ dVm1;
            
end    

return