Cell-Type-Specific Synaptic Scaling Mechanisms Differentially Contribute to Associative Learning.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches

[1] § Methods › Multi-compartment model ↔ model.py, lines 439–480 · score 0.61 · apical dendrite, basal dendrite, SST populations, soma, firing rate, model
[2] § Results › Memory undergoes transient generalization caused by Hebbian plasticity before gradually achieving specificity ↔ model_analysis.py, lines 1134–1200 · score 0.56 · steady state, associative memory, active, synaptic scaling, pre, onset
[3] § Results › Characterization of the temporal evolution of memory representations ↔ model_analysis.py, lines 193–241 · score 0.54 · excitatory firing rate, factor Hebbian learning, longer, timescale, synaptic scaling, post

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,331 lines · 81 KB · no license · 2 matches

import numpy as np
import matplotlib.pyplot as plt
from util import *
import sys
from model import *
from plotting_functions import *
import os
# from parameters import *
import pickle
def analyze_model(hour_sim, flags_list, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\',
K=0.25, flag_only_S_on=False, run_simulation=True, save_results = False, plot_results=False,modulation_SST=0):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list: contains a list of tuples. Each tuple is a collection of all the flags (e.g. synaptic scaling, hebbian learning, ...)
:param flags_theta: used to study the behaviour of the model (no longer useful). Theta1 for population1 and Theta2 for population2
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param g_top_down_to_S: Represents the top-down signal to SST neurons triggered by the hyperexcitation. It reaches
SST neurons at the offset of the conditioning
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially (thermalization)
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 10) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200_000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim) #time points for the first 15s
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2) ##time points for the seoncd phase 4/24/48h
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(1, 0), (0, 1)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(0.5, 0), (0, 0.5)])
elif modulation_SST < 0:
g_stim_S = np.array([(-0.5, 0), (0, -0.5)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_hebb = 240 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta = 24 * (60 * 60) # time constant of target activity in seconds(24h)
tau_beta = 28 * (60 * 60) # time constant of target activity regulator in seconds(28h)
tau_scaling_E = 8 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 8 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 8 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S = 1.5, 1.5, 1.5
rheobases = (rheobase_E, rheobase_P, rheobase_S)
# Background inputs
g_E = 4.5
g_P = 3.2
g_S = 3
back_inputs = (g_E, g_P, g_S)
# Initial conditions for plastic weights
# w_EP_within = 0.81; w_EP_cross = 0.41
# w_ES_within = 0.81; w_ES_cross = 0.31
# w_EE_within = 0.71; w_EE_cross = 0.41
w_EP_within = 0.91; w_EP_cross = 0.41
w_ES_within = 0.51; w_ES_cross = 0.31
w_EE_within = 0.51; w_EE_cross = 0.51
# # Initial conditions for plastic weights
# w_EP_within = 0.7; w_EP_cross = w_EP_within*0.3
# w_ES_within = 0.7; w_ES_cross = w_ES_within*0.3
# w_EE_within = 0.5; w_EE_cross = 0.4
# Weights
w_PE_within = 0.3; w_PE_cross = 0.1
w_PP_within = 0.2; w_PP_cross = 0.1
w_PS_within = 0.3; w_PS_cross = 0.1
w_SE_within = 0.4; w_SE_cross = 0.1
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
J_EE_phase1 = np.zeros((4, n_time_points_stim), dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((10, n_time_points_phase2),
dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
J_phase2 = np.zeros((12, n_time_points_phase2),
dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, r_phase2, r_phase3, max_E)
l_res_weights = (J_EE_phase1, J_phase2)
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_' + str(hour_sim) + 'h' + '_k' + str(K).replace(".","")
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
#All flags = 0 and simulation is 30 seconds long. It is used to evaluate what happens when activating E1 what's the response of E2. Afterwards it is evaluating the av_threshold given the result
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=(0,0,0,0,0,0), flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=flags, flags_theta=flags_theta)
if save_results:
l_results = [l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates,
l_res_weights,
av_threshold, stim_times, stim_duration, sim_duration]
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates, l_res_weights,
av_threshold, stim_times, stim_duration, sim_duration] = l_results
if plot_results:
print('Plotting the results.')
time_plots([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights, av_threshold,
stim_times, dir_plot + name, hour_sim,modulation_SST, flag_only_S_on=flag_only_S_on, format='.png')
time_plots([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights, av_threshold,
stim_times, dir_plot + name, hour_sim,modulation_SST, flag_only_S_on=flag_only_S_on, format='.pdf')
#this function is a generalized version of the one above. this one, with the right flags, is the only one necessary. For clarity, they are separated
def analyze_model_timescales(hour_sim, flags_list, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\',
K=0.25, flag_only_S_on=False, run_simulation=True, save_results = False, plot_results=False,modulation_SST=0,timescales_exploration=False):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list: contains a list of tuples. Each tuple is a collection of all the flags (e.g. synaptic scaling, hebbian learning, ...)
:param flags_theta: used to study the behaviour of the model (no longer useful). Theta1 for population1 and Theta2 for population2
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially (thermalization)
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 10) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200_000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim) #time points for the first 15s
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2) ##time points for the seoncd phase 4/24/48h
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(1, 0), (0, 1)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(0.5, 0), (0, 0.5)])
elif modulation_SST < 0:
g_stim_S = np.array([(-0.5, 0), (0, -0.5)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_hebb = 240 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta_list = [12*3600]
tau_beta_list = [0.01*3600]
# tau_theta_list = (np.arange(24-12, 24+12) * 3600).tolist()
# tau_beta_list = (np.arange(28-12, 28+12) * 3600).tolist()
# tau_theta_list = [6*3600, 26*3600, 260*3600]
# tau_beta_list = [1e-20*3600, 1e20*3600]
tau_scaling_E = 8 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 8 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 8 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S = 1.5, 1.5, 1.5
rheobases = (rheobase_E, rheobase_P, rheobase_S)
# Background inputs
g_E = 4.5
g_P = 3.2
g_S = 3
back_inputs = (g_E, g_P, g_S)
for tau_beta in tau_beta_list:
for tau_theta in tau_theta_list:
taus = (tau_E, tau_P, tau_S, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Initial conditions for plastic weights
# w_EP_within = 0.81; w_EP_cross = 0.41
# w_ES_within = 0.81; w_ES_cross = 0.31
# w_EE_within = 0.71; w_EE_cross = 0.41
w_EP_within = 0.91; w_EP_cross = 0.41
w_ES_within = 0.51; w_ES_cross = 0.31
w_EE_within = 0.51; w_EE_cross = 0.51
# # Initial conditions for plastic weights
# w_EP_within = 0.7; w_EP_cross = w_EP_within*0.3
# w_ES_within = 0.7; w_ES_cross = w_ES_within*0.3
# w_EE_within = 0.5; w_EE_cross = 0.4
# Weights
w_PE_within = 0.3; w_PE_cross = 0.1
w_PP_within = 0.2; w_PP_cross = 0.1
w_PS_within = 0.3; w_PS_cross = 0.1
w_SE_within = 0.4; w_SE_cross = 0.1
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
J_EE_phase1 = np.zeros((4, n_time_points_stim), dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((10, n_time_points_phase2),
dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
J_phase2 = np.zeros((12, n_time_points_phase2),
dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, r_phase2, r_phase3, max_E)
l_res_weights = (J_EE_phase1, J_phase2)
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_' + str(hour_sim) + 'h' + '_k' + str(K).replace(".","") + '_theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","")
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
#All flags = 0 and simulation is 30 seconds long. It is used to evaluate what happens when activating E1 what's the response of E2. Afterwards it is evaluating the av_threshold given the result
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=(0,0,0,0,0,0), flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=flags, flags_theta=flags_theta)
if save_results:
l_results = [l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates,
l_res_weights,
av_threshold, stim_times, stim_duration, sim_duration]
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates, l_res_weights,
av_threshold, stim_times, stim_duration, sim_duration] = l_results
if plot_results:
print('Plotting the results.')
os.makedirs(dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","") + "/", exist_ok=True)
time_plots([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights, av_threshold,
stim_times, dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","") + "/" + name, hour_sim,modulation_SST, flag_only_S_on=flag_only_S_on, format='.png')
time_plots([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights, av_threshold,
stim_times, dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","")+ "/" + name, hour_sim,modulation_SST, flag_only_S_on=flag_only_S_on, format='.pdf')
def analyze_model_3_compartmental_v3(hour_sim, flags_list, dir_data=r'\figures\data\\', dir_plot=r'\figures\\', modulation_SST=0, run_simulation=True, save_results=False, plot_results=False):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
delta_t = 0.0001 # time step in seconds (0.1 ms)
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially
sim_duration = int(((hour_sim) * 60 * 60 + (stim_duration + 10) * 2 + 2)*(1/delta_t))
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200_000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim)
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2)
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_dend = tau_E # time constant of E population firing rate in seconds(20ms)
tau_hebb = 120 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta = 24 * (60 * 60) # time constant of target activity in seconds(24h)
tau_beta = 28 * (60 * 60) # time constant of target activity regulator in seconds(28h)
tau_scaling_E = 2.5 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 6.5 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 2.5 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_dend, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# g_stim, rheobases, g, K, lambdas = get_model_parameters(coupling='strong', excitation='high')
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(2.5, 0), (0, 2.5)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(1.25, 0), (0, 1.25)])
elif modulation_SST < 0:
g_stim_S = np.array([(-1.25, 0), (0, -1.25)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S, rheobase_A,rheobase_B = 1, 1.5, 1.5, 3, 9
rheobases = (rheobase_E, rheobase_P, rheobase_S, rheobase_A,rheobase_B)
# Background inputs
g_AD = 4
g_BD = 6
g_E = 0
g_P = 4
g_S = 3.2
g = (g_AD, g_BD, g_E, g_P, g_S)
# K parameter in target regulator equation, it tunes the steady state value of target activity and its regulator
K = 0.3
# Constant that define the contribution of each current
lambda_AD = 0.4 # in stronger lambda config it is 0.5
lambda_BD = 0.3 # in stronger lambda config it is 0.3
lambdas = (lambda_AD, lambda_BD)
# Initial conditions for plastic weights
w_EP_within = 0.6; w_EP_cross = 0.18 #PV -> soma
w_DS_within = 0.4; w_DS_cross = 0.18# SST -> apioal
w_DE_within = 0.5; w_DE_cross = 0.3 #E -> apical
w_EE_within = 0.5; w_EE_cross = 0.3 #E -> basal
# Weights
w_PE_within = 0.35; w_PE_cross = 0.10 #E -> PV
w_PP_within = 0.20; w_PP_cross = 0.10 #PV -> PV
w_PS_within = 0.30; w_PS_cross = 0.10 #SST -> PV
w_SE_within = 0.15; w_SE_cross = 0.10 #E -> SST
weights = (w_DE_within, w_EE_within, w_EP_within, w_DS_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_DE_cross, w_EE_cross, w_EP_cross, w_DS_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
I_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # IAD1, IAD2, IBD1, IBD2, IE1, IE2
J_exc_phase1 = np.zeros((8, n_time_points_stim), dtype=np.float32) # WDE11,WDE12,WDE21,WDE22,WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((6, n_time_points_phase2), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
set_phase2 = np.zeros((12, n_time_points_phase2), dtype=np.float32) # thetaDD1,thetaDD2,thetaBD1,thetaBD2,thetaE1,thetaE2,betaAD1,betaAD2,betaBD1,betaBD2,betaE1,betaE2
I_phase2 = np.zeros((6, n_time_points_phase2), dtype=np.float32) # IAD1, IAD2, IBD1, IBD2, IE1, IE2
J_phase2 = np.zeros((20, n_time_points_phase2),
dtype=np.float32) # WDE11,WDE12,WDE21,WDE22,WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, I_phase1, r_phase2, I_phase2, set_phase2, r_phase3, max_E)
l_res_weights = (J_exc_phase1, J_phase2)
# The flags for activating the following plasticity mechanisms in the given order: Hebbian learning, three-factor Hebbian learning,
# adaptive set point, E-to-E scaling, P-to-E scaling, S-to-E scaling
flags_theta = (1, 1)
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_' + str(hour_sim) + 'h' # + '_k' + str(K).replace(".","") + '_td' + str(g_top_down_to_S)
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
model_3_compartmental_v3(delta_t, sampling_rate, l_res_rates, l_res_weights, sim_duration, weights, g, g_stim,
stim_times, taus, K, rheobases, lambdas, flags=flags,flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
if save_results:
l_results = [l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates,
l_res_weights, av_threshold, stim_times, stim_duration, sim_duration]
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[l_time_points_stim, l_time_points_phase2, delta_t, sampling_rate, l_res_rates, l_res_weights,
av_threshold, stim_times, stim_duration, sim_duration] = l_results
if plot_results:
print('Plotting the results.')
plot_all_3_compartmental([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights,
av_threshold, stim_times,modulation_SST, dir_plot + name, hour_sim, format='.png', scale_y=False)
plot_all_3_compartmental([l_time_points_stim, l_time_points_phase2], l_res_rates, l_res_weights,
av_threshold, stim_times, modulation_SST, dir_plot + name, hour_sim, format='.pdf')
def plot_testing_at_regular_intervals(flags_list, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\',
K=0.25, flag_only_S_on=False, run_simulation=True,
save_results = False, plot_results=False,modulation_SST=0):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list
:param flags_theta
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param g_top_down_to_S: Represents the top-down signal to SST neurons triggered by the hyperexcitation. It reaches
SST neurons at the offset of the conditioning
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
hour_sims = np.arange(48) + 1
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_test_every_h' + '_k' + str(K).replace(".","")
l_delta_rE1 = []
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
for hour_sim in hour_sims:
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 5 + 5) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
# l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim)
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2)
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(1, 0), (0, 1)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(0.5, 0), (0, 0.5)])
elif modulation_SST < 0:
g_stim_S = np.array([(-0.5, 0), (0, -0.5)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_hebb = 240 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta = 24 * (60 * 60) # time constant of target activity in seconds(24h)
tau_beta = 28 * (60 * 60) # time constant of target activity regulator in seconds(28h)
tau_scaling_E = 8 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 8 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 8 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S = 1.5, 1.5, 1.5
rheobases = (rheobase_E, rheobase_P, rheobase_S)
# Background inputs
g_E = 4.5
g_P = 3.2
g_S = 3
back_inputs = (g_E, g_P, g_S)
# Initial conditions for plastic weights
# w_EP_within = 0.81; w_EP_cross = 0.41
# w_ES_within = 0.81; w_ES_cross = 0.31
# w_EE_within = 0.71; w_EE_cross = 0.41
w_EP_within = 0.91; w_EP_cross = 0.41
w_ES_within = 0.51; w_ES_cross = 0.31
w_EE_within = 0.51; w_EE_cross = 0.51
# # Initial conditions for plastic weights
# w_EP_within = 0.7; w_EP_cross = w_EP_within*0.3
# w_ES_within = 0.7; w_ES_cross = w_ES_within*0.3
# w_EE_within = 0.5; w_EE_cross = 0.4
# Weights
w_PE_within = 0.3; w_PE_cross = 0.1
w_PP_within = 0.2; w_PP_cross = 0.1
w_PS_within = 0.3; w_PS_cross = 0.1
w_SE_within = 0.4; w_SE_cross = 0.1
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
J_EE_phase1 = np.zeros((4, n_time_points_stim), dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((10, n_time_points_phase2), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
J_phase2 = np.zeros((12, n_time_points_phase2),dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, r_phase2, r_phase3, max_E)
l_res_weights = (J_EE_phase1, J_phase2)
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=(0,0,0,0,0,0), flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=flags, flags_theta=flags_theta)
l_delta_rE1.append(np.max(r_phase3[0][int(stim_times[0][0] * (1 / (delta_t * sampling_rate_stim))):
int(stim_times[0][1] * (1 / (delta_t * sampling_rate_stim)))]).copy())
print('Simulation of ' + str(hour_sim) + ' hours is completed')
if save_results:
l_results = [r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] #added weights for analysis with Kris
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] = l_results
if plot_results:
print('Plotting the results.')
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.png')
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.pdf')
#this function is a generalized version of the one above. this one, with the right flags, is the only one necessary. For clarity, they are separated
def plot_testing_at_regular_intervals_timescales(flags_list, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\',
K=0.25, flag_only_S_on=False, run_simulation=True,
save_results = False, plot_results=False,modulation_SST=0):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list
:param flags_theta
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param g_top_down_to_S: Represents the top-down signal to SST neurons triggered by the hyperexcitation. It reaches
SST neurons at the offset of the conditioning
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
hour_sims = np.arange(48) + 1
tau_theta_list = (np.arange(22, 31,2) * 3600).tolist()
# tau_beta_list = (np.arange(16, 37) * 3600).tolist()
tau_beta_list = [30*3600]
# tau_theta_list = [12*3600]
# tau_beta_list = [0.01*3600]
# tau_theta_list = (np.arange(24-12, 24+12) * 3600).tolist()
# tau_beta_list = (np.arange(28-12, 28+12) * 3600).tolist()
# tau_theta_list = [6*3600, 26*3600, 260*3600]
# tau_beta_list = [1e-20*3600, 1e20*3600]
# tau_theta_list = [6*3600, 26*3600, 260*3600]
# tau_beta_list = [1e-20*3600, 1e20*3600]
for tau_beta in tau_beta_list:
for tau_theta in tau_theta_list:
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_test_every_h' + '_k' + str(K).replace(".","") + '_theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","")
l_delta_rE1 = []
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
for hour_sim in hour_sims:
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 5 + 5) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
# l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim)
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2)
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(1, 0), (0, 1)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(0.5, 0), (0, 0.5)])
elif modulation_SST < 0:
g_stim_S = np.array([(-0.5, 0), (0, -0.5)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_hebb = 240 # time constant of three-factor Hebbian learning in seconds(2min)
tau_scaling_E = 8 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 8 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 8 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S = 1.5, 1.5, 1.5
rheobases = (rheobase_E, rheobase_P, rheobase_S)
# Background inputs
g_E = 4.5
g_P = 3.2
g_S = 3
back_inputs = (g_E, g_P, g_S)
# Initial conditions for plastic weights
# w_EP_within = 0.81; w_EP_cross = 0.41
# w_ES_within = 0.81; w_ES_cross = 0.31
# w_EE_within = 0.71; w_EE_cross = 0.41
w_EP_within = 0.91; w_EP_cross = 0.41
w_ES_within = 0.51; w_ES_cross = 0.31
w_EE_within = 0.51; w_EE_cross = 0.51
# # Initial conditions for plastic weights
# w_EP_within = 0.7; w_EP_cross = w_EP_within*0.3
# w_ES_within = 0.7; w_ES_cross = w_ES_within*0.3
# w_EE_within = 0.5; w_EE_cross = 0.4
# Weights
w_PE_within = 0.3; w_PE_cross = 0.1
w_PP_within = 0.2; w_PP_cross = 0.1
w_PS_within = 0.3; w_PS_cross = 0.1
w_SE_within = 0.4; w_SE_cross = 0.1
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
J_EE_phase1 = np.zeros((4, n_time_points_stim), dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((10, n_time_points_phase2), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
J_phase2 = np.zeros((12, n_time_points_phase2),dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, r_phase2, r_phase3, max_E)
l_res_weights = (J_EE_phase1, J_phase2)
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=(0,0,0,0,0,0), flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=flags, flags_theta=flags_theta)
l_delta_rE1.append(np.max(r_phase3[0][int(stim_times[0][0] * (1 / (delta_t * sampling_rate_stim))):
int(stim_times[0][1] * (1 / (delta_t * sampling_rate_stim)))]).copy())
print('Simulation of ' + str(hour_sim) + ' hours is completed')
if save_results:
l_results = [r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] #added weights for analysis with Kris
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] = l_results
if plot_results:
print('Plotting the results.')
os.makedirs(dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","") + "/", exist_ok=True)
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","") + "/" + name, flag_only_S_on=flag_only_S_on, format='.png')
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + 'theta_' + str(int(tau_theta/3600)).replace(".","") + '_beta_' + str(int(tau_beta/3600)).replace(".","") + "/" + name, flag_only_S_on=flag_only_S_on, format='.pdf')
def plot_testing_at_regular_intervals_dendrites_v3(flags_list, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\', flag_only_S_on=False, modulation_SST=0, run_simulation=True,
save_results = False, plot_results=False):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list
:param flags_theta
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param g_top_down_to_S: Represents the top-down signal to SST neurons triggered by the hyperexcitation. It reaches
SST neurons at the offset of the conditioning
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
hour_sims = np.arange(48) + 1
# hour_sims = np.array([1, 4, 24, 48], dtype=int)
# K parameter in target regulator equation, it tunes the steady state value of target activity and its regulator
K = 0.25
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_test_every_h' + '_k' + str(K).replace(".","") + '_td'
l_delta_rE1 = []
print('*****', title, '*****')
if run_simulation:
print('Simulation started.')
print('\n')
for hour_sim in hour_sims:
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 5 + 5) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
# l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim)
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2)
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(2.5, 0), (0, 2.5)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(1.25, 0), (0, 1.25)])
elif modulation_SST < 0:
g_stim_S = np.array([(-1.25, 0), (0, -1.25)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_dend = tau_E
tau_hebb = 120 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta = 24 * (60 * 60) # time constant of target activity in seconds(24h)
tau_beta = 28 * (60 * 60) # time constant of target activity regulator in seconds(28h)
tau_scaling_E = 2.5 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 6.5 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 2.5 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_dend, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S, rheobase_A,rheobase_B = 1, 1.5, 1.5, 3, 9
rheobases = (rheobase_E, rheobase_P, rheobase_S, rheobase_A,rheobase_B)
# Background inputs
g_AD = 4
g_BD = 6
g_E = 0
g_P = 4
g_S = 3.2
back_inputs = (g_AD, g_BD, g_E, g_P, g_S)
# Constant that define the contribution of each current
lambda_AD = 0.4 # in stronger lambda config it is 0.5
lambda_BD = 0.3 # in stronger lambda config it is 0.3
lambdas = (lambda_AD, lambda_BD)
# Initial conditions for plastic weights
w_EP_within = 0.6; w_EP_cross = 0.18
w_DS_within = 0.4; w_DS_cross = 0.18
w_DE_within = 0.5; w_DE_cross = 0.3
w_EE_within = 0.5; w_EE_cross = 0.3
# Weights
w_PE_within = 0.35; w_PE_cross = 0.10
w_PP_within = 0.20; w_PP_cross = 0.10
w_PS_within = 0.30; w_PS_cross = 0.10
w_SE_within = 0.15; w_SE_cross = 0.10
weights = (w_DE_within, w_EE_within, w_EP_within, w_DS_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_DE_cross, w_EE_cross, w_EP_cross, w_DS_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
I_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # IAD1, IAD2, IBD1, IBD2, IE1, IE2
J_exc_phase1 = np.zeros((8, n_time_points_stim), dtype=np.float32) # WDE11,WDE12,WDE21,WDE22,WEE11,WEE12,WEE21,WEE22
r_phase2 = np.zeros((6, n_time_points_phase2), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
set_phase2 = np.zeros((12, n_time_points_phase2), dtype=np.float32) # thetaDD1,thetaDD2,thetaBD1,thetaBD2,thetaE1,thetaE2,betaAD1,betaAD2,betaBD1,betaBD2,betaE1,betaE2
I_phase2 = np.zeros((6, n_time_points_phase2), dtype=np.float32) # IAD1, IAD2, IBD1, IBD2, IE1, IE2
J_phase2 = np.zeros((20, n_time_points_phase2),
dtype=np.float32) # WDE11,WDE12,WDE21,WDE22,WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
# Lists to hold data arrays
l_res_rates = (r_phase1, I_phase1, r_phase2, I_phase2, set_phase2, r_phase3, max_E)
l_res_weights = (J_exc_phase1, J_phase2)
model_3_compartmental_v3(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights, back_inputs, g_stim,
stim_times, taus, K, rheobases, lambdas, flags=(0,0,0,0,0,0),flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
# av_threshold = r_phase1[1][idx_av_threshold] * 1.15 #it is defined with an extra 15% for old reason. not required anymore
av_threshold = r_phase1[1][idx_av_threshold]
model_3_compartmental_v3(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights, back_inputs, g_stim,
stim_times, taus, K, rheobases, lambdas, flags=flags,flags_theta=flags_theta)
# print(r_phase3)
l_delta_rE1.append(np.max(r_phase3[0][int(stim_times[0][0] * (1 / (delta_t * sampling_rate_stim))):
int(stim_times[0][1] * (1 / (delta_t * sampling_rate_stim)))]).copy())
print('Simulation of ' + str(hour_sim) + ' hours is completed')
if save_results:
l_results = [r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] #added weights for analysis with Kris
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
print('Data is saved.')
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] = l_results
if plot_results:
print('Plotting the results.')
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.png')
change_in_reactivation_every_h_vslides(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.pdf')
def plot_testing_at_regular_intervals_weights(ww_weights,flags_list, plastic_flag, flags_theta=(1,1), dir_data=r'\figures\data\\', dir_plot=r'\figures\\',
K=0.25, flag_only_S_on=False, run_simulation=True,
save_results = False, plot_results=False,modulation_SST=0):
"""
:param hour_sim: Defines how many hours does the simulation lasts
:param flags_list
:param flags_theta
:param dir_data
:param dir_plot
:param K: Tunes the steady state value of target activity and its regulator
:param g_top_down_to_S: Represents the top-down signal to SST neurons triggered by the hyperexcitation. It reaches
SST neurons at the offset of the conditioning
:param run_simulation: True to run the numerical simulation, False to read the already saved data
:param save_results: True to save the results
:param plot_results: True to plot the results
Multi-purpose function to analyze the model. Here we run (if run_simulation is True) our computational model to
investigate the role of cell-type dependent synaptic scaling mechanisms in associative learning. We replicate the
experimental procedure in [1] in model() in model.py. The model has two subnetworks, each consisted of a canonical
circuit of excitatory pyramidal neurons (E), parvalbumin-positive interneurons (P), somatostatin-positive
interneurons (S). The simulation procedure is divided into three phases:
Phase 1 - Conditioning: The first subnetwork receives extra input representing the conditioned stimulus in [1].
The parameters describing the stimulation (when and how much stimulation) is described in analyze_model()
function. The onset response of the excitatory firing rate of the first subnetwork is defined as the aversion
threshold of this network. Three-factor Hebbian learning is active during this period. Also, synaptic scaling
mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 2: In the experiment [1], the novel stimulus is presented to the mice at 4h/24h/48h after conditioning.
This phase corresponds to the waiting time after conditioning and before testing. During this phase, synaptic
scaling mechanisms, adaptive target activity (theta) and target activity regulator (beta) are active.
Phase 3 - Testing: In this phase, the second subnetwork receives extra input corresponds to the novel stimulus
in [1]. The memory specificity/overgeneralization is determined whether the excitatory rate in the second
subnetwork is below/above the aversion threshold, respectively.
During simulation model() writes data to the data arrays. This data can be saved (if save_results is set to True)
and the results can be plotted (if plot_results is set to True).
[1] Wu, C. H., Ramos, R., Katz, D. B., & Turrigiano, G. G. (2021). Homeostatic synaptic scaling establishes the
specificity of an associative memory. Current biology, 31(11), 2274-2285.
"""
hour_sims = np.arange(48) + 1
os.makedirs(dir_data, exist_ok=True)
os.makedirs(dir_plot, exist_ok=True)
# print('Simulation started.')
for flags in flags_list:
id, title = determine_name(flags)
name = 'Case' + id + '_test_every_h' + '_k' + str(K).replace(".","") + "_"
name += '_'.join(str(weight).replace(".","") for weight in ww_weights)
l_delta_rE1 = []
# print('*****', title, '*****')
if run_simulation:
# print('Simulation started.')
# print('\n')
for hour_sim in hour_sims:
stim_duration = 15 # stimulation duration in seconds
# Simulation duration in seconds, 5 extra seconds for pre- and post-stimulation each, 2 extra seconds to reach steady state initially
sim_duration = int((hour_sim) * 60 * 60 + (stim_duration + 5 + 5) * 2 + 2)
delta_t = 0.0001 # time step in seconds (0.1 ms)
sampling_rate_stim = 20 # register data at every 20 step during phase 1 and 3 (conditioning and testing)
sampling_rate_sim = 200000 # register data at every 2e5 time step (20 seconds) during phase 2 (in between conditioning and testing)
sampling_rate = (sampling_rate_stim, sampling_rate_sim)
# Total number of timepoints for stimulation and simulation
n_time_points_stim = int((stim_duration + 10) * (1 / delta_t) * (1 / sampling_rate_stim))
n_time_points_phase2 = int((hour_sim * 60 * 60 - 20) * (1 / delta_t) * (1 / sampling_rate_sim)) + 1 # total no the rest
# l_time_points_stim = np.linspace(0, stim_duration + 10, n_time_points_stim)
l_time_points_phase2 = np.linspace(0, hour_sim, n_time_points_phase2)
# Timepoints of the onset (first column) and offset (second column) of the first (first row) and second (second) stimuli.
stim_times = np.array([[5, 5 + stim_duration],
[int(hour_sim * 60 * 60) + 5, int(hour_sim * 60 * 60) + 5 + stim_duration]]).reshape(2, 2)
# The stimuli are given as inputs to the populations.
g_stim_E = np.array([(1, 0), (0, 1)])
g_stim_P = np.array([(0.5, 0), (0, 0.5)])
if modulation_SST == 0:
g_stim_S = np.array([(0, 0), (0, 0)])
elif modulation_SST > 0:
g_stim_S = np.array([(0.5, 0), (0, 0.5)])
elif modulation_SST < 0:
g_stim_S = np.array([(-0.5, 0), (0, -0.5)])
g_stim = (g_stim_E, g_stim_P, g_stim_S)
# Time constants
tau_E = 0.02 # time constant of E population firing rate in seconds(20ms)
tau_P = 0.005 # time constant of P population firing rate in seconds(5ms)
tau_S = 0.01 # time constant of S population firing rate in seconds(10ms)
tau_hebb = 240 # time constant of three-factor Hebbian learning in seconds(2min)
tau_theta = 24 * (60 * 60) # time constant of target activity in seconds(24h)
tau_beta = 28 * (60 * 60) # time constant of target activity regulator in seconds(28h)
tau_scaling_E = 8 * (60 * 60) # time constant of E-to-E scaling in seconds (15h)
tau_scaling_P = 8 * (60 * 60) # time constant of P-to-E scaling in seconds (15h)
tau_scaling_S = 8 * (60 * 60) # time constant of S-to-E scaling in seconds (15h)
taus = (tau_E, tau_P, tau_S, tau_hebb, tau_scaling_E, tau_scaling_P, tau_scaling_S, tau_theta, tau_beta)
# Rheobases (minimum input needed for firing rates to be above zero)
rheobase_E, rheobase_P, rheobase_S = 1.5, 1.5, 1.5
rheobases = (rheobase_E, rheobase_P, rheobase_S)
# Background inputs
g_E = 4.5
g_P = 3.2
g_S = 3
back_inputs = (g_E, g_P, g_S)
if plastic_flag == True: #I want to modify plastic weights
# Initial conditions for plastic weights
(w_EP_within, w_EP_cross, w_ES_within, w_ES_cross, w_EE_within, w_EE_cross) = ww_weights
# Weights
w_PE_within = 0.3; w_PE_cross = 0.1
w_PP_within = 0.2; w_PP_cross = 0.1
# w_PS_within = 0.3; w_PS_cross = 0.1
# w_SE_within = 0.4; w_SE_cross = 0.1
#strong_connection version
w_PS_within = 0.95; w_PS_cross = 0.1
w_SE_within = 0.1; w_SE_cross = 0.1
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
else:
# Initial conditions for plastic weights
(w_PE_within,w_PP_within,w_PS_within, w_SE_within) = ww_weights
w_PE_cross = 0.1
w_PP_cross = 0.1
w_PS_cross = 0.1
w_SE_cross = 0.1
w_EP_within = 0.91; w_EP_cross = 0.41
w_ES_within = 0.51; w_ES_cross = 0.31
w_EE_within = 0.51; w_EE_cross = 0.51
weights = (w_EE_within, w_EP_within, w_ES_within, w_PE_within, w_PP_within, w_PS_within, w_SE_within,
w_EE_cross, w_EP_cross, w_ES_cross, w_PE_cross, w_PP_cross, w_PS_cross, w_SE_cross)
# Arrays created to hold data
# r_phase1 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
# J_EE_phase1 = np.zeros((4, n_time_points_stim), dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
# r_phase2 = np.zeros((10, n_time_points_phase2), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
# J_phase2 = np.zeros((12, n_time_points_phase2),dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
# r_phase3 = np.zeros((6, n_time_points_stim), dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
max_E = np.zeros(1, dtype=np.float32)
r_phase1 = np.full((6, n_time_points_stim), np.nan, dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
J_EE_phase1 = np.full((4, n_time_points_stim), np.nan, dtype=np.float32) # WEE11,WEE12,WEE21,WEE22
r_phase2 = np.full((10, n_time_points_phase2), np.nan, dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2,theta1,theta2,beta1,beta2
J_phase2 = np.full((12, n_time_points_phase2), np.nan, dtype=np.float32) # WEE11,WEE12,WEE21,WEE22,WEP11,WEP12,WEP21,WEP22,WES11,WES12,WES21,WES22
r_phase3 = np.full((6, n_time_points_stim), np.nan, dtype=np.float32) # rE1,rE2,rP1,rP2,rS1,rS2
# Lists to hold data arrays
l_res_rates = (r_phase1, r_phase2, r_phase3, max_E)
l_res_weights = (J_EE_phase1, J_phase2)
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(30 * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=(0,0,0,0,0,0), flags_theta=flags_theta)
idx_av_threshold = int(15 * (1 / delta_t) * (1 / sampling_rate_stim))
av_threshold = r_phase1[1][idx_av_threshold] * 1.15
model(delta_t, sampling_rate, l_res_rates, l_res_weights, int(sim_duration * (1 / delta_t)), weights,
back_inputs, g_stim, stim_times, taus, K, rheobases, flags=flags, flags_theta=flags_theta)
l_delta_rE1.append(np.max(r_phase3[0][int(stim_times[0][0] * (1 / (delta_t * sampling_rate_stim))):
int(stim_times[0][1] * (1 / (delta_t * sampling_rate_stim)))]).copy())
# print('Simulation of ' + str(hour_sim) + ' hours is completed')
if save_results:
l_results = [r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights]
# Open a file and save
with open(dir_data + name + '.pkl', 'wb') as file:
# A new file will be created
pickle.dump(l_results, file)
else:
# Open the file and read
with open(dir_data + name + '.pkl', 'rb') as file:
l_results = pickle.load(file)
print('Data is read.')
[r_phase1, l_time_points_phase2, r_phase2, l_delta_rE1, av_threshold, delta_t, sampling_rate_sim,l_res_weights] = l_results
#it doesn't go through here
if plot_results:
print('Plotting the results.')
change_in_reactivation_every_h(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.png')
change_in_reactivation_every_h(l_time_points_phase2, hour_sims, l_delta_rE1, av_threshold,
dir_plot + name, flag_only_S_on=flag_only_S_on, format='.pdf')
print("Data for", '_'.join(str(weight).replace(".", "") for weight in ww_weights), "is saved\n")

model_analysis.py at commit 3955ed2, no license · at the source

Overview

Authors: Fabio Veneto^1,2, Ayça Kepçe^2,3, Yue Kris Wu^2,4, Julijana Gjorgjieva^1,2,5

ORCID iDs: Fabio Veneto, Yue Kris Wu, Julijana Gjorgjieva

School of Medicine and Health, Institute for Neuroscience, Technical University of Munich, Munich 81675, Germany
School of Life Sciences, Technical University of Munich, Freising 85354, Germany
Elite Master Program in Neuroengineering, Technical University of Munich, Munich 80333, Germany
Zuckerman Mind Brain Behavior Institute, Columbia University, New York, New York 10027
Munich Cluster for Systems Neurology (SyNergy), Munich 81377, Germany

Institutions: Technical University of Munich (Germany); Mortimer B. Zuckerman Mind Brain Behavior Institute (United States); Columbia University (United States)

Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience, volume 46, issue 27, article e0987252026

Dates: received 20 May 2025; accepted 3 May 2026; published online 28 May 2026; in print 8 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1523/jneurosci.0987-25.2026 · PMID 42209264 · PMCID PMC13421141 · OpenAlex W4410351841

Open access: hybrid, a free copy (OpenAlex)

Preprint: doi:10.1101/2025.05.14.654005

Status: code verified

Categories: cellular / molecular (subfield)

Methods: Single-unit activity, calcium imaging

Keywords: associative learning, circuit model, interneurons, synaptic scaling, top-down modulation

MeSH: Association Learning*, Neuronal Plasticity*, Neurons*, Synapses*, Animals, Memory, Models, Neurological, Parvalbumins, Pyramidal Cells, Somatostatin (* major topic)

Topic: Neuroscience and Neuropharmacology Research (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex

Funding: Deutsche Forschungsgemeinschaft (CRC 1080); EC | Horizon Europe | Excellent Science | HORIZON EUROPE European Research Council (101170267); European Research Council (101170267); Joachim Herz Stiftung

Citations: not cited yet (Europe PMC); 82 references in the paper

Abstract

Excitatory synaptic scaling regulates network dynamics by proportionally adjusting excitatory synaptic strengths after sensory perturbations. During associative learning, blocking excitatory scaling in conditioned taste aversion paradigms prolongs generalized aversive responses and delays memory specificity. Recent evidence also implicates inhibitory synaptic scaling in the regulation of network dynamics. Specifically, parvalbumin (PV)-expressing inhibitory neurons, targeting perisomatic regions of excitatory (E) pyramidal neurons, and somatostatin (SST)-expressing neurons, targeting distal dendrites, exhibit distinct scaling responses. This leaves open the question of how complex plasticity mechanisms regulate recurrent excitatory-inhibitory circuit dynamics in associative learning. Using computational approaches, we demonstrate that Hebbian plasticity drives memory generalization to novel stimuli not presented during conditioning. Following conditioning, diverse synaptic scaling mechanisms progressively induce memory specificity, which can be regulated by top-down inputs. Our results reveal that, in the absence of excitatory scaling, PV-to-E scaling can effectively compensate and rescue memory specificity, highlighting the presence of degenerate mechanisms in the brain. Notably, in the process of establishing memory specificity, excitatory scaling and PV-to-E scaling function synergistically, while concurrently opposing SST-to-E scaling. The synergistic and antagonistic plasticity mechanisms are orchestrated to shape the temporal evolution of memory representations, from generalized to precise.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

comp-neural-circuits/cell-type-specific-synaptic-scaling

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 3955ed23c29e65a6bc2bfe238beb380ce63675e1, 14 November 2025

Languages: Python (7), Jupyter (3), Shell (1)

Size: 63 files, 11 scripts

Software Heritage: not archived

Found in: “Data Availability”

Holds: README, 3 notebooks

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: NumPy (8 files), Matplotlib (6 files), seaborn (4 files), Numba (1 file), pandas (1 file), SciPy (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

12 files

analytics.ipynb, Jupyter, 1,600 lines
main_for_paper.py, Python, 261 lines
model.py, Python, 560 lines, 1 match
model_analysis.py, Python, 1,331 lines, 2 matches
parameter_exploration.py, Python, 166 lines
parameter_generator_syn_scaling.ipynb, Jupyter, 170 lines
plotting_functions.py, Python, 1,705 lines
plotting_functions_analytics.py, Python, 541 lines
robustness.ipynb, Jupyter, 705 lines
run_simulation_efficient_syn_scal.sh, Shell, 69 lines
util.py, Python, 2,170 lines
README.md, Text, 5 lines

At the source: github.com/comp-neural-circuits/cell-type-specific-synaptic-scaling

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
11 scripts, each with its path and the digest of its content;
3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

The simulation code is publicly available at https://github.com/comp-neural-circuits/cell-type-specific-synaptic-scaling.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 5 keywords, 10 MeSH terms, 4 funders, 76 references.

Cite

This paper

Veneto, F., Kepçe, A., Wu, Y. K., & Gjorgjieva, J. (2026). Cell-Type-Specific Synaptic Scaling Mechanisms Differentially Contribute to Associative Learning. The Journal of neuroscience : the official journal of the Society for Neuroscience, 46(27), e0987252026. https://doi.org/10.1523/jneurosci.0987-25.2026

BibTeX

@article{veneto2026cell, author = {Veneto, Fabio and Kepçe, Ayça and Wu, Yue Kris and Gjorgjieva, Julijana}, title = {{Cell-Type-Specific Synaptic Scaling Mechanisms Differentially Contribute to Associative Learning}}, journal = {The Journal of neuroscience : the official journal of the Society for Neuroscience}, year = {2026}, month = jul, volume = {46}, number = {27}, pages = {e0987252026}, publisher = {Society for Neuroscience}, issn = {0270-6474}, doi = {10.1523/jneurosci.0987-25.2026}, url = {https://doi.org/10.1523/jneurosci.0987-25.2026}, pmid = {42209264}, pmcid = {PMC13421141} }

RIS

TY - JOUR AU - Veneto, Fabio AU - Kepçe, Ayça AU - Wu, Yue Kris AU - Gjorgjieva, Julijana TI - Cell-Type-Specific Synaptic Scaling Mechanisms Differentially Contribute to Associative Learning T2 - The Journal of neuroscience : the official journal of the Society for Neuroscience J2 - J Neurosci PY - 2026 DA - 2026/07/08 VL - 46 IS - 27 SP - e0987252026 SN - 0270-6474 PB - Society for Neuroscience DO - 10.1523/jneurosci.0987-25.2026 UR - https://doi.org/10.1523/jneurosci.0987-25.2026 LA - en ER -

CSL-JSON

{ "id": "10.1523/jneurosci.0987-25.2026", "type": "article-journal", "title": "Cell-Type-Specific Synaptic Scaling Mechanisms Differentially Contribute to Associative Learning", "container-title": "The Journal of neuroscience : the official journal of the Society for Neuroscience", "author": [ { "family": "Veneto", "given": "Fabio" }, { "family": "Kepçe", "given": "Ayça" }, { "family": "Wu", "given": "Yue Kris" }, { "family": "Gjorgjieva", "given": "Julijana" } ], "container-title-short": "J Neurosci", "volume": "46", "issue": "27", "page": "e0987252026", "DOI": "10.1523/jneurosci.0987-25.2026", "PMID": "42209264", "PMCID": "PMC13421141", "ISSN": "0270-6474", "publisher": "Society for Neuroscience", "URL": "https://doi.org/10.1523/jneurosci.0987-25.2026", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 8 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.