WM     Working Models

The model is described below in roughly the order in which operations occur.

Visual stimulus
Noise in the input
Cone photoreceptors
Converting RGB stimuli to Cone excitation
Horizontal grid
BP (bipolar cell)
Long term adaptation
Gain control
Ganglion cell
Responses and output
Automatic stimulus centering
Stimulus override
Models of this type


Latest development

21-Sep-2008  ~wyeth/w/model/wm/spop/ret/devel/
   Connected mod_mesh to mod_pop
17-Sep-2008  ~wyeth/w/model/wm/mesh/m_ad/
   Added long-term adaptation, tested with stat1_anti 

Visual stimulus

The visual stimulus S(x,y,t) is defined on a cartesian grid, regardless of the arrangement of the photoreceptors. S(x,y,t) takes a value from 0 to 1, and is xn by yn by tn pixels (values are set in the .moo file). Here, tn can be any sufficiently large positive integer, and does not have to be a power of two.

Optical blurring of the stimulus is controlled in the .stm file. Note, this might currently require the space-time dimensions being set to powers of two if the FFT is used for the blurring.

Color, if present, is controlled in the .stm file with the parameters amp_r, amp_g, amp_b, which are intended to vary from 0 to 1. Please note that using -1 will not simply invert the phase of a grating, because the stimulus is non-negative.


In the WM modeling framework, the entire visual stimulus for one trial is generated before the model processes that trial; therefore, memory requirements will limit the maximum single-trial simulation length. As of Sept. 2008, the reasonable time limit for a 32x32 spatial-pixel (xn by yn) model is 131,072 time units (tn). In total, this is 2^27 space-time pixels, or 1/8 Gigapixels. Thus, a simple 32x32 model will run for at most about a couple minutes of wall clock time. To overcome this limit, this retinal model is able to hold its state so that it can pick up processing on a new trial where it left off on the previous trial. In this case, the following apply,
  • To cause the model to hold its state, in the .moo file, use:
    mod_trial_reset  0  # Do not reset state between trials.
  • The simulation must be run with WM (do not use MM).
  • MM is problematic because there is essentially no meaning to holding on to a state from one trial to the next if the trials are executed in parallel.
  • Stimuli should have a period of constant luminance around the boundary, because for simplicity the model currently assumes that the value of the stimulus at time zero on any trial is also the value that applies back in time as far as is required for the model to initialize any processing for the trial. For example, this is assumed when the temporal filtering of the cones, described below, is applied to the visual stimulus.
  • Future. A more robust implementation could allow for a period of overlap by holding in memory the stimuli two consecutive trials at any given time.

Noise in the input

Noise is added to the visual stimulus S(x,y,t) before it reaches the cones. The added noise signal, N0(x,y,t), at each pixel is Gaussian distributed with mean zero and SD proportional to the square root of S(x,y,t). Relevant parameters are:
  • mesh_qno_f - [1.0] noise amplitude. Set to zero to suppress stimulus photon noise.

Cone photoreceptor (C) signal

Non-cartesian cone mosaic

Each cone sees only the stimulus values at the four points that surround its position with the cartesian stimulus grid. The weight of the four locations is proportional to the horiztonal and vertical position of the cone within the square grid box. For example, if the cone sits exactly on a cartesian grid point, it will see only the input from that point. If it sits in the middle of a grid box, it will weight the four corners equally.

For simple color stimuli, the input to the cone from a location in the stimulus is determined by applying the LMS matrix (in the .moo file) to the stimulus, which has already been weighted by the color amplitudes specified in the .stm file.

Let C_raw be the signal that is derived by the spatial and color weighting scheme as described.

Noise. At this point (or probably before the color signals are added), it would be appropriate to transform C_raw into a noisy signal that accounts for the random arrival of photons. We are skipping this step for now.

Temporal filtering. C_raw is convolved with a temporal filter that is defined as the difference of two Maxwell distributions: M1 - a*M2. Three parameters control the temporal filter:

  • mod_mesh_photo_m1 - spread of the M1 distribution
  • mod_mesh_photo_m2 - spread of the M2 distribution
  • mod_mesh_photo_a2 - scale factor for M2 ('a' above)

  • photo_delta - [0] set to 1 to override and use delta-function
The resulting filter is then scaled so that its integral is 1. Thus, the response to a constant stimulus of value A is A.

To estimate appropriate temporal filtering, see the responses of horizontal cells reported in Smith, Pokorny, Lee, and Dacey (2001). The function given in Schneeweiss and Schnapf, 1999, seems too slow.

Customizing cone colors. In the <mosaic> object, a filename (or the reserved word "null") can be specified to read custom cone IDs from a text file:

  custom_cid_file  cone_file_01.txt 
where the text file is a list of pairs (ConeIndex, ConeID):
0 S
10 M
103 L

Converting RGB stimuli to Cone excitation levels

C_ex is the resulting cone excitation signal that will be applied to the horizontal-cell resistive mesh. This signal is computed from a matrix of values:
    lr lg lb 
    mr mg mb 
    sr sg sb 
such that the excition of any L cone, for example, to a color with components (r,g,b) would be:
    EL = lr*r + lg*g + lb*b          (Eqn 1)

The matrix values are computed as the Frobenious normalized form of the values specified in <cone_mosaic>:
    L_rgb  0.0085  0.0271  0.0059   # L cone sensitivity to R G B
    M_rgb  0.0035  0.0281  0.0086   # M cone sensitivity to R G B
    S_rgb  0.0009  0.0024  0.0237   # S cone sensitivity to R G B

Isoluminant stimulus

To make a stimulus that is isoluminant with respect to the M+L signal (the summed response of the M and L cones), we can define two colors that are to have the same luminance:
    C1 = (r1, g1, bo)
    C2 = (r2, g2, bo)
where the blue component is constant. We can equate the summed responses of the L and M cones for the two colors,
    lr*r1 + lg*g1 + mr*r1 + mg*g1  =  lr*r2 + lg*g2 + mr*r2 + mg*g2.          (Eqn 2)
If we now assume that the red components are modulated by an amplitude, Ar, evenly but oppositely around a mean value, a, (for example, a = 0.5), and likewise for the green components, then
    r1 = a + Ar
    r2 = a - Ar
    g1 = a - Ag
    g2 = a + Ag
where the modulation favors red in C1 and green in C2, then the Eqn 2 above becomes:
         (lr + mr)
    Ag = ---------  Ar.
         (lg + mg)
Notice that the mean value, a, has canceled out of the equation. Using the values specified in the color matrix above,
    Ag = 0.2174 Ar.
If we assume,
    a = 0.5,
then the largest value that Ar can take is 0.5 under the assumption that RGB color values are limited to the range [0-1], in which case,
    Ag = 0.1087.
For this example, the final RGB colors can be written to two decimal places as,
    C1 = (1.00, 0.39, bo)
    C2 = (0.00, 0.61, bo).
The smaller modulation of green reflects the fact that, for our CRT, the cones are more sensitive to the green gun than to the red gun.

Horizontal (H) mesh signal

The smoothing operation of the horizontal cells is implemented by solving the diffusion equation on a cartesian resistive grid (or mesh). This grid has the same spatial dimensions as the stimulus, S(x,y,t), which is xn by yn. For each grid point:
  • A weighted average of nearby cone exitations is applied through a conductance, Gp, mesh_gp.
  • Connections to the four neighboring grid points have conductance, Gh, mesh_gh.
  • Capaticance, C, is mesh_c.
The ratio of Gp to Gh controls the space constant. For example, for a 1D chain, L = sqrt(Gh/Gp), where L is the characteristic length, i.e., the distance over which an input applied (and held constant) at one point decays to 1/e of its value.

The differential equations are solved with a time step, DT, given by mesh_dt, which is in units specified by tscale (in the .moo file). Increasing DT should make the model run faster, but with less accuracy. One should try varying DT to see how the output changes, to be sure to pick a reasonable value. You may be able to speed the model up by 10 times or more and sacrifice almost nothing.

BP (bipolar) signal

The BP signal, or photoreceptor minus horizontal signal, PH(t) is the difference: C_ex - c*H, where c is a constant set by mesh_whc, which controls the strength of the surround relative to the center. Each BP signal is centered on a cone, thus its location is taken to be that of the cone. The H signal input to the BP cell is computed as a weighted average of the four grid points that surround the cone/BP location.

Bipolar cells are supposedly highly nonlinear. This is not implemented here. They are very sensitive to small incr. and decr. in contrast (Burkhardt & Fahey, 1998).

Long-term Adaptation

A long-term adaptation signal ad1(t) is computed by low-pass filtering (ideal LP, i.e., decaying exponential) the PH(t) signal. This offset is subtracted from PH(t) to form PHAD(t), the adapted difference signal.

To control long-term adaptation, set these parameters in the <retina0> object:

  • ad1_tau - time constant of adaptation (s) [0.0]
  • ad1_amp - amplitude of adaptation, 1=fully adapt [0.0]
Set both of these values to 0 to turn off this form of adaptation.

Gain control

Two time-varying control signals are computed from global (across the entire stimulus grid) spatial averages:
  • Gain1(t) - the average of the absolute value of PHAD(t)
  • Gain2(t) - the average of the square of PHAD(t), LP-filtered

Gain-controlled HP-filter, GCHP

This form of gain control applies a high-pass filter to the PHAD(t) signal by subtracting a scaled, low-pass signal (LP(t), where the low-pass is performed by causal convolution with an alpha-function). The scaling is controlled by Gain1S(t), which is a sigmoid-squashed version of Gain1 that can only take values from 0 to 1. The resulting signal is known as PHADG(t). Relevant parameters are:
  • mesh_gain_hp_flag - [0] 1-apply a time constant of adaptation (s)
  • mesh_lp_tau - [0.0] 1/alpha for low-pass filtering with an alpha-function (s). Must be > 0.0 to use this gain control.
  • mesh_gsig_mean - the gain1 value that will map to 0.5 * lp
  • mesh_gsig_slope - the slope of the sigmoid at 'mesh_gsig_mean'

Inhibitory gain-control, GCI

This form of gain control applies an inhibitory conductance to the integrate-and-fire spike generator. The conductance is set to the Gain2(t) signal. Parameters that control the application of this form of gain control:

  • mesh_gain_inhib_g2 - (0) 1-apply the gain control. If 0, the inhibitory conductance is set to zero.
  • mesh_lp2_tau - (0.0) 1/alpha for alpha-function LP-filter. Must be set if mesh_gain_inhib_g2 is 1.

Ganglion cell (GC) signal

The ganglion cell takes a single PHADG(t) signal or a pool of PHADG(t) signals as input. The input pool is controlled by the following parameters:

  • sum_bp_distrib - [0 default] Distribution of synaptic weights. 0-Gaussian, 1-uniform.

  • sum_bp_dist - [required, no default] Distance of spread of spatial summation (stimulus pixels). For a Gaussian, this specifies the SD; for a uniform distribution, the radius. Set to zero to use a single BP signal.

  • sum_bp_minw - [0.02 default] For a Gaussian distribution, set any weight to zero that is below this fraction of the peak.

  • sum_bp_ccode - [111 default] Cone code. 111-SML, 001-L only, 110-S and M, etc. This controls the types of cones that will be sampled.
The pooled signal can be taken as the output, or used to drive the excitatory conductance of an integrate-and-fire spiking model. The GCI gain control, if present, is implemented by setting the inhibitory conductance to be Gain2(t). This signal can be regulated using the standard parameters for the IFC model, e.g., gi_scale, gi_bias, and gi_noise.

Responses and output

Specific examples of output for the retina0 mesh model are given here. For a more general description of output and response files,
click here.

Standard output

To save output to an ndata (.nd) file, use lines like the following in the .rsp file:
    save_pop_unit_as_unit0 s 1000.0 rgc 0 0 1 gc_spikes
where the x-index (the first '0') specifies the cone index for ganglion cell to be saved, and the z-index ('1') specifies the ON sublayer ('0' being the OFF sublayer). The following line will save a set of responses around the center of the mosaic,

    save_pop_grid_as_on   s 1000.0 rgc 0 0 1 100 1 1 gc_spikes
where the x-index (the first '0') specifies the starting cone index and the x-length ('100' above) specifies the number of contiguous cones to save. The z-index (the final '1') specifies the ON sublayer. In the examples above, the following response components can be requested:

Dumping output

To dump intermediate traces for testing and development, the following options can be used:

Automatic stimulus centering

As part of the model preparation (in the prep routine), the stimulus can be centered on a particular mosaic index automatically by using the following command in the <retina0> object:
  stim_auto_center_mi  0
where the number, 0, specifies any mosaic index, or -1 to deactivate this featuer. By default, this feature is turned off.

The centering works by changing the values for cy and cx, which must appear as non-varying parameters in the .stm file. These changes will be reflected in the output .nd file.

Stimulus Override

The visual stimulus can be replaced with an artificial stimulus that provides exact values to the cones in the mosaic. This allows testing and computing, for example, cone purity measures of the surround.

This feature is controlled by lines such as the following in the <retina0> object:

  stim_override 0           # 0-Default, do not override
                            # 1-The stimulus will be ignored and replaced

  stim_override_binary  L    # 0,L,M,S,all,none - set to 1, all others to 0
  stim_override_delay   1    # 0-full override, 1-override at tn/2, zero before
  stim_override_zero_ci 0    # Set this ci to zero always, -1 to deactivate

The feature is activated by setting stim_override to 1. When activated, all cone excitation levels are set to zero, except those specified by stim_override_binary, which will take the value 1.0. The stim_override_binary parameter may be one of the following:
  • L - all L cones will have excitation level 1.0
  • M - all M cones will have excitation level 1.0
  • S - all S cones will have excitation level 1.0
  • none - no cones will have excitation level 1.0
  • [ci] - an integer cone index will have level 1.0
  • all - all cones will have excitation level 1.0
The stim_override_delay parameter can be set to 1 to cause all values to be zero before time tn/2, and the override conditions will hold thereafter.

The stim_override_zero_ci parameter allows one particular cone to be held at 0, regardless of all other conditions that might have otherwise caused it to be set to 1. This can be used to prevent a center cone from influencing the measurement of surround purity.

Models of this type