Biological tissue-inspired tunable photonic fluid

Model for Epithelial Cell Packing in 2D.

When epithelial and endothelial cells pack densely in 2D to form a confluent monolayer, the structure of the resulting tissue can be described by a polygonal tiling (18). A great variety of cell shape structures have been observed in tissue monolayers, ranging from near-regular tiling of cells that resembles a dry foam or honeycomb lattice (19) to highly irregular tilings of elongated cells (20). To better understand how cell shapes arise from cell-level interactions, a framework called the Self-Propelled Voronoi (SPV) model has been developed recently (21, 22). In the SPV model, the basic degrees of freedom are the set of 2D cell centers {ri}, and cell shapes are given by the resulting Voronoi tessellation. The complex biomechanics that govern intracellular and intercellular interactions can be coarse-grained (18, 19, 23⇓⇓⇓–27) and expressed in terms of a mechanical energy functional for individual cell shapes.E=∑i=1NKA(Ai−A0)2+KP(Pi−P0)2.[1]The SPV energy functional is quadratic in both cell areas ({Ai}) with modulus KA and cell perimeters ({Pi}) with modulus KP. The parameters A0 and P0 set the preferred values for area and perimeter, respectively. Changes to cell perimeters are directly related to the deformation of the acto-myosin cortex concentrated near the cell membrane. After expanding Eq. 1, the term KPPi2 corresponds to the elastic energy associated with deforming the cortex. The linear term in cell perimeter, −2KPP0Pi, represents the effective line tension in the cortex and gives rise to a “preferred perimeter” P0. The value of P0 can be decreased by up-regulating the contractile tension in the cortex (19, 24, 27), and it can be increased by up-regulating cell–cell adhesion. We simulate tissues containing N cells under periodic boundary conditions; the value of N has been varied from N=64 to 1600 to check for finite-size effects (SI Appendix, Fig. S6). A0 is set to be equal to the average area per cell, and A0 is used as the unit of length. After nondimensionalizing Eq. 1 by KAA02 as the unit energy scale, we choose KP/(KAA0)=1 such that the perimeter and area terms contribute equally to the cell shapes. The choice of KP does not affect the results presented. The preferred cell perimeter is rescaled p0=P0/A0 and varied between 3.7 (corresponding to the perimeter of a regular hexagon with unit area) and 4.6 (corresponding to the perimeter of an equilateral triangle with unit area) (27). We obtain disordered ground states of the SPV model by minimizing E, using the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm (28) starting from a random Poisson point pattern. We also test the finite temperature behavior of the SPV model by performing Brownian dynamics (21).

The ground states of Eq. 1 are amorphous tilings where the cells have approximately equal area but varying perimeters as dictated by the preferred cell perimeter p0. It has been shown that at a critical value of p0*≈3.81, the tissue collectively undergoes a solid–fluid transition (27). When p0<p0*, cells must overcome finite energy barriers to rearrange and the tissue behaves as a solid, while above p0*, the tissue becomes a fluid with a vanishing shear modulus as well as vanishing energy barriers for rearrangements (27). Coupled to these mechanical changes, there is a clear signature in cell shapes at the transition (21, 29, 30): The shape-based order parameter calculated by averaging the observed cell perimeter-to-area ratio s=⟨P/A⟩ grows linearly with p0 when p0>p0* in the fluid phase but remains at a constant (s≈p0*) in the solid phase. In Fig. 1A, we show three representative snapshots of the ground states at various values of p0. We take advantage of the diversity and tunability of the point patterns and cell structures (Fig. 1A) produced by the SPV model and use them as templates to engineer photonic materials.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Tissue structure in the SPV model. (A) Simulation snapshots at three different values of the preferred cell perimeters p0. Cell centers are indicated by points, and cell shapes are given by their Voronoi tessellation (red outlines). (B) Contour plot of the structure factor S(qx,qy) corresponding to the states shown in A. Scale bar has length 2π/D in reciprocal space, where D is the average spacing between cell centers. (C) Pair-correlation function g(r) at different values of p0. (D) Structure factor S(q) at different values of p0.

Characterization of 2D Structure.

To better understand the ground states of the SPV model, we first probe short-range order by analyzing the pair-correlation function g(r) (see details in Materials and Methods) of cell centers. In Fig. 1C, g(r) for the solid phase (p0<3.81) shows mutual exclusion between nearest neighbors and becomes constant at large distances. These features are similar to those observed in other amorphous materials with short-range repulsion and a lack of long-range positional order, such as jammed granular packings (31) and dense colloidal arrays (32). When p0 is increased, the tissue enters into the fluid phase at p0>3.81. To satisfy the higher preferred perimeters, cells must become more elongated. And when two neighboring elongated cells align locally, their centers can be near each other, whereas cells not aligned will have centers that are further apart. As a result, the first peak of g(r) broadens. Hence, as p0 is increased, short-range order is reduced. Deeper in the fluid regime, when p0>4.2, g(r) starts to develop a peak close to r=0, which means that cell centers can come arbitrarily close. Interestingly, the loss of short-range order does not coincide with the solid–fluid transition, which yields an intermediate fluid state that retains short-range order.

Next we focus on the structural order at long length scales. While the SPV ground states are aperiodic by construction, they show interesting long-range density correlations. The structure factor S(q) (see details in Materials and Methods) is plotted for various p0 values in Fig. 1 B and D. Strikingly, for all p0 values tested, the structure factor vanishes as q→0, corresponding to a persistent density correlation at long distances. This type of “hidden” long-range order is characteristic of patterns that are known as hyperuniformity (33, 34). In real space, hyperuniformity is equivalent to when the variance of the number of points σR2 in an observation window of radius R grows as function of the surface area of the window—that is, σR2∝Rd−1, where d is the space dimension (33). This is in contrast to the σR2∝Rd scaling that holds for uncorrelated random patterns. Indeed, real space measurements in the SPV model also confirm that the distribution of cell centers is strongly hyperuniform for all values of p0 and system sizes tested that include both solid and fluid states (SI Appendix, Fig. S1).

It has been suggested that hyperuniform amorphous patterns can be used to design photonic materials that yield PBGs (9, 10). Florescu et al. (9) further conjectured that hyperuniformity is necessary for the creation of PBGs. However, recent work by Froufe-Pérez et al. (14) have demonstrated that short-range order rather than hyperuniformity may be more important for PBGs. The SPV model provides a unique example of a hyperuniform point pattern with a short-range order that can be turned on and off. This tunability will allow a direct test of the ideas proposed in refs. 9, 10, and 14.

2D Photonic Material Design and Properties.

For any point pattern (crystalline or amorphous), the first step in the engineering of PBGs is to decorate it with a high dielectric contrast material. The simplest protocol is to place cylinders centered at each point {ri=(xi,yi)} that are infinitely tall in the z−direction. Such design typically yields bandgaps in the transverse magnetic (TM) polarization (the magnetic field is parallel to the xy-plane) (35). Based on this design, we first construct a material using SPV point patterns. To maximize the size of the bandgap, the cylinders are endowed with dielectric constant ϵ=11.56 and radius r/D=0.189 (14). We will refer to this construction as “TM-optimized.” We also use a second decoration method (9, 14) that has been shown to yield complete PBGs—that is, gaps in both TM and Transverse Electric (TE) polarizations. We use a design based on the Delaunay triangulation of a point pattern. Cylinders with ϵ=11.56 and radius r/D=0.18 are placed at the nodes of the Delaunay triangulation while walls with ϵ=11.56 and thickness w/D=0.05 are placed on the bonds of this trivalent network. We refer to this construction as “TM + TE-optimized.” Photonic properties are numerically calculated using the plane wave expansion method (36) implemented in the MIT Photonic Bands program. We use the supercell approximation in which a finite sample of N cells is repeated periodically. The photonic band structure is calculated by following the path of k∥=(0,0)→(12,12)→(−13,13)→(0,0) in reciprocal space. SI Appendix includes a sample script used for photonic band calculations.

The TM-optimized band structure based on a SPV ground state with N=64 cells and p0=3.85 is shown in Fig. 2A. Due to the aperiodic nature of the structure, the PBG is isotropic in k∥. We also calculate the ODOS (Fig. 2B) by binning eigenfrequencies from 10 samples with the same p0 but different initial seeds. The relative size of the PBG can be characterized by the gap–midgap ratio Δω/ω0, which is plotted as a function of p0 in Fig. 3A for both TM- and TM + TE-optimized structures. We find that the size of the PBG is constant in the solid phase of the SPV model (p0<p0*) with width Δω/ω0≈0.36 for the TM-optimized structure and Δω/ω0≈0.1 in the TM + TE-optimized structure. In the fluid phase, Δω/ω0 decreases as p0 increases yet stays finite in the range of 3.81<p0≲4.0. At even larger p0 values, the PBG vanishes. The location at which the PBG vanishes appears to coincide with the loss of short-range order in the structure. To quantify this, we plot Δω/ω0 versus the first peak height of the pair-correlation g1 in Fig. 3C. This is in agreement with the findings of Yang et al. (11) and suggests that short-range positional order is essential to obtaining a PBG that allows for collective Bragg backscattering of the dielectric material. This also shows that PBGs are absent in states that are hyperuniform but missing short-range order.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Characterizing photonic properties in 2D. (A) Photonic band structure of Transverse Magnetic (TM) for the material constructed by placing dielectric cylinders at the cell centers exhibits a PBG. Design is based on a ground state of the SPV at p0=3.85. k∥ is the in-plane wave vector. (B) The optical density of states (ODOS) of TM at different values of p0. The width of the bandgap Δω has a strong dependence on p0, while the midgap frequency ω0 remains constant (SI Appendix, Fig. S4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Tunability of photonic properties in 2D. (A) The structure of the PBG as a function of p0. The TM-optimized structure corresponds to the size of the TM bandgap for a material constructed by placing dielectric cylinders at cell centers. The TM + TE-Optimized structure corresponds to complete PBGs for a material constructed using a trivalent network design. Phases are colored according to the mechanical property of the material. At p0≈3.81, the material undergoes a solid–fluid transition where the shear modulus vanishes. The dependence of ω0 on p0 is weak as shown in SI Appendix, Fig. S4A. (B) Effect of heating on the PBG in the TM-optimized structure. At fixed p0=3.7, temperature T is gradually increased. At T≈0.25, the material begins to fluidize through melting, while the bandgap still persists into the mechanical fluid phase. (C) The relationship between the bandgap size in the TM-optimized structure and the short-range order of the system. The explicit dependence of short-range order on p0 and T is shown in SI Appendix, Figs. S2 and S3.

To test for finite-size dependence of the ODOS and the PBG, we carry out the photonic band calculations at p0=3.7 for various system sizes ranging from N=64 to N=1600. At each system size, the ODOS was generated by tabulating TM frequencies along the k∥=(0,0) and k∥=(0.5,0.5) directions in reciprocal space for 10 randomly generated states. As shown in SI Appendix, Fig. S6, while low-frequency modes may depend on the system size and bigger fluctuations exist for smaller systems, the PBG is always located between mode numbers N and N+1 and has a width that does not change as a function of N.

Since the model behaves as a fluid-like state above p0*≈3.81 (21, 27) and the PBG does not vanish in the fluid-like regime (Fig. 3A), this gives rise to a photonic fluid where a PBG can exist without a static and rigid structure. To show this explicitly, we test the effects of fluid-like cell rearrangements by analyzing the photonic properties of the dynamical fluid phase at finite temperature. At a fixed value of p0=3.7, we simulate Brownian dynamics in the SPV model at different temperatures T. At each T, we take 10 steady-state samples and construct the TM-optimized structures to calculate their ODOS. In Fig. 3B, we plot the bandgap size as a function of increasing temperature. Note that even past the melting temperature of T≈0.25 (21), the PBG does not vanish, again giving rise to a robust photonic fluid phase. Inside the fluid phase, we also find that the PBG is not affected by the positional changes due to cell rearrangements (SI Appendix, Fig. S5). Finally, we analyze the effect of heating on the short-range order. In Fig. 3C, increasing T results in another “path” in manipulating the short-range order.

Extension to 3D Photonic Material.

To further demonstrate the viability and versatility of tissue-inspired structures as design templates for photonic materials, we extend this study to 3D. Recently, Merkel and Manning (37, 38) generalized the 2D SPV model to simulate cell shapes in 3D tissue aggregates by replacing the cell area and perimeter with the cell volume and surface area, respectively. This results in a quadratic energy functional that is a direct analog of Eq. 1:E=∑i=1NKS(Si−S0)2+KV(Vi−V0)2,[2]where Si and Vi are the surface area and volume of the i-th cell in 3D, with the preferred surface area and volume being S0 and V0, respectively. Similar to the 2D version of the model, we have two moduli—KS for the surface area and KV for volume. Following ref. 37, we make our model dimensionless by setting V01/3 as the unit of length and KSV04/3 as the unit of energy. This gives us a dimensionless energy and dimensionless shape factor s0=S0/V02/3 as the single tunable parameter in 3D, which is the anaolog of the parameter p0 in the 2D scenario. For the 3D model, we perform Monte Carlo simulations (39) with the cell centers at a scaled temperature T=1, expressed in the unit KSV04/3/kB where kB is the Boltzmann constant. During the simulation, each randomly chosen cell center is given by a displacement following the Metropolis algorithm (40), where the Boltzmann factor is calculated using the energy function in Eq. 2. For this purpose, we extract the cell surface areas and volumes from 3D Voronoi tessellations generated using the Voro++ library (41). N such moves constitute a Monte Carlo step, and we perform 105 such steps to ensure that the cells have reached a steady state. Then we run for another 105 steps to compute the g(r) and S(q). The average volume per cell is held constant for this procedure. Periodic boundary conditions were applied in all directions. For all of the simulations, we keep the scaled modulus KVV02/3/KS=1. It was previously found (37) that the T=0 ground states undergo a similar solid-to-fluid transition at s0∼5.4—that is, solid for s0<5.4 and fluid when s0≥5.4. Our Monte Carlo simulations also recover this transition near s0=5.4 when we plot an effective diffusivity, extracted from the mean-square displacement (MSD) of the cells, at different s0 values (SI Appendix, Fig. S9). States below the transition (s0=5) and above the transition (s0=5.82) are shown in Fig. 4 for a system of N=100 cells.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Structure characterization in the 3D Voronoi cell model. (A) Cell centers and their corresponding Voronoi tessellation as well as the decorated 3D photonic structure are shown for a state at s0=5 and N=100 cells. The cell centers are drawn with a finite size and different colors only to aid visualization. (B) Cell centers and their corresponding Voronoi tessellation as well as the decorated 3D photonic structure are shown for a state at s0=5.82 and N=100 cells. (C) Pair-correlation function g(r) at different values of s0. (D) Structure factor S(q) at different values of s0.

We characterize the structure of the 3D cell model by calculating their pair-correlation function g(r). In Fig. 4C, the short-range order behaves similar to the 2D case. At low values of s0 that correspond to rigid solids, the first peak in g(r) is pronounced, suggesting an effective repulsion between nearest neighbors. As s0 is increased, this short-range order gets eroded as a preference for larger surface area-to-volume ratio allows nearest neighbors to be close. The s0=6.1 state in Fig. 4C is an extreme example where it is possible for two cell centers to be arbitrarily near. However, the long-range order remains throughout all values of s0 tested as shown by the S(q) in Fig. 4D. The small value of S(q→0) is a consequence of suppressed density fluctuations across the system when cells all prefer the same volume dictated by Eq. 2. Whether they are truly hyperuniform would require sampling at higher system sizes, which is beyond the scope of this study.

Next, to make a photonic material, we decorate the Voronoi tessellations to create a connected dielectric network. While a Voronoi tessellation is already a connected network made of vertices and edges, it possesses a large dispersion of edge lengths. This can result in vertices that are arbitrarily close to each other and could hinder the creation of PBGs (42). To overcome this, we adopt a method described in refs. 9, 15, and 43 to make the structure locally more uniform. In a 3D Voronoi tessellation, each vertex is calculated from the circumcenter of the four neighboring cell centers and edges are formed by connecting adjacent vertices. In this design protocol, the connectivity of the dielectric network is the same as the network of the vertices and edges in the Voronoi tessellation. However, the vertex positions of the dielectric network are replaced by the center-of-mass (barycenters) of the four neighboring Voronoi cell centers. The resulting structure is a tetrahedrally connected network where the edges are more uniform in length, two representative samples are shown in Fig. 5A. Next we decorate the network with dielectric rods of width W running along the edges. For the dielectric rods, we again use ϵ=11.56, and their width W is chosen such that the volume filling fraction of the network is Vrod/Vbox=20% (42).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Characterizing photonic properties in the 3D design. (A) The optical density of states is shown at various values of s0. Here frequency is in units of 2πc/L, where L is the average edge length in the photonic network. (B) The size of the gap–midgap ratio Δω/ω0 as a function of s0. (C) Midgap frequency ω0 as a function of s0.

The photonic properties of the 3D dielectric network are calculated using the MIT Photonic Bands program (36). We calculate the ODOS for structures based on different values of s0. Here we have chosen structures containing N=100 cells. Due to the isotropic nature of the photonic band structure (SI Appendix, Fig. S7), we generate the ODOS based on two reciprocal vectors k=(0,0,0) & k=(0.5,0.5,0). At each value of s0, we average over 10 different random samples. We find the first complete PBG between mode number nV and nV+1, where nV is the number of vertices. In the solid phase of the model (s0<5.4), we find an average gap–midgap ratio of ∼6%; this decreases as s0 is increased and becomes vanishingly small when the model is in its fluid phase (s0>5.4) (Fig. 5B). Interestingly, the midgap frequency ω0 also shifts slightly toward lower values with increasing s0 (Fig. 5C). The appearance of the PBG here also coincides with the presence of short-range order. The PBG vanishes when there is no longer a pronounced first peak in the g(r).