A minimal statistical-mechanical model for multihyperuniform patterns in avian retina
Enrique Lomba, Jean-Jacques Weis, Leandro Guisández, Salvatore Torquato
AA minimal statistical-mechanical model for multihyperuniform patterns in avian retina
Enrique Lomba , Jean-Jacques Weis , Leandro Guis´andez , , and Salvatore Torquato , Instituto de Qu´ımica F´ısica Rocasolano, CSIC, Calle Serrano 119, E-28006 Madrid, Spain IFLYSIB (UNLP, CONICET), 59 No. 789, B1900BTE La Plata, Argentina Universit´e de Paris-Saclay, Laboratoire de Physique Th´eorique, Bˆatiment 210, 91405 Orsay Cedex, France Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA Princeton Institute for the Science and Technology of Materials,Princeton University, Princeton, New Jersey 08544, USA
Birds are known for their extremely acute sense of vision. The very peculiar structural distribu-tion of five different types of cones in the retina underlies this exquisite ability to sample light. Itwas recently found that each cone population as well as their total population display a disorderedpattern in which long wave-length density fluctuations vanish [Jiao et al., Phys. Rev. E, , 022721(2014)]. This property, known as hyperuniformity is also present in perfect crystals. In situationslike the avian retina in which both the global structure and that of each component display hy-peruniformity, the system is said to be multi-hyperuniform. In this work, we aim at devising aminimal statistical-mechanical model that can reproduce the main features of the spatial distri-bution of photoreceptors in avian retina, namely the presence of disorder, multi-hyperuniformityand local hetero-coordination. This last feature is key to avoid local clustering of the same typeof photoreceptors, an undesirable feature for the efficient sampling of light. For this purpose weformulate a simple model that definitively exhibits the required structural properties, namely anequimolar three-component mixture (one component to sample each primary color, red, green, andblue) of non-additive hard disks to which a long-range logarithmic repulsion is added between likeparticles. A Voronoi analysis of our idealized system of photoreceptors shows that the space-fillingVoronoi polygons interestingly display a rather uniform area distribution, symmetrically centeredaround that of a regular lattice, a structural property also found in human retina. Disordered multi-hyperuniformity offers an alternative to generate photoreceptor patterns with minimal long-rangeconcentration and density fluctuations. This is the key to overcome the difficulties in devising anefficient visual system in which crystal-like order is absent. I. INTRODUCTION
Sampling light is one of the essential activities thatenables the interaction of living organisms with the sur-rounding environment. From simple devices such asthe stigma that provides “vision” in certain classes ofmicroalgae [1], to the sophisticated compound eyes ofarthropods [2, 3] living organisms have developed increas-ingly efficient ways to map visual information from theexternal world onto signals that can be processed by theircognitive systems. The case of arthropod eyes is particu-larly interesting. It is known from classical sampling the-ory [4] that an optimal sampling of light can be achievedby an hexagonal array of photodetectors. This is actu-ally the pattern adopted by ommatidia (the optical unitsforming a compound eye) in arthropods. Compound eyesare imaging systems with low aberration, wide-angle fieldof view and infinite field depth[5]. These properties havemotivated intense research into the development of bioniccompound eyes intended for small robots[6] or sensors fordigital cameras[5].When it comes to vertebrates, with the exception ofsome teleost fish [7, 8] and some reptiles [8, 9], the situa-tion is different and structural disorder in photoreceptorpatterns is the general trend. In this connection, birdsare in a class of their own. They possess one of the mostelaborate visual systems among vertebrates. In avianretina, one can find five different types of cones, [10, 11],one type for luminance detection and the remaining four building a tetrachromatic color sensing device coveringwavelengths from red to ultraviolet [12]. In contrast withthe regular shape of ommatidia in insects, photorecep-tors in bird retina are polydisperse, in size and number[13, 14]. This variation provides an adaptative advan-tage: changing the relative numbers and even pigmenta-tion of the cones bird species can have visual capabilitiesadapted to different habitats (sea birds have high den-sity of red/yellow cones for hazy conditions, nocturnalbirds have a extremely high density of luminance cones,...). However, polydispersity is known to frustrate crys-tallization [15], so an alternative to the regular hexagonalpattern of arthropod eyes is needed if we want to pre-serve a good sampling of light. In this connection, Jiaoand coworkers [8] found that the spatial distribution ofphotoreceptors in chicken retina retained some “hiddenorder” reminiscent of crystalline patterns. Namely, theyfound that long-range density and concentration fluctu-ations were vanishingly small. This feature can be quan-tified by means of two intimately connected structuralproperties. In two dimensions, we have first the numbervariance of cones over a sampling area of radius R, de-fined as σ N ( R ) = (cid:104) N (cid:105) R − (cid:104) N (cid:105) R , where N is the numberof cones contained in the sample area and (cid:104) . . . (cid:105) R denotesthe average over a certain number of sampling areas. InRef. [8], it was found that this quantity obeys the follow-ing large- R asymptotic scaling σ N ( R ) ∝ R (1) a r X i v : . [ phy s i c s . b i o - ph ] M a y in the plane. This is one of the possible scalings of hy-peruniform systems, also characteristic of crystalline-likeorder in two dimensions (class I following Ref. [16]). Sec-ondly, it is known that density fluctuations in Fourierspace are directly related to the structure factor. This isdefined for a set of points/particles with number density ρ by S ( Q ) = 1 + ρ ˜ h ( Q ) , (2)where Q is the wave vector, ˜ h ( Q ) is the spatial Fouriertransform of h ( r ) = g ( r ) −
1, being g ( r ) the pair distri-bution function of the point/particle configuration. It ispossible to show[16] that for a system satisfying Eq. (1)in two dimensions then S ( Q ) ∝ Q α ( Q →
0) (3)with α >
0. Since Eq. (1) holds for each cone distribu-tion, then we will have a relation like (3) for the structurefactor computed from each cone pattern, as was found byJiao and coworkers [8], i.e.lim Q → S ii ( Q ) = 0 (4)for each cone type i . This implies that density fluctua-tions of the corresponding point patterns will vanish forlong wavelengths, i.e. when Q →
0. The same applies tothe overall point pattern. This property was termed inRef. [8] as “multi-hyperuniformity”.Interestingly, since Torquato and Stillinger [17] intro-duced the concept of hyperuniformity and stressed itssignificance in structurally disordered materials, such ex-otic “states of matter” have been found in a wide vari-ety of systems. A partial list of examples include amor-phous dielectric networks with large and complete pho-tonic band gaps [18, 19], dense transparent disorderedmedia [20], the enhanced pinning of vortices in arraysin superconductors[21], certain composites with desirabletransport, dielectric and fracture properties [22–25], sandpiles and other avalanche models [26, 27], driven nonequi-librium granular and colloidal systems [28–30] and evenimmune system receptors [31] all have in common thepresence of hyperuniformity.One might ask why hyperuniformity plays such a cru-cial role in the quality of vision in birds ?. As men-tioned, the optimal sampling configuration of photore-ceptors corresponds to a fully regular hexagonal arrange-ment. Hyperuniformity prevents long-wavelength fluc-tuations in the photoreceptor density (or concentrationof different species) that would be otherwise be presentin a structurally disordered configuration of photorecep-tors. The presence of such fluctuations is certainly not adesirable property for an accurate image representation.Fully regular arrangements such as the hexagonal pat-terns of ommatidia are hyperuniform, but in the case ofbird retina, crystal-like order is preempted by polydisper-sity. Thus hyperuniform patterns might well be a good compromise solution. Multi-hyperuniformity will guar-antee the same sampling quality for each type of photore-ceptors and aids in ensuring local hetero-coordination,which is key to prevent the unwanted clustering of samecolor photoreceptors.After these considerations, it is our aim to build aminimal statistical mechanical model that can reproducethe main characteristics of the photoreceptor distribu-tion. These are, in addition to disorder, on one handmulti-hyperuniformity, and on the other local hetero-coordination. By this we mean that photoreceptors ofthe same type should not be allowed to cluster togetherif color sensitivity is to be uniformly distributed on space.In fact, in Ref. [34] it was shown that a system can bemulti-hyperuniform and display a strong degree of clus-tering (chain formation). From pictures of actual chickencone distributions (see Figure 1 in Ref. [35]) it is read-ily apparent that cones of different types tend to clustertogether, i.e. their spatial distribution displays hetero-coordination.The findings of Ref. [34] suggest that a mixture withlogarithmic long-range repulsions and non-additive hard-core volume exclusions can display the sought character-istics. Strictly speaking the model in question was a two-dimensional Coulomb plasma. Interestingly, in Ref. [8] itwas found that the structure factor derived from pho-toreceptor patterns displays a small wave number de-cay consistent with ∼ Q (or ∼ Q when fitted into amultiscale packing model [8]). Strictly two-dimensionalCoulomb plasmas are known to have structure factorsthat decay quadratically with the wavenumber as Q → II. MODEL AND METHODS
As mentioned, our minimal model of “retina” consistsof three classes of photoreceptors, (red-green-blue=RGB)in which, following Ref. [34], interactions will be definedin terms of a purely repulsive logarithmic potential. Inaddition, in order to guarantee hetero-coordination frommoderate to high densities, the particles will have a hard-core volume exclusion defined by a hard-disk diameter σ ,with unlike particles having a distance of minimum ap-proach (1 + ∆) σ , with ∆ <
0. From Ref. [34], we knowthat ∆ > i and j can be explicitly written as βu ij ( r ) = (cid:26) ∞ if r < (1 + ∆(1 − δ ij )) σ − γ ij log r/σ if r ≥ (1 + ∆(1 − δ ij )) σ (5)where γ ij is an effective coupling parameter, and δ ij isKronecker’s symbol. Our minimal model is fully sym-metric, with u ii = u ∀ i , and u ij = u ∀ i (cid:54) = j . Forthe logarithmic repulsion the coupling parameter is ex-pressed as γ ij = ( λ + (1 − λ ) δ ij )Γ (6)with 0 ≤ λ ≤
1. The parameter λ controls the non-additivity of the long-range interactions, and we willsee it determines whether the system displays multi-hyperuniformity or not.Now, from our study on binary mixtures in Refs.[34,37] we know that disordered systems with long-rangedrepulsive interactions, whose small wavenumber scalingin Fourier space followslim Q → β ˜ u ij ( Q ) ∝ Q − α (7)with α > Q → | ˜u ( Q ) | (cid:54) = 0 (8)where | . . . | denotes a matrix determinant, and the ele-ments of the matrix ˜u ( Q ) are the Fourier transform ofthe species-species interactions, u ij ( r ). It can be shownthat a sufficient condition for Eq. (8) to be fulfilled isthat lim Q → (cid:2) ˜ u ii ( Q )˜ u jj ( Q ) − ˜ u ij ( Q ) (cid:3) (cid:54) = 0 , (9) which actually means that cross interactions must not comply with the Lorentz-Berthelot mixing rules in thelong wavelength limit. This we had already found forbinary mixtures in Ref. [34]. In practice, for our modelsystem this means that λ <
1. Here we will simply set λ = 0 which reduces cross interactions to bare hard-corerepulsions.The low- Q asymptotics of the structure factor whenall densities are identical ( ρ i = ρ/ ∀ i ) simplifies consid-erably. A detailed derivation can be found in the sup-plementary information based on the low- Q expansion ofthe OZ equation. In our particular case, given the sym-metry of the interactions and compositions and setting λ = 0, from Eq. (S.11) in the SI the limiting behavior ofthe partial structure factors reduces tolim Q → S ii ( Q ) = Q / (2 πρ Γ) , ∀ i lim Q → S ij ( Q ) = ρ ˜ c Rij (0) Q / (2 πρ Γ) , ∀ i (cid:54) = j. (10)with ˜ c Rij ( Q ) being the Fourier transform of the shortrange component of the direct correlation function (cf SIfor further details), which is non-zero and finite as Q → S NN ( Q ) = (cid:88) i,j S ij ( Q ) . (11)From Eq. (10) we then havelim Q → S NN ( Q ) = 3 Q / (2 πρ Γ) + bQ . (12)where b = 3 ρ ˜ c Rij (0) / (2 πρ Γ) .The systems studied in this work have been analyzedusing an integral equation approach based on the OZequation, Eq. (S.1), with a Reference Hypernetted Chain(RHNC) closure (Eq. (11) in Ref. [37]). We refer thereader to [37] for further details on the numerical ap-proach to solve this equation. We have also performedextensive canonical (NvT) Monte Carlo simulations, inwhich the energy of the periodic system is evaluatedusing the Ewald technique with conducting boundaryconditions[37, 38]. Computational details of the simu-lations are identical to those of Ref. [37]. III. RESULTS
We will first consider two instances of photoreceptorpatterns for low density ( ρσ = 0 . ρσ = 0 . − . σ ij = 0 . σ ii ). The cou-pling factor of the long range interaction is set to Γ = 5, (a) ρσ = 0 . ρσ = 0 . FIG. 1. Snapshots of Monte Carlo configurations of our three component minimal model of retina with rgb receptors (as shownin the Figure). The interaction is defined by a coupling γ = 5 and a non-additivity parameter ∆ = − and the long-range cross interactions are set to zero (i.e. λ = 0 in (5)). This means that unlike particles will onlyinteract via a pure hard core exclusion. For comparisonwe will also show results for λ = 1 which will only dis-play global hyperuniformity. Partial densities for eachphotoreceptor type are ρσ / ρσ = 0 . ≈ − σ . When com-paring this illustration with real representations of birdcone distributions (see Figure 1 in Ref. [35]) the similar-ity is evident. At higher densities ( ρσ = 0 .
8) packingeffects become dominant and clustering is not so appar-ent, but hetero-coordination is still clearly seen in thesnapshot of Figure 1(b). The cluster size distribution(not shown) is monotonously decreasing, with no domi-nant cluster size. This is a consequence of the lack of acompeting short range attraction that would counteractthe long-range repulsion and would thus stabilize finitesize clusters, as it was the case for ∆ > A. Structure factor analysis
In Figure 2 we plot the partial and total structure fac-tors corresponding to the systems described above. Themulti-hyperuniform character of the system is clearly il-lustrated by their vanishing behavior for low- Q . In theinsets one can observe that they closely follow the asymp-totic behavior described by Eqs. (10) and (12). Theory and simulation agree to a very large extent.For comparison we also plot the theoretical results for ρσ = 0 .
2, and λ = 1. Now, this choice of the long-range cross interactions leads to a globally hyperuniformconfiguration, as confirmed by the behavior of S NN ( Q )as Q →
0. In contrast, the partial structure factor doesnot vanish for Q → S ii ( Q ) ≈ x i ∀ Q . Reducing λ , which actually implies decreasing (or in our presentcase, eliminate) the unlike long-range repulsive interac-tions, induces a certain degree of clustering between un-like particles. This effect is visible when comparing thetotal structure factor at low density (lower graph, redcurve in Figure 2) for λ = 0 and λ = 1. Only in thecase of λ = 0, S αα ( Q ) exhibits a prepeak at Qσ ≈ . ≈ . σ which we have already qualitatively de-tected in the snapshot of Figure 1(a). In summary, thecombination of very long-ranged repulsions between likeparticles with non-additive unlike interactions both in theshort and long range reproduces the features sought forin our minimal statistical mechanical model of retina. B. Mimicking avian retina
How do our model results compare with a real struc-ture factor obtained from a distribution of avian photore-ceptors ? One must first bear in mind that in bird retinafive different types of photoreceptors [8] are present inunequal numbers, so in principle our model departs sig-nificantly from the real situation. Nonetheless, a sim-ple inspection of the experimental structure factors pre- S aa ( Q ) rs =0.8rs =0.2 s S NN ( Q ) rs =0.2,l=1, RHNC /(2 prG) /(2 prG )+bQ G=5l=0
FIG. 2. Total and partial structure factors of our model pho-toreceptor system displaying multi-hyperuniformity for lowand moderate number densities (see legend). Solid and dash-dotted curves correspond to theoretical (RHNC) calculations,symbols denote Monte Carlo data. Dashed curves in the in-sets represent the low- Q regime derived from Eq.(10). Partialstructure factors correspond to correlations between like par-ticles. For comparison we show on blue dash-dotted curvesthe structure factors for a system displaying only global hy-peruniformity ( λ = 1 in Eq. (5). sented in Figure 9 of Ref. [8] indicates that basically allphotoreceptor species qualitatively display similar partialstructure factors. The global structure factor is qualita-tively different, with very little structure at low Q values.Therefore, it seems reasonable to compare our model sys-tem with actual experimental results from a qualitativestandpoint. To that aim, we have adjusted the couplingconstant Γ of our effective potential to match the resultsof [8]. Density is basically coupled to Γ (except for subtlehard core effects not visible in the low Q behavior of thestructure factor), so we have set ρσ = 0 . Q is scaled with the position of thestructure factor maximum, which sets the length scale tothe appropriate value in order to ease the comparison.This is equivalent to rescaling the data so as to accountfor the appropriate sizes of the photoreceptors. We ob-serve that the behavior of our simple model depicted inFigure 3 agrees qualitatively with the experimental data.As a matter of fact, even if in Ref. [8] the experimentallow- Q behavior seems to follow a linear decay instead of S rr ( Q ) Red photoreceptorsmodel system,
G=30 m S NN ( Q ) All photoreceptors
G=5l=0
FIG. 3. Total and partial structure factor of the symmet-ric three component plasma with negative non-additivity for λ = 0, ρσ = 0 .
8, Γ = 30 compared with those from avianphotoreceptors from Ref. [8]. In the upper curve only thered photoreceptors are presented. Solid curves correspond toRHNC calculations, symbols to experimental data[8]. Thepartial structure factor is normalized to one and the Q axisis scaled with the position of the maximum which is equiv-alent to adjust σ to the effective experimental value in thephotoreceptor correlations. the Q dependence of our model, the quadratic depen-dence appears acceptable. Interestingly, the multiscalepacking model also proposed by Jiao and coworkers[8]displays the same quadratic decay. The other salient fea-ture that is observed in Figure 3 is the lack of structureof S NN ( Q ) for low Q values. This feature is visible bothin our simple model (although somewhat enhanced) andin the experimental data. It is apparent that our mini-mal model is capable of reproducing key features of thespatial patterns display by photoreceptors in actual birdretina. C. Voronoi analysis
When thinking of photoreceptors, one must also takeinto account that their ability to reproduce an image isdirectly related with the area they sample. This suggeststhat a Voronoi analysis of our point configurations willprovide information as to the sampling area correspond- (a)Photoreceptor model, red cones with ρ red σ = 0 . / (b)Photoreceptor model, rgb cones, ρσ = 0 . FIG. 4. Voronoi tessellations corresponding to the red cones (left) of the photoreceptor model and all the photoreceptors (right)for total density ρσ = 0 . ing to each particle. We have therefore performed a char-acterization of the spatial configurations of our modelsystem using Voronoi tessellations. We have studied thecorresponding area distribution of the Voronoi polygons.In order to put these results in perspective, we have alsoperformed a corresponding analysis for purely randomtwo-dimensional point configurations, as well as config-urations obtained from Molecular Dynamics simulationsfor 2D fluids of Lennard-Jones (LJ) particles and LJ par-ticles with added Coulomb repulsion. In the last instance,the competition between short-range attractive and long-range repulsive forces leads to the formation of stableclusters that nonetheless display hyperuniformity. Forthese cases, we have used similar density conditions andsupercritical temperatures ( k B T /(cid:15) = 2 .
0, where k B isBoltzmann’s constant and T the absolute temperature).When referring to LJ results, (cid:15) and σ correspond to thewell depth and particle size respectively. These resultsare presented in detail in the Supplementary Information.Upon examination of Figure 4, one can clearly appreci-ate that in our model system the Voronoi tessellation ex-hibits a fairly regular distribution of the polygon areas.A similar observation was made by Legras et al. [44]when analyzing the cone distribution in human retina.However, when comparing with tessellations for LJ flu-ids or random distributions (see Figure 1 of the SI) atsimilar density, one finds that these have a much largerdispersion in their areas. This can be more quantitativelyanalyzed by examining the normalized area distributions.These are plotted in Figure 5 vs the area scaled with thecorresponding particle densities, ρA . The correspondingfigure for random and LJ can be found in the SI.Figure 5 reveals that our model leads to the area distri-butions that are symmetrized with respect to the regular lattice result, ρA = 1. It is interesting to note that inall cases the curves apparently follow a Gaussian distri-bution. In contrast, area distributions for random con-figurations and LJ particles at low density follow highlyasymmetric Γ-distributions, and denser packings of LJparticles can be fit to log-normal distributions (see SI). Itis important to note that the symmetrization of the areadistributions is not a consequence of hyperuniformity. InRef. [33], it was found that certain stealthy hyperuniformpatterns led to asymmetric distributions similar to thoseof random configurations. This is also illustrated by theanalysis of LJ particles with added Coulombic repulsions,which form hyperuniform patterns with a strong degreeof clustering. The area distribution of the Voronoi poly-gons is a short/medium range property, unlike hyperuni-formity which is a large-scale property.For our retina model, the marked symmetry of theVoronoi area distribution is due to the fact that the in-teractions are monotonic and repulsive, and thus tendto produce very regular local environments. This is illus-trated by the uniform linear behavior of number variance σ N ( R ) for both the global and the single species pho-toreceptor configurations in our model, as can be seenin Figure 6. The same linear dependence is found inthe experimental photoreceptor distributions (see Figure4 in Ref. [8]). Conversely, hyperuniform configurationsthat display clustering (and asymmetric area distribu-tions of their Voronoi tessellations) present a numbervariance with a clear non-monotonic behavior for smallsampling windows. Such is the case of particles withLJ+Coulombic interactions as illustrated in Fig. 3 of theSI.One can then interpret interpret the symmetrizationof the area distribution in a disordered media as the con- D ( A ) rs =0.8 l=1 (all) r s =0.2667 l=0 (r/g/b species) r a A0123 D ( A ) rs =0.2 l=0 (all) r s =0.0667 l=0 (r/g/b species) FIG. 5. Scaled area distribution of the Voronoi polygons forhyperuniform states: high density (upper graph) and low den-sity (lower graph) single species and global configurations ofthe three-component plasma. Hyperuniformity symmetrizesthe area distributions around the value of the square regularlattice ( δ ( ρA − sequence of the minimization of repulsive interactions,maximizing the area around each point in the config-uration, and reflecting random deviations from the thecrystalline (ordered hyperuniform) state. This, togetherwith the strong suppression of long-wavelength densityand concentration fluctuations leads to what could pos-sibly be optimal photoreceptor patterns. IV. CONCLUDING REMARKS
In summary, we have shown that two key featuresof the experimental patterns of photoreceptors in birdretina, multi-hyperuniformity and hetero-coordination,can be captured by a simple model with logarithmic re-pulsions between like particles and hard core exclusionswith negative non-additivity between the unlike ones.The fact that disordered hyperuniform systems represent s < N > R - < N > R Red global
FIG. 6. Number variance σ N ( R ) dependence of sample win-dow radius, R , for both the global and the single species pho-toreceptor model configurations. topological states of matter sharing fluid and crystal-likeproperties makes them the solution of choice when reg-ular arrangements such as those of arthropod eyes arehampered by the variability of the photoreceptors (e.g.unequal sizes and numbers). Present day bio-inspiredoptical devices rely on regular arrangements [5, 41]. Incertain instances, the combination of different types of re-ceptors and in different number might be required com-promising the feasibility of regular arrays of receptors.Disordered multi-hyperuniformity might then offer an al-ternative to overcome these difficulties.A natural extension of this work should be the exten-sion to non-Euclidean geometries. Steps in this directioncan be found in the works of Meyra et al. [42] and Boˇziˇcand ˇCopar [43] for spherical surfaces. Actually designson curved surfaces have been already proposed for reg-ular arrays in Ref. [5]. Disorder hyperuniform systemson curved surfaces might well have a potentially largerimpact on technological applications. An analysis alongthese lines of photoreceptor patterns in humans[44] mightalso be of interest to further our understanding of ourcomplex visual system. In fact, in Ref. [44] it was foundthat human cones tend to preserve locally hexagonal ar-rangements. A preliminary Voronoi analysis of the areadistributions taken from Ref. [44] show that these arealso relatively symmetric. On the other hand, photore-ceptor patterns in human retina are also dependent onthe eccentricity of the sampling area, and this consider-ably complicates the analysis. ACKNOWLEDGMENTS
EL acknowledges the support from the Agencia Es-tatal de Investigaci´on and Fondo Europeo de DesarrolloRegional (FEDER) under grant No. FIS2017-89361-C3-2-P. EL and LG were also supported by the EuropeanUnion’s Horizon 2020 Research and Innovation Staff Ex-change programme under the Marie Sk(cid:32)lodowska-Curiegrant agreement No 734276. S.T. was supported by the National Science Foundation under Award No. DMR-1714722. The authors are grateful to Prof. Y. Yiao forkindly providing the structure factor data from Ref. [8]. [1] P. Hegemann, Planta , 265 (1997).[2] D. F. Ready, T. E. Hanson, and S. Benzer, Dev. Biol. , 217 (1976).[3] D. K. Lubensky, M. W. Pennington, B. I. Shraiman, andN. E. Baker, Proc. Natl. Acad. Sci. U.S.A. , 11145(2011).[4] D. P. Petersen and D. Middleton, Inf. Control , 279(1962).[5] Y. M. Song, Y. Xie, V. Malyarchuk, J. Xiao, I. Jung, K.-J. Choi, Z. Liu, H. Park, C. Lu, R.-H. Kim, R. Li, K. B.Crozier, Y. Huang, and J. A. Rogers, Nature , 95(2013).[6] Y. Zheng, L. Song, J. Huang, H. Zhang, and F. Fang,Optics Letters , 4143 (2019).[7] P. A. Raymond and L. K. Barthel, Int. J. Dev. Biol. ,935 (2004).[8] Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C.Corbo, and S. Torquato, Phys. Rev. E , 022721 (2014).[9] R. F. Dunn, J. Ultrastruct. Res. , 672 (1966).[10] M. Ruggeri, J. C. Major, C. McKeown, R. W. Knighton,C. A. Puliafito, and S. Jiao, Invest. Opthalmol. Vis. Sci. , 5789 (2010).[11] M. B. Toomey and J. C. Corbo, Front. Neural Circuit. , 97 (2017).[12] J. Withgott, BioScience , 854 (2000).[13] T. H. Goldsmith, J. S. Collins, and S. Licht, Vision Res. , 1661 (1984).[14] B. A. Moore, P. Baumhardt, M. Doppler, J. Rando-let, B. F. Blackwell, T. L. DeVault, E. R. Loew, andE. Fernandez-Juricic, J. Exp. Biol. , 3442 (2012).[15] D. Frenkel, “Introduction to colloidal systems,” (Taylorand Francis, 2006).[16] S. Torquato, J. Phys. : Condens. Matter , 414012(2016).[17] S. Torquato and F. H. Stillinger, Phys. Rev. E , 041113(2003).[18] M. Florescu, S. Torquato, and P. J. Steinhardta, Proc.Natl. Acad. Sci. U.S.A. , 20658 (2009).[19] L. S. Froufe-P´erez, M. Engel, J. J. S´aenz, and F. Schef-fold, Proc. Natl. Acad. Sci. U.S.A. , 9570 (2017).[20] O. Leseur, R. Pierrat, and R. Carminati, Optica , 763(2016).[21] Q. Le Thien, D. McDermott, C. J. O. Reichhardt, andC. Reichhardt, Phys. Rev. B , 094516 (2017).[22] G. Zhang, F. H. Stillinger, and S. Torquato, J. Chem.Phys , 244109 (2016).[23] D. Chen and S. Torquato, Acta Mater. , 152 (2018). [24] Y. Xu, S. Chen, P.-E. Chen, W. Xu, and Y. Jiao, Phys.Rev. E , 043301 (2017).[25] B.-Y. Wu, X.-Q. Sheng, and Y. Hao, PloS one ,e0185921 (2017).[26] R. Dickman and S. D. da Cunha, Phys. Rev. E ,020104R (2015).[27] R. Garcia-Millan, G. Pruessner, L. Pickering, andK. Christensen, Europhys. Lett. , 50003 (2018).[28] D. Hexner and D. Levine, Phys. Rev. Lett. , 110602(2015).[29] J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo,Phys. Rev. Lett. , 108301 (2015).[30] E. Tjhung and L. Berthier, Phys. Rev. Lett. , 148301(2015).[31] A. Mayer, V. Balasubramanian, T. Mora, and A. M.Walczak, Proc. Nat. Acad. Sci. , 5950 (2015).[32] O. U. Uche, F. H. Stillinger, and S. Torquato, Phys. Rev.E , 046122 (2004).[33] G. Zhang, F. H. Stillinger, and S. Torquato, Phys. Rev.E , 022119 (2015).[34] E. Lomba, J.-J. Weis, and S. Torquato, Phys. Rev. E , 010102(R) (2018).[35] Y. A. Kram, S. Mantey, and J. C. Corbo, PLoS ONE ,e8992 (2010).[36] J. M. Caillol, D. Levesque, J. J. Weis, and J. P. Hansen,J. Stat. Phys. , 325 (1982).[37] E. Lomba, J. J. Weis, and S. Torquato, Phys. Rev. E , 062126 (2017).[38] S. Leeuw and J.W.Perram, Physica A , 546 (1982).[39] M. Ferraro and L. Zaninetti, Physica A: Statistical Me-chanics and its Applications , 4575 (2012).[40] M. A. Klatt, J. Lovri´c, D. Chen, S. C. Kapfer, F. M.Schaller, P. W. A. Sch¨onh¨ofer, B. S. Gardiner, A.-S.Smith, G. E. Schr¨oder-Turk, and S. Torquato, NatureComm. , 811 (2019).[41] M. Garcia, C. Edmiston, T. York, R. Marinov, S. Mon-dal, N. Zhu, G. P. Sudlow, W. J. Akers, J. Margenthaler,S. Achilefu, R. Liang, M. A. Zayed, M. Y. Pepino, andV. Gruev, Optica , 413 (2018).[42] A. G. Meyra, G. J. Zarragoicoechea, A. L. Maltz,E. Lomba, and S. Torquato, Phys. Rev. E ,022104(R) (2019).[43] A. L. Boˇziˇc and S. ˇCopar, Phys. Rev. E , 032601(2019).[44] R. Legras, A. Gaudric, and K. Woog, PLOS ONE13