% # iNMDA_PYdr_PYso_JB12:
%
% Normalized synaptic NMDAergic excitatory current, with synaptic
% depression and minis, FROM the pyramidal axo-soma TO pyramidal dendrite
% PYdr<-PYso connection used in the DynaSim implementation of (Benita et
% al., 2012).
%
% - Dependencies:
%     - netconNearestNeighbors.m
%
% - References:
%     - Benita, J. M., Guillamon, A., Deco, G., & Sanchez-Vives, M. V. (2012).
%     Synaptic depression and slow oscillatory activity in a biophysical
%     network model of the cerebral cortex. Frontiers in Computational
%     Neuroscience, 6. https://doi.org/10.3389/fncom.2012.00064
%
% - Tags: synapse, connection, excitation, ampa
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
% For DynaSim, we need to convert synaptic maximal conductance from
% absolute to area-relative terms, while also converting the units to
% mS/cm^2 to keep pace with the rest of the equations:
%
% 0.9 nS (gNMDA EE) / 0.035 mm^2 (dendrite area) * ...
% 1 mS / 1000000 nS * 100 mm^2 / 1 cm^2 = ...
% 0.00257 mS/cm^2
%
gNMDA = 0.00257       % mS/cm^2
ENMDA = 0 % mV

alphaS = 0.5
tauS = 100 % ms
alphaX = 3.48
tauX = 2 % ms

deprFactor = 0.9
tauRes = 400 % ms

resIC = 0.8
resNoiseIC = 0.1
sNMDAIC = 0
sNMDANoiseIC = 0.1
xNMDAIC = 0
xNMDANoiseIC = 0.1

% Connectivity
% Connective radius, aka how many target cells each source cell connects
% to, from the source's perspective.
radius = 10

% Remove autapses to the dendrite corresponding to this soma
removeRecurrentBool = 1
% We also need to normalize the conductance in mS/cm^2 by the number of
% connections each target cell is receiving on average, so that the TOTAL
% sum of their conductive inputs adds to our overall maximal conductance
% above.
normalizingFactor = min(((2*radius + (1-removeRecurrentBool)) / (N_post/N_pre)), N_pre)

% Note that what is passed is 2x the radius
netcon = netconNearestNeighbors(2*radius, N_pre, N_post, removeRecurrentBool)


% Functions
INMDA_PYdr_PYso_JB12(X,sNMDA,rec) = -gNMDA/normalizingFactor.*((res.*sNMDA)*netcon).*(X-ENMDA)

monitor INMDA_PYdr_PYso_JB12

% This is the synaptic activity state variable
sNMDA' = alphaS.*xNMDA.*(1-sNMDA) - sNMDA./tauS
sNMDA(0) = sNMDAIC+sNMDANoiseIC.*rand(1,N_pre)
xNMDA' = alphaX./(1 + exp(-(X_pre-20)./2)) - xNMDA./tauX
xNMDA(0) = xNMDAIC+xNMDANoiseIC.*rand(1,N_pre)

% This is the 'resources' state variable
res' = (1 - res)./tauRes + max((t-tspike_pre)<=(2*dt)).*(-(1 - res)./tauRes + (deprFactor.*res - res)./dt)
res(0) = resIC-resNoiseIC.*rand(1,N_pre)

% Linker
@current += INMDA_PYdr_PYso_JB12(X_post,sNMDA,res)