1# -*- coding: utf-8 -*- 2# 3# urbanczik_synapse_example.py 4# 5# This file is part of NEST. 6# 7# Copyright (C) 2004 The NEST Initiative 8# 9# NEST is free software: you can redistribute it and/or modify 10# it under the terms of the GNU General Public License as published by 11# the Free Software Foundation, either version 2 of the License, or 12# (at your option) any later version. 13# 14# NEST is distributed in the hope that it will be useful, 15# but WITHOUT ANY WARRANTY; without even the implied warranty of 16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17# GNU General Public License for more details. 18# 19# You should have received a copy of the GNU General Public License 20# along with NEST. If not, see <http://www.gnu.org/licenses/>. 21 22""" 23Weight adaptation according to the Urbanczik-Senn plasticity 24------------------------------------------------------------ 25 26This script demonstrates the learning in a compartmental neuron where the 27dendritic synapses adapt their weight according to the plasticity rule by 28Urbanczik and Senn [1]_. In this simple setup, a spike pattern of 200 poisson 29spike trains is repeatedly presented to a neuron that is composed of one 30somatic and one dendritic compartment. At the same time, the somatic 31conductances are activated to produce a time-varying matching potential. 32After the learning, this signal is then reproreproduced by the membrane 33potential of the neuron. This script produces Fig. 1B in [1]_ but uses standard 34units instead of the unitless quantities used in the paper. 35 36References 37~~~~~~~~~~ 38 39.. [1] R. Urbanczik, W. Senn (2014): Learning by the Dendritic Prediction of 40 Somatic Spiking. Neuron, 81, 521-528. 41""" 42import numpy as np 43from matplotlib import pyplot as plt 44import nest 45 46 47def g_inh(amplitude, t_start, t_end): 48 """ 49 returns weights for the spike generator that drives the inhibitory 50 somatic conductance. 51 """ 52 return lambda t: np.piecewise(t, [(t >= t_start) & (t < t_end)], 53 [amplitude, 0.0]) 54 55 56def g_exc(amplitude, freq, offset, t_start, t_end): 57 """ 58 returns weights for the spike generator that drives the excitatory 59 somatic conductance. 60 """ 61 return lambda t: np.piecewise(t, [(t >= t_start) & (t < t_end)], 62 [lambda t: amplitude * np.sin(freq * t) + offset, 0.0]) 63 64 65def matching_potential(g_E, g_I, nrn_params): 66 """ 67 returns the matching potential as a function of the somatic conductances. 68 """ 69 E_E = nrn_params['soma']['E_ex'] 70 E_I = nrn_params['soma']['E_in'] 71 return (g_E * E_E + g_I * E_I) / (g_E + g_I) 72 73 74def V_w_star(V_w, nrn_params): 75 """ 76 returns the dendritic prediction of the somatic membrane potential. 77 """ 78 g_D = nrn_params['g_sp'] 79 g_L = nrn_params['soma']['g_L'] 80 E_L = nrn_params['soma']['E_L'] 81 return (g_L * E_L + g_D * V_w) / (g_L + g_D) 82 83 84def phi(U, nrn_params): 85 """ 86 rate function of the soma 87 """ 88 phi_max = nrn_params['phi_max'] 89 k = nrn_params['rate_slope'] 90 beta = nrn_params['beta'] 91 theta = nrn_params['theta'] 92 return phi_max / (1.0 + k * np.exp(beta * (theta - U))) 93 94 95def h(U, nrn_params): 96 """ 97 derivative of the rate function phi 98 """ 99 k = nrn_params['rate_slope'] 100 beta = nrn_params['beta'] 101 theta = nrn_params['theta'] 102 return 15.0 * beta / (1.0 + np.exp(-beta * (theta - U)) / k) 103 104 105# simulation params 106n_pattern_rep = 100 # number of repetitions of the spike pattern 107pattern_duration = 200.0 108t_start = 2.0 * pattern_duration 109t_end = n_pattern_rep * pattern_duration + t_start 110simulation_time = t_end + 2.0 * pattern_duration 111n_rep_total = int(np.around(simulation_time / pattern_duration)) 112resolution = 0.1 113nest.resolution = resolution 114 115# neuron parameters 116nrn_model = 'pp_cond_exp_mc_urbanczik' 117nrn_params = { 118 't_ref': 3.0, # refractory period 119 'g_sp': 600.0, # soma-to-dendritic coupling conductance 120 'soma': { 121 'V_m': -70.0, # initial value of V_m 122 'C_m': 300.0, # capacitance of membrane 123 'E_L': -70.0, # resting potential 124 'g_L': 30.0, # somatic leak conductance 125 'E_ex': 0.0, # resting potential for exc input 126 'E_in': -75.0, # resting potential for inh input 127 'tau_syn_ex': 3.0, # time constant of exc conductance 128 'tau_syn_in': 3.0, # time constant of inh conductance 129 }, 130 'dendritic': { 131 'V_m': -70.0, # initial value of V_m 132 'C_m': 300.0, # capacitance of membrane 133 'E_L': -70.0, # resting potential 134 'g_L': 30.0, # dendritic leak conductance 135 'tau_syn_ex': 3.0, # time constant of exc input current 136 'tau_syn_in': 3.0, # time constant of inh input current 137 }, 138 # parameters of rate function 139 'phi_max': 0.15, # max rate 140 'rate_slope': 0.5, # called 'k' in the paper 141 'beta': 1.0 / 3.0, 142 'theta': -55.0, 143} 144 145# synapse params 146syns = nest.GetDefaults(nrn_model)['receptor_types'] 147init_w = 0.3 * nrn_params['dendritic']['C_m'] 148syn_params = { 149 'synapse_model': 'urbanczik_synapse_wr', 150 'receptor_type': syns['dendritic_exc'], 151 'tau_Delta': 100.0, # time constant of low pass filtering of the weight change 152 'eta': 0.17, # learning rate 153 'weight': init_w, 154 'Wmax': 4.5 * nrn_params['dendritic']['C_m'], 155 'delay': resolution, 156} 157 158""" 159# in case you want to use the unitless quantities as in [1]: 160 161# neuron params: 162 163nrn_model = 'pp_cond_exp_mc_urbanczik' 164nrn_params = { 165 't_ref': 3.0, 166 'g_sp': 2.0, 167 'soma': { 168 'V_m': 0.0, 169 'C_m': 1.0, 170 'E_L': 0.0, 171 'g_L': 0.1, 172 'E_ex': 14.0 / 3.0, 173 'E_in': -1.0 / 3.0, 174 'tau_syn_ex': 3.0, 175 'tau_syn_in': 3.0, 176 }, 177 'dendritic': { 178 'V_m': 0.0, 179 'C_m': 1.0, 180 'E_L': 0.0, 181 'g_L': 0.1, 182 'tau_syn_ex': 3.0, 183 'tau_syn_in': 3.0, 184 }, 185 # parameters of rate function 186 'phi_max': 0.15, 187 'rate_slope': 0.5, 188 'beta': 5.0, 189 'theta': 1.0, 190} 191 192# synapse params: 193 194syns = nest.GetDefaults(nrn_model)['receptor_types'] 195init_w = 0.2*nrn_params['dendritic']['g_L'] 196syn_params = { 197 'synapse_model': 'urbanczik_synapse_wr', 198 'receptor_type': syns['dendritic_exc'], 199 'tau_Delta': 100.0, 200 'eta': 0.0003 / (15.0*15.0*nrn_params['dendritic']['C_m']), 201 'weight': init_w, 202 'Wmax': 3.0*nrn_params['dendritic']['g_L'], 203 'delay': resolution, 204} 205""" 206 207# somatic input 208ampl_exc = 0.016 * nrn_params['dendritic']['C_m'] 209offset = 0.018 * nrn_params['dendritic']['C_m'] 210ampl_inh = 0.06 * nrn_params['dendritic']['C_m'] 211freq = 2.0 / pattern_duration 212soma_exc_inp = g_exc(ampl_exc, 2.0 * np.pi * freq, offset, t_start, t_end) 213soma_inh_inp = g_inh(ampl_inh, t_start, t_end) 214 215# dendritic input 216# create spike pattern by recording the spikes of a simulation of n_pg 217# poisson generators. The recorded spike times are then given to spike 218# generators. 219n_pg = 200 # number of poisson generators 220p_rate = 10.0 # rate in Hz 221 222pgs = nest.Create('poisson_generator', n=n_pg, params={'rate': p_rate}) 223 224prrt_nrns_pg = nest.Create('parrot_neuron', n_pg) 225 226nest.Connect(pgs, prrt_nrns_pg, {'rule': 'one_to_one'}) 227 228sr = nest.Create('spike_recorder', n_pg) 229nest.Connect(prrt_nrns_pg, sr, {'rule': 'one_to_one'}) 230 231nest.Simulate(pattern_duration) 232t_srs = [ssr.get('events', 'times') for ssr in sr] 233 234nest.ResetKernel() 235nest.resolution = resolution 236 237""" 238neuron and devices 239""" 240nrn = nest.Create(nrn_model, params=nrn_params) 241 242# poisson generators are connected to parrot neurons which are 243# connected to the mc neuron 244prrt_nrns = nest.Create('parrot_neuron', n_pg) 245 246# excitatory input to the soma 247spike_times_soma_inp = np.arange(resolution, simulation_time, resolution) 248sg_soma_exc = nest.Create('spike_generator', 249 params={'spike_times': spike_times_soma_inp, 250 'spike_weights': soma_exc_inp(spike_times_soma_inp)}) 251# inhibitory input to the soma 252sg_soma_inh = nest.Create('spike_generator', 253 params={'spike_times': spike_times_soma_inp, 254 'spike_weights': soma_inh_inp(spike_times_soma_inp)}) 255 256# excitatory input to the dendrite 257sg_prox = nest.Create('spike_generator', n=n_pg) 258 259# for recording all parameters of the Urbanczik neuron 260rqs = nest.GetDefaults(nrn_model)['recordables'] 261mm = nest.Create('multimeter', params={'record_from': rqs, 'interval': 0.1}) 262 263# for recoding the synaptic weights of the Urbanczik synapses 264wr = nest.Create('weight_recorder') 265 266# for recording the spiking of the soma 267sr_soma = nest.Create('spike_recorder') 268 269 270# create connections 271nest.Connect(sg_prox, prrt_nrns, {'rule': 'one_to_one'}) 272nest.CopyModel('urbanczik_synapse', 'urbanczik_synapse_wr', 273 {'weight_recorder': wr[0]}) 274nest.Connect(prrt_nrns, nrn, syn_spec=syn_params) 275nest.Connect(mm, nrn, syn_spec={'delay': 0.1}) 276nest.Connect(sg_soma_exc, nrn, 277 syn_spec={'receptor_type': syns['soma_exc'], 'weight': 10.0 * resolution, 'delay': resolution}) 278nest.Connect(sg_soma_inh, nrn, 279 syn_spec={'receptor_type': syns['soma_inh'], 'weight': 10.0 * resolution, 'delay': resolution}) 280nest.Connect(nrn, sr_soma) 281 282# simulation divided into intervals of the pattern duration 283for i in np.arange(n_rep_total): 284 # Set the spike times of the pattern for each spike generator 285 for (sg, t_sp) in zip(sg_prox, t_srs): 286 nest.SetStatus( 287 sg, {'spike_times': np.array(t_sp) + i * pattern_duration}) 288 289 nest.Simulate(pattern_duration) 290 291 292# read out devices 293 294# multimeter 295mm_events = mm.events 296t = mm_events['times'] 297V_s = mm_events['V_m.s'] 298V_d = mm_events['V_m.p'] 299V_d_star = V_w_star(V_d, nrn_params) 300g_in = mm_events['g_in.s'] 301g_ex = mm_events['g_ex.s'] 302I_ex = mm_events['I_ex.p'] 303I_in = mm_events['I_in.p'] 304U_M = matching_potential(g_ex, g_in, nrn_params) 305 306# weight recorder 307wr_events = wr.events 308senders = wr_events['senders'] 309targets = wr_events['targets'] 310weights = wr_events['weights'] 311times = wr_events['times'] 312 313# spike recorder 314spike_times_soma = sr_soma.get('events', 'times') 315 316 317# plot results 318fs = 22 319lw = 2.5 320fig1, (axA, axB, axC, axD) = plt.subplots(4, 1, sharex=True) 321 322# membrane potentials and matching potential 323axA.plot(t, V_s, lw=lw, label=r'$U$ (soma)', color='darkblue') 324axA.plot(t, V_d, lw=lw, label=r'$V_W$ (dendrit)', color='deepskyblue') 325axA.plot(t, V_d_star, lw=lw, label=r'$V_W^\ast$ (dendrit)', color='b', ls='--') 326axA.plot(t, U_M, lw=lw, label=r'$U_M$ (soma)', color='r', ls='-') 327axA.set_ylabel('membrane pot [mV]', fontsize=fs) 328axA.legend(fontsize=fs) 329 330# somatic conductances 331axB.plot(t, g_in, lw=lw, label=r'$g_I$', color='r') 332axB.plot(t, g_ex, lw=lw, label=r'$g_E$', color='coral') 333axB.set_ylabel('somatic cond', fontsize=fs) 334axB.legend(fontsize=fs) 335 336# dendritic currents 337axC.plot(t, I_ex, lw=lw, label=r'$I_ex$', color='r') 338axC.plot(t, I_in, lw=lw, label=r'$I_in$', color='coral') 339axC.set_ylabel('dend current', fontsize=fs) 340axC.legend(fontsize=fs) 341 342# rates 343axD.plot(t, phi(V_s, nrn_params), lw=lw, label=r'$\phi(U)$', color='darkblue') 344axD.plot(t, phi(V_d, nrn_params), lw=lw, 345 label=r'$\phi(V_W)$', color='deepskyblue') 346axD.plot(t, phi(V_d_star, nrn_params), lw=lw, 347 label=r'$\phi(V_W^\ast)$', color='b', ls='--') 348axD.plot(t, h(V_d_star, nrn_params), lw=lw, 349 label=r'$h(V_W^\ast)$', color='g', ls='--') 350axD.plot(t, phi(V_s, nrn_params) - phi(V_d_star, nrn_params), lw=lw, 351 label=r'$\phi(U) - \phi(V_W^\ast)$', color='r', ls='-') 352axD.plot(spike_times_soma, 0.0 * np.ones(len(spike_times_soma)), 353 's', color='k', markersize=2) 354axD.legend(fontsize=fs) 355 356# synaptic weights 357fig2, axA = plt.subplots(1, 1) 358for i in np.arange(2, 200, 10): 359 index = np.intersect1d(np.where(senders == i), np.where(targets == 1)) 360 if not len(index) == 0: 361 axA.step(times[index], weights[index], label='pg_{}'.format(i - 2), 362 lw=lw) 363 364axA.set_title('Synaptic weights of Urbanczik synapses') 365axA.set_xlabel('time [ms]', fontsize=fs) 366axA.set_ylabel('weight', fontsize=fs) 367axA.legend(fontsize=fs - 4) 368plt.show() 369