We account for the increasing size and complexity of the receptive fields as one moves up the cortical hierarchy by assuming that in area ℓ the units receive input from two neighboring units in area ℓ − 1. For example, units at position 1 in area 1 receive inputs from the units at position 1 and 2 of the input, so layer 0, while units at position 2 receive feedforward inputs from area 0 units at position 3 and 4, etc. At each position in layer 1 two units with different receptive field types. In the first the input is proportional to the sum of the outputs of the two units in layer 0 which project to it, in the other the input is proportional to their difference. In layer 2 units at position i receive inputs from units at position 2i −1 and 2i in layer 1, and there are 4 kinds of receptive fields. The input into units with type 1 receptive fields take as argument the sum of output of the two units with type 1 receptive fields in layer 1. For type 2 units the argument is the difference between these. Type 3 receptive fields have as input the sum of the outputs of the type 2 units in layer 1, while type 4 units have the difference of these two as input. This algorithm is repeated for higher layers. This results in a system in which for layer ℓ there 2^{L − ℓ} positions, at each of which there are 2^ℓ different types of receptive fields (Figure 2).

We denote the rate of the unit of layer ℓ at position i with receptive field type k with rk,iℓ and its feedforward input by Ik,iℓ(FF) for 1 ≤ ℓ ≤ L − 1 the feedforward input is given by

for i = 1, …, 2^{L − ℓ} and k = 1, …, 2^ℓ. Here W_FF is the strength of the feedforward connections.

For ℓ = 1 two types of receptive field exists, k = 1, 2 and the feedforward input is, for i = 1, …, 2^L − 1, given by I1,i1(FF)=WFF(r2i-10+r2i0) and I2,i1(FF)=WFF(r2i-10-r2i0), where ri0 is the output of the Ith LGN unit.

Units in area ℓ receive reciprocal input from those units in area ℓ + 1 onto which they project. This feedback input has the same sign as the feedforward input but is modulated by connection strength W_FB. The cortical feedback input, Ik,iℓ(FB), in the unit with receptive field type k at position i in layer ℓ, is given by:

Units in the L_th cortical area have no feedback inputs. To compensate for this we assume that the strength of the feedforward connection to the last area is W_FF + W_FB rather than W_FF. Note also that in the layer L the receptive fields span the whole input range and there is only 1 position (i = 1). The feedforward input into layer L units is given by

It is well known that the synaptic connections from one cortical area to another emanate from pyramidal neurons (Rodney et al., 2004). So, both feedforward and the feedback pathways are are excitatory. However, in the model here we are considering effective inputs. The effective connection are positive if the excitatory population in the presynaptic “On” unit project to the excitatory “On” population and the inhibitory “Off” population of the postsynaptic unit while the presynaptic excitatory “Off” population projects to the excitatory “Off” and inhibitory “On” populations in the postsynaptic unit. The effective connection is negative if the presynaptic “On” cells project to the inhibitory “On” and excitatory “Off” cells in the postsynaptic unit and similar for the excitatory presynaptic “Off” cells.

The Pul is the largest thalamic nucleus in primates and it presents anatomical and physiological properties that involve with visual cortical transmission. The Pul has two topographic maps that traverse retinotopically the lateral (PL) and inferior (PI) subdivisions of the Pul. These two subdivisions connect directly with the ventral stream of the cortex. The other two subdivision of the Pul, medial (PM) and anterior divisions (PA), connect partially to cortical areas of the dorsal stream (Stepniewska, 2003; Kaas and Lyon, 2007). RFs of pulvinar neurons has simple visual features that correspond with the cortical areas that they target. RFs of cat and monkey pulvinar neurons have small and large diameter sizes, are driven by orientated bars with a broad tuning responses as well as motion to textured patterns, and color-sensitive attributes (Casanova, 2003).

The major source of visual inputs to Pul come from the visual cortex. Lesions in striate cortex of macaque eliminated the visual response of pulvinar neuron. This input from the cortex is represented as a gradient inside the Pul. Shipp (2003), based on cortico-thalamic and thalamo-cortical connections, postulates the existence of a “cortical gradient” in the Pul. While injections with dual tracer in V1 and V4 label preferentially respective medio-caudal and latero-rostral pulvinar areas, injection in V2 target lateral within Pul, and inferior temporal cortical areas, medial within Pul. On the other hand, injections in area V1, that represent retinotopic position of either the upper and lower contralateral hemivisual field, label neurons in respective hemifield of both PL and PI. Thus, a fronto-occipital axis in the cortex is reproduced as a medio-lateral gradient in the pulvinar (medio-lateral cortical axis rotates to a rostro-caudal gradient in the thalamus). In addition to the gradient observed in the Pul, the cortico-pulvino-cortical loop has two types of connections. In one subgroup the cortical layer VI project to a pulvinar region, which in turn send back to same cortical area but to layer IV (“Reciprocal connections,” similar to connections between V1 and the LGN). In the other, connections arise from cortical layer V and end in a non-reciprocal pulvinar region. In turn, this pulvinar region sends back orthogonally to the cortex. These latter are known as “non-reciprocal connections” and it is considered here as an open loop (Sherman, 2007).

In addition to projections from and to the cortex, the Pul also has a local circuitry. Recently works show at least four intrinsic interactions: (i) axons type I are branched and highly divergent (1.0–3.0 mm), to the extent that they can easily be shown to cross over subdivisions (Rockland, 1996, 1998); (ii) Long range inhibitory interneurons traverse areas in 1.0 mm of length (Imura and Rockland, 2006); (iii) The existence of “bridges” between PI subdivisions that stain to calcium binding protein calbindin and to substance P (Stepniewska, 2003); (iv) Inhibitory inputs from the reticular nucleus which receives excitatory branches from the cortico-thalamic and thalamo-cortical axons (Sherman and Guillery, 2000).

For the Pul architecture the previous attributes are considered. The Pul is modeled similarly to the cortex. However, each pulvinar area has at most 4 types of RFs, the patterns corresponding to k = 1, 2, 2^{ℓ − 1}, 2^ℓ, for l ≥ 2. Reciprocal cortico-pulvino-cortical interactions mainly have the effect of changing the effective if modifying the effective cortical transfer functions, so that, in the interest of simplicity only the non-reciprocal cortico-pulvino-cortical pathway is assumed and the gradient inside the Pul is modeled as a feedforward pathway with long range connections. Pulvinar units in area ℓ receive input from cortical units in area ℓ and from pulvinar units in areas 1 to ℓ – 1. The input Jk,iℓ(PC) from cortex to unit i of type k in pulvinar region ℓ is given by

for l = 3, …, L − 1, i = 1, …, 2^{L − ℓ − 1} and k = 1, 2^{ℓ − 1}.

The long range interactions in the pulvinar are mediated through large GABAergic interneurons (Imura and Rockland, 2006). Thus the connections between units in different pulvinar layers are through inhibitory synapses. Nevertheless, as for cortico-cortical interactions, the effective coupling can be positive or negative, depending on whether the postsynaptic target neurons are excitatory or inhibitory. The input Jk,iℓ(PP) from the rest of the pulvinar satisfies:

Here we have used sk,iℓ denote the rate of the pulvinar units. The units in pulvinar layer 1 do not receive input from the rest of the pulvinar, Jk,i1(PP)=0, while for pulvinar layer 2 we assume that J(PP) is given by

This specifies how the pulvinar input to pulvinar units depends on the activity in the previous areas. For example in pulvinar layer 4, the input Jk,i4(PP) for k = 1 is given by the combination of the pulvinar-feedforward and the long range connections

A direct expression of Jk,iℓ in the activity of units in the previous layers for different values of k and ℓ is straightforward.

Finally, the input from pulvinar to cortex, I(CP), is given by:

where sk,iℓ is the output of pulvinar unit for k = 1, 2, 2^{ℓ − 1}, 2^ℓ, and sk,iℓ=0 otherwise.

The activity of the last cortical area is the summary of the previous ones. When a discontinuous or very small change of its response occurs for an increase in contrast of the input, the output is not very useful to estimate the stimulus. The optimal output of layer L is one that spans the whole range of outputs and varies more or less linearly with the input, r₀. However, these two requirements are in conflict. For example, when a homogeneous input from −1, 1 is applied to the feedforward model, a large W_FF assures utilization of the whole dynamic range between −1 and 1, but yields an extremely non-linear curve. Instead, for small W_FF the curve that plots r_L against r₀ appears much more linear, but only covers a small part of the output range. An optimization criterion that penalizes both these extreme cases and measures how good the network is able to transmit information about stimulus contrast is the entropy of the output distribution if the r₀ is distributed homogeneously between −1 and 1. The entropy is low both in the case where the input-output relation is close to a step-function and also when the output range is small. Thus, to optimize the output of the unit r1,1L of different networks when a homogeneous input is applied in layer 1, we determine parameters for which the entropy, H, given by

is maximal. Here, P_L(r) is the probability density distribution of the Lth layer.

We now derive an expression for H, where the relation r_L = F(r₀) is known and r₀ is drawn from a homogeneous distribution between −1 and 1. For a small Δr this probability will be

Since r_L = F(r₀), this is equal to

Here F ′⁻¹ is the derivative of F^–1, the inverse of F.

Since r₀ is drawn from a uniform distribution between −1 and 1, Prob (F⁻¹(r) < r₀ < F⁻¹ (r) + F ′⁻¹(r) Δr) is equal to F ′⁻¹(r)Δr/2. Together with F ′⁻¹(r_L) = 1/F ′(r₀) this yields

Inserting this into equation (14), the entropy is given by

where we have used dr_L = F ′(r₀)dr₀.

That was the optimization procedure for networks in which a homogeneous input is used, ri0=r0. We use similar optimization principle for the analysis of the visual input. In this case, a random input is applied in ri0=σxi, where σ codes for contrast and x_i is independently drawn from a Gaussian with zero mean and unit variance. Here, changes in contrast are changes in σ without changing x_i. To ensure that the response is sensitive to contrast and the tuning is maximally contrast invariant, we want that the amplitude of the 2^L dimensional vector V→ given by Vk=rk,1L varies smoothly with σ², while its direction changes as little as possible as σ is varied.

We define the output amplitude, F, as the average length of the output for LGN inputs with contrast σ, F(σ)=⟨|V→|⟩, where the average is over the variables x_i. Similar to the analysis of the homogeneous input condition, we aim for an amplitude function F that is as linear as possible and uses the dynamic range maximally. This we ensure by imposing a cost function, H_L, defined as HL=∫01dσlog(F′(σ)), for this property. If H_L = L log2, the output scales linearly and exploits the whole dynamic ranges. It decreases if less of the dynamic range is used or the response scales non-linearly.

To explore whether the network can maintain contrast invariant tuning, we calculate the mean of separation distance S between normalizes output vectors e→(σ)=V→(σ)∕|V→(σ)| and e→(σ′)=V→(σ′)∕|V→(σ′)| for LGN inputs ri0=σxk and ri0=σ′xk respectively, S=∫01dσ∫01dσ′e→(σ)⋅e→(σ′). If S = 1, the vectors are in the same direction for all contrasts. As the direction changes more with contrast, S decreases.

For the optimization of the network parameters define a total error E=-2S-eHL that takes both these factors into account. In both optimization criteria, the homogeneous and the visual input, we attempt to minimize the error of the functions. Thus, we want to maximize both S and H_L. So, given the parameter space of models we use the Powell’s method (Press et al., 1992).

In the purely feedforward model (W_FB = W_CP = W_PC = 0), the activity is propagated sequentially through the cortical areas until reach l = L. The input is varied in magnitude to mimic changes in contrast of the stimulus.

If ri0 is the same for all units, r1,j1 will also be identical for all j. As a result, r1,i2 will also be independent of i, etc. At the same time, because ri0 is the same for all i, r2,j1 will be zero for all j. Extending this logic to larger ℓ, we see that rk,iℓ=0 for K = 2, 3, …, 2^{ℓ − 1}, while r1,iℓ=rℓ is the same for all i. The equilibrium rates are given by r₁ = f(W_FFr₀), and

for l ≥ 1.

Figure 3 shows the equilibrium rate as a plotted against r₀ for layers l = 1, 5, and 10, for different values of the feedforward strength and threshold. When W_FF is small the rate in r₁ increases smoothly from approximately −1 to approximately 1 as the input r₀ varies from −2 to 2. For larger ℓ the response is progressively smaller, until for l = 10 the response almost stays at 0. On the other hand, for large W_FF with Ith = 0, r₀ shows a clear sigmoidal response. For larger ℓ the steepness of the sigmoid increases so that for l = 10 the response in almost a step-function. For Ith ≠ 0, the response evolves to a sum of two sigmoids, with thresholds at −Ith/W_FF and Ith/W_FF. While these sigmoids are not as steep as the corresponding sigmoid for Ith, still in layer 10 the response takes values near ±1 or 0, for most of the input range.

To better understand this behavior we analyze the system near r₀ = 0. Since f (0) = 0 we have that if r₀ = 0, r_ℓ = 0 for all ℓ. For 0<r0=δr0≪1, r_ℓ = δr_ℓ is also much less than 1, and, to leading order, satisfies δr_ℓ = 2W_FFf ′(0)δr_{ℓ − 1} for l ≥ 1. Here f′ is the derivative of effective transfer function, f. This means that we can write δr_ℓ as δr_ℓ = Λ^ℓδ₀, where Λ = 2W_FFf ′(0). Thus, when |Λ| < 1 the response gets progressively smaller as ℓ is increased, while to |Λ| > 1 the size of the response increases with ℓ. Note that this decrease/increase is geometric, so that even for Λ relatively close to 1, δr_ℓ will deviate a lot from δr₀ for large ℓ. As a result, for Λ ≠ 1, at r₀ = 0, the slope of the function that plots r_ℓ against r₀ will be either very large or very small when ℓ is large.

Using equation (4) we have that f′(0)=[1+cosh(Ith)]-1. If Ith = γW_FF, Λ is given by Λ = 2W_FF/[1 + cosh(γW_FF)]. For γ = 0, Λ = W_FF and, near r₀ = 0, the slope of the transfer is very low for W_FF < 1 for large ℓ, while for W_FF > 1 it becomes very steep. If γ > 0, Λ increases with W_FF for small W_FF, but it decreases asymptotically to 0 as W_FF is increased further and further. This is because for large W_FF, cosh(γW_FF) increases faster than W_FF. So as W_FF is increased from 0, Λ first increases from 0, until it reaches its maximum, then it decreases again to 0. The maximum value Λ takes depends on γ. If γ is to large, maximum value of Λ is less then 1, so that for large γ the slope of the transfer function near r₀ = 0 is small for any value of W_FF. This is demonstrated in Figure 4, where Λ is plotted against W_FF for different values of γ.

We will no examine the implications of these findings in the limit where the total number of layers L becomes infinite.

If we have a long chain of layers, the rates r_ℓ for large ℓ will approach a constant value. The values it can approach are the stable solutions of the equation r_∞ = f (2W_FFr_∞). With Ith = 0 one solution of this equation is r_∞ = 0, for any value of W_FF. For W_FF < 1 the curve of f will have a slope of less than 1 at r = 0, and r_∞ = 0 is a stable solution. It is also the only solution. For larger W_FF the slope at r = 0 is larger than 1 and there are two extra solutions, one with r_∞ < 0 and one with r_∞ > 0 (see Figure 3A). In this case the solution r_∞ = 0 is unstable, that is, if r₀ deviates slightly from 0, the deviation from this value increases as ℓ is increased. For an infinite chain of hierarchical levels, r_ℓ approaches one of the other two solutions with increasing ℓ. r_ℓ Goes to the smaller value if r₀ < 0, while it goes to the larger value for r₀ > 0. As show in Figure 5A, as W_FF in increased, the two stable non-zero solutions approach −1 and 1 respectively.

When Ith = γ W_FF ≠ 0, the solution r_∞ = 0 also exists for all values of W_FF. However, the stability of this solution depends on γ. For sufficiently small γ there is a transition from a stable to an unstable solution, followed by a second transition from unstable to stable, as W_FF is increased. The first transition is the same as that for the case where Ith = 0. Below this transition r_∞ = 0 is a unique solution which is also stable. Above the transition two new stable solutions appear. These solutions approach −1 and 1 as W_FF is increased. The second transition is also a pitchfork bifurcation of the r_∞ = 0 solution. The r_∞ = 0 switches from unstable to stable and two more solutions appear. These are unstable and asymptotically approach the values ±γ/2. The value of W_FF at which the first bifurcation occurs increases with γ, from W_FF = 1, the bifurcations point of the Ith = 0 solution. The point where the second bifurcation occurs starts at W_FF → ∞ for γ → 0 and decreases with increasing γ. See Figure 5B.

These two critical values of W_FF converge in one bifurcation point as γ increases. We have observed that at the first transition point we move from 1 solution to 3 solutions. At the second, we move from 3 to 5 solutions. As γ is increased the separation between both bifurcation points becomes narrower and narrower. When distance between these two points reaches 0, both points converge turning out only one critical value. At this value we have a transition from 1 to 5 solutions. For still higher values of γ there no longer is a bifurcation from the solution r_∞ = 0, this solution stays stable. Instead at a critical value of W_FF two new pairs of solutions emerge. For r > 0 a stable solution with large r_∞ appears, together with an unstable solution that lies between this solution and the r_∞ = 0 solution. There is also a corresponding pair of solutions with r < 0. These solutions are shown in Figure 5C.

To obtain the regions that 1, 3, or 5 fixed-point solutions we determine the number of solutions in the plane (W_FF, γ). The transition from 1 solution to 3 solutions and from 3 to 5 solutions are the solutions of Λ = 1 with the smaller and larger solution respectively of W_FF at fixed γ. These are obtained by solving

Taking the parametrization x = γ W_FF, this can be separated in a set of two equations:

Using this, we plot 1/W_FF against γ in relationship with the parameter x. Figure 6 shows the solutions of the purely feedforward model when a spatial homogeneous input is applied in layer 1. In the plane, the region for 1 solution is the only stable solution for W_FF small. As W_FF increases and γ is small, 3 solutions appear while two of them are stables. Holding the previous threshold condition and increasing W_FF even further, 5 solutions show up three stables fixed-points and 2 unstable ones. As γ is progressively increased to 1, the range over which there are 3 solutions shrinks. When γ ≈ 0.89 both boundaries are merged. For γ > 0.89 there is a transition from 1 solution to 5 solutions. This transition was determined from solving the fixed-point equations directly.

We now explore the effects of feedback connections when a homogeneous input is applied. As before rk,iℓ=0 for k ≥ 2 and we can write r1,iℓ=rℓ. In this case, the input for equation (1) have W_FB ≠ 0 and W_CP = 0. The new system is described by the equations

and

Effective W_FF values of the new system – The objective here is to analyze if the feedback connections improve the sequential transmission through the cortical hierarchy. As we know from the analytical results of the feedforward model there are two possibles final states. Here, we look for the influence of the factor W_FB on the behavior of this final states. We consider the system equations (22) and (23) in the steady state for a constant input, r₀

and

As before, for r₀ = 0 the solution is r_ℓ = 0. For small input r₀ = δr₀ we can expand the solution, and write

and

Making the Ansatz δr_ℓ = Λ^ℓδr₀, we obtain from equation (26) that Λ has to satisfy

which has the solutions Λ₊ and Λ₋. These solutions are given by

The general solution to the system can be written as

where ψ₊ and ψ₋ are two constants.

From boundary condition, equation (27) we obtain that

while from l = 0 we have ψ₊ + ψ₋ = 1.

After some algebra, one obtains that δr_ℓ is given by

where κ = [1 – f ′(0)(2W_FF + W_FB)/Λ_]/[1 – f ′(0)(2W_FF + W_FB)/Λ₊]. If 8[f ′(0)]²W_FFW_FB < 1, |Λ₋/Λ₊| < 1, otherwise |Λ₋/Λ₊| = 1. So, in the large L limit δrL=Λ-L[1-κ]∕[1-κ(Λ-∕Λ+)L]δr0 goes to zero in the large L limit, if |Λ₋| < 1, while for |Λ₋| > 1 δr_L becomes very large. The transition occurs when Λ₋ = 1, or f ′(0)[2W_FF + W_FB] = 1.

Let us now consider the behavior high up in the hierarchy when r₀ is not approximately 0. As in the network with only feedforward connections, for arbitrary r₀, r_ℓ will approach a set of fixed value as ℓ is increased. Here, these fixed values are the stable solutions of

Thus we can conclude that a network with feedback connections with strength W_FB and feedforward connections with strength W_FF, behaves very similarly to a purely feedforward network in which the strength is W˜FF=WFF+WFB/2. Both networks have the same fixed values to which the rates evolve as ℓ is increased, and if r₀ is small, for both the rate geometrically increases (decreases) when W˜FF<1(>1). Thus adding feedback connections does not qualitatively improve the ability of the network to have an output in the higher areas of the hierarchy that vary smoothly with r₀.

We now consider the effect of adding a Pulvinar like structure (W_CP ≠ 0 and W_PC ≠ 0) on the response of the network to constant input, ri0=r0. To simplify the calculations we assume that there are no feedback connections (W_FB = 0). The model with feedback connections behaves qualitatively similar.

It is straightforward to verify that, as before, rk,iℓ=0 for k ≠ 1 and that r1,iℓ=rℓ. Likewise, for the Pulvinar we have that sk,iℓ=0 for k ≠ 1 and s1,iℓ=sℓ. Taking this into account, the equilibrium rates, r_ℓ and s_ℓ, are, for l = 1, 2, …, L, given by

Here f_ctx and f_pul are the transfer functions of the cortical and Pulvinar units respectively, while s₀ = 0, J₁ = 0, J₂ = 2W_FPs₁, and J_ℓ = (2W_FPs_{ℓ − 1} + W_LPJ_{ℓ − 1})/(1 + W_LP).

Due to the parallel processing and long range connections, a full analysis of this system is much more involved that the analysis of the system without Pulvinar that we have analyzed above. This analysis is beyond the scope of this paper, here we will concentrate on the response of the system to small input, r₀ = δr₀ and indicate how the system behave differently from the network without Pulvinar.

If r₀ = δr₀ is small, the response of all layers of the cortex and Pulvinar will be small, r_ℓ = δr_ℓ and s_ℓ = δs_ℓ, with

where δJ_ℓ = (2W_PFδs_{ℓ − 1} + W_LPδJ_{ℓ − 1})/(1 + W_LP), Fctx=fctx′(0) and Fpul=fpul′(0).

Analogous to what happened in the cortical network with feedback we can here, for l ≥ 2, write δrℓ=(ψ-Λ-ℓ+ψ+Λ+ℓ)δr0. For δs_ℓ we have δsℓ=(ϕ-Λ-ℓ+ϕ+Λ+ℓ)δr0. After some tedious algebra one can show that if the largest of eigenvalues, Λ₊, is larger than 1, δr_L is much larger than δr₀, while if it is smaller, δr_L will be much smaller than δr₀. Thus for a gradual increase of r_L with r₀, Λ₊ should be close to 1.

The eigenvalues Λ₊ and Λ₋ can be found by making the Ansatz δx_ℓ = Λ^ℓδx₀, where x is r, s or J. Inserting this into equations (35) we obtain

while from δJ_ℓ = (2W_PFδs_{ℓ − 1} +  W_LPδJ_{ℓ − 1}) we obtain

Solving these equations under the assumption δr₀ ≠ 0, we find that Λ satisfies the quadratic equation

where B = W_FFW_PCW_CPF_ctxF_pul. Thus Λ₊ and Λ₋ satisfy the equation

It might seem that we have not gained much by adding the Pulvinar. As before we have 2 eigenvalues whose values depend on the network parameters, so it would seem that we again need fine-tuning of the parameters to get an output δr_L of the same order as δr₀. Like before, for these fine-tuned parameters we may obtain that r_L is comparable to r₀ for small r₀, but not if r₀ becomes larger.

However, this is not the case: If we take that WLP≫1 the two eigenvalues satisfy Λ₊ ≈ 2F_ctxW_FF(1 + W_CPW_PCF_pul) and Λ₋ ≈ 1. Thus in this case δr_L is comparable to δr₀, provided that |Λ₊| < 1.

When the activity in area L is plotted as a function of the input, the best relation is one that show a smooth linear increment of firing rate when the contrast gradually increase. This relation in the feedforward-pulvinar model is observed when the relation δr_ℓ ≈ δr₀ is satisfied, so Λ = 1. Therefore, we require 1 = 2F_ctxW_FF(1 + W_CPW_PCF_pul) when W_LP is large with the constraint of an arbitrary W_FP < W_LP. Without loss of generality, we consider when Ithcx=0, so F_ctx = 1/2, and we plug this result in Λ₁ to obtain

where [W_FF]_cr is the value of the cortical gain of the effective transfer function equation (2).

Analytical and Numerical results of the models – The analytical treatment has shown that the cortico-cortical and the cortico-pulvinar-cortical network behave differently when a spatially constant input is applied. In the case of the purely feedforward model, considering L large, the rate will approach r_∞ to fixed-point solutions depending of the value of the threshold and the gain, W_FF. In the unimodal transmission, Ith = 0, the solutions go to one stable rate to two stable and one unstable as one moves from small to large W_FF. This bifurcation also exists at the bimodal transmission, Ith ≠ 0, but here occurring twice as the gain is progressively increased given at the end 5 solutions, three of them stables and two unstable. If we plot r_∞ against r₀ at large W_FF in the bimodal transmission, a double step response appears while both inflections points will be at ±Ith_.

The meaning of these results is that for small W_FF when Ith = 0, the response, r_ℓ, always approaches 0 with increasing ℓ for any input r₀, while for large W_FF it approaches upper stable solution for r₀ > 0 and the lower one for r₀ < 0. For Ith ≠ 0, solutions approach similar as before, but for W_FF large 2 unstable fixed-points appears in r₀ = ±Ith/2 and the 2 stables responses move forward to r₀ = ±1. As Ith > 0.9, r_ℓ = 0 becomes stable and the other 4 solutions maintain the same previous stability. Thus, information about r₀ is lost in the higher areas for any value of W_FF either for Ith = 0 or otherwise.

Adding feedback connections only, W_FB ≠ 0 and W_CP = 0, does not qualitatively improve the situation. The bifurcation point is adjusted to [W_FF]_cr = 2/W_FF(2 + W_FB), but for a sufficiently large L we still have an almost constant output for small W_FF and a step or double step response for larger W_FF when Ith = 0 or Ith ≠ 0, respectively.

The response of the network to spatially constant input with the pulvinar included, W_CP ≠ 0 and W_PC ≠ 0, could also be treated analytically. One solution that will satisfy a smooth linear increment of firing rate when the contrast input gradually increases will show up. This relation in the feedforward-pulvinar model is observed when WLP≫1 with the constraint W_FP < W_LP. At Ithcx=0, this process is satisfied as W_FF = 1/(1 + F_pulW_CPW_PC), where W_FF is the cortical gain and F_pul the derivative of the effective transfer function for the pulvinar. This improvement in the last area’s activity is because pulvinar area ℓ receives input from all lower areas and passes directly to higher areas. Because of these long range interactions, responses in the higher pulvinar regions may not tend to a bimodal output distribution. This will be confirmed by numerical simulations.

In this part we investigate what are the values of parameters that better explain a smooth linear increment of firing rate at r_L when an input is gradually varied in contrast. The optimal output of layer L is one that spans largely the dynamic range of outputs and conserves as much as possible a linear relation with the input, r₀. In the simulations, when a homogeneous input −1, 1 is applied to the feedforward model, a large W_FF assures utilization of the whole dynamic range between –1 and 1, but yields an extremely non-linear curve. Instead, for small W_FF the curve that plots r_L against r₀ appears much more linear, but only covers a small part of the output range. So, to combine properly these two requirements we work out the entropy of the output distribution by assuming that the r₀ is distributed homogeneously between −1 and 1. The entropy is high both in the case where the input-output relation is linear and also when the output range is large (see Materials and Methods for the mathematical description).

We look for the parameters values that maximize the entropy of the rate distribution of the last layer. Because, most of the analytical treatment considered both unimodal and bimodal transmissions, we analyze separately the cases where Ith = 0 and Ith ≠ 0. Tables 1 and 2 recapitulate the results for the different models. We consider as networks the purely feedforward (FF), feedforward-feedback (FF-FB), feedforward-pulvinar (FF-Pul), and feedforward-feedback-pulvinar (FF-FB-Pul). The highest entropy is obtained for models that have Pul. The lowest entropy is observed for the purely feedforward model. However, qualitatively FF and FF-FB are almost the same for both Ith = 0 and Ith ≠ 0. In spite of the several results, for Ith = 0, FF-FB-Pul always was more informative than the feedforward and feedforward-feedback model. For Ith = 0, FF-Pul has this role. Figure 7 shows the output of the four models with the optimal parameters.

Compare with the cortico-cortical networks, a smooth linear increase of firing rate is present when the Pul is included and almost the whole output range is used. For the FF-Pul model this solution is obtained when WLP≫1 with W_FP < W_LP. For Ith = 0, W_FF was worked out with equation (40) and it was the best optimization value of the entire possibles values of variables analyzed (* in Table 1). To get these values for the networks with Pul, we also fit the reciprocal connections as W_CP > W_PC. Moreover, in the case of Ithcx≠0, W_CP were much stronger than that for the W_PC. However, without the reciprocal connectivity from the Pul the model present a low entropy finishing a constant rate in r^L. So that, the output for Pul models encodes optimally better an input than that for cortico-cortical networks.

In the network with a non-linear transfer function these requirements cannot be met exactly. In a purely feedforward model we can make the transfer function effectively linear by taking a small value of W_FF. This will ensure contrast invariance of the tuning, but in higher layers the response will by extremely small, so that, in the presence of noise a readout of the activity will give very little information about the stimulus. A larger W_FF will exploit the dynamic range of the system better, but introduces non-linearity in the response which may destroy the contrast invariance of the tuning and may make the contrast response function of the neurons in areas with large ℓ less smooth. A compromise between these two extremes needs to be made. We will now investigate how bad this compromise is in the purely feedforward cortical model and whether adding feedback connections and interactions with the pulvinar improves the network response.

By assuming visual propagation in the feedforward network with a non-linear transfer function, f, the response of the system either decreases or increases with ℓ as the threshold, Ith, and the feedforward strength, W_FF, are varied. We investigate this sequential propagation quantifying the firing rate from layer 1 to L, using histograms of distribution at two values of contrast, σ = 0.1 and σ = 1.0. As a first approximation, we analyze the response of each layer with L = 10 keeping Ith constant and varying W_FF.

Independent of the contrast, the visual input passes through the layers like in the case of linear transmission. The response of layer L first decreases to zero, then blows as the strength, W_FF, is increased gradually from small to large values. For small W_FF, any given value of contrast produces firing rate distributions that progressively evolve from a broad to narrow distribution as l → L. The activity moving higher through the system stays in small region of the firing rate range surrounding the mean (rk,iℓ)=0. As in the model with linear transfer function, the network integrates the visual input to a narrow distribution keeping only a small representation of the contrast. On the other hand, when W_FF is increased the network changes the modality of transmission. As W_FF becomes large, the rate distribution moves widens the range of activity and the distribution becomes multimodal as we move from layer to layer. As can be seen from Figure 8A at W_FF = 2, the distribution in the last layer is broad and unimodal for both low (σ = 0.1) and high (σ = 1.0) contrast. For larger W_FF, the distribution of the activity of units gradually becomes bimodal with peaks at the extrema of the response. At W_FF = 3, the visual input applied to l = 1 evolves sequentially through the network ending in net firing rates equal to –1 and 1 for most units. This tendency of the response to saturate at the borders is enhanced when W_FF becomes even larger. In each layer the activity moves progressively to r = −1 and r = 1 showing also a small peak around the value 0. Three clear solutions appear at W_FF ≤ 5. As we can see in Figure 8A at W_FF = 5, r¹ either starts with a broad distribution or already has net firing rates when a stimuli of low or high contrast is applied to the network. As the activity move from low to higher areas, large peaks are observed is at the maximum and minimum activity as well as a small cluster surrounding zero.

The distribution of the activity in layer L can be one of types. For small W_FF the distribution is clustered around 0 for all contrast levels. For an intermediate strength the distribution is broad for large σ and narrower for smaller σ. Finally for large W_FF the distribution is multimodal for all contrasts with peaks near −1, 0, and 1. These results can be readily understood from the transfer function of the units shown in Figure 8B. Small W_FF corresponds to a transfer function f whose maximal slope is less then 1∕2, resulting in responses that get progressively small as ℓ is increased. For intermediate values the average slope is close to the critical value for a significant range of inputs. For large W_FF the gain is much larger than 1∕2 for a significant range around 0 and saturates at ±1 for larger values. This results in the outputs being pushed to the extrema as ℓ increases. The small peak near 0 in the output distribution can be understood from the fact that in the input we are summing two inputs from the previous layer. If these are at the extrema, but have opposite sign, the total input will be close to 0, resulting in a response of almost 0.

We now explore how the response of the network depends on the threshold. We determine the distribution of the activity in layer L. These distributions are shown in Figure 9. Here, we test the effect of varying W_FF for five different levels of Ith. In each case the distribution of the activity has been analyzed at two values of the contrast, σ = 0.1 and σ = 1.0. To make the presentation more clear, we discuss separately the results for Ith = 0 and for Ith ≠ 0.

From the linear feedforward network, the know that the activity collapses to 0 if the transfer function has a derivative at 0, which is less than 1∕2. which means, for Ith = 0 that W_FF has to be larger than 2 to get a distribution with a finite width. This is confirmed in Figure 9. For W_FF less than approximately 2 the distribution is concentrated around 0. For larger values it is it is spread out and eventually becomes multimodal as explained above.

Now we analyze the layer L output when Ith ≠ 0. As in the previous case, the behavior of the system is pretty much the same on the interval 0 < Ith ≤ W_FF/2. The activity moves from one to three peaks at any value of contrast with a gradual transition between both extremes. As the threshold increases, qualitatively different behavior is observed. For σ = 0.1 the response is peaked around 0 for any value of W_FF. For σ = 1 the response distribution is broadened for W_FF > 1.4 when Ith = 3 W_FF/4, but this distribution narrows again for even larger values of W_FF. If we take Ith = W_FF, the output distribution is always narrow.

Why does the activity stays near 0 when Ith = W_FF? The assumption of a linear transfer function f produces a narrow distribution of rates in layer L for a small gain. For Ith ≠ 0, the function f always has a bimodal derivative as W_FF is sufficiently large, and the derivative near 0 becomes very small at 0. Thus for small inputs the response decreases as ℓ is increased. Whether there is a significant response for broader input distribution depends on whether a large enough fraction of the inputs is beyond the thresholds of the transfer function. For larger Ith are further apart so that is less likely that enough of the input exceeds the thresholds. Thus for Ith = 3W_FF/4 one can have a relatively broad output distribution for σ = 1, while for Ith = W_FF this is not the case.

In the analysis of the feedforward network we observed that the firing rate of units in layer L, rk,1L, stays near 0 when W_FF is small and approximately takes one of the three fixed solutions for large W_FF. In these cases the distribution of output rates is also nearly independent of σ. Only in a small range of W_FF the model shows a gradual variation of the distribution to changes in stimulus contrast. Adjusting the threshold does not produce any clear improvement. The respond of the feedforward model to variations in contrast is far from the desired result and we consider alternative network architectures to see whether they give an improvement. To do that we consider response rk,iL of the last cortical area, when feedback connections and interaction with the pulvinar are included. Moreover, we will not only consider the amplitude of the response over the network, but will also consider contrast invariance of the tuning and the contrast response function of the units. Here we take contrast invariance of the tuning to mean that the input elicits in the last layer an output vector rk,1L whose direction is independent of contrast. Smooth increase of the contrast response means that the length of this vector increases gradually with contrast. To account both these properties, we varied the optimization procedure used before for the spatially constant input.

To test whether the network response is sensitive to changes in contrast, we measure the response rk,1L when the input standard deviation, σ, is now gradually changed. We estimate two properties of the output firing rate: The output amplitude is a function of the contrast and the mean angle between different values of contrast. The output amplitude, F, is defined as average length of the response of layer L for LGN inputs with widths σ, F(σ)=|V→|, where the average is over input patterns with standard deviation σ. We aim for an amplitude function F that is as linear as possible and uses the dynamic range maximally. We use H_L defined as HL=∫01dσlog(F′(σ)) as the cost-function for this property. If H_L → L log2, the output scales linearly and exploits the whole dynamic ranges. It decreases if less of the dynamic range is used or the response is non-linearly. To explore whether the network can maintain contrast invariant tuning, we calculate the mean of distance S between normalizes output vectors e→(σ) and e→(σ′) for LGN inputs ri0=σxk and ri0=σ′xk respectively, S=∫01dσ∫01dσ′e→(σ)⋅e→(σ′). Here e→=V→∕||V→||. If S = 1 then vectors are in the same direction for all contrasts. As the direction changes more with contrast, S decreases. We define an error E=-2S-eHL that takes both these factors into account.

For networks with different architectures we determine the parameters which minimize the cost function E. We analyze separately when Ith = 0 and Ith ≠ 0. To optimize each network we have used Powell’s method in multiple dimensions (Press et al., 1992). Tables 3 and 4 show the results for the different models. We consider same systems from the homogeneous spatial input analysis, but considering more cases for Ith = 0. First in models which Ith = 0, we observe that the best minimization is produced by the network FF-FB-Pul with threshold for both cortex and Pul Ithctx,pul=0. The less improvement is for FF, nevertheless this network presents the best S value. Between these values, models that have both Ith = 0 are better than those that have only the Ithctx=0. Thus, for the cost function E, FFIth=0-FB-PulIth=0<FFIth=0-PulIth=0<FFIth=0-FB-Pul<FFIth=0-FB-Pul<FF-FB<FF. It is surprising that the network with Ithctx,pul=0 has an optimal enough minimization of E that it approaches the value for networks with Ith ≠ 0. The model with Ithctx,pul=0 does better in S but the range of magnitudes is smaller. It seems that the action of a shortcut between low and high cortical levels overcomes partly the problem of non-linearity regardless the presence of a threshold. For cases with Ith ≠ 0, the model that minimize the cost function better is the FF-FB-Pul network, while the purely feedforward system has the least optimization. However, models with feedback do not qualitatively change the cost, compared to those without feedback. For example, the system FF-FB-Pul has a small negative value of W_FB and its presence produces only a slight improvement in S. Indeed, networks with feedback input tend to converge toward small negative values, except for FF-FB. Here, the feedback input produces a slight improvement in the H_L and keeps almost the same value of S. For these networks, the cost function E is FF-FB-Pul < FF-Pul < FF-FB < FF.

To graphically observe how the networks optimize this error we plot the response of a neuron at contrast σ varied against the response of the same cell at contrast σ = 1. Only feedforward (Figure 10 Left) and feedforward-pulvinar networks (Figure 10 Right) are plotted. For both architectures the contrast invariance of the tuning reasonable good, as reflected by the fact that the points fall nearly on a straight line. However, for the feedforward the slope of this line does not change much as σ is varied, reflecting the fact that the response amplitude only changes weakly with the contrast. Furthermore the dynamic range is not fully used here.

In the cortico-pulvinar-cortical model the dynamical range is almost fully exploited and the response amplitude increases by almost a factor of 5 as the contrast is increased from σ = 0.1 to σ = 1. This is further illustrated in Figure 11. In Figure 11A we plot the average response amplitude against the contrast for both architectures. In Figure 11B the separation between the normalized response vector for contrast σ and the normalized response vector averaged over contrasts, is plotted against σ with Ith = 0 and Ith ≠ 0. For both models S varies over the range 0.99 to 0.93, but the response amplitude clearly increases more linearly and uses more of the dynamic range for the cortico-pulvinar-cortical model. When Ith = 0, the separation response is better for the FF network followed for the FF-FB system. The inclusion of Pul to those system produces a sharp tuning of the separation response while for both low and high contrast the amplitude decreases. For the case Ith ≠ 0, is clearly that the small negative feedback input produces a shift of the average response curve to the right, producing a more linear output in the FF-FB-Pul system. However, compare to the FF-Pul network, the separation response decreases in amplitude as a function of the contrast. Systems without Pul present a wider and higher tuning (separation) response to contrast, but differences in magnitude are qualitatively similar. So that, systems with Pul included always show a less cost function while the improvement is produced overall for an enlargement of the length response as a function of contrast. In almost all the cases, this improvement occurs as |WLP|≫1 and |W_FP| <  |W_LP|.

A recent work of Theyel et al. (2010) has shown that higher-order thalamic nuclei can drive the activity of cortex. With the optimization procedure we show that the best values are for the model with Pul-Cortex while |W_CP| > |W_PC|. A problem with the optimization analysis however is that we regard of qualitatively other good solutions that satisfy for a sufficient transmission. Then, models that included Pulvinar input can display another optimization when |W_CP| < |W_PC| or |W_CP| > |W_PC|. Thus, we investigate a simple case for the feedforward-pulvinar network when W_FF = W_FP and Ithcx=IthPul, and vary gradually W_CP and W_PC to observe whether exist more than a solution. As we can see in Figure 12 solutions that minimize the cost function of the network appear as a |W_CP|/|W_PC| ratio. We observe that an improvement is present when connections from cortex to pulvinar are negative (positive) while connections from pulvinar to cortex are positive (negative). Surprisingly there are at least two almost equally good solutions: in one |W_CP| is large and |W_PC| small. In the other |W_PC| is large and |W_CP| small. Our model does well when the cortex modulates the pulvinar while the pulvinar drives the cortex, but it does equally well when it is the other way around.