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