Virtual Tissue Engineering and Optic Pathways: Plotting the Course of the Axons in the Retinal Nerve Fiber Layer

Abstract

Purpose.: As part of a larger project on virtual tissue engineering of the optic pathways, we describe the conditions that guide axons extending from the retina to the optic nerve head and formulate algorithms that meet such conditions. To find the entrance site on the optic nerve head of each axon, we challenge the fibers to comply with current models of axonal pathfinding.

Methods.: First, we build a retinal map using a single type of retinal ganglion cell (RGC) using density functions from the literature. Dendritic arbors are equated to receptive fields. Shape and size of retinal surface and optic nerve head (ONH) are defined. A computer model relates each soma to the corresponding entry point of its axon into the optic disc. Weights are given to the heuristics that guide the preference entry order in the nerve.

Results.: Retinal ganglion cells from the area centralis saturate the temporal section of the disc. Retinal ganglion cells temporal to the area centralis curve their paths surrounding the fovea; some of these cells enter the disc centrally rather than peripherally. Nasal regions of the disc receive mixed axons from the far periphery of the temporal hemiretina, together with axons from the nasal half. The model plots the course of the axon using Bezier curves and compares them with clinical data, for a coincidence level of 86% or higher.

Conclusions.: Our model is able to simulate basic data of the early optic pathways including certain singularities and to mimic mechanisms operating during development, such as timing and fasciculation.

Introduction

Knowledge of nerve-fiber-bundle trajectories in the retina is considered a prerequisite for any form of computer-assisted perimetry, including new diagnostic techniques such as fundus-oriented perimetry (FOP) ¹ and scotoma-oriented perimetry (SCOPE; Paetzold J, et al. IOVS 2005;46:ARVO E-Abstract 636) as well as for visual-field modeling. ^2,3 Comparing standard automatic perimetry (SAP) sensitivity loss with structural damage—measured by any state-of-the-art clinical device—requires a map that relates local regions on the retina to local regions on the optic disc. Maps of the spatial correspondence of retinal locations to peripapillary retinal nerve fiber layer (RNFL) sectors are an active area of structure-function research. ³ Accurate analysis of visual-field defects resulting from localized optic-nerve damage rests in precise topographical localization of the retinal ganglion cell (RGC) somas, the size of the receptive fields, the course of their axons on the retinal surface, and the point of entrance into the optic-nerve head. 

Ever since the mid-19th century, with works on ophthalmoscopy by Hemholtz and campimetry by Von Graefe, ⁵ structure–function relationship has been an active area of intensive clinical research. ^6–9 Experimental studies on the course of axons on the optic nerve also predate the 20th century, as in the case of Usher and Dean. ¹⁰ Knowledge of the exact course of axons on the retinal surface has been gained through considerable debate. Traquair's mention that segregation of crossed and uncrossed chiasmal bundles took place just behind the eyeball moved Posner and Schlossman ¹¹ to assume that fibers coming from the area of the retina between the fovea and the optic disc had to enter the optic nerve on the nasal side of the disc, an error refuted by Hoyt and Tudor ¹² by experimental photocoagulation in a cynomolgus monkey. Also in a photocoagulation-based study in owls and rhesus monkeys, Radius and Anderson ¹³ found that, at any point of the nerve fiber layer (NFL), axons from the periphery lie deeply, close to the somas of RGCs, and that axons from more centrally situated RGCs lie more superficially, near the vitreous-retinal interface. Similar results were reported by Minckler ¹⁴ injecting horseradish peroxidase into the retina of Macaca fasicularis . Additionally, they mentioned that the locations of labeled ganglion cells could be accurately predicted by following the course of the arcuate NFL reflexes in fundus photographs. Figure 1 is based on Minckler's diagram, showing the speculative lamination of axons within the retina. 

Figure 1

 

A consequence of the retinotopic order of the axons in the retina simplified in an idealized diagram. The axons coming from the peripheral cells send axons to the periphery of the optic stalk, while axons coming from cells closer to the future optic nerve enter through the center of the optic stalk. An exception to this rule, discussed in the text, are midperipheral RGCs, from the temporal region of the raphe close to the macula. Axons from those RGCs travel on top of the retinal nerve fiber layer, and enter through the center of the optic stalk, close to the rim of the cup.

One fundamental working assumption in this field has been that of retinotopy—that is, that the topographic order of retinal ganglion cell axons along the visual pathway is determinant for the formation of topographic maps within the lateral geniculate nucleus and the visual cortex. ¹⁵ How the retinal surface is ordered, how this order relates to central visual structures, and how this order is acquired during development has been greatly debated over recent decades (see Ref. 16 for review). The notion of retinotopy has, nevertheless, been challenged by some histologic reports. Early autoradiographic studies by retinal injection of radionucleotides ¹⁷ have shown no apparent stratification of peripheral versus central axons. Similar results were reported by Hoyt and Luis ¹⁸ and Naito. ¹⁹ Fitzgibbon and Taylor, ¹⁵ working with donor human eyes and incrustation of crystals of dyes at different retinal locations confirmed some of the tenets of retinotopy, but also defended the idea that peripheral retinal axons were scattered throughout the thickness of the nerve fiber layer. As results depend heavily on the techniques employed, a criticism of those studies based on the spreading of markers (up to 24-week incubation periods in the case of Fitzgibbon and Taylor) is that diffusion of the intraretinally inserted label may reduce the resolution of the technique. ²⁰  

In clinical terms, ever since the histologic studies of Radius, ^20,21 which established the relationship between glaucomatous injury to axon bundles within the lamina and paracentral scotomas, attempts have been made to locate the portions of the neuroretinal rim affected by damage in the RGCs and the corresponding visual-field defects. Reportedly, a simple linear relationship describes the data relating SAP to a particular structural test such as optical coherence tomography (OCT) ²² or spectral-domain (SD)-OCT. ^23,24 As the location of nerve fiber–layer defects and field defects appear to correspond closely, ²⁵ an important task has been to relate delimited areas of the retina, which projects in clinically defined units of the visual field, to sectors in the optic disc or sectors around the disc, as given in the RNFL profile. ^3,26–28 Axon bundles are delimited by Müller cells and are very thin, compared to the resolution of visual-field testing. According to Radius and the results of Anderson, ¹³ the dimensions of fiber bundles in the glia tunnels delimited by Müller cell processes approximate those of the fiber-bundle striations seen ophthalmoscopically. According to the micrographs of Radius, this would approximate 30 μm per bundle. However, probably due to the inherent imprecision of the measuring devices, and despite using very large regions of the visual field, Garway-Heath et al. ²⁸ noted that the uncertainty of measurements following the course of the retinal-nerve-fiber bundles and bundle defects up to the disc margin was on the order of 30°. 

Jansonius et al. ^29,30 developed a mathematical model able to describe the average course and variability of retinal nerve fiber–bundle trajectories. These researchers determined empirical nerve fiber–bundle trajectories from fundus photographs, as previously done by Garway-Heath et al., ²⁸ and fitted the mathematical model to them. A recent work by Denniss et al. ³¹ corroborates this modeling approach with the study of variability in nerve fiber tracing. This mathematical approach, together with the wealth of empirical measurements it is based on, can be used as a tool to guide attempts to generate the course of the axons in a computerized model. The advantage of a model, in our case a computerized model, is that it can be improved in a stepwise fashion until it is clinically useful or, otherwise, discarded. Additionally, such a model might offer insight into the embryological development of the optic stalk. 

An ideal objective would be for each RGC on the retinal surface to be able to follow the course of its axon on the retinal fiber layer and to locate the entrance point in the optic disc, but clinical devices as well as experimental techniques are still far from this objective. Our approach is to build a computerized map of the optic pathways that can be honed on the basis of pertinent biological research and clinical data from health and disease. In the present work, our purpose is to formulate an algorithm able to establish a relationship between each individual ganglion cell's soma on the retinal surface and its axon's point of entrance into the optic nerve. This is a fundamental step in the building of a realistic model of the primary visual pathways and their corresponding visual fields. The reproduction of the course of axons on the retinal surface—conforming to empirical data—by a single algorithm would mean an improvement of virtual tissue engineering over our previous model, which is intended to assist in the clinical evaluation of visual field damage in glaucoma and other pathologies. 

Methods

We developed a computer model of a section of the anterior visual pathways comprising the retinal ganglion cells and the axons up to the optic nerve head. ³² Here we detail the simulation steps necessary to establish the correct relationship between the virtual soma of the RGC and its entrance point into the optic disc, which is vital in order to determine how the nerve damage would project in the visual field. 

Our simulation process includes several steps: 

1.  
   Distribution of ganglion cells according to user-defined density functions;
2.  
   Distribution of entrance points into the optic disc;
3.  
   Establishing the priority order of the entrance of RGC axons into the optic disc. This will determine the preference when choosing a free slot in the optic disc; and
4.  
   Linking each RGC soma to its chosen entrance point. This step allows us to minimize the set of requirements in the process of finding an entrance point for each axon.

As in our previous model ³² (refined here), in the first step, we created a retinal-cell distribution, or mosaic. The following empirical data were used to generate the mosaic: the total number of RGCs, density of RGCs in each location of the retina, entry site of the axons on the optic nerve head (ONH), density of axons on the ONH, size of the ONH, size and position of the optic cup, and area of the neuroretinal rim. The model requires two parameters to generate the mosaic: position in the retina and density of cells at each spot on the retina. We based our modeling on the assumption that a majority of RGCs in the human retina form part of the geniculostriate pathways. ^33,34 We purposefully ignore the many morphologically distinct RGCs and represent a single idealized RGC that participates in the most basic light perception or contrast sensitivity. 

Ganglion Cell Distribution

Data on the densities of the two most abundant RGCs in human (i.e., the midget and parasol ganglion cells), were taken from Curcio and Allen ³⁴ for the peripheral retina and Sjöstrand et al. ³⁵ for the area centralis. Following Denniss et al., ⁵⁶ we have integrated a distribution function based on Figure 5 in Curcio and Allen, ³⁴ to establish the range density and have applied a linear interpolation to the data between (0,1). Each time the program is run, our model generates a random variable retina from this normalized function, using the maximal value in the fovea, to determine F(x,y). For the present report, we generated a retinal matrix of 1,200,000 RGCs plotted with the same color scale as in Curcio and Allen ³⁴ for easy comparison (Fig. 2). 

Figure 2

 

Density and distribution of virtual retinal ganglion cells. (a) Data of the densities of virtual RGCs in the simulated retina, with data taken from Curcio and Allen, ³⁵ for the peripheral retina, using the same color code, and Sjöstrand et al. ³⁴ for the area centralis. (b) Plot of a sample distribution of the RGC somas on a generated retinal surface, extended to 90° in the nasal meridian. Each point represents the center of the dendritic arbor of the RGC, which is made equivalent to the corresponding receptive field.

Depending of the purpose of the simulation, the number of RGCs represented is varied accordingly. For instance, if the purpose is to reproduce the arrangement of the axons in the optic nerve, 1,200,000 fibers is the preferred option (Fig. 3). If a relationship is to be sought between the number of RGC and visual-field sensitivity, the number of RGCs should not be higher than 600,000 because of the histological structure of the receptive fields. In the simulation, we equate dendritic field with the center of the receptive field, which can be considered the perceptive field. Each receptive field is served by two RGCs, ON and OFF, which function as a coupled unity and can be represented by a single functional ON-OFF RGC. If only one subtype of RGCs is to be explored, for instance to design new visual stimuli, the number of cells can be adapted accordingly (Swanson et al. ³⁶ ). As detailed below, the course of the axons is not negatively affected by the number of RGCs. To generate random cell locations according to the chosen distribution functions, we integrate a mathematical process explained in the Appendix. 

Figure 3

 

The course of axons in the NFL plotted by a single algorithm that relates each soma in the retina to an entrance point in the optic disc. The arching of axons surrounding the area centralis is automatically traced in a process that simulates fasciculation of the growing axons during development. T, temporal; N, nasal.

Distribution of Entrance Points Into the Optic Disc

Once the RGC somas are distributed, we need to simulate the entry point of their axons into the optic disc. As we detail in following sections, RGC density is higher in the temporal hemiretina, axon density in the temporal region of the disc is also correspondingly higher. The entry slots are located in a circle representing the neuroretinal rim with the exclusion of a small inner circle representing the optic cup. The relative sizes of the circles and the position of the optic cup are parameters adjustable by the user. 

For the simulation of the retinal trajectories of the axons, the entry points in the disc can be distributed with different densities (i.e., number of points per area) in different regions of the nerve. This is because initially the disc is symmetrical, with the cup located in the geometrical center of the disc. 

Due to the high number of fibers coming from the retina temporal to the optic nerve, we use a progressive linear function of decreasing density from the temporal boundary of the disc to the nasal boundary. The nerve distribution has been simulated with random entry points fit to the linear function, where the probability in the temporal boundary is 1 decreasing to 1 − k in the temporal boundary (see Appendix for details, Equation A1). This parameter k can be adjusted to suit different retinal typologies. 

Establishing the Priority Order of the Entrance of RGC Axons in the Optic Disc

The entry order of RGCs is critical for determining how and when a precise location in the disc is saturated and a new location needs to be found. In retinas with an area centralis, the competition of the axons to find a free slot is determined by the fact that the maculopapillary bundle has priority access to the disc. Accordingly, the remaining fibers have to avoid covering the maculopapillary bundle and have to deviate, surrounding the extant fibers, in order to find the next most proximal entrance point available in the disc (Fig. 4). Our model assumes that the guiding process determinant in the layering of the axons is the intimate contact that growing axons establish with preceding axons, a process known as fasciculation (see Appendix for details, Equation F2). In order to make fasciculation responsible for the arcuate path of the axons, a sequential order of growth, also known as timing, is a fundamental parameter. 

Figure 4

 

Trajectory of the axons on the retinal surface and corresponding entry point on the disc. (a) The path is arcuated in higher density areas and straight in lower density areas. (b) A close view of the disc showing exact correspondences.

Figure 5

 

Rules for the order of assigning free slots as entrance points in the optic stalk to axons from RGC somas: (a) The narrower angle formed by the line that connects the RGC soma and the entry point in the disc and the line that connects the fovea with the center of the optic stalk is selected. (b) If two angles have the same value, the slot situated at the nearest point to the RGC soma is selected. (c) The decision order is inversely proportional to the distance from the RGC soma to the optic stalk.

Linking the Retinal Ganglion Cells

Once the priority of retina cells has been ordered, they compete to find a free nerve slot and occupy it. Thus, in the minimization process, each ganglion cell will choose a nerve-axon location occupying it (see Appendix for details, Equation A3). The leftover nerve axon locations are removed after the minimization process. A strength of the present model is that the connections between retinal somas and nerve entry points is not affected by the number elements in the retina and nerve because of the use of density functions independents of the number of elements. 

Generating the Course of the Fibers

The nerve location is chosen using a heuristic minimization process, which represents the biochemical signals and the adherence of fasciculation. This heuristic is represented by a weighted function that includes three parameters (Fig. 5) that determine the course of the axons bending toward the optic nerve head (see Appendix for details, Equation A4). 

We assayed Bézier curves (Figs. 3, 4, 6) to connect the soma of the RGC and their corresponding entry points in the optic stalk. These curves are traced from one starting to an end point using control points that curve the path gradually. 

Figure 6

 

Each axon enters the disc at a precise location determined by an angle. (a) A single path is represented by a continuous line. This case shows the singularity for some of the RGCs that, despite having the soma peripheral to the fovea, enter the disc near the center of the ONH, not following the general arrangement shown in Figure 1. (b) Shows the angle formed as the incoming fiber intersects the perimeter of the disc. The angle of the intersections for each tested axon are compared with the range of values given by Jansonius et al. ³⁰

A sample retina, with a limited number of points to allow easy visualization of the individual trajectories, is shown in Figure 3. See Appendix for mathematical and computational details. 

Results

The order of entry of each axon into the optic nerve is a critical step, and has been tentatively resolved by means of an elliptical pattern. A sequence of the wave that generates the nerve is shown in Figure 7. All the geometric parameters of the wave can be configured by the user of the program and the output can be adapted to the desired situations. 

Figure 7

 

Sequential order of addition of fibers during the virtual nerve development. (a–c) Samples from the propagating wave that generates the entry order of axons into the optic stalk.

The course of the axons complies with the demands of a virtual retina, which mimics important events known to take place during retinal development. One notable feature is that its course deviates near the macula and avoids passing above the area centralis. The resulting arcuate course of the axons, a landmark of the retina in humans and other primates, is reproduced by our model using a single and configurable algorithm. The order of entrance of fibers into the optic disc appears in the modeling of the early visual pathways as one of the critical features of the system. 

As a consequence of the correct relationship between soma location, axon path, and entry point into the optic disc, the correspondence between areas of the retina and sector of the optic disc is a straightforward output of the model. An example of a generated retina and nerve showing the relationship between areas of the retina and the corresponding portions of the optic disc is displayed in Figure 8. As commented above, the entrance point for each axon is critical in generating the early optic pathway and has consequences for the course of the axons on the retinal surface. Precise timing regarding the entry order of axons results automatically in ordered arcuate pathways on the retinal surface. 

Figure 8

 

A sample retinal optic disc generated by the model. Left: Retina. Right: Optic disc. In alternate shades of gray, corresponding regions in the retina and disc. 1, maculopapillary bundle; 2 (not marked), supramacular bundle; 3, superotemporal bundle; 4 and 5 (both not marked), upper quadrant.

A singularity affecting the axons coming from RGCs close to the raphe, temporal to the fovea, previously modeled in our 2011 paper, is that they enter the disc centrally rather than peripherally, as would be deduced from Figure 1. This singularity would help explain the early presence of the nasal step in the glaucomatous visual fields. The contribution of the present model is that the singularity appears spontaneously generated by the algorithm that simulates fasciculation and timing and persists when the model is validated with the data of Jansonius. ^29,30  

To validate our model, we compare the course of a sample of axons finding the location and angle of entrance in the perimeter of the disc in our model and compare it with the corresponding fiber in Figure 5B from Jansonius et al., ³⁰ which also provides the location, average entry angle, and 95% range angle. For a comparative evaluation, we have set the nerve in the same location (15°, 2°) as in the work by Jansonius et al., ³⁰ and we have set the same disc size (1.6, 1.8 mm) as in Denniss et al. ⁵⁶  

A noteworthy contribution of our approach is the relation between ganglion cells and the entry point location in the disc. We can estimate the angle formed by a line connecting the entry point and the center of the disc with the line connecting the fovea and the center of the disc. This is the angle referred to by Jansonius. ³⁰ The point where this line intersects the perimeter of the disc is, for all practical purposes, identical to the point where the Bézier curve connects the RGC soma, and the entry point in the disc intersects the perimeter of the disc. The value of the angle at this last point is used for our validation (nerve entry angle in Fig. 6B). To start the validation process, we fixed an initial value of k in the density function of the nerve at k = 0.2, in the nasal half of the disc [−60°, +60°] where the Bézier trajectories are straight, as shown in Figure 6A of Jansonius et al. ³⁰ Also, as in the aforementioned work, we explored different weights for the heuristic values in the upper and lower hemiretinas. The weights of the heuristics have been tested several times because the retinal nerve configuration is very sensitive to changes in the density distribution and in the guided random relation. As differences between an accurate or imprecise configuration could be the result of chance, we execute the process several times with the same weights, getting the average values of Table A1. On the other hand, the model improves the results in single distributions, reaching the percentages in Table A2. For each value of the heuristics in Table A1, the random process was iterated 10 times. We measured the percentage of matches with respect to the ranges, and mean absolute error with respect to the mean angle value provided by Jansonius et al. ³⁰ The results show a high level of coincidence, 85% on average (see Appendix, Table A1) and up to 90% in certain cases (see Appendix, Table A2). Figure 9 presents a graph of the relationship of parameter k and the saturation angle. The curve has an area of interest in which the slope changes. This area is the border between the overpopulated temporal central retina and the less populated nasal region. It is useful to relate the manipulation of k values to the possibility of reproducing from a normal retina up to a situs inversus. 

Figure 9

 

Graph of the relationship of k factor and percentage of matching with the ranges provided by Jansonius et al. ³¹ The change in density between temporal and nasal halves of the nerve is directly proportional to k. This means that, for our simulation, with a slight variation in density, the percentage of matching reaches almost 90%.

Discussion

A computational approach toward understanding how neural systems may function involves collecting the most relevant anatomical, physiological, and pathological information available at the moment concerning the structure of the nervous system, as well as using feedback data to direct further experimental research. Our computer model is in line with current retinal optic nerve data that reproduces and helps explain some critical features of the progressive damage that develops in the visual fields of glaucoma patients. Glaucoma is a pathology that can be ideally explored with a computer-based model. ³⁷ The well-founded suspicion that the damage site is confined to the intraocular portion of the optic nerve head allows the testing of different patterns of damage impinging upon the optic nerve. Damage to optic nerve fibers can be translated both to the retinal and visual-field models, including mechanisms for simulating intra- as well as interindividual variability. ³² Here, we discuss aspects related to an important debated feature of early visual pathways, the course of axons on the retinal surface and their entrance point into the optic nerve, and leave the structure/function relationship between the model and the visual fields to be dealt with in a later work. This first implies questioning whether the relevant features of the retina are structured enough for a simulation to be meaningful; and second, an analysis of the correspondence between the processes simulated and their biological counterpart. 

Relevance of Retinotopy

High levels of order and coherence in the pattern organization of the central nervous system (CNS) is evident from the rearrangement of blastoderm cells during gastrulation. ³⁸ In the CNS, topographic maps are a basic architectural method of organizing neural circuits. Each map permits the comparison and combination of the information carried by various specialized neuronal subpopulations. ³⁹ Retinal map architecture is one of the best understood. Visual field maps preserve the spatial structure of the scene itself: nearby scene points are represented in the responses of nearby neurons. ⁴⁰ This orderly arrangement is maintained throughout the different strata in the visual pathways, including the successive visual cortical regions, and is encompassed in the concept of retinotopy. 

Retinal nerve cell classes are distributed as regular retinal arrays, commonly referred to as retinal mosaics. We simulate a mosaic of RGCs according to empirical data, in which the RGC soma is the central point of an imaginary dendritic arbor. It has been shown by live imaging of labeled retinal cells that early during retinal development the somas of the retinal cells move and that their movements contribute to the improvement in mosaic patterning. Dendritic interactions help to establish a domain within which the cell soma will occupy a central location, ¹⁶ a model carefully recorded by Dacey. ⁴¹  

Regarding the axons of the RGCs, it has been argued that retinotopic fiber organization is mediated by a number of processes that affect the axonal outgrowth during embryonic development. ⁴² In nonmammal vertebrates, this organization has been deemed responsible for the arrival of fibers in orderly arrays to the central nuclei during normal development. However, in mammals, and especially among primates, the retro-ocular optic nerve undergoes deviations, sometimes bizarre, from the ideal retinotopic ordering. In primates, the retro-ocular portion of the optic nerve is quite disorderly, especially near the chiasm. ⁴³ Some experimental studies seeking to ascertain the path of axons have highlighted the lack of evident order, rather than the presence of it, as a substantial feature of the optic pathways. ^21,44 Nevertheless, it is important to distinguish conceptually between the notion of retinotopy and the presence or lack of apparent order in a section of the optic pathways. In our model the pertinent fact is that both timing and spatial positioning of fibers as they grow along the pathway are critically preserved in the intraocular optic pathways (i.e., NFL and prelaminar optic nerve), and allow RGC axons all to project correctly toward specific points, ⁴⁵ although such retinotopic fiber order becomes degraded in the retro-ocular and prechiasmatic portion of the pathway. During developmental axon pathfinding in the retina, the deeper retinal lamina exerts a repellant effect on axons to prevent their escaping from the RGC axon layer. ⁴⁶ Multiple axon guidance mechanisms concentrically organized around the optic nerve have been identified, and those include both growth-promoting and growth-inhibiting guidance molecules. ⁴⁵ Additionally, it is assumed that ample dendritic arbors in the lateral geniculate nucleus (LGN) and visual cortex (V1) compensate for the intrinsic disorder. The terminal field gradually emerges from axonal arbors that are initially diffuse, particularly in mammals. ⁴⁷ In this process, the formation or exclusion of stable mature synapses is the result of an orderly series of hierarchical processes. ⁴⁸ Other guidance mechanisms are attractive-repulsive chemogradients. Additionally, the participation of chemospecific receptors in the mapping of optic pathways is mediated by gradients and counter-gradients of EphA receptor expression for the membrane-linked molecules, Ephrin-A2 and Ephrin-A5. ⁴⁹ Also in play are attractive interactions mediating the effect of some of those gradients. ⁵⁰ Cell death—eliminating useless contacts—plays a sculptor's role in this process, ⁵¹ although this aspect has also been minimized by some researchers. ⁵² There is enough evidence to state that complementary gradients across the retina participate in retinotopic map formation by providing retinal growth cones with the positional information to locate each cell in the general arrangement. 

Consequently, our simulation is based largely on assumptions of retinotopy. How this retinotopic order is being investigated has been recently reviewed by Reese. ¹⁶ Retinotopy is a functionally adequate concept, although the morphologic elements of the pathways can show different degrees of apparent disorganization, depending on the area considered. We have seen that the RFL is ordered, although it includes some fiber scattering. ²⁰ In relation to this order: 

1.  
   Foveal fibers occupy a large proportion of the temporal aspect of the optic nerve head and nerve, whereas
2.  
   Fibers from areas temporal to the fovea appear to be displaced to more superior and inferior regions, maintaining a flat profile (without invading other tissue or stepping over (covering) the area centralis.
3.  
   A general topographic representation within the optic nerve head agrees with the retinal topology, so that both peripheral and central fibers are correspondingly projected radially into the peripheral and central portions of the nerve (Fig. 10).
4.  
   This basic ordering conforms to a system of signals that contribute to axonal pathfinding during retinogenesis. An attractive/repulsive system would tentatively include an attraction signal with a gradient of concentration from the optic stalk that would direct axons to the future nerve in the center of the optic cup. Figure 11 shows two hypothetical attractive and repulsive foci that would give rise to an arciform course of axons in the simulated retina. A neuroattractant focus is situated in the optic stalk, and a neurorepellant one in the fovea, following the hypothesis of Airaksinen, ⁵³ both with a concentration gradient thinning toward the periphery. Our model shows that arcuate paths are autogenerated by reproducing two well-established developmental processes, namely a neuroattractant at the optic stalk and axonal fasciculation (i.e., adhesion of one axon to the immediately preceding one). The result of both processes combined is an arcuate pattern that conforms to the empirical data from the human retina.
5.  
   A crucial step is the timing of the axons reaching the optic stalk. The first axons to reach the optic stalk are those in the horizontal midline of the maculopapillary bundle and closest to the disc. Figure 7 shows samples from the propagation wave determining the entry order.

Figure 10

 

The direction of the growing axons of the RGCs toward the optic stalk depends on the diffusion of biochemical messengers. In the absence of an area centralis (a, b), the track toward the optic stalk would be direct, and would project in a plane as a straight line, as in the case of the much-studied zebrafish. In case of an area centralis, the axons of mammals and some birds, avoid entering the macular area and follow a curved trajectory bordering the maculopapillary bundle (gray elongated dot in [c, d]) before reaching the optic stalk. Arching would attenuate as axons are located away from the maculopapillary bundle ([c, d]; based on a drawing by Oster et al. ⁴⁶ ).

Figure 11

 

Growth of axons in simple and complex retinas. (a) In a retina without area centralis, axons are attracted toward the optic nerve by the increasing concentration gradient of a biochemical signal. Fasciculation assures minimal or no axon wandering. (b) A hypothetical factor of rejection spreading outward from the fovea, acting in conjunction with the attraction from the optic stalk, could theoretically justify the arcuate trajectories. However, our model resolves the problem of the circular course of axons with a simpler approach, exclusively by means of fasciculation and order of entry into the optic stalk (see text).

To sum up this part, we question whether the relevant features of the retina are structured enough for a simulation to be meaningful. We conclude that it is consistent with current empirical data that the retinal distribution of RGC somas and axons are orderly enough to allow the use of maps generated by algorithms that simulate aspects of the pathfinding mechanisms. Additionally, the results can be completed with a variable degree of randomness to account for interindividual variations or tailored to the clinical images of real patients. Computer simulation of retinal sensitivity already has a place in the interpretation of visual function. ^36,54  

Simulation and Biological Processes

We are reproducing here a chronotopic model of axonal addition seen in the retinal genesis of mammals and nonmammals alike. In mammals, retinal neurogenesis proceeds in overlapping centroperipheral waves. ^20,55 Our model mimics the centroperipheral sequence of RGCs neurogenesis. In it, axons from somas situated at similar distances from the optic stalk would be forced to compete for an entry slot. The order of entry into the optic nerve can be resolved by means of an elliptical pattern of variable flatness (Fig. 7). 

To validate the model, we have used data from Jansonius et al. ^29,30 In fact, one of the singularities in Jansonius' work seems to relate to the parameter k. The number of available entry points in the optic disc is a parameter that has great impact on the course of the axons in the retina when retinal somas and entry points are connected by a plotted axon. Small variations in the density of free slots available at the disc are needed if a predetermined course of the axons is sought. In this work, we present the density values needed to reproduce the average course of axons reported by Jansonius ³⁰ based on fundus photographs (Tables A1, A2). 

In our model, when the axons do not have to compete for the empty slots in the nerve, when their most preferred entry points are available, then the axons stop curving and follow a rectilinear path. In the model developed here, the nerve density is either homogeneous or slightly less populated in the nasal half compared to the temporal half (Fig. 4). The fuzzy parameter k helps create a progressive gradient of density in the nerve between the temporal and nasal areas. 

Our model reproduces the singularity pointed out by Jansonius et al., ³⁰ in the sense that fibers originating from the far temporal periphery have to enter the nasal half of the nerve mixed with fibers from the nasal hemiretina (if we divide the retina in two halves bisecting the nerve). 

In our model, the most important singularity is the central location of axons coming from the temporal periphery, a fact that, although reported in our previous study, ³² is generated spontaneously by the model. In relation to Figure 1, our model follows the general array of the axons on the z axis depicted in the figure with one important exception, due to the abundance of RGC in the center of the retina and the deflection of the axons in an arcuate pattern encircling the area centralis. Axons coming from the periphery temporal to the macula may enter the disc in the central regions close to the optic cup, and not in the periphery of the disc, as the figure would indicate. This singularity that was simulated in our first model ³² to account for the early defects of the glaucomatous visual field loss known as nasal steps is spontaneously reproduced by the ordered fasciculation of the retina modeled here. If we assume that the axons coming from the areas adjacent to the temporal raphe assume a superficial position on the nerve fiber layer, then the singularity would affect only the peripheral fibers (on the retina) entering centrally (in the disc), while the common arrangement of superficial fibers (on the retinal fiber layer) entering centrally (in the disc) and deep (in the retinal fiber layer) fibers entering peripherally (in the disc) would be preserved. 

Concerning the order of entrance of nerve fibers at the optic stalk, a remarkable fact in the simulation to represent the optic nerve is that if axons enter sequentially, starting from the closest location of the temporal horizontal midline and follow the principle of fasciculation (i.e., each axon “adheres” to the preceding one), then all the fibers distribute automatically with a single algorithm. Furthermore, when both soma and entrance point are then linked with the most suitable arc, the resulting path of fibers on the surface of the retina is strikingly reproduced. This process may partially simulate the correct entry order of fibers into the optic stalk. Depending on the relative elliptical shape, what is indicated is the order of entrance of the axons. The more flattened the ellipse, the more privilege of entrance the axons coming attached to the previous axons have. If the pattern is circular, all the axons in the circumference arrive at the same time, and a free slot can be disputed by more than one axon. Together with timing, the distance to the optic stalk is another determining factor. Thus, without the need of additional parameters, depending of the size of the disc, when the first fibers have occupied the central and peripheral portions of the disc near the horizontal meridian, fibers temporal to the macula will tend to find an entrance point again near the center of the disc, as soon as the periphery of the nerve is saturated in the lower row. Fibers coming from RGC somas temporally situated at a further distance from the macula will compete with fibers coming from the supero- and inferotemporal regions to find a free slot in the optic stalk. Fitzgibbon and Taylor ¹⁵ defended the contention that axons of the arcuate bundles show fiber scatter across the nerve fiber layer, a situation that is not so evident in other retinal regions. This is compatible with our approach, mentioned before, that the bundle coming from the temporal raphe has to mingle with fibers from the radial periphery and, consequently, some fibers may enter the central optic disc. This would explain the early damage of those fibers in glaucoma, giving rise to the nasal-step type of field defect, when the optic cup increases in size due to optic nerve atrophy. 

A compelling issue raised by our modeling process is the relevance of the position of the entry points of the axons in the disc. Previous models emphasize the retinal course of the axons due to the clinical relevance of the relationship between visual field defects, which point to the retinal location of deranged RGC somas, and defects in the peripapillary nerve fiber layer. One advanced model of this latter type is that of Denniss et al. ³ Although the objective of Denniss et al. ³ is not the axonal distribution of the axons in the optic disc, the algorithms employed in our simulation, developed in parallel, appear in retrospect to be similar to those used by them. The first four stages are very similar in the use of angles and distances to the main retinal features (fovea, raphe ONH). As our model avoids the coarse use of sectors in the disc, it renders as an output a finer distribution of cells and axons, including the position of every axon in the ONH, as well as the plotted course of the axons on the retina that match current clinical data. The curves that link the points representing RGC somas and entrance points in the optic disc very closely reproduce the paths of the retinal nerve fiber bundles taken from fundus photographs of real subjects. It would be possible to fit exactly the same type of curve to the clinical data of each specific patient. 

Conclusions

A computational approach to reproduce in great detail the optic pathways in the retina and optic disc is feasible due to the accumulated knowledge about the distribution of RGC density in the retina and the predictable course of the axons in the RFL, based on the concept of retinotopy and the biological processes governing it, notably fasciculation. The importance of timing in this context is emphasized in our model with the autogeneration of the course of axons in the NFL based on the competition for an entry point in the optic stalk. The notable similarity of samples generated by our model with sample photographs of the normal fundus make this simulation a configurable tool to explore structure/function relationships in normal and pathological situations, such as glaucoma. The Table succinctly summarizes the backbone of the model, including assumptions (model input), findings (model output), agreement/disagreement with the literature, and new predictions that might feed new studies. 

Table

 

Modeling the Prelaminar Optic Pathways

1. Assumptions, model input

– Density functions of retina and ONH in order to locate the RGCs and entry points in the optic disc.

– Dendritic arbor and receptive fields coincide.

– Shape and size of retinal surface and ONH. Optic nerve head is superiorly displaced in relation to fovea.

– Weights are given to the heuristics that guide the preference entry order in the nerve.

– The inputs of model have been adjusted using data from the literature (Jansonius et al.).29,30

2. Findings, including sensitivity to assumptions, model output

– From this configuration, we simulate a distribution and competition process that relates each RGC to an entry into a nerve, giving nerve obtaining a relation that offers the following findings.

– The order in which RGCs enter the nerve. Retinal ganglion cells closer to the optic disc and nearer the horizontal meridian have preference.

– The entry point where each RGC projects in the nerve. Orderly disposition avoids crisscrossing of the axons.

– The previous two steps determine the course of the axons that merge at both points. The course is plotted using Bézier curves.

– Retinal ganglion cells from the area centralis saturate the temporal section of the disc. Retinal ganglion cells temporal to the area centralis curve their paths surrounding the fovea. Retinal ganglion cells from the nasal side follow a straight course.

– Nasal regions of the disc receive mixed axons from the far periphery of the temporal hemiretina along with axons from the nasal half, reproducing the singularity referred to by Jansonius et al.30

– Sensitivity to the weighting process can be deduced from Tables A1 and A2 and the graph of Figure 9. As ONH is superiorly displaced in relation to the fovea, this creates an asymmetry than warrants a tailored approach in the validation process.

3. Agreement with the literature

– Our model shares several elements with Denniss et al.,56 in the relation to the course of the axons in the retina, refining the precision and introducing the exact entry point for each axon in the optic disc.

– When compared with the empirical analysis of Jansonius et al.,29,30 our model reaches a coincidence level of 86% on average and up to 90% in certain cases.

4. New predictions

– The finding that RGCs close to the peripheral temporal raphe enter the disc centrally is reproduced spontaneously by the model. This singularity may explain early nasal step defects in glaucoma.

– Relating each soma in the retina to an entry point in the nerve simulates the migration of the axons into the optic stalk during development. Following a strict consecutive sequence reproduces the biological phenomenon of fasciculation, in which timing is considered a decisive factor. This may imply that only orderly growth (timing) and fasciculation may be sufficient to give rise to the ordered retina ONH during embryological development.

Acknowledgments

Supported by a CEI-Biotic grant as part of the program “Compromise with Research and Development, 2012,” of the University of Granada, Spain. The authors do not have any commercial interest in the software presented in this article. A version of the software, and the data on which it is based, is already freely available to all researchers consulting our website, at http://ic.ugr.es/vifilab/index.php, with the only condition that its use will be properly acknowledged. 

Disclosure: F.J. Carreras, None; J. Medina, None; M. Ruiz-Lozano, None; I. Carreras, None; J.L. Castro, None 

Appendix

Methods

Ganglion Cell Distribution.

Given a function F(x,y), which represents the distribution density, we can determine random N values based on this process: 

1.  
   Generate a normal distribution for random value x in the function range.
2.  
   Generate a normal distribution for random value y in the domain range.
3.  
   If y < F(x,y) return y;
4.  
   Otherwise go to 1.

This process is executed N times to determine the RGC soma locations (Fig. 2). 

Distribution of Entrance Points Into the Optic Disc.

To simulate the entry point location, we distribute N points in the disc according to a linear progressive function to locate more entry points in nasal areas in the nerve. The linear function evaluates the X location of the entry points in order to evaluate their probability. At nasal entry points the functions are closer to 1 and at temporal points decrease to 1 − k has minimal probability in limit. Given a k value (0,1) and a point n(x,y) in nerve the normalized function distribution is  

 

where  

• –  
  n(x, y) is the nerve-axon location;
• –  
  k is the minimal probability in the temporal margin of the disc;
• –  
  nerve.center(x,y) is the nerve-center location; and
• –  
  nerve.size(x,y) is the nerve size.

Table A1

 

Best Result for the Median of 10 Execution Times

Entry Order Upper Hemiretina			Entry Order Lower Hemiretina			Entry Preference Optic Stalk			Range, %	Median Error, %
wα+	wF+	wN+	wα−	wF−	wN−	wα	wD	wN
0.5	0.1	0.3	0.3	0.2	0.8	0.7	0.3	0.3	86	4.00
0.5	0.1	0.3	0.3	0.3	0.9	0.8	0.3	0.4	86	5.00
0.6	0.2	0.4	0.4	0.3	0.9	0.7	0.3	0.3	86	4.53

Table A2

 

Best Result for a Single Execution

Entry Order Upper Hemiretina			Entry Order Lower Hemiretina			Entry Preference Optic Stalk			Range, %	Median Error, %
wα+	wF+	wN+	wα−	wF−	wN−	wα	wD	wN
0.5	0.1	0.3	0.3	0.3	0.8	0.8	0.3	0.4	94	3.45
0.5	0.1	0.3	0.3	0.1	0.8	0.7	0.3	0.4	93	3.90
0.5	0.1	0.3	0.3	0.2	0.9	0.7	0.3	0.4	93	3.29
0.7	0.2	0.4	0.4	0.2	0.9	0.8	0.3	0.4	93	4.13
0.7	0.2	0.5	0.4	0.2	0.8	0.8	0.3	0.4	93	4.56
0.5	0.1	0.3	0.3	0.3	0.9	0.8	0.3	0.4	92	3.22
0.5	0.1	0.3	0.3	0.1	0.9	0.7	0.3	0.4	92	3.22

Establishing the Priority Order of the Entrance of RGC Axons in the Optic Disc.

To represent the priority of the cells, we define a weight function that represents the value of priority:  

 

where:  

• –  
  r(x,y) is the retina-cell location;
• –  
  fove is the fovea-center location;
• –  
  nerve is the disc-center location;
• –  
  |p_i − p_j | is the Euclidean distance of the two points p_i , p_j ;
• –  
  p_i × p_j is the line between p_i , p_j points;
• –  
  α(r_i , r_j ) is the angle between r_i , r_j lines;
• –  
  w_F is a weight [0,1], which increases with fovea proximity;
• –  
  w_N is a weight [0,1], which increases with nerve proximity; and
• –  
  w_α is the weight [0,1], of the line minimization between the nerve and the retina locations.

In result, we discuss the values for comparative works. 

Linking the Retinal Ganglion Cells.

The process for linking the RGC in the nerve is iterative. We apply the heuristic function for each RGC, evaluating all free entry points into the nerve. The minimum value of H represents the best-matching nerve location. To establish a stochastic relation, we link the RGC to a good lower candidate using a random approximation. The random value is generated using a normal probability distribution, which represents the nature of the nonstochastic method. 

 

 

 

Where: 

• –  
  r(x, y) is the retina-cell location;
• –  
  n(x, y) is the nerve-axon location;
• –  
  p_i × p_j is the line between p_i , p_j points;
• –  
  α(r_i , r_j ) is the angle between r_i , r_j lines;
• –  
  |p_i − p_j | is the Euclidean distance of the two points p_i , p_j ;
• –  
  w_α is the weight of the line minimization between the nerve and the retina locations;
• –  
  w_D is the weight of the distance between the nerve and the retina locations; and
• –  
  w_N is the weight of the distance between the nerve-axon location and the optic-nerve center.

In this way, this function contains three factors. The first one, w _α , represents how the axons find the path closest to the fovea and the nerve location. When two angles coincide, the second parameter, w _D , chooses the entry point closest to the RGC soma. The third parameter, w _N , is the distance from the RGC soma to the center of the disc, in order to establish an entry order, or timing. 

Finally, the steps in the iterative process are: 

• 1.  
  Order RGC according to O(r(x,y)); and
• 2.  
  For each retina cell r(x,y) in RGC
   
  •    2.1.  
    bestCandidate = null
  •    2.2.  
    For each location n(x,y) in Nerve where n(x,y).free
  •  
     If (H[r(x,y),n(x,y)] + random < bestCandidate)
     
    •         2.2.1.  
      bestCandidate = n(x,y)
  •    2.3.  
    bestCandidate.free = false

Generating the Course of the Fibers.

In our case, we have determined two control points, which generate a quadratic Bézier curve. The Bézier points are defined by the ganglion cell location r(x, y) and the nerve entry point n(x, y): 

• –  
  Init point. It is r(x, y), the ganglion cell location in the retina that indicates the init of the axon path.
• –  
  End point. This is n(x, y), the nerve entry point, where the axon enters the optic disk.
• –  
  First control point. This is the center of the optic disk, c ₁ = nerve, representing the convergence of the path toward the nerve.
• –  
  Second control point. This represents the angle saturation in the nerve, that is, the curve of the path to find the entry point. This control point is estimated by the angle between the ganglion cell and the entry point and the distance between the ganglion cell and the entry point.

 

 

 

Where 

• –  
  |p_i − p_j | is the Euclidean distance of the two points p_i , p_j
• –  
  p_i × p_j is the line between p_i , p_j points.
• –  
  α(r_i , r_j ) is the angle between r_i , r_j lines.

Thus, the second control point represents the curve due to the accumulative density in the nerve. When there is no stress between the ganglion cell and the nerve entry in nasal areas, the entry point is located on the straight line of r(x, y) and n(x, y). On the other hand, when the entry point is not in this straight line, it is moved to another location due to the competition of temporal areas, widening the angle of the Bézier curve. In Figure 6, we plot the Bezier curve, the control points in circular dots, and the nerve entry point, RGC, and entry angle in square dots.