Hello Brian2 community,
I am trying to implement a thalamocortical model for propofol induced anesthesia from this paper, " Rapid thalamocortical network switching mediated by cortical synchronization underlies propofol-induced EEG signatures: a biophysical model" by Soplata et al. published in JNP in 2023.
Keywords: PY(Pyramidal neurons), IN(Interneurons), TC(Thalamic relay cells), TRN(Thalamic reticular nucleus)
The authors have implemented their model and study in Dynasim (MATLAB based toolbox for neural simulations). When I ran their code in Dynasim then the results were almost similar to as given in their paper. My goal is to implement their proposed model in Brian2 for more openness and scalability and extend it to a mesoscale whole-brain model for anesthesia with the help of experimentally observed data. I translated their code and implemented their complete study in Brian2. In their study they have included a short-term depression of intracortical excitatory connections for AMPAergic and NMDAergic synapses (PY>PY, PY>IN) and for short-term depression of inhibitory GABAergic synaptic connections (IN-PY, IN>IN). I guess that it is much similar to STDP (spike time dependent plasticity) or maybe it is STDP itself because the authors haven’t mentioned about it much in the paper. However, I am not able to implement this mechanism in Brian2 exactly similar to Dynasim. The following equations from the paper indicate this mechanism:
The complete equations for the model can be seen from the supplementary material available on the paper link mentioned earlier.
The Dynasim implementation of this mechanism for PY>PY AMPAergic synapse is as follows:
% # iAMPA_PYdr_PYso_JB12:
%
% Normalized synaptic AMPAergic 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:
%
% 5.4 nS (gAMPA EE) / 0.035 mm^2 (dendrite area) * ...
% 1 mS / 1000000 nS * 100 mm^2 / 1 cm^2 = ...
% 0.0154 mS/cm^2
%
gAMPAMult = 1 % for experimenting with gAMPA
gAMPA = gAMPAMult * 0.0154 % mS/cm^2
EAMPA = 0 % mV
alpha = 3.48
tauS = 2 % ms
deprFactor = 0.9
tauRes = 400 % ms
sAMPAIC = 0.0
sAMPANoiseIC = 0.1
resIC = 0.8
resNoiseIC = 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
IAMPA_PYdr_PYso_JB12(X,sAMPA,rec) = -gAMPA/normalizingFactor.*((res.*sAMPA)*netcon).*(X-EAMPA)
monitor IAMPA_PYdr_PYso_JB12
% This is the synaptic activity state variable
sAMPA' = alpha./(1 + exp(-(X_pre-20)./2)) - sAMPA./tauS
sAMPA(0) = sAMPAIC+sAMPANoiseIC.*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 += IAMPA_PYdr_PYso_JB12(X_post,sAMPA,res)
So what I am facing the issue is that I am not getting the same results as Dynasim in Brian2 because of some implementation issue in Brian2 or possibly the internal environment differences between the two which I am not able to figure out. From the equations and Dynasim code we can see that there is a resources variable which is depressed by some factor whenever a presynaptic spike occurs. I tried to implement exactly the same logic as given in the paper in Brian2 but what I observe is that the plots for the resources variable are quite different between Brian2 and Dynasim.
I am also confused about whether to include the refractory period (I have kept around 2ms) for all the cells in Brian2 or not. Because Dynasim based model does not include it in the code or mention anywhere in the paper. The model is based on Hodgkin-Huxley dynamics. If I don’t include that then the cells fire at a much higher rate than as compared to Dynasim.
So in summary, any help in solving this issue or the exact implementation of model like Dynasim in Brian2 would be really appreaciated. Apologies for such a long post :'). Looking forward to hearing from you all if you are interested in it and I am eager to learn from this wonderful community. It would be great help for my project as well as it might inform others to translate Dynasim based models to Brian2 or possibly create such large models on their own in Brian2.
Thanks in advance!
). This is about individual synapses that get depressed after repeated activation, and does not depend on the spiking activity of the post-synaptic cell (as STDP does)
Now I get it, thanks for clearing confusion between STDP and STP.
for upgrading me to basis user 
.