corticothalamic_spatial_integration/figs/sim_edog.py

# 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)