Mechanisms for stable, robust, and adaptive development of orientation maps (Stevens et al. 2013)

 Download zip file 
Help downloading and running models
Accession:151951
GCAL (Gain Control, Adaptation, Laterally connected). Simple but robust single-population V1 orientation map model.
Reference:
1 . Stevens JL, Law JS, AntolĂ­k J, Bednar JA (2013) Mechanisms for stable, robust, and adaptive development of orientation maps in the primary visual cortex. J Neurosci 33:15747-66 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Connectionist Network;
Brain Region(s)/Organism: Neocortex; Thalamus; Retina;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: Topographica;
Model Concept(s): Development; Homeostasis; Orientation selectivity; Contrast-gain control; Vision;
Implementer(s): Bednar, James [jbednar at alumni.utexas.net]; Stevens, Jean-Luc R ; Law, Judith ; Antolik, Jan ;
/
StevensEtAl2013
README.txt
gcal.ty
                            
# -*- mode: python;-*-
"""
GCAL: Gain Control, Adaptation, Laterally connected

Simple but robust single-population V1 model orientation map from:

   Jean-Luc R. Stevens, Judith S. Law, Jan Antolik, and James A. Bednar.
   Mechanisms for stable, robust, and adaptive development of orientation
   maps in the primary visual cortex.
   Journal of Neuroscience, 33:15747-15766, 2013.

   Development of orientation maps in ferret and cat primary visual
   cortex (V1) has been shown to be stable, in that the earliest
   measurable maps are similar in form to the eventual adult map,
   robust, in that similar maps develop in both dark rearing and in a
   variety of normal visual environments, and yet adaptive, in that
   the final map pattern reflects the statistics of the specific
   visual environment. How can these three properties be reconciled?
   Using mechanistic models of the development of neural connectivity
   in V1, we show for the first time that realistic stable, robust,
   and adaptive map development can be achieved by including two
   low-level mechanisms originally motivated from single-neuron
   results. Specifically, contrast-gain control in the retinal
   ganglion cells and the lateral geniculate nucleus reduces variation
   in the presynaptic drive due to differences in input patterns,
   while homeostatic plasticity of V1 neuron excitability reduces the
   postsyn- aptic variability in firing rates. Together these two
   mechanisms, thought to be applicable across sensory systems in
   general, lead to biological maps that develop stably and robustly,
   yet adapt to the visual environment. The modeling results suggest
   that topographic map stability is a natural outcome of low-level
   processes of adaptation and normalization. The resulting model is
   more realistic, simpler, and far more robust, and is thus a good
   starting point for future studies of cortical map development.


   @article{Stevens02102013,
      author = {Stevens, Jean-Luc R. and Law, Judith S. and Antolik,
      Jan and Bednar, James A.},
      title = {Mechanisms for Stable, Robust, and Adaptive Development
      of Orientation Maps in the Primary Visual Cortex},
      journal = {The Journal of Neuroscience}
      volume = {33},
      number = {40},
      pages = {15747-15766},
      year = {2013},
      doi = {10.1523/JNEUROSCI.1037-13.2013},
      url = {http://www.jneurosci.org/content/33/40/15747.full}
   }

=================
Running the model
=================

This model may be run with the Topographica simulator as follow:

./topographica -g gcal.ty

A suitable version of the Topographica simulator may be obtained using
git:

git clone https://github.com/ioam/topographica.git
cd topographica
git submodule update --init
git checkout GCALModelDB

After 10000 iterations (~ 15 minutes on a quad-core 3Ghz machine), the
orientation map corresponding to the condition shown in Figure 9F can
be measured (GCAL model at 100% contrast). In the GUI, this may be
viewed by opening the Orientation Preference window (Plots ->
Preference Maps -> Orientation Preference).

The default settings are appropriate for quick simulations, suitable
for regular exploratory work. An exact match to Figure 9E requires a
slower simulation:

./topographica -p cortex_density=98 -p area=1.5 -g gcal.ty

The linear density of cortical units is doubled for a cortical area of
1.5 x 1.5. Making the cortical area 2.25 times larger allows border
effects to be eliminated from the central 1.0 x 1.0 area. With four
times more cortical units due to the higher density, the overall
network is nine times larger.


=================================
Reproducing all published figures
=================================

An IPython Notebook detailing how to reproduce the entire paper
(including all figures and simulations) may be viewed or run from
here:

http://topographica.org/_static/stevens_jn13.html
"""

from math import pi
import os

import numpy
import param

from topo import learningfn,numbergen,transferfn,pattern,projection,responsefn,sheet

import topo.learningfn.optimized
import topo.learningfn.projfn
import topo.transferfn.optimized
import topo.pattern.random
import topo.pattern.image
import topo.responsefn.optimized
import topo.sheet.lissom
import topo.sheet.optimized

import topo.transferfn.misc
from topo.base.arrayutil import DivideWithConstant
from topo.transferfn.misc import HomeostaticResponse
from topo.misc.commandline import global_params as p

p.add(

    input_seed = param.Integer(default=102, doc="""
         The random seed to set the input patterns"""),

    dataset=param.ObjectSelector(default='Gaussian',objects=
        ['Gaussian','Nature', 'VGR'],doc="""
        Set of input patterns to use::
          :'Gaussian': Two-dimensional Gaussians
          :'Nature':   Shouval's 1999 monochrome 256x256 images.
          :'VGR':       Simulated vertical goggle rearing
                       (anisotropically blurred Shouval)"""),

    num_inputs=param.Integer(default=2,bounds=(1,None),doc="""
        How many input patterns to present per unit area at each
        iteration, when using discrete patterns (e.g. Gaussians)."""),

    aff_strength=param.Number(default=1.5,bounds=(0.0,None),doc="""
        Overall strength of the afferent projection to V1."""),

    exc_strength=param.Number(default=1.7,bounds=(0.0,None),doc="""
        Overall strength of the lateral excitatory projection to V1."""),

    inh_strength=param.Number(default=1.4,bounds=(0.0,None),doc="""
        Overall strength of the lateral inhibitory projection to V1."""),

    aff_lr=param.Number(default=0.1,bounds=(0.0,None),doc="""
        Learning rate for the afferent projection to V1."""),

    exc_lr=param.Number(default=0.0,bounds=(0.0,None),doc="""
        Learning rate for the lateral excitatory projection to V1."""),

    inh_lr=param.Number(default=0.3,bounds=(0.0,None),doc="""
        Learning rate for the lateral inhibitory projection to V1."""),

    area=param.Number(default=1.0,bounds=(0,None),
        inclusive_bounds=(False,True),doc="""
        Linear size of cortical area to simulate.
        2.0 gives a 2.0x2.0 Sheet area in V1."""),

    retina_density=param.Number(default=24.0,bounds=(0,None),
        inclusive_bounds=(False,True),doc="""
        The nominal_density to use for the retina."""),

    lgn_density=param.Number(default=24.0,bounds=(0,None),
        inclusive_bounds=(False,True),doc="""
        The nominal_density to use for the LGN."""),

    cortex_density=param.Number(default=49.0,bounds=(0,None),
        inclusive_bounds=(False,True),doc="""
        The nominal_density to use for V1."""),

    # Noisy disk parameters

    retinal_waves=param.Integer(default=0,bounds=(0,None),doc="""
        How many retinal wave (noisy disk) presentations before the
        chosen dataset is displayed. Zero for no noisy disk
        presentations otherwise 6000 has been typically used.."""),

    percent_noise = param.Number(default=50, doc="""
       The percentage of the noise strength to the disk strength for noisy disks"""),

    contrast=param.Number(default=100, bounds=(0,100),
        inclusive_bounds=(True,True),doc="""
        Brightness of the input patterns as a contrast (percent)."""),

    # Toggle between the L, GCL, AL and GCAL models

    gain_control = param.Boolean(default=True, doc="""
        Whether or not gain-control (divisive lateral inhibitory) is to be
        applied in the LGN"""),

    homeostasis = param.Boolean(default=True, doc="""
        Whether or not the homeostatic adaption should be applied in V1"""),

    t_init = param.Number(default=0.20, doc="""
        The initial V1 threshold value. This value is static in the L and GCL models
        and adaptive in the AL and GCAL models.""")
)
### Input patterns

# Scale of 1.0 is equivalent to 100% contrast.
contrast_scale = p.contrast / 100.0
if p.dataset=="Gaussian":
    total_num_inputs=int(p.num_inputs*p.area**2)
    inputs=[pattern.Gaussian(x=numbergen.UniformRandom(lbound=-(p.area/2.0+0.25),
                                                 ubound= (p.area/2.0+0.25),seed=p.input_seed+12+i),
                       y=numbergen.UniformRandom(lbound=-(p.area/2.0+0.25),
                                                 ubound= (p.area/2.0+0.25),seed=p.input_seed+35+i),
                       orientation=numbergen.UniformRandom(lbound=-pi,ubound=pi,seed=p.input_seed+21+i),
                       size=0.088388, aspect_ratio=4.66667, scale=contrast_scale)
            for i in xrange(total_num_inputs)]
    combined_inputs = pattern.SeparatedComposite(min_separation=0,generators=inputs)
if p.dataset in ["Nature", "VGR"]:
    # Do not randomly rotate GR patches (otherwise horizontal bias is meaningless)
    patch_orientation = 0.0 if p.dataset=="VGR" else numbergen.UniformRandom(lbound=-pi,ubound=pi,seed=p.input_seed+65)

    if p.dataset=="Nature":
        image_filenames= [param.resolve_path("images/shouval/combined%02d.png"%(i+1)) for i in xrange(25)]
    else:
        # Actual ordering used in the paper for historical reasons,
        # makes little difference if ordering=range(1,26) is used instead.
        ordering = [1,25,8,24,22,23,4,15,17,13,5,20,18,14,3,6,12,10,21,7,9,2,16,11,19]
        image_filenames  = [param.resolve_path("images/VGR/combined%02d_VGR.png" % i) for i in ordering]

    inputs=[pattern.image.FileImage(filename=f,
                scale=contrast_scale,
                size=10.0,
                x=numbergen.UniformRandom(lbound=-0.75,ubound=0.75,seed=p.input_seed+12),
                y=numbergen.UniformRandom(lbound=-0.75,ubound=0.75,seed=p.input_seed+36),
                orientation=patch_orientation)  for f in image_filenames]

    combined_inputs = pattern.Selector(generators=inputs)

noise_ratio = (p.percent_noise / 100.0)
disk_scale= 1.0 / (1.0 + noise_ratio)
rand_scale= noise_ratio / (1.0 + noise_ratio)

disks_inputs=[topo.pattern.Composite(operator=numpy.add,
                   scale=contrast_scale,
                   generators=[topo.pattern.Disk(
                        x=numbergen.UniformRandom(lbound=-2.125,ubound=2.125, seed=p.input_seed+12),
                        y=numbergen.UniformRandom(lbound=-2.125,ubound=2.125, seed=p.input_seed+36),
                                                 size=2.0, aspect_ratio=1.0, scale = disk_scale,
                                                 bounds=sheet.BoundingBox(radius=1.125),smoothing=0.1),

                             topo.pattern.random.UniformRandom(scale=rand_scale)])]

retina_inputs = topo.pattern.Selector(generators=disks_inputs)

if p.retinal_waves == 0:
    retina_inputs = combined_inputs
else:
    topo.sim.schedule_command(p.retinal_waves, 'topo.sim["Retina"].set_input_generator(combined_inputs, push_existing=False)')
### Specify weight initialization, response function, and learning function
projection.CFProjection.cf_shape=pattern.Disk(smoothing=0.0)
projection.CFProjection.response_fn=responsefn.optimized.CFPRF_DotProduct_opt()
projection.CFProjection.learning_fn=learningfn.optimized.CFPLF_Hebbian_opt()
projection.CFProjection.weights_output_fns=[transferfn.optimized.CFPOF_DivisiveNormalizeL1_opt()]
projection.SharedWeightCFProjection.response_fn=responsefn.optimized.CFPRF_DotProduct_opt()

topo.sim['Retina']=sheet.GeneratorSheet(nominal_density=p.retina_density,
    input_generator=retina_inputs, period=1.0, phase=0.05,
    nominal_bounds=sheet.BoundingBox(radius=p.area/2.0+0.25+0.375+0.5))


lgn_surroundg = pattern.Gaussian(size=0.25,aspect_ratio=1.0, output_fns=[transferfn.DivisiveNormalizeL1()])

# LGN has lateral connections for divisive normalization for GCL and GCAL models
for s in ['LGNOn','LGNOff']:

    extra_kwargs = dict(tsettle=2,strict_tsettle=1) if p.gain_control else dict(tsettle=0,strict_tsettle=0)

    topo.sim[s]=sheet.optimized.LISSOM_Opt(
                    nominal_density=p.lgn_density,
                    nominal_bounds=sheet.BoundingBox(radius=p.area/2.0+0.25+0.5),
                    output_fns=[transferfn.misc.HalfRectify()],
                    measure_maps=False, **extra_kwargs)

    if p.gain_control:

        topo.sim.connect(s,s,delay=0.05, name='LateralGC',
                    dest_port=("Activity"),activity_group=(0.6, DivideWithConstant(c=0.11)),
                    connection_type=projection.SharedWeightCFProjection,
                    strength=0.6,weights_generator=lgn_surroundg,
                    nominal_bounds_template=sheet.BoundingBox(radius=0.25))

learning_rate = 0.01 if p.homeostasis else 0.0

topo.sim["V1"] = sheet.lissom.LISSOM(nominal_density=p.cortex_density,
    tsettle=16, plastic=True,
    nominal_bounds=sheet.BoundingBox(radius=p.area/2.0),
    output_fns = [HomeostaticResponse(t_init=p.t_init, learning_rate=learning_rate)])

topo.sim["V1"].joint_norm_fn=topo.sheet.optimized.compute_joint_norm_totals_opt

centerg   = pattern.Gaussian(size=0.07385,aspect_ratio=1.0,
                             output_fns=[transferfn.DivisiveNormalizeL1()])

surroundg = pattern.Gaussian(size=0.29540,aspect_ratio=1.0,
                             output_fns=[transferfn.DivisiveNormalizeL1()])

on_weights = pattern.Composite(generators=[centerg,surroundg],operator=numpy.subtract)

off_weights = pattern.Composite(generators=[surroundg,centerg],operator=numpy.subtract)


strength_factor = 6.0
topo.sim.connect(
    'Retina','LGNOn',delay=0.05,strength=2.33*strength_factor, name='Afferent',
    connection_type=projection.SharedWeightCFProjection,
    nominal_bounds_template=sheet.BoundingBox(radius=0.375),
    weights_generator=on_weights)

topo.sim.connect(
    'Retina','LGNOff',delay=0.05,strength=2.33*strength_factor, name='Afferent',
    connection_type=projection.SharedWeightCFProjection,
    nominal_bounds_template=sheet.BoundingBox(radius=0.375),
    weights_generator=off_weights)

"Center surround (difference-of-Gaussian) weights successfully generated"

# Adjust feedforward delays to allow a common measurement protocol with and without gain control.

LGN_V1_delay = 0.05 if p.gain_control else 0.10

topo.sim.connect(
    'LGNOn','V1',delay=LGN_V1_delay,strength=p.aff_strength,name='LGNOnAfferent',
    dest_port=('Activity','JointNormalize','Afferent'),
    connection_type=projection.CFProjection,learning_rate=p.aff_lr,
    nominal_bounds_template=sheet.BoundingBox(radius=0.27083),
    weights_generator= pattern.random.GaussianCloud(gaussian_size=2*0.27083),
    learning_fn=learningfn.optimized.CFPLF_Hebbian_opt())

topo.sim.connect(
    'LGNOff','V1',delay=LGN_V1_delay,strength=p.aff_strength,name='LGNOffAfferent',
    dest_port=('Activity','JointNormalize','Afferent'),
    connection_type=projection.CFProjection,learning_rate=p.aff_lr,
    nominal_bounds_template=sheet.BoundingBox(radius=0.27083),
    weights_generator= pattern.random.GaussianCloud(gaussian_size=2*0.27083),
    learning_fn=learningfn.optimized.CFPLF_Hebbian_opt())

"Afferent GaussianCloud weights successfully generated."

lateral_excitatory_weights = pattern.Gaussian(aspect_ratio=1.0, size=0.05)

topo.sim.connect(
    'V1','V1',delay=0.05,strength=p.exc_strength,name='LateralExcitatory',
    connection_type=projection.CFProjection,learning_rate=p.exc_lr,
    nominal_bounds_template=sheet.BoundingBox(radius=0.104),
    weights_generator=lateral_excitatory_weights)

"Lateral excitatory weights successfully generated"

lateral_inhibitory_weights = pattern.random.GaussianCloud(gaussian_size=0.15)

topo.sim.connect(
    'V1','V1',delay=0.05,strength=-1.0*p.inh_strength,name='LateralInhibitory',
    connection_type=projection.CFProjection,learning_rate=p.inh_lr,
    nominal_bounds_template=sheet.BoundingBox(radius=0.22917),
    weights_generator=lateral_inhibitory_weights)

"Lateral inhibitory weights successfully generated"

### Default locations for model editor
topo.sim.grid_layout([[None,    'V1',     None],
                      ['LGNOn', None,     'LGNOff'],
                      [None,    'Retina', None]], xstart=150,item_scale=0.8)


import topo.analysis.featureresponses
topo.analysis.featureresponses.FeatureMaps.selectivity_multiplier=2.0
contrasts =  [{'contrast':c} for c in [100, 80, 60, 40, 20, 10]]
topo.analysis.featureresponses.FeatureCurveCommand.curve_parameters= contrasts

Loading data, please wait...