Retinal self-organization: a model of RGC and SAC mosaic formation

Individual retinal cell types exhibit semi-regular spatial patterns called retinal mosaics. These mosaics enable uniform sampling of visual information and are formed to varying degrees across cell types. Retinal ganglion cells (RGC) and amacrine cells (including starburst amacrine cells (SAC)) are notably known to exhibit such layouts. Mechanisms responsible for the formation of such organised structures and their requirements are still not well understood. Mosaic formation follows three main principles: (1) homotypic cells prevent nearby cells from adopting the same type, (2) cell tangential migration, with homotypic cell repulsion, (3) cell death (with RGCs exhibiting high rates of apoptosis). Here, we use BioDynaMo, an agent-based simulation framework, to build a detailed and mechanistic model of mosaic formation. In particular, we investigate the implications of the three theories for RGC’s mosaic formation. We report that the cell migration mechanism yields the most regular mosaics and that cell death can create regular mosaics only if the death rate is kept below 30%, after which cell death has a negative impact on mosaic regularity. In addition, and in accordance with recent studies, we propose here that low density RGC type mosaics exhibit on average low regularities, and thus we question the relevance of regular spacing as a criterion for a group of RGCs to form a RGC type. We also investigate SAC mosaics formation and possible interactions between the ganglion cell layer (GCL) and inner nuclear layer (INL) populations. Investigations are conducted both experimentally and by applying our simulation model to the SAC population. We report that homotypic interactions between the GCL and INL populations during mosaics creation are required to reproduce the observed SAC mosaics’ characteristics. This suggests that the GCL and INL populations of SACs might not be independent during retinal development. Author Summary Retinal function depends on cells self-organisation during early development. Understanding the mechanisms underlying this self-organisation could improve not only our comprehension of the retina and its development but also of the cortex. Ultimately, this could lead to novel therapeutic approaches for developmental diseases. Computational models can be of precious help to study this process of self-organisation, given that they are biologically plausible. In this sense, it is important that implemented developmental mechanisms follow the principle of locally available information, without any global knowledge or external supervisor. Here, we follow this principle to investigate mosaic formation during retinal development. In this work, we demonstrate that tangential migration is the only mechanism able to form regular mosaics and that the GCL/INL SAC populations might not be independent during their mosaic formation. More, we question the relevance of regular spacing for RGC types classification.


Introduction 44
The mammalian retina is composed of six main types of neuronal cells (cones, rods, horizontal, combination. 157 or in combination with the CF mechanism. As shown in Figure 1C, the average RI value increases 159 from random (around 1.8) to 3.31 (±0.33) at the end of CD. This death rate amounts to around 65% 160 when it reaches a steady state at the end of the simulation. These death rate dynamics are very 161 similar to rates observed in-vitro (see Figure 2A). Moreover, and unlike for the case of the CF 162 mechanism, CD is able to generate mosaics of medium regularities (RI > 3). (blue) and in simulations (orange). In-vitro population at day 1 is an estimation based on a final CD 166 of 65%. Error bars represent standard deviation. B: RI score depending on final CD rate in a 167 simulation implementing only the CD mechanism, for selected RGC populations of high density 168 (blue curve, initial density = 571 cells/mm 2) and low density (orange curve, initial density = 114 169 cells/mm 2 ). The vertical dotted line indicate a death rate of 65%. Error bars represent standard 170 deviation. 171 Interestingly, and as shown by Figure 2B, the death rate measured in-vitro and selected for our 172 simulations (grey vertical dashed line) is not the one generating the highest regularity. Indeed, the 173 highest scores of RI are achieved for death rates between 5 and 30% of the RGC population, 174 regardless of the initial density of the considered population. After 30% of cell death, RI decreases 175 until it reaches a random distribution once around 90% of cell death is achieved. This is observed 176 under both high and low-density conditions. This is observed in simulations with CD alone or in 177 differences between populations of high and low initial densities. As shown in Figure 1B, the high-179 density population is able to generate more regular mosaics than the population of low density, both 180 for their maximum value (RI > 9 and RI > 5, respectively, when cell death is below 30%) and at 181 65% of cell death (RI = 4.49 ±0.36 and RI = 3.61 ±0.59, respectively). Therefore, a positive 182 correlation between cell density and the final regularity is observed when the death rate is set to 183 65%. Mosaics of low density exhibit low RI values, while those of cell density higher than 65 184 cells/mm 2 (vertical dashed line of Figure 1D) exhibit a higher average RI score of 3.35 (horizontal 185 dashed line of Figure 1D). 186 While no differences are observed in RI scores between simulation of the CD mechanism and a 187 combination of the CF and CD mechanisms if all mosaics are considered (3.31 ±0.33 and 3.48 188 ±0.44 respectively, p = 0.76), a positive impact on dense mosaics' regularity (for cell densities 189 higher than 125 cells/mm 2 ) is to be noted. Thereby, RI values in the case of CF and CD combination 190 plateaus around 4.1 instead of 3.6 if CD is the only implemented mechanism. As the CD and CM mechanisms require cells to be differentiated, CF is simulated beforehand. A 202 first RI increase corresponding to the effect of CF is observed (see Figure 1E). After the CD and 203 CM mechanisms are triggered (first dashed line), they give rise to a significant second increase, 204 until the RI value stagnates toward the end of CD (simulation day 4 to 5.5 depending on the cell 205 type). Finally, a third RI increase is observed after CD is over (second dashed line) due to the CM 206 mechanism, leading to an average RI score of 4.01 (±0.75) at the end of the simulation. Unlike any 207 other mechanism alone, and thanks to tangential migration, this simulation condition is able to 208 generate highly regular mosaics (RI > 5). Moreover, a strong correlation appears between cell 209 density and RI values (linear correlation magnitude of 0.87) as shown by Figure 1F. Thereby, only 210 RGC types exhibiting a cell density higher than 125 cell/mm 2 are able to generate mosaics with a RI 211 value higher than 5. Thus, as illustrated by the blue and orange lines in Figure 1E  When all three mechanisms are implemented, surviving cells migrate tangentially with an average 215 distance of 8.72 μm (±0.11, n = 8 simulations), which is in accordance with in-vivo measurements 216 reporting migration distance below 30µm (22). Important disparities in migration distance between 217 cells are to be noted, as shown in Figure 4A, with an average migration distance standard deviation 218 of 9.44 (±0.18). No correlation between final RI and migration distance can be seen. Likewise, no 219 correlation appears between final density and migration distance if the whole population is 220 considered. However, if only populations with a final density higher than 100 cells/mm 2 are 221 considered, a strong correlation can be observed (correlation coefficient r = 0.92, see Figure 4B red 222 line). Hence, the denser the cell type the larger the distance cells migrate. 223 Migration distance distribution. B: Relation between cell type density and migrating distance. The 226 red line represents the correlation between migration distance and cell type density for densities 227 higher than 100 cells/mm 2 (correlation coefficient r=0.92). 228

SAC mosaic development 229
The SAC population is divided between two different cellular layers, the GCL and the INL, forming 230 two separate populations (see Figure 5 for an illustration). Our In-vitro results exhibit no significant 231 differences in the GCL and INL populations densities from P4 to P10 (p=0.27 and p=0.32 232 respectively, see Figure 6A). These two SAC populations exhibit regular pattern organisation, and 233 no significant difference over time of their RI is measured from P4 to P10, as shown by Figure 6B.  the simulation. This is observed in both developmental conditions, using either one common or two 256 separate (one for the GCL population, one for the INL population) chemical substances for mosaic 257 formation (RI of 3.57±0.12 and 4.11±0.12 respectively when one substance is used, RI of 3.36±0.07 258 and 4.37±0.15 respectively when two substances are used, n=8 for each group, p<0.0001). This 259 mosaic regularity disparity is in accordance with observations in mouse (see Figure 6B) where the 260 INL population has been reported to be more regular than in the GCL population. In our 261 simulations, this disparity can be explained by the cell density difference between these two layers. 262 Indeed, and as previously demonstrated in our simulations, the denser a cell population is, the more 263 regular its mosaic can be. Hence, our model provides a mechanistic explanation for this observed 264 difference in RIs between the two SAC populations. what has been measured in-vitro (0.74±0.09, n=5, see Figure 6). This indicates that the GCL and 269 INL populations' mosaics tend not to overlap, and so are not fully independent of each other. 270 However, if two distinct developmental cues are used, the exclusion factor is lower, at 0.31 (±0.1, 271 n=8), denoting independent mosaics that tend to overlap. In this second condition, the measured 272 exclusion factor is significantly lower than the one observed in mouse (p<0.0001 using a T-test for 273 two independent samples, n=8 and 5 respectively). Thus, only the first condition is able to Masland (2015) speculated that known RGC types represent only about 60% of the total RGC 284 population, corresponding to around 1740 cells/mm 2 from the total 3000 cells/mm 2 observed in the 285 mouse retina. In addition, it is important to note that from these known RGC populations, only 286 12.4% are On type. As On, Off and On-Off are equally numerous (30% to 35% each), a great 287 number of On cells still needs to be discovered in order to reach the theoretical percentage of On 288 RGC in the total RGC population. Thus, we can hypothesize that either: 1. Several high density On 289 types have not yet been discovered. 2. There are more On types than Off or On-Off types. 290 The first hypothesis appears unlikely as RGCs are widely studied, especially with the emergence of 291 large-scale and high density MEA recordings, but also using morphological and molecular 292 characterisations. Thus, it is unlikely that the existence of dense On RGC types (representing the 293 majority of the On population, and so being the most common On type) has not been captured by at 294 least one of these techniques. The second hypothesis appears to be supported by experimental 295 evidences because mice, similarly to other nocturnal animals, have rod-dominated vision. Indeed, 296 rods are known to project their dendrites and to establish synaptic connections only to On bipolar 297 cells, that in turn establish synaptic connections to On RGCs. In order to extract as many features as 298 possible from a visual scene using mainly rod vision, a great diversity of specialized RGCs can be 299 justified. The hypothesis of a great diversity of low density On types is also in agreement with 300 Masland et al., (2015), who speculate that around 30 low density RGC types exist and are yet to be 301 discovered. Baden et al. (2016) also estimate the total number of RGC types to be over 40, 302 supporting the hypothesis of numerous low density RGC types, including On types. As it is still 303 possible that On types of mid density has not been discovered, we chose to allow the possibility for 304 this hypothesis in our simulations, in addition to adding multiple low density On RGCs. 305 One major basis of our simulations is the presence of chemical cues supporting cellular self-306 organization mechanisms. Evidences of such chemical cues have been previously reported 307 (20,29,30). 308

RGC mosaic formation 309
CF implication on RGC mosaics' regularity is particularly difficult to study in-vitro or in-vivo as 310 RGC progenitor cells do not express RGC type-specific markers cells will differentiate into. Despite 311 experimental studies on RGC progenitors, no evidence has been found for RGC type-specific 312 progenitors (31). Hence, RGC types are probably not pre-determined early on and so are likely to 313 depend on extrinsic factors, such as the presence of chemical cues (7). Thereby, it allows for the 314 contribution of a mechanism such as CF for RGC type differentiation, and its potential implication 315 in mosaic formation. One major conclusion from our simulations is that highly regular mosaics 316 (RI > 2.5) cannot be explained only through the CF mechanism. Likewise, the CF mechanism does 317 not significantly increase the power of other mechanisms (CD and/or CM) neither. This suggests 318 that RGC types are unlikely to be defined by cell body mosaics. They may instead be dictated by 319 intrinsic factors (that still remain to be discovered), functional determination (dictated by the input 320 from other cells), or a combination of intrinsic factors interacting with extrinsic factors. 321 In our simulations, the CD mechanism (alone, or in combination with the CF mechanism) is able to 322 create regular mosaics (RI > 3.5) with a death rate of 65%. As this mechanism is based on a locally 323 diffused chemical substance, homotypic cellular spacing (and therefore cell type initial density) has 324 an important impact on the CD mechanism. For this reason, only populations with a high initial cell 325 density exhibit regular mosaics. Importantly, our CD implementation is able to match measured 326 RGC death dynamics during development, thus strengthening its plausibility. It should be pointed 327 out that CD serves additional purposes in the retinal maturation process and is not only geared 328 towards mosaic creation. Indeed, some cell types which do exhibit mosaic regularity do not undergo 329 any significant levels of CD (such as horizontal cells or photoreceptors). In addition, as 330 demonstrated here, the maximum positive impact of CD upon RI is reached at death rate lower than 331 30%, which is below the 60-70% death rate observed in mouse. This implies that even if CD can be 332 involved in mosaic formation at early stages, cell death at levels above 30% is likely to be driven by 333 other mechanisms and for other purposes than mosaic regularity. For instance, CD could be 334 implicated in refining retinal functional connectivity and activity. CD could also have evolutionary 335 advantage with regard to generating an optimised neural architecture (32). 336 Finally, CM is the only mechanism able to explain the formation of highly regular mosaics (RI > 5). 337 As for the CD mechanism, the efficacy of the CM mechanism is dependent on cell density as it is 338 based on local interactions. The shorter homotypic cellular distances are, the more they can sense 339 and repulse each other. Thereby, a strong correlation emerges between RGC type populations 340 densities and the regularity of their mosaics. Therefore, we propose here that low density RGC type 341 mosaics exhibit on average significantly lower regularities than high density RGC types mosaics. It 342 would be very informative to experimentally verify this prediction. To this date, this question 343 remains unanswered. This hypothesis is in accordance with recent studies showing that some low 344 density RGCs do not exhibit regular spacing (33). Moreover, we question here the relevance of 345 regular spacing as a criterion for a group of RGC to form a RGC type. Indeed, if all low density 346 RGC types do not exhibit highly regular spacing as predicted here, this criterion does not 347 discriminate RGC types. 348 We show here that high mosaic regularity can be achieved with limited migration distance (8.72μm 349 ±0.11 in average, n = 8). This average migration distance is in accordance with in-vivo 350 measurements, reporting that RGCs and SACs tangential migration does not exceed 30μm (22). 351 However, the average migration distance measured in our simulations is notably lower than the 352 average migration distance experimentally measured at around 20μm (19) and could be explained 353 by the absence of retinal surface expansion implementation in our simulations. The CM mechanism 354 implemented here is based only on local cues and short-distance interactions, and thereby follows 355 the description of tangential dispersion in mouse, reported as a local, short-distance, phenomenon 356 (21). Our results are consistent with previous studies showing that a tangential cell dispersion does 357 not appear to be directly related to the cell time of birth, but rather to its cell type (21). 358 The cellular migration and RI dynamics resulting from the CM mechanism are in agreement with 359 the literature, where it is reported that RI increases mostly between P1 and P5, with the spacing 360 between cells still increasing after that period, until P10 (2). After reaching the correct cell layers, a 361 slower and finer tangential positioning phase of RGC within the GCL has been reported (19,34). 362 Cellular movement during this period has been described as random but important for exact cellular 363 positioning (34). In accordance with our results and as stated by other studies (27), these highly 364 varied movements are likely to be related to mosaic formation and refinement. Indeed, these 365 movements appear random as the whole RGC population (On, Off, On-Off population) is 366 considered, while RGC populations should themselves be divided into types in order to 367 meaningfully investigate RGCs lateral migration. If it were possible to examine each type 368 independently, our model suggests that these movements, reported as random, would appear as 369 coherent, as illustrated by Figure 7. populations. We found no RI variation from P3-4, indicating that these two SACs populations have 379 already created their mosaics by P3-4, hence shortly after SACs migration into their respective 380 cellular layer. In addition, we observed that the calculated exclusion factor does not vary, also 381 supporting this assumption. Moreover, the observed complementarity of GCL and INL mosaics 382 perhaps indicates interactions between these two SAC populations during their cellular 383 organisation, before they migrate to their respective layer. 384 Here, we investigated this GCL/INL population interaction hypothesis further by building a 385 simulation of SACs mosaics development. These simulations clearly show that our modelling 386 procedure can successfully be applied to another cell population, without changing any simulation 387 parameters. Indeed, we have been able to explain differences in GCL and INL mosaic regularities 388 (the RI of the INL population being higher than that of the GCL population) by using only local 389 interactions between SACs. This is the case if SACs constitute a unique population, or if GCL and 390 migrating to the INL compared to the GCL. Precisely, the percentage of a population characterised 393 by a highly regular mosaic dictates the regularity of the resulting sub-population. Hence, the bigger 394 the sub-population, the closer the obtained RI will be to the RI of the initial population, if cells 395 constituting this sub-population are chosen randomly. For instance, if a population with a high RI is 396 randomly divided into two sub-populations of 80% and 20% of the initial population (denoted 397 respectively sub-population A and B), the sub-population A will have a RI closer to the initial 398 population than the sub-population B. Thus, in the latter case, this observed RI difference between The average RGC and SAC density for each developmental day was measured by performing a 433 manual cell count from P2 to P10 for RGCs and from P4 to P10 for SACs. By accounting for the 434 surface expansion observed during retinal development, we estimated changes in populations 435 through development. The estimated total RGC and SAC populations for a given retina are 436 calculated by multiplying the averaged cell density (obtained from 3-6 sample areas per retina) by 437 its corresponding retinal surface. These individual measurements are then averaged for each 438 developmental day to give an estimation of the total population from P2 to P10. Cell population 439 death rate during development is then calculated. In detail, the population of each retina measured 440 on day D+1 is subtracted from the population of each retina measured on day D to calculate the 441 amount of CD between day D and D+1. The amount of apoptosis measured between two 442 consecutive days is then averaged to calculate the daily death rate of RGC and SAC populations. Simulations were conducted using the agent-based simulation framework BioDynaMo (35). 461 such as its 3D geometry, mass, adherence and position in space. Individual neurons are represented 463 by a sphere. Diffusion in 3D of chemical substances in the extracellular space has also been 464 implemented, with the discrete central difference method. This diffusion is supported by grids 465 representing substances concentration and gradients. Mechanical forces are also taken into account 466 between all simulation objects such that they cannot overlap, but mechanically repulse each other. 467 Each simulation object can have a biology module attached to it, that describes its behaviour at each 468 simulation time step, such as substance secretion, cell migration or cell growth. 469 As an AB simulation framework, each simulation object is independent, without a central

Simulations 477
All simulations took place in a cubic space of 1,300μm 3 , with cells of 7 to 8μm diameter randomly 478 distributed (uniform distribution) in a space of 1,000μm×1,000μm×22μm. The initial cell density 479 has been set to 8600 cells/mm 2 , in order to reach the RGC density once programmed CD 480 mechanism is over -around 3000cells/mm 2 reported in literature (3)  The CD mechanism corresponds to the cells removing themselves from the simulation if their 513 corresponding substance concentration is higher than a defined threshold. In this way, the clusters of 514 homotypic cells exhibit high death rates and become sparser. As the cell density decreases, the 515 initial multilayer collapses into a RGC monolayer. This is implemented as cells moving along the z 516 axis toward the centre of the RGC layer, using their chemical cue. CD is triggered after completion 517 of CF, and continues until a steady-state is reached, around a death rate of 65%. This steady-state is 518 reached without global controllers but depends on the chosen concentration threshold triggering cell 519 death. If this threshold is low the steady-state will be reached with a high death rate, while if this 520 threshold is high the steady-state will be reached with a low death rate. 521

RGC mosaic formation: CM 522
CM is implemented such that the homotypic substances act as a repulsive factor. Thereby, cells 523 exhibit short distance avoidance, moving tangentially against their substance gradient, distancing 524 themselves from homotypic neighbours. We assume that CM is triggered after completion of CF, at 525 the same time as CD, and continues either until a steady state or day 13 is reached. 526 Development conditions incorporating all combinations of these three mechanisms have been 527 investigated. As the mechanisms influence each other, parameters vary depending on the 528 implemented mechanisms. The CD mechanism parameters were chosen for each RGC type such 529 that its final death rate is about 65%. The CM parameters were chosen depending on the CD 530 parameter value and such that the interaction is kept to close range distance. Pseudocode 531 corresponding to these three biology modules can be found in Pseudocode 1. Table 1 summarises 532 the parameters used for RGC mosaics formation mechanisms. triggering homotypic avoidance (tangential migration mechanism). Once mosaics are formed, the 539 two populations migrate to their respective layers, along the Z axis. Two different developmental 540 conditions have been implemented, using either one common or two separated (one for GCL 541 population, one for INL population) chemical substances for mosaics formation. Importantly, 542 concentration thresholds are identical for the GCL and INL populations. Parameters have been set 543 such that the mosaic RIs match the measured RIs in mouse SACs mosaics. 544

Data analysis 545
The RI was used to assess the regularity of the mosaics. It is computed as the average value of the 546 closest neighbour distribution (distribution of the closest neighbour measured for each cell) divided 547 by its standard deviation (36). The RI offers a single score that is able to discriminate regularity 548 differences between mosaics of low regularities. In addition, and as previously reported by Reese 549 and Keeley (2015), the RI offers a scale-invariant measure of mosaic regularity and thus more direct 550 evidence of any change in the mosaic spatial organisation during development. It is not only the 551 absolute RI value that carries information, but also its evolution across development, related to the 552 contribution of each mosaic developmental mechanism (CF, CD, CM). However, it should be 553 pointed out that RI is sensitive to a low sampling rate, leading to significant variability in RI scores 554 for mosaics constituted of few cells. 555 Comparisons between two RI values have been conducted using T-tests for two independent 556 samples. 557