BACKGROUND
topographic ordering is an important feature of the visual system, which is conserved among many visual areas  <cit> .
thus, the projection from retina to superior colliculus  is established in a way, which retains neighbourhood relationships between neurons  <cit> .
this implies that two axons of retinal ganglion cells , which originate from neighbouring points in retina, terminate proximally in sc.
it is assumed that this facilitates visual processing, which involves wiring local to the termination zone  <cit> .
the mechanisms responsible for topographic ordering have been lately under thorough examination.
following the original suggestion by sperry  <cit> , it was shown that chemical labels play an essential role in formation of the map .
for the projection from retina to sc the eph family of receptor tyrosine kinases and their ligands ephrins were shown to be necessary for establishing correct topographic maps  <cit> .
the coordinate system is encoded chemically in retina through graded expression of the eph receptors by the rgcs.
thus, in mouse retina, two receptors of the family, epha <dig> and a <dig>  are expressed in the low nasal – high temporal gradient  <cit> .
the recipient coordinate system in the sc is established through high caudal – low rostral gradient of ephrin-a <dig> and a <dig> ligands  <cit> .
since rgc axons expressing epha receptors are repelled by high levels of ephrin-a ligands this system of reciprocal gradients allows sorting of the projecting axons in the order of increasing density of receptors, whereby contributing to the formation of topographic map  <cit>  .
thus, the system of reciprocal gradients is involved in formation of topographic representation along the nasal-temporal axis, albeit some additional fine-tuning is provided by activity-dependent mechanisms  <cit> .
in this study we address the results of gain-of-function experiments, in which the retinocollicular maps were modified by genetic manipulations  <cit> .
rgcs of the wild-type mouse express the lim homeobox gene islet <dig>   <cit> .
retina of a single animal is composed of two types of cells with regard to their expression of isl <dig> gene, isl2+ and isl2-, which are intermixed in roughly equal proportion throughout the rgc layer .
to test the mechanisms of the retinocollicular map formation brown et al.
<cit>  generated "knock-in" mice, in which the isl <dig> and epha <dig> genes are coexpressed.
this implies that each isl2+ rgc and its axons, in addition to epha <dig> and a <dig>  also expresses epha <dig>  not found in the wild-type rgcs.
the isl2- cells remain epha3-, as the wild-type cells.
by doing so brown et al.
<cit>  increased the total level of epha receptors in a given fraction of retinal cells.
since the overall level of ephas is increased in isl2+/epha3+ cells, axons of two neighboring cells, knock-in and wild-type, should terminate in quite different places in sc .
the knock-in cells, interacting more strongly with the repellent should terminate at the position of decreased density of ephrins, i.e.
more rostrally with respect to the wild-type cells.
the neighborhood relationships between axons should be lost, the new map should lose its continuous nature, and it should split into two maps: one for wild-type rgcs, one for knock-in cells.
this prediction was confirmed by experiments of brown et al.
<cit>  .
in addition to the observation of the overall map doubling in homozygous knock-ins , brown et al.
discovered a curious behavior of the map in heterozygous animals.
in these animals the exogenous levels of epha <dig> were reduced roughly by a factor of two with respect to the homozygous knock-ins .
in terms of the expression density of epha <dig> these animals stand between the wild-type and knock-in animals.
accordingly, the structure of the map resembles a hybrid of the wild-type and homozygous maps.
the more rostral part of the map is single-valued, similarly to the wild-type, whereas about 60% of the caudal-most part is double-valued, like in the homozygous animals.
this observation suggests that the map bifurcates somewhere between double-and single-valued regions.
although overall doubling of the map in homozygotes is easy to understand, any true model for the retinocollicular map formation should be able to account for the bifurcating behavior of map in heterozygotes.
therefore, experiments in heterozygotes represent a powerful tool to falsify various theoretical models.
brown et al.
<cit>  suggest that the bifurcating behavior of the map is consistent with the importance of relative rather than absolute values of the expression levels.
indeed, the relative difference of exogenous epha <dig> to endogenous epha5/ <dig> is maximal in nasal retina , where the doubled map is observed .
in the temporal retina  the epha <dig> to epha5/ <dig> ratio is not so large, which may account for the fact that the map is single-valued there.
thus a model for the topographic map from retina to sc should rely on the relative but not absolute levels of epha signaling.
the point, which we make in this study, is that more experimental tests are needed to justify the suggestion about relative expression levels.
to make our point clear we present a model for the retinocollicular map formation, which is based upon differences in the absolute values of eph/ephrin expression levels, rather than relative differences.
our model manages to reproduce all the essential features of experiments described in brown et al.
<cit> , including bifurcation of the map in heterozygotes.
in the model presented here the map is single-valued in rostral part of heterozygous sc due to inhomogeneous gradients of ligand and receptor, rather than reduced relative difference of epha receptors.
below we suggest experimental tests, which may distinguish these two classes of models.
this model was presented previously at society for neuroscience meeting  <dig> and on the arxiv preprint server  <cit> .
to test various hypotheses we use a model for retinocollicular map formation employing stochastic markov chain process.
our model is based upon three principles: chemoaffinity, axonal competition, and stochasticity.
some features of our model are similar to arrow model of hope, hammond, and gaze  <cit> .
the implementation of the model used here is available in  <cit> .
RESULTS
markov chain model
let us first describe the 1d version of the model.
we consider a linear chain of  <dig> rgc, each expressing individual level of epha receptors given by ra, where i =  <dig> .. <dig> is the rgc index, which also determines a discrete position of the cell in the retina.
we have verified that results presented below do not depend on the number of cells, as long as this number is large enough.
each rgc is attached by an axon to one and only terminal cell in sc, which has an expression level of ligand given by la, where k =  <dig> .. <dig> is the index in sc, also describing the terminal's topographic position.
the receptor density ra is an overall increasing function of its index i, while the ligand density la is decreasing, when going from k =  <dig>  to k =  <dig>  positions .
this determines the layout of chemical "tags" used to set up map's "topography".
an additional feature is that no two cells can project to the same spot in sc, which is meant to mimic axonal repulsion/competition for positive factors in sc, described in detail by ref.
<cit> .
we start from a random map, in which the terminal positions of all rgc axons in sc are chosen randomly.
we then modify the map probabilistically, using the following rule.
we consider two axons projecting to the neighboring points in sc .
we attempt to exchange these axons in sc with probability
here α >  <dig> is the parameter of our model.
the probability of the axons to stay unchanged pretain is determined from pexchange + pretain =  <dig> and is therefore given by
since the only difference between these probabilities is the sign in front of α, it is important to describe the nature of this sign.
assume that the product of gradients in eq.
is negative, i.e.
the gradients run in the opposite directions, which corresponds to the correct order of axonal terminals in sc.
then pexchange < 1/ <dig> and pexchange <pretain, i.e.
the probability or retaining the current ordering of the axonal pair is larger than changing it.
this is consistent with the chemorepellent interactions of receptors and ligands.
in the opposite case of the wrong order, i.e.
when the product of gradients in eq.
is positive and gradients run in the same directions, pexchange >pretain by the same reasoning.
the described process will tend to exchange the order of gradients and, therefore establish the correct order of topographic projections.
by using probabilities described by eqs.
and  we incorporate the chemoaffinity principle into our stochastic model.
this step is then repeated for another nearest neighbor couple, chosen randomly, and so on, until a stationary distribution of projections is reached.
such a process belongs to the class of markov chain processes, since transformations to the map are determined only by the present state of the mapping and are not otherwise affected by development history  <cit> .
let us first consider cases in which final distribution can be understood without the use of computer.
the model described by  can be solved exactly for at least two limiting cases: when α =  <dig> and when α is very large.
in the former case  the information about chemical labels cannot affect the solution, since it is multiplied by  <dig> in eq.
.
hence, the map is completely random .
in the latter case  the molecular cues are very strong.
they eventually produce solution in which the axons are perfectly sorted in sc in the order of increasing density of receptor .
an intermediate situation with certain finite value of parameter α is described by a compromise between noise and chemical cues, with former randomizing the map on the finer scale, while the latter inducing the overall correct ordering .
we conclude that mean position of projections is controlled by the chemical signal, while the spread of projections or the size of tz is determined by noise .
it should be noted that in the case of large α  our model is equivalent to the arrow model, introduced by hope, hammond, and gaze in  <cit> .
the arrow model uses exchanges between nearest axons if they terminate the wrong way in tectum/sc.
it can also include stochastic steps, as described in  <cit> .
the stochastic behaviors of the model described here  and the arrow model are not the same .
topographic maps in knock-in mice
maps in both wild-type mice and homozygote conditions  can be understood on the basis of axonal sorting in sc in the order of monotonously increasing levels of ephas.
results of such sorting are displayed in the bottom row in figure  <dig>  it is clear that maps resulting from simple sorting reproduce all essential features of mapping observed in wild-type and homozygotes.
at the same time, bifurcation, observed in heterozygotes , is not captured by simple sorting procedure.
this observation led us to develop the version of noisy sorting, based on chemical cues, described here.
why does the blending of two maps occur in rostral sc?
this question is addressed below in discussion section.
here we display manipulations with 'chemical tags', which can shift the bifurcation point in our model.
these manipulations point to two factors, which contribute to the location of the bifurcation point.
the first factor is inhomogeneous endogenous epha5/ <dig> density.
it results in a smaller separation between two branches of the map observed in rostral sc .
the second factor is a smaller gradient of ligand in rostral sc  than in caudal sc.
both these factors are addressed below.
let us now demonstrate the impact of inhomogeneous epha gradient on bifurcation.
figure  <dig>  shows that if the endogenous gradient of receptor is made more inhomogeneous, the point of bifurcation is shifted caudally.
this means that the single-valued part of the map becomes larger.
to see this, notice that the density of receptor in figure  <dig>  has a minimum value of about  <dig> .
the minimum value of the receptor density in figure  <dig>  is about twice as small.
hence the endogenous receptor density in figure 6a changes faster than in figure  <dig>  the results of simple sorting of axons according to increasing level of receptor are also shown in figure 6a .
two branches of the map approach each other closer than in figure 5b.
this is consistent with the expanded single-valued part of the map observed in figure 6a.
but is receptor distribution the sole determinant of the position of the bifurcation point?
to demonstrate that the latter is also controlled by the gradient of ligand in sc we reduce the density of ligand uniformly by 25% .
this may mimic experiments in which increment in rgc receptor is combined with reduction in ephrin-a ligand density.
as a result, the point of transition between single-valued and double valued parts of the map is located more caudally in figure 6b than in figure 5a.
hence, we can affect the point of transition by changing the densities of both receptor and ligand to a similar degree.
these results are explained below in the discussion section.
finally, we verify that increasing of the gradient of ligand leads to a small expansion of double-valued part of the map.
this result is demonstrated in figure 6c.
this example shows once again that ligand concentration can affect the position of bifurcation point and that increasing the ligand profile inhomogeneity  leads to a more pronounced bifurcation effect .
results for 2d model
we simulated 2d development using the hypothesis that another pair of chemical tags, ephb family of receptors and their ligands, ephrins-b, are responsible for establishing topographic projection from dorsal-ventral  axis on retina to lateral-medial axis in sc  <cit> .
ephb2/3/ <dig> are expressed in high-ventral-to-low-dorsal gradient by rgcs  <cit> , while ephrins-b are expressed in high-medial-to-low-lateral gradient in tectum/sc  <cit> .
since dorsal/ventral axons project to lateral/medial sc this implies attractive interactions between ephb+ axons and ephrin-b rich environment  <cit>  .
in our model the attractive interactions are modeled by the following exchange probability of two axonal terminals in the dv direction:
here rb, rb, lb, and lb are ephb receptor and ephrin-b ligand densities at neighboring points  <dig> and  <dig> in sc.
this probability is similar to eq.
.
notice a sign change compared to eq.
, which insures that pexchange >pretain if the order of gradients is wrong, i.e.
if the gradients of receptor and ligand are antiparallel.
by choosing this sign we therefore ensure attraction between axons and ligands.
the details of our simulations are described in methods.
our model allows not only exploration of two-dimensional maps  but also observing and modeling temporal development .
videos with detailed evolution of the map are available in  <cit> .
discussion
why does the map in heterozygotes bifurcate?
in our model the map is formed through interaction of three factors: eph/ephrin-based chemorepulsion/attraction, competition between axons for space, and noise.
it is the latter that fuses two maps together in rostral sc  in this model.
therefore, to understand position of the bifurcation point one has to consider the interplay between signal and noise at different positions in the map.
as we have shown above both ligand and receptor distributions influence the range of single-valued portion of the map independently .
let us first address the impact of ligand distribution.
figures  <dig> and 12a show that the gradient of ligand is the smallest in rostral sc.
this leads to a larger impact of noise there.
since noise drives blending of two branches of the map, this blending first occurs in rostral sc, in agreement with brown et al.
<cit> .
interestingly, brown et al.
also shows larger diameters of axonal tz in the single-valued part of the map, which is consistent with the larger impact of noise there.
the second factor contributing to bifurcation in heterozygotes is inhomogeneity of epha gradient in retina.
consider the case of no noise.
mapping in this case is obtained by sorting axonal terminals in the order of increasing density of epha .
separation between two maps is the smallest in rostral part .
this is because of inhomogeneous gradient of receptor in 'retinal' cells .
therefore, even if noise were the same in all parts of the map, rostral part has the smallest signal in terms of separation between two maps, and the largest potential to be blended by noise.
we conclude that two factors, increased noise and reduced signal, cooperate in rostral sc in fusing the wild-type and knock-in maps.
this leads to the formation of single-valued map there.
in caudal part both noise is reduced and distance between maps is larger.
hence, the map is double-valued in caudal sc.
mapping in isl2/ephb knock-ins has two bifurcations
it is possible to spatially separate these two blending factors, increased noise and decreased signal, if one applies the same logic to the dv axis of the map.
in our model this mapping is implemented by attractive interactions between ephb+ axons and ephrin-b rich environment.
hence, dv mapping is "flipped" with respect to the tn one in the sense that high ephb gradient region of retina maps to a high ephrin-b gradient region in sc.
two blending effects described above  are therefore spatially separated for the dv axis.
to observe mapping in these conditions we performed a numerical "experiment" on the isl2/ephb knock-in conditions.
this may have relevance to mapping in dv direction.
the results are shown in figure  <dig>
two bifurcations observed in figure  <dig> confirm the hypothesis about two factors operating in the numerical model.
the ventral bifurcation is associated with receptor, since separation between two maps in perfectly ordered conditions is the smallest in medial sc.
the second bifurcation, dorsal, occurs due to noise, since noise is maximal where the gradient of ligand is the smallest, i.e.
in lateral sc.
thus, we suggest that experiments on isl2/ephb knock-ins should make clear if inhomogeneity in receptor density or noise is more important.
it is also possible that activity-dependent mechanisms drive blending of two maps.
activity leads to focusing of projections whereby axons with close locations in retina are effectively attracted to each other in sc.
activity-dependent attraction will blend axons positioned proximally in sc, therefore ventral bifurcation, described above, may be robust with respect to these factors.
the dorsal bifurcation, on the other hand, may or may not be observable if activity-dependent focusing of projections takes place.
these questions will be addressed in future studies.
absolute versus relative
brown et al.
<cit>  demonstrates that retinocollicular mapping is based on relative levels of epha/ephrin-a expression in the broad meaning of this term.
indeed, the absolute value of epha density does not determine where an axon terminates in colliculus.
this is because axonal tz may shift in the presence of axons with altered expression of chemical tags.
for example, wild-type axons terminate more caudally in the presence of isl2+/epha3+ axons.
hence, an important factor is the presence of other axons, relative to which given axon establishes its termination point.
this idea is also evident from retinal and collicular/tectal ablation experiments in rodents  <cit>  and other species  <cit> .
can we take this idea to the next level and hypothesize that relative differences between neighboring retinal cells represent the chemical signal?
this suggestion was used  <cit>  to explain blending of the two maps in heterozygous rostral sc, since relative differences in receptor levels are the smallest in the corresponding part of retina .
in this study we present a model, which uses differences in absolute values of chemical label, as seen from eq.
.
indeed, in our model adding a constant value to all densities does not change the resulting mapping, since  depends only on differences in expression levels.
but this manipulation decreases the relative differences in the expression of ephas between neighboring knock-in and wild-type cells.
hence, our model is not based on relative differences between receptor densities.
yet, we demonstrate that it can account for experimental results in detail.
thus, we suggest that existing experimental evidence is not sufficient to distinguish relative and absolute labeling in the narrower sense.
of course, our model also accounts for the caudal displacement of wild-type tzs, thus resulting in a relative labeling system in the broad sense.
in the first approximation, this model performs a sorting procedure, understood mathematically, of the fibers based on the expression levels of epha.
our procedure uses differences in absolute values of epha densities rather than relative differences.
we suggest that more quantitative evidence is needed to distinguish these two "relativity principles".
relative labeling in the narrow sense can be incorporated in our model too, if coefficient α is a function of label densities.
thus, the condition α ∝ 1/ ensures the weber's law for axonal "perceptual thresholds", since chemical signal is proportional to the relative differences.
comparison to other theoretical models
theories based on chemoaffinity principle are reviewed in  <cit> .
some features of our approach are similar to the arrow models described in ref.
<cit> , which employs exchanges between neighbouring axons to establish ordered retinotectal/collicular maps.
at the same time some features of our model are different from the arrow model.
first, we employ information about chemical labels, ephs and ephrins.
at the heart of our model are equations , which rely on the known distributions of chemical labels.
these equations are unique to our approach.
as it was noted above, in the absence of stochasticity , when perfectly ordered map is formed, our 1d model with nearest neighbor exchanges is equivalent to the 1d version of the arrow model.
however, in the stochastic regime, description of developmental noise is different here.
in particular, we relate features of the map, to the distribution of chemical labels.
we argue that this feature is important in understanding experiments  <cit> , since distribution of labels determines where tzs fuse to form bifurcations.
second, we consider both nearest neighbour and distant neighbour exchanges .
indeed, eq.
can be applied to determine exchange probability for a pair of distant axons too.
this feature may be crucial, since development of map in the rc direction is determined by original primary axonal overshoot with subsequent retraction of inappropriate projections  <cit> .
in this process the interstitial branches of the same primary axon are eliminated and subsequently added non-locally.
we show below in method section that in many cases local and global exchanges produce the same results in terms of final distribution of projections.
but, the process of development is different in the local and global exchange cases.
prestige and willshaw  <cit>  suggested to divide developmental mechanisms into two groups.
in group i mechanisms each rgc axon has maximum affinity to certain unique point in the target, even without other axons.
in group ii mechanisms, the position of tz results from competitive interactions with other axons.
our model definitely belongs to the second group, since we assume that all axons experience maximum affinity to rostral medial sc and are spread over entire sc by competition.
our approach is similar to described in prestige and willshaw in the way how graded distributions of molecular tags are represented.
the details of map modifications are somewhat different in this study and are precisely defined by eqs.
.
in a recent study honda  <cit>  considered results of experiments  <cit> .
he used servomechanism model to explain the overall structure of the maps in mutants.
servomechanism model is a hybrid between group i and ii models in terminology of prestige and willshaw, since it assumes that axons have equilibrium points in sc and they are subject to competition with each other.
although ref.
<cit>  reproduces doubling of the map in homozygotes it does not succeed in obtaining the bifurcation observed in heterozygotes, which is one of the purposes of the present study.
on the biological realism
when dealing with numerical simulations one always faces the question of the degree of realism with which to model the data.
does one have to model behaviors of individual atoms, or description on the level of axons is sufficient?
in this work we choose the level of description on the basis of what is known about this system.
we realize that our model does not capture many behaviors, but we argue that the mechanisms involved are unclear at the moment to be incorporated into a more detailed model.
our approach also fulfils its original goal, which is to reproduce the results of experiments  <cit>  and to generate experimentally testable predictions, thus satisfying the requirement of parsimony.
model presented here does not describe the difference between development along tn and dv axes.
the former mapping is controlled by original axonal overshoot along the rc direction in sc, with subsequent elimination of topographically inappropriate projections  <cit> .
in contrast, primary axons from the same dv retinal position enter sc in a broad distribution along ml axis.
topographically precise termination is provided by producing additional interstitial branches in the ml direction  <cit> .
these findings cannot be reproduced by our model, since no distinction is made between the primary rgc axon and its branches.
instead, our model deals with terminal points of interstitial branches produced by rgc axons.
CONCLUSIONS
we present a model for retinocollicular map development, which can account intriguing behaviors observed in gain-of-function experiments by brown et al.
<cit> , including bifurcation in heterozygous isl2/epha <dig> knock-ins.
the model is based on chemoaffinity, axonal repulsion/competition, and stochasticity.
we discuss possible mappings in ephrina-/isl2+/epha3+ knock-out/ins and isl2/ephb knock-ins.

 topographic ordering is an important feature of the visual system, which is conserved among many visual areas. thus, the projection from retina to superior colliculus is established in a way, which retains neighbourhood relationships between neurons. this implies that two axons of retinal ganglion cells, which originate from neighbouring points in retina, terminate proximally in sc. it is assumed that this facilitates visual processing, which involves wiring local to the termination zone. the mechanisms responsible for topographic ordering have been lately under thorough examination. following the original suggestion by sperry, it was shown that chemical labels play an essential role in formation of the map. for the projection from retina to sc the eph family of receptor tyrosine kinases and their ligands ephrins were shown to be necessary for establishing correct topographic maps. the coordinate system is encoded chemically in retina through graded expression of the eph receptors by the rgcs. thus, in mouse retina, two receptors of the family, epha and a are expressed in the low nasal – high temporal gradient. the recipient coordinate system in the sc is established through high caudal – low rostral gradient of ephrin-a and a ligands. since rgc axons expressing epha receptors are repelled by high levels of ephrin-a ligands this system of reciprocal gradients allows sorting of the projecting axons in the order of increasing density of receptors, whereby contributing to the formation of topographic map. thus, the system of reciprocal gradients is involved in formation of topographic representation along the nasal-temporal axis, albeit some additional fine-tuning is provided by activity-dependent mechanisms. in this study we address the results of gain-of-function experiments, in which the retinocollicular maps were modified by genetic manipulations. rgcs of the wild-type mouse express the lim homeobox gene islet. retina of a single animal is composed of two types of cells with regard to their expression of isl gene, isl2+ and isl2-, which are intermixed in roughly equal proportion throughout the rgc layer. to test the mechanisms of the retinocollicular map formation brown et al. generated "knock-in" mice, in which the isl and epha genes are coexpressed. this implies that each isl2+ rgc and its axons, in addition to epha and a also expresses epha not found in the wild-type rgcs. the isl2- cells remain epha3-, as the wild-type cells. by doing so brown et al. increased the total level of epha receptors in a given fraction of retinal cells. since the overall level of ephas is increased in isl2+/epha3+ cells, axons of two neighboring cells, knock-in and wild-type, should terminate in quite different places in sc. the knock-in cells, interacting more strongly with the repellent should terminate at the position of decreased density of ephrins, i.e. more rostrally with respect to the wild-type cells. the neighborhood relationships between axons should be lost, the new map should lose its continuous nature, and it should split into two maps: one for wild-type rgcs, one for knock-in cells. this prediction was confirmed by experiments of brown et al.. in addition to the observation of the overall map doubling in homozygous knock-ins, brown et al. discovered a curious behavior of the map in heterozygous animals. in these animals the exogenous levels of epha were reduced roughly by a factor of two with respect to the homozygous knock-ins. in terms of the expression density of epha these animals stand between the wild-type and knock-in animals. accordingly, the structure of the map resembles a hybrid of the wild-type and homozygous maps. the more rostral part of the map is single-valued, similarly to the wild-type, whereas about 60% of the caudal-most part is double-valued, like in the homozygous animals. this observation suggests that the map bifurcates somewhere between double-and single-valued regions. although overall doubling of the map in homozygotes is easy to understand, any true model for the retinocollicular map formation should be able to account for the bifurcating behavior of map in heterozygotes. therefore, experiments in heterozygotes represent a powerful tool to falsify various theoretical models. brown et al. suggest that the bifurcating behavior of the map is consistent with the importance of relative rather than absolute values of the expression levels. indeed, the relative difference of exogenous epha to endogenous epha5/ is maximal in nasal retina, where the doubled map is observed. in the temporal retina the epha to epha5/ ratio is not so large, which may account for the fact that the map is single-valued there. thus a model for the topographic map from retina to sc should rely on the relative but not absolute levels of epha signaling. the point, which we make in this study, is that more experimental tests are needed to justify the suggestion about relative expression levels. to make our point clear we present a model for the retinocollicular map formation, which is based upon differences in the absolute values of eph/ephrin expression levels, rather than relative differences. our model manages to reproduce all the essential features of experiments described in brown et al., including bifurcation of the map in heterozygotes. in the model presented here the map is single-valued in rostral part of heterozygous sc due to inhomogeneous gradients of ligand and receptor, rather than reduced relative difference of epha receptors. below we suggest experimental tests, which may distinguish these two classes of models. this model was presented previously at society for neuroscience meeting and on the arxiv preprint server. to test various hypotheses we use a model for retinocollicular map formation employing stochastic markov chain process. our model is based upon three principles: chemoaffinity, axonal competition, and stochasticity. some features of our model are similar to arrow model of hope, hammond, and gaze. the implementation of the model used here is available in.  markov chain model let us first describe the 1d version of the model. we consider a linear chain of rgc, each expressing individual level of epha receptors given by ra, where i =.. is the rgc index, which also determines a discrete position of the cell in the retina. we have verified that results presented below do not depend on the number of cells, as long as this number is large enough. each rgc is attached by an axon to one and only terminal cell in sc, which has an expression level of ligand given by la, where k =.. is the index in sc, also describing the terminal's topographic position. the receptor density ra is an overall increasing function of its index i, while the ligand density la is decreasing, when going from k = to k = positions. this determines the layout of chemical "tags" used to set up map's "topography". an additional feature is that no two cells can project to the same spot in sc, which is meant to mimic axonal repulsion/competition for positive factors in sc, described in detail by ref.. we start from a random map, in which the terminal positions of all rgc axons in sc are chosen randomly. we then modify the map probabilistically, using the following rule. we consider two axons projecting to the neighboring points in sc. we attempt to exchange these axons in sc with probability here α > is the parameter of our model. the probability of the axons to stay unchanged pretain is determined from pexchange + pretain = and is therefore given by since the only difference between these probabilities is the sign in front of α, it is important to describe the nature of this sign. assume that the product of gradients in eq. is negative, i.e. the gradients run in the opposite directions, which corresponds to the correct order of axonal terminals in sc. then pexchange < 1/ and pexchange <pretain, i.e. the probability or retaining the current ordering of the axonal pair is larger than changing it. this is consistent with the chemorepellent interactions of receptors and ligands. in the opposite case of the wrong order, i.e. when the product of gradients in eq. is positive and gradients run in the same directions, pexchange >pretain by the same reasoning. the described process will tend to exchange the order of gradients and, therefore establish the correct order of topographic projections. by using probabilities described by eqs. and we incorporate the chemoaffinity principle into our stochastic model. this step is then repeated for another nearest neighbor couple, chosen randomly, and so on, until a stationary distribution of projections is reached. such a process belongs to the class of markov chain processes, since transformations to the map are determined only by the present state of the mapping and are not otherwise affected by development history. let us first consider cases in which final distribution can be understood without the use of computer. the model described by can be solved exactly for at least two limiting cases: when α = and when α is very large. in the former case the information about chemical labels cannot affect the solution, since it is multiplied by in eq.. hence, the map is completely random. in the latter case the molecular cues are very strong. they eventually produce solution in which the axons are perfectly sorted in sc in the order of increasing density of receptor. an intermediate situation with certain finite value of parameter α is described by a compromise between noise and chemical cues, with former randomizing the map on the finer scale, while the latter inducing the overall correct ordering. we conclude that mean position of projections is controlled by the chemical signal, while the spread of projections or the size of tz is determined by noise. it should be noted that in the case of large α our model is equivalent to the arrow model, introduced by hope, hammond, and gaze in. the arrow model uses exchanges between nearest axons if they terminate the wrong way in tectum/sc. it can also include stochastic steps, as described in. the stochastic behaviors of the model described here and the arrow model are not the same. topographic maps in knock-in mice maps in both wild-type mice and homozygote conditions can be understood on the basis of axonal sorting in sc in the order of monotonously increasing levels of ephas. results of such sorting are displayed in the bottom row in figure it is clear that maps resulting from simple sorting reproduce all essential features of mapping observed in wild-type and homozygotes. at the same time, bifurcation, observed in heterozygotes, is not captured by simple sorting procedure. this observation led us to develop the version of noisy sorting, based on chemical cues, described here. why does the blending of two maps occur in rostral sc? this question is addressed below in discussion section. here we display manipulations with'chemical tags', which can shift the bifurcation point in our model. these manipulations point to two factors, which contribute to the location of the bifurcation point. the first factor is inhomogeneous endogenous epha5/ density. it results in a smaller separation between two branches of the map observed in rostral sc. the second factor is a smaller gradient of ligand in rostral sc than in caudal sc. both these factors are addressed below. let us now demonstrate the impact of inhomogeneous epha gradient on bifurcation. figure shows that if the endogenous gradient of receptor is made more inhomogeneous, the point of bifurcation is shifted caudally. this means that the single-valued part of the map becomes larger. to see this, notice that the density of receptor in figure has a minimum value of about. the minimum value of the receptor density in figure is about twice as small. hence the endogenous receptor density in figure 6a changes faster than in figure the results of simple sorting of axons according to increasing level of receptor are also shown in figure 6a. two branches of the map approach each other closer than in figure 5b. this is consistent with the expanded single-valued part of the map observed in figure 6a. but is receptor distribution the sole determinant of the position of the bifurcation point? to demonstrate that the latter is also controlled by the gradient of ligand in sc we reduce the density of ligand uniformly by 25%. this may mimic experiments in which increment in rgc receptor is combined with reduction in ephrin-a ligand density. as a result, the point of transition between single-valued and double valued parts of the map is located more caudally in figure 6b than in figure 5a. hence, we can affect the point of transition by changing the densities of both receptor and ligand to a similar degree. these results are explained below in the discussion section. finally, we verify that increasing of the gradient of ligand leads to a small expansion of double-valued part of the map. this result is demonstrated in figure 6c. this example shows once again that ligand concentration can affect the position of bifurcation point and that increasing the ligand profile inhomogeneity leads to a more pronounced bifurcation effect. results for 2d model we simulated 2d development using the hypothesis that another pair of chemical tags, ephb family of receptors and their ligands, ephrins-b, are responsible for establishing topographic projection from dorsal-ventral axis on retina to lateral-medial axis in sc. ephb2/3/ are expressed in high-ventral-to-low-dorsal gradient by rgcs, while ephrins-b are expressed in high-medial-to-low-lateral gradient in tectum/sc. since dorsal/ventral axons project to lateral/medial sc this implies attractive interactions between ephb+ axons and ephrin-b rich environment. in our model the attractive interactions are modeled by the following exchange probability of two axonal terminals in the dv direction: here rb, rb, lb, and lb are ephb receptor and ephrin-b ligand densities at neighboring points and in sc. this probability is similar to eq.. notice a sign change compared to eq., which insures that pexchange >pretain if the order of gradients is wrong, i.e. if the gradients of receptor and ligand are antiparallel. by choosing this sign we therefore ensure attraction between axons and ligands. the details of our simulations are described in methods. our model allows not only exploration of two-dimensional maps but also observing and modeling temporal development. videos with detailed evolution of the map are available in. discussion why does the map in heterozygotes bifurcate? in our model the map is formed through interaction of three factors: eph/ephrin-based chemorepulsion/attraction, competition between axons for space, and noise. it is the latter that fuses two maps together in rostral sc in this model. therefore, to understand position of the bifurcation point one has to consider the interplay between signal and noise at different positions in the map. as we have shown above both ligand and receptor distributions influence the range of single-valued portion of the map independently. let us first address the impact of ligand distribution. figures and 12a show that the gradient of ligand is the smallest in rostral sc. this leads to a larger impact of noise there. since noise drives blending of two branches of the map, this blending first occurs in rostral sc, in agreement with brown et al.. interestingly, brown et al. also shows larger diameters of axonal tz in the single-valued part of the map, which is consistent with the larger impact of noise there. the second factor contributing to bifurcation in heterozygotes is inhomogeneity of epha gradient in retina. consider the case of no noise. mapping in this case is obtained by sorting axonal terminals in the order of increasing density of epha. separation between two maps is the smallest in rostral part. this is because of inhomogeneous gradient of receptor in'retinal' cells. therefore, even if noise were the same in all parts of the map, rostral part has the smallest signal in terms of separation between two maps, and the largest potential to be blended by noise. we conclude that two factors, increased noise and reduced signal, cooperate in rostral sc in fusing the wild-type and knock-in maps. this leads to the formation of single-valued map there. in caudal part both noise is reduced and distance between maps is larger. hence, the map is double-valued in caudal sc. mapping in isl2/ephb knock-ins has two bifurcations it is possible to spatially separate these two blending factors, increased noise and decreased signal, if one applies the same logic to the dv axis of the map. in our model this mapping is implemented by attractive interactions between ephb+ axons and ephrin-b rich environment. hence, dv mapping is "flipped" with respect to the tn one in the sense that high ephb gradient region of retina maps to a high ephrin-b gradient region in sc. two blending effects described above are therefore spatially separated for the dv axis. to observe mapping in these conditions we performed a numerical "experiment" on the isl2/ephb knock-in conditions. this may have relevance to mapping in dv direction. the results are shown in figure two bifurcations observed in figure confirm the hypothesis about two factors operating in the numerical model. the ventral bifurcation is associated with receptor, since separation between two maps in perfectly ordered conditions is the smallest in medial sc. the second bifurcation, dorsal, occurs due to noise, since noise is maximal where the gradient of ligand is the smallest, i.e. in lateral sc. thus, we suggest that experiments on isl2/ephb knock-ins should make clear if inhomogeneity in receptor density or noise is more important. it is also possible that activity-dependent mechanisms drive blending of two maps. activity leads to focusing of projections whereby axons with close locations in retina are effectively attracted to each other in sc. activity-dependent attraction will blend axons positioned proximally in sc, therefore ventral bifurcation, described above, may be robust with respect to these factors. the dorsal bifurcation, on the other hand, may or may not be observable if activity-dependent focusing of projections takes place. these questions will be addressed in future studies. absolute versus relative brown et al. demonstrates that retinocollicular mapping is based on relative levels of epha/ephrin-a expression in the broad meaning of this term. indeed, the absolute value of epha density does not determine where an axon terminates in colliculus. this is because axonal tz may shift in the presence of axons with altered expression of chemical tags. for example, wild-type axons terminate more caudally in the presence of isl2+/epha3+ axons. hence, an important factor is the presence of other axons, relative to which given axon establishes its termination point. this idea is also evident from retinal and collicular/tectal ablation experiments in rodents and other species. can we take this idea to the next level and hypothesize that relative differences between neighboring retinal cells represent the chemical signal? this suggestion was used to explain blending of the two maps in heterozygous rostral sc, since relative differences in receptor levels are the smallest in the corresponding part of retina. in this study we present a model, which uses differences in absolute values of chemical label, as seen from eq.. indeed, in our model adding a constant value to all densities does not change the resulting mapping, since depends only on differences in expression levels. but this manipulation decreases the relative differences in the expression of ephas between neighboring knock-in and wild-type cells. hence, our model is not based on relative differences between receptor densities. yet, we demonstrate that it can account for experimental results in detail. thus, we suggest that existing experimental evidence is not sufficient to distinguish relative and absolute labeling in the narrower sense. of course, our model also accounts for the caudal displacement of wild-type tzs, thus resulting in a relative labeling system in the broad sense. in the first approximation, this model performs a sorting procedure, understood mathematically, of the fibers based on the expression levels of epha. our procedure uses differences in absolute values of epha densities rather than relative differences. we suggest that more quantitative evidence is needed to distinguish these two "relativity principles". relative labeling in the narrow sense can be incorporated in our model too, if coefficient α is a function of label densities. thus, the condition α ∝ 1/ ensures the weber's law for axonal "perceptual thresholds", since chemical signal is proportional to the relative differences. comparison to other theoretical models theories based on chemoaffinity principle are reviewed in. some features of our approach are similar to the arrow models described in ref., which employs exchanges between neighbouring axons to establish ordered retinotectal/collicular maps. at the same time some features of our model are different from the arrow model. first, we employ information about chemical labels, ephs and ephrins. at the heart of our model are equations, which rely on the known distributions of chemical labels. these equations are unique to our approach. as it was noted above, in the absence of stochasticity, when perfectly ordered map is formed, our 1d model with nearest neighbor exchanges is equivalent to the 1d version of the arrow model. however, in the stochastic regime, description of developmental noise is different here. in particular, we relate features of the map, to the distribution of chemical labels. we argue that this feature is important in understanding experiments, since distribution of labels determines where tzs fuse to form bifurcations. second, we consider both nearest neighbour and distant neighbour exchanges. indeed, eq. can be applied to determine exchange probability for a pair of distant axons too. this feature may be crucial, since development of map in the rc direction is determined by original primary axonal overshoot with subsequent retraction of inappropriate projections. in this process the interstitial branches of the same primary axon are eliminated and subsequently added non-locally. we show below in method section that in many cases local and global exchanges produce the same results in terms of final distribution of projections. but, the process of development is different in the local and global exchange cases. prestige and willshaw suggested to divide developmental mechanisms into two groups. in group i mechanisms each rgc axon has maximum affinity to certain unique point in the target, even without other axons. in group ii mechanisms, the position of tz results from competitive interactions with other axons. our model definitely belongs to the second group, since we assume that all axons experience maximum affinity to rostral medial sc and are spread over entire sc by competition. our approach is similar to described in prestige and willshaw in the way how graded distributions of molecular tags are represented. the details of map modifications are somewhat different in this study and are precisely defined by eqs.. in a recent study honda considered results of experiments. he used servomechanism model to explain the overall structure of the maps in mutants. servomechanism model is a hybrid between group i and ii models in terminology of prestige and willshaw, since it assumes that axons have equilibrium points in sc and they are subject to competition with each other. although ref. reproduces doubling of the map in homozygotes it does not succeed in obtaining the bifurcation observed in heterozygotes, which is one of the purposes of the present study. on the biological realism when dealing with numerical simulations one always faces the question of the degree of realism with which to model the data. does one have to model behaviors of individual atoms, or description on the level of axons is sufficient? in this work we choose the level of description on the basis of what is known about this system. we realize that our model does not capture many behaviors, but we argue that the mechanisms involved are unclear at the moment to be incorporated into a more detailed model. our approach also fulfils its original goal, which is to reproduce the results of experiments and to generate experimentally testable predictions, thus satisfying the requirement of parsimony. model presented here does not describe the difference between development along tn and dv axes. the former mapping is controlled by original axonal overshoot along the rc direction in sc, with subsequent elimination of topographically inappropriate projections. in contrast, primary axons from the same dv retinal position enter sc in a broad distribution along ml axis. topographically precise termination is provided by producing additional interstitial branches in the ml direction. these findings cannot be reproduced by our model, since no distinction is made between the primary rgc axon and its branches. instead, our model deals with terminal points of interstitial branches produced by rgc axons.  we present a model for retinocollicular map development, which can account intriguing behaviors observed in gain-of-function experiments by brown et al., including bifurcation in heterozygous isl2/epha knock-ins. the model is based on chemoaffinity, axonal repulsion/competition, and stochasticity. we discuss possible mappings in ephrina-/isl2+/epha3+ knock-out/ins and isl2/ephb knock-ins. 
we present a theoretical model for retinocollicular map development, which can account for intriguing behaviors observed in gain-of-function experiments by brown et al., including bifurcation in heterozygous isl2/epha knock-ins. in our model the map undergoes a series of reconstructions with probabilities dependent upon a set of chemical cues. we study the map formation using stochastic model based on markov chains. first, the inhomogeneous gradient of endogenous receptor in retina makes the impact of exogenous receptor less significant in temporal retina. we examine results of gain-of-function experiments on retinocollicular maps in knock-in mice. the model is based on known chemical labels, axonal repulsion/competition, and stochasticity. possible mapping in isl2/ephb knock-ins is also discussed. the established map is single-valued, which implies that each point in retina maps to a unique termination zone in superior colliculus. in homozygous isl2/epha knock-in mice the map is double-valued, which means that each point on retina maps to two termination zones in superior colliculus. in wild-type mice the temporal-nasal axis of retina is mapped to the rostral-caudal axis of superior colliculus. second, the gradient of ephrin in the corresponding region of superior colliculus is smaller, which reduces the chemical signal-to-noise ratio. we predict that if gradient of ephrin is reduced by a genetic manipulation, the single-valued region of the map should extend to a larger portion of temporal retina, i.e.
 we examine results of gain-of-function experiments on retinocollicular maps in knock-in mice. in wild-type mice the temporal-nasal axis of retina is mapped to the rostral-caudal axis of superior colliculus. the established map is single-valued, which implies that each point in retina maps to a unique termination zone in superior colliculus. in homozygous isl2/epha knock-in mice the map is double-valued, which means that each point on retina maps to two termination zones in superior colliculus. this is because about percent of cells in retina express isl and two types of projections, wild-type and isl2/epha positive, form two branches of the map. in heterozygous isl2/epha knock-ins the map is intermediate between the homozygous and wild-type: it is single-valued in temporal and double-valued in the nasal parts of retina. in this study we address possible reasons for such a bifurcation of the map.  we study the map formation using stochastic model based on markov chains. in our model the map undergoes a series of reconstructions with probabilities dependent upon a set of chemical cues. our model suggests that the map in heterozygotes is single-valued in temporal region of retina for two reasons. first, the inhomogeneous gradient of endogenous receptor in retina makes the impact of exogenous receptor less significant in temporal retina. second, the gradient of ephrin in the corresponding region of superior colliculus is smaller, which reduces the chemical signal-to-noise ratio. we predict that if gradient of ephrin is reduced by a genetic manipulation, the single-valued region of the map should extend to a larger portion of temporal retina, i.e. the point of transition between single-and doulble-valued maps should move to a more nasal position in isl2-epha heterozygotes.  we present a theoretical model for retinocollicular map development, which can account for intriguing behaviors observed in gain-of-function experiments by brown et al., including bifurcation in heterozygous isl2/epha knock-ins. the model is based on known chemical labels, axonal repulsion/competition, and stochasticity. possible mapping in isl2/ephb knock-ins is also discussed. 