123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383 |
- # code for simulating size tuning curves using the edog model
- # import
- import matplotlib.pyplot as plt
- from operator import itemgetter
- from edog.tools import*
- import spatint_utils
- from importlib import reload
- import os
- # reload modules
- reload(spatint_utils)
- # load params
- spatint_utils.plot_params()
- _, _, optocolor = spatint_utils.get_colors()
- class SimEdog:
- def __init__(self, start=1, stop=41):
- """Init class
- Parameters
- -----
- start: int
- starting width for inhibitory feedback kernel
- stop: int
- end width for inhibitory feedback kernel width + 1
- """
- self.start = start
- self.stop = stop
- self.step = 1
- self.lowbounds = None
- self.highbounds = None
- self.sids = None
- self.rfcsds = None
- self.smallrateds = None
- self.largerateds = None
- self.siss = None
- self.rfcsss = None
- self.wavenumber = None
- self.patch_diameter = None
- self.w_rc_mix = None
- def describe(self):
- """Get parameter boundaries in which experimental results are qualitatively replicated.
- Creates two dictionaries, one for the low boundaries, one for the high boundaries. Each
- dict contains one key for each parameter (smallrateds: response to preferred size;
- largerateds: response to largest stimulus; sids: surround suppression; rfcsds:
- preferred size) plus one key containing the intersection of all parameter boundaries
- (all)."""
- if self.sids is None:
- # run simulations if this has not been done
- self.iterate()
- # init boundary dicts
- lowbounds = {}
- highbounds = {}
- # check for which kernel size CT feedback increases responses to optimal stimulus size
- smallrateds_increaseis = (np.where(np.array(self.smallrateds) > 0))[0] + self.start
- lowbounds['smallrateds'] = np.min(smallrateds_increaseis)
- highbounds['smallrateds'] = np.max(smallrateds_increaseis)
- print('Responses to optimal stimulus diameter inreases with feedback for kernel '
- 'widths between %d and %d degree.'
- % (lowbounds['smallrateds'], highbounds['smallrateds']))
- # check for which kernel size CT feedback decreases responses to largest stimulus size
- largerateds_decreaseis = (np.where(np.array(self.largerateds) < 0))[0] + self.start
- lowbounds['largerateds'] = np.min(largerateds_decreaseis)
- highbounds['largerateds'] = np.max(largerateds_decreaseis)
- print('Responses to largest stimulus diameter decreases with feedback for kernel '
- 'widths between %d and %d degree.'
- % (lowbounds['largerateds'], highbounds['largerateds']))
- # check for which kernel size CT feedback increases surround suppression
- sids_increaseis = (np.where(np.array(self.sids) > 0))[0] + self.start
- lowbounds['sids'] = np.min(sids_increaseis)
- highbounds['sids'] = np.max(sids_increaseis)
- print('SI increases with feedback for kernel widths between %d and %d degree.'
- % (lowbounds['sids'], highbounds['sids']))
- # check for which kernel size CT feedback decreases receptive field center size
- rfcsds_decreaseis = (np.where(np.array(self.rfcsds) < 0))[0] + self.start
- lowbounds['rfcsds'] = np.min(rfcsds_decreaseis)
- highbounds['rfcsds'] = np.max(rfcsds_decreaseis)
- print('RFCS decreases with feedback for kernel widths between %d and %d degree.'
- % (lowbounds['rfcsds'], highbounds['rfcsds']))
- # check range for all paramters
- all_lowi = max(lowbounds.items(), key=itemgetter(1))[0]
- all_low = lowbounds[all_lowi]
- all_highi = min(highbounds.items(), key=itemgetter(1))[0]
- all_high = highbounds[all_highi]
- lowbounds['all'] = all_low
- highbounds['all'] = all_high
- print('Qualitatively matching range between %d and %d degree.'
- % (lowbounds['all'], highbounds['all']))
- # store dicts as attributes
- self.lowbounds = lowbounds
- self.highbounds = highbounds
- def iterate(self):
- """Run size tuning simulation for range of inhibitory FB kernel widths specified for
- conditions in which feedback is intact (weight: 1) and feedback is abolished (
- weight: 0). For each stimulation store percent change in surround suppression (
- sids); preferred size (rfcsds), response to preferred size (smallrateds),
- and response to the largest stimulus size (largerateds), as well as absolute values
- for suppression index (siss) and preferred size (rfcsss) in both conditions"""
- # init lists to store values
- sids = []
- rfcsds = []
- smallrateds = []
- largerateds = []
- siss = []
- rfcsss = []
- # create vector of inhibitory feedback widths
- inhib_fbs = np.arange(self.start, self.stop, self.step)
- # iterate over widths
- for inhib_fb in inhib_fbs:
- # run simulation
- curves = self.simulate(inhib_fb=inhib_fb)
- # compute target values
- sid, rfcsd, smallrated, largerated, sis, rfcss = self.characterize(curves=curves)
- # store values in list
- sids.append(sid)
- rfcsds.append(rfcsd)
- smallrateds.append(smallrated)
- largerateds.append(largerated)
- siss.append(sis)
- rfcsss.append(rfcss)
- # store as attributes
- self.sids = sids
- self.rfcsds = rfcsds
- self.smallrateds = smallrateds
- self.largerateds = largerateds
- self.siss = siss
- self.rfcsss = rfcsss
- def simulate(self, inhib_fb=None):
- """Run the edog model and return the two tuning curves
- Parameters
- ----
- inhib_fb: int
- width of inhibitory feedback kernel
- """
- # get parameters
- parentdir = os. path. dirname(os. getcwd())
- filename = parentdir + '/data/params_mouse.yaml'
- params = parse_parameters(filename)
- nt, nr, dt, dr = itemgetter("nt", "nr", "dt", "dr")(params["grid"])
- k_id, w_id, patch_diameter = itemgetter("k_id", "w_id", "patch_diameter")(
- params["stimulus"])
- A_g, a_g, B_g, b_g = itemgetter("A", "a", "B", "b")(params["ganglion"])
- w_rg, A_rg, a_rg = itemgetter("w", "A", "a")(params["relay"]["Krg"])
- w_rig, A_rig, a_rig = itemgetter("w", "A", "a")(params["relay"]["Krig"])
- w_rc_mix = itemgetter("w")(params["relay"]["Krc_mix"])
- A_rc_mix_in, a_rc_mix_in = itemgetter("A", "a")(params["relay"]["Krc_mix"]["Krc_in"])
- A_rc_mix_ex, a_rc_mix_ex = itemgetter("A", "a")(params["relay"]["Krc_mix"]["Krc_ex"])
- if inhib_fb is not None:
- # if input to inhib_fb is not none replace default value of inhibitory FB kernel
- # width
- a_rc_mix_in = inhib_fb * pq.deg
- # initiate variables
- tuning_curve = np.zeros([len(w_rc_mix), len(patch_diameter)])
- cen_size = np.zeros(len(w_rc_mix))
- supp_index = np.zeros(len(w_rc_mix))
- tuning_curve[:] = np.nan
- cen_size[:] = np.nan
- supp_index[:] = np.nan
- # run model
- for i, w in enumerate(w_rc_mix):
- network = create_spatial_network(nt=nt, nr=nr, dt=dt, dr=dr,
- A_g=A_g, a_g=a_g, B_g=B_g, b_g=b_g,
- w_rg=w_rg, A_rg=A_rg, a_rg=a_rg,
- w_rig=w_rig, A_rig=A_rig, a_rig=a_rig,
- w_rc_in=w, A_rc_in=A_rc_mix_in,
- a_rc_in=a_rc_mix_in, w_rc_ex=w,
- A_rc_ex=A_rc_mix_ex, a_rc_ex=a_rc_mix_ex)
- angular_freq = network.integrator.temporal_angular_freqs[int(w_id)]
- wavenumber = network.integrator.spatial_angular_freqs[int(k_id)]
- spatiotemporal_tuning = spatiotemporal_size_tuning(network=network,
- angular_freq=angular_freq,
- wavenumber=wavenumber,
- patch_diameter=patch_diameter)
- # store tuning curve, preferred size and suppression index
- tuning_curve[i, :] = spatiotemporal_tuning[0, :]
- # normalize by no feedback condition
- curves = tuning_curve / tuning_curve[0].max()
- # store anglular frequencies, FB weighs and patch diameters
- self.wavenumber = network.integrator.spatial_angular_freqs[int(k_id)]
- self.patch_diameter = patch_diameter
- self.w_rc_mix = w_rc_mix
- return curves
- def characterize(self, curves):
- """Compute differences between size tuning curves with and without feedback. Return
- percent change in surround suppression (sid); preferred size (rfcsd), response to
- preferred size (smallrated), and response to the largest stimulus size (largerated),
- as well as absolute values for suppression index (sis) and preferred size (rfcss)
- Parameters
- ----
- curves: ndarray
- 2D array containing size tuning curves in both conditions
- Returns
- ----
- sid: float
- percent change in suppression index
- rfcsd: float
- percent change in preferred size
- smallrated: float
- percent change in response to small stimulus
- largerated: float
- percent change in response to large stimulus
- sis: list
- suppression indices under both condtitions
- rfcss: list
- preferred sizes under both conditions
- """
- # get patch_diameters
- patch_diameter = self.patch_diameter
- # init lists to store target values
- sis = []
- rfcss = []
- smallrates = []
- largerates = []
- # preferred size is defined by stimulus where an increase in one degree does not
- # increase response by 0.05%. This is analogous to how visTRN size tuning is analyzed
- perc_thres = 0.05
- # iterate over the two conditions
- for icurve in range(curves.shape[0]):
- # get x and y
- y = curves[1-icurve, :] # start with control condition
- x = patch_diameter
- # compute RF center size: check where 1deg in increase of visual stim fails to
- # increase resp by perc_thres (output 2 to deal with divide by 0 -> will set
- # first increase to 100%)
- rel_change = np.divide(y[1:], y[:-1], out=np.ones_like(y[1:])*2, where=y[:-1] != 0)
- # compute percent change
- perc_change = (rel_change - 1) * 100
- # check where percent change is lower than perc_thres
- max_resp = np.where(perc_change < perc_thres)[0]
- # receptive field center size is either first max_resp or last stimulus size
- if any(max_resp):
- rfcs = max_resp[0]
- else:
- rfcs = x[-1]
- if icurve == 0:
- # store rfcs of control curve
- cont_rfcs = rfcs
- # compute surround size
- rfss = int(x[-1].item())
- # determine modeled firing rates
- smallrate_cont = y[cont_rfcs]
- smallrate = y[rfcs]
- largerate = y[rfss]
- # compute suppression index, if si is smaller 0 because of small increase for large
- # stimuli, set it to 0
- si = max(((smallrate - largerate) / smallrate), 0)
- # store values
- sis.append(si)
- rfcss.append(rfcs)
- smallrates.append(smallrate_cont)
- largerates.append(largerate)
- # compute percent change
- sid = ((sis[0] / sis[1]) - 1) * 100
- rfcsd = ((rfcss[0] / rfcss[1]) - 1) * 100
- smallrated = ((smallrates[0] / smallrates[1]) - 1) * 100
- largerated = ((largerates[0] / largerates[1]) - 1) * 100
- return sid, rfcsd, smallrated, largerated, sis, rfcss
- def plot(self, figsize=(6, 6), inhib_fb=None, quant=False, ax=None):
- """Plot the two tuning curves for inhibitory feedback kernel width
- Parameters
- ----
- figsize: tuple
- figure size (width, height)
- inhib_fb: int
- inhibitory feedback kernel width
- quant: bool
- if true print info regarding
- ax: mpl axis
- axis for plotting
- """
- # simulate data
- curves = self.simulate(inhib_fb=inhib_fb)
- # get xvalues
- patch_diameter = self.patch_diameter
- # set colors
- colors = [optocolor, 'k']
- if ax is None:
- # init fig if doesnt exist
- fig, ax = plt.subplots(figsize=figsize)
- # plot curves
- for curve, color in zip(curves, colors):
- ax.plot(patch_diameter, curve, '-', color=color)
- # layout
- ax.set_ylabel("Normalized response")
- ax.set_xlabel("Diameter ($\degree$)")
- ax.set_xticks((0, 25, 50, 75))
- ax.set_xticklabels((0, 25, 50, 75))
- ax.set_xlim([0, 75])
- ax.set_ylim([0, 1.5])
- ax.spines['bottom'].set_bounds(0, 75)
- if quant:
- # plot information about differences
- # init variables
- title = "Mixed feedback"
- wavenumber = self.wavenumber
- # get quant measures
- sid, rfcsd, smallrated, largerated, sis, rfcss = self.characterize(curves=curves)
- # insert model results
- ax.text(50, 0.1, 'opto si = %0.2f' % sis[1])
- ax.text(50, 0.14, 'control si = %0.2f' % sis[0])
- ax.text(50, 0.18, 'perc change si = %0.2f' % sid)
- ax.text(50, 0.22, 'opto rfcs = %0.2f' % rfcss[1])
- ax.text(50, 0.26, 'control rfcs = %0.2f' % rfcss[0])
- ax.text(50, 0.3, 'perc change rfcs = %0.2f' % rfcsd)
- ax.text(50, 0.34, 'perc change smallrate = %0.2f' % smallrated)
- ax.text(50, 0.38, 'perc change largerate = %0.2f' % largerated)
- ax.text(50, 0.42, 'model')
- # add legends and titles
- ax.set_title(title)
- fig.text(0.39, 0.97, "Patch grating ({})".format(round(wavenumber.item(), 2)),
- fontsize=12)
- # plot preferred size
- for rfcs, color in zip(rfcss[::-1], colors):
- ax.axvline(rfcs, color=color)
|