We use statistical mechanics to model and characterize aggregates of cells and their dynamics.

Self-Similar Dynamics of Nuclear Packing in the Early Drosophila Embryo

## Self-Similar Dynamics of Nuclear Packing in the Early Drosophila Embryo

Embryonic development starts with cleavages, a rapid sequence of reductive divisions that result in an exponential increase of cell number without changing the overall size of the embryo. In Drosophila, the final four rounds of cleavages occur at the surface of the embryo and give rise to $$\sim6000$$ nuclei under a common plasma membrane. In Ref. [1], we use live imaging to study the dynamics of this process and to characterize the emergent nuclear packing in this system. We show that the characteristic length scale of the internuclear interaction scales with the density, which allows the densifying embryo to sustain the level of structural order at progressively smaller length scales. This is different from nonliving materials, which typically undergo disorder-order transition upon compression. To explain this dynamics, we use a particle-based model that accounts for density-dependent nuclear interactions and synchronous divisions (see Videos 1 and 2). We reproduce the pair statistics of the disordered packings observed in embryos and recover the scaling relation between the characteristic length scale and the density both in real and reciprocal space. This result reveals how the embryo can robustly preserve the nuclear-packing structure while being densified. In addition to providing quantitative description of self-similar dynamics of nuclear packings, this model generates dynamic meshes for the computational analysis of pattern formation and tissue morphogenesis.

Video 1. Simulations of dynamics of packing of nuclei by mitotic division on the surface of a prolate spheroid, similar to early fruit-fly (Drosophila) embryos. A density-adaptive force-field model ensures the packings are self-similar from cycle to cycle as observed in real embryo.

Video 2. Dynamics of packing of spheres by $$10$$ rounds of mitotic divisions starting from a single particle in a periodic box. The final packing fraction is $$0.63$$.

1. S. Dutta, N. J.-V. Djabrayan, S. Torquato, S. Y. Shvartsman, and M. Krajnc, Self-Similar Dynamics of Nuclear Packing in the Early Drosophila Embryo, Biophysical Journal, 117(4) 743 (2019).
Disordered Multihyperuniformity in Avian Photoreceptor Cell Patterns and Models

## Disordered Multihyperuniformity in Avian Photoreceptor Cell Patterns and Models

The evolution of animal eyes has been an intense subject of research since Darwin. The purpose of a visual system is to sample light in such a way as to provide an animal with actionable knowledge of its surroundings that will permit it to survive and reproduce. Cone photoreceptor cells in the retina are responsible for detecting colors and they are often spatially arranged in a regular array (e.g., insects, some fish and reptiles), which is often a superior arrangement to sample light. In the absence of any other constraints, classical sampling theory tells us that the triangular lattice (i.e., a hexagonal array) is the best arrangement.

Diurnal birds have one of the most sophisticated cone visual systems of any vertebrate, consisting of four types of single cone (violet, blue, green and red) which mediate color vision and double cones involved in luminance detection; see Figure 1. Given the utility of the perfect triangular-lattice arrangement of photoreceptors for vision, the presence of disorder in the spatial arrangement of avian cone patterns was puzzling.

Our investigation in collaboration with Joseph Corbo at Washington University presents a stunning example of how fundamental physical principles can constrain and limit optimization in a biological system [1]. By analyzing the chicken cone photoreceptor system consisting of five different cell types using a variety of sensitive microstructural descriptors, we found that the disordered photoreceptor patterns are “hyperuniform” (as defined above), a property that had heretofore been identified in a unique subset of physical systems, but had never been observed in any living organism.

Disordered hyperuniform structures can be regarded as a new exotic state of disordered matter in that they behave more like crystals or quasicrystals in the manner in which they suppress density fluctuations on large length scales, and yet are also like liquids and glasses in that they are statistically isotropic structures with no Bragg peaks. Thus, hyperuniform disordered materials can be regarded to possess a “hidden order” that is not apparent on short length scales.

Remarkably, the patterns of both the total population and the individual cell types are simultaneously hyperuniform, which has never been observed in any system before, physical or not. We term such patterns “multi-hyperuniform” because multiple distinct subsets of the overall point pattern are themselves hyperuniform. We devised a unique multiscale cell packing model in two dimensions that suggests that photoreceptor types interact with both short- and long-ranged repulsive forces and that the resultant competition between the types gives rise to the aforementioned singular spatial features characterizing the system, including multi-hyperuniformity (see Figure 1). These findings suggest that a disordered hyperuniform pattern may represent the most uniform sampling arrangement attainable in the avian system, given intrinsic packing constraints within the photoreceptor epithelium. In addition, they show how fundamental physical constraints can change the course of a biological optimization process. Our results suggest that multi-hyperuniform disordered structures have implications for the design of materials with novel physical properties and therefore may represent a fruitful area for future research. For more information, see Ref. [1].

Figure 1: (left) Experimentally obtained configurations representing the spatial arrangements of centers of the chicken cone photoreceptors (violet, blue, green, red and black cones). (right) Simulated point configurations representing the spatial arrangements of chicken cone photoreceptors. The photoreceptor types interact with both short- and long-ranged repulsive forces such that the resultant competition between the types gives rise to the aforementioned singular spatial features characterizing the system, including multi-hyperuniformity. The simulated patterns for individual photoreceptor species are virtually indistinguishable from the actual patterns obtained from experimental measurements. Figures taken from Ref. [1].

We have developed two additional theoretical models of these fascinating disordered multihyperuniform systems. In Ref. [2], we present a statistical-mechanical model that rigorously achieves disordered multihyperunifom many-body systems by tuning interactions in binary mixtures of nonadditive plasmas. This model provides a powerful method to generate disordered multihyperuniform systems via equilibrium mixtures and will facilitate the future study of their potentially unique photonic, phononic, electronic, and transport properties.

In Ref. [3], we further model the avian retina via an equimolar three-component mixture (one component to sample each primary color: red, green, and blue) of nonadditive hard disks to which a long-range logarithmic repulsion is added between like particles. Notably, this model captures the salient features of the spatial distribution of photoreceptors in avian retina; namely, the presence of disorder, multihyperuniformity, and local heterocoordination. The latter feature being critical for the efficient sampling of light.

1. Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato, Avian Photoreceptor Patterns Represent a Disordered Hyperuniform Solution to a Multiscale Packing Problem, Physical Review E, 89, 022721 (2014).
2. E. Lomba, J.-J. Weis, and S. Torquato, Disordered Multihyperuniformity Derived from Binary Plasmas, Physical Review E, 97, 010102(R) (2018).
3. E. Lomba, J.-J. Weis, L. Guisández, and S. Torquato, Minimal Statistical-Mechanical Model for Multihyperuniform Patterns in Avian Retina, Physical Review E, 102 012134 (2020).
Interdisciplinary Development of a Computational Tumor Model

## Interdisciplinary Development of a Computational Tumor Model

This work began as a result of a small grant from the National Institutes of Health Grant CA 84509 which was a joint grant between Harvard Medical School and Princeton University.

Currently, the dynamics of malignant brain tumor growth is a medical mystery. The incidence of primary malignant brain tumors is already 8/100,000 persons per year and is still increasing. The vast majority consists of high-grade malignant tumors such as glioblastoma multiforme (GBM) (see Figure 1). In spite of aggressive conventional and advanced treatments, the prognosis remains uniformly fatal with a median survival time for patients with GBM of 8 months.

The rapid growth and resilience of tumors make it difficult to believe that they behave as random, disorganized and diffuse cell masses and suggest that they are emerging opportunistic systems, that not only adapt to their environment but also change their environment for survival purposes. If this hypothesis holds true, a growing tumor must be investigated and treated as a self-organizing complex dynamical system. This cannot be done with currently available in vitro/in vivo models or common modeling approaches.

Figure 1. (left) An MRI scan of a brain. The tumor is located in the upper left portion of the image. The white rim around the tumor is composed of highly active cells, corresponding to the red rim in the right panel. (right) The outer, red rim corresponds to cells that are dividing rapidly. The middle, yellow region are cells that are alive but are not dividing. The innermost, black region are cells that are necrotic (dead).

The aim of the project is to show that we can model the growth (proliferation and invasion) of brain tumors using concepts from statistical physics, materials science and dynamical systems, as well as data from novel oncological experiments. We expect the work not only to increase our understanding of tumorigenesis, but to provide insight into novel ways to treat malignant brain tumors.

In Ref. [1], we develop a novel and versatile three-dimensional cellular automaton model of brain tumor growth. We show that macroscopic tumor behavior can be realistically modeled using microscopic parameters. Using only four parameters, this model simulates Gompertzian growth for a tumor growing over nearly three orders of magnitude in radius. It also predicts the composition and dynamics of the tumor at selected time points in agreement with medical literature. We also demonstrate the flexibility of the model by showing the emergence, and eventual dominance, of a second tumor clone with a different genotype. The model incorporates several important and novel features, both in the rules governing the model and in the underlying structure of the model. Among these are a new definition of how to model proliferative and non-proliferative cells, an isotropic lattice, and an adaptive grid lattice.

Malignant brain tumors consist of a number of distinct subclonal populations. Each of these subpopulations may be characterized by its own behaviors and properties. These subpopulations arise from the constant genetic and epigenetic alteration of existing cells in the rapidly growing tumor. However, since each single cell mutation only leads to a small number of offspring initially, very few newly arisen subpopulations survive more than a short time. In Ref. [2], we quantify “emergence,” i.e., the likelihood of an isolated subpopulation surviving for an extended period of time. Our model only considers competition between clones; no cooperative effects are considered. We find that the probability that a subpopulation emerges under these conditions to be a sigmoidal function of the degree of change in cell division rates. This function has a non-zero value for mutations that confer no advantage in growth-rate, which represents the emergence of a distinct subpopulation with an advantage that has yet to be selected for, such as hypoxia tolerance or treatment resistance. We also observe a logarithmic dependence on the size of the mutated population in time. A significant probability of emergence is found for subpopulations with any growth advantage that comprise even 0.1% of the proliferative cells in a tumor. The impact of even two clonal populations within a tumor is shown to be sufficient such that a prognosis based on the assumption of a monoclonal tumor can be markedly inaccurate.

In Ref. [3], we extend the cellular automaton model to study the effects of treatment. By varying three treatment parameters, we simulate tumors that display clinically plausible survival times. Much of this work is dedicated to heterogeneous tumors with both treatment-sensitive and treatment-resistant cells. First, we investigate two-strain systems in which resistant cells are initialized within predominantly sensitive tumors. We find that when resistant cells are not confined to a particular location, they compete more effectively with the sensitive population. Moreover, in this case, the fraction of resistant cells within the tumor is a less important indicator of patient prognosis. In additional simulations, we study tumors that are initially monoclonal and treatment-sensitive, but that undergo resistance-mutations in response to treatment. Here, the tumors with both very frequent and very infrequent mutations develop with more spherical geometries. Tumors with intermediate mutational responses exhibited multi-lobed geometries, as mutant strains develop at localized points on the tumor’s surfaces.

1. A. R. Kansal, S. Torquato, G. R. Harsh, E. A. Chiocca, and T. S. Deisboeck, Simulated Brain Tumor Growth using a Three-Dimensional Cellular Automaton, Journal of Theoretical Biology, 203, 367 (2000).
2. A. R. Kansal, S. Torquato, E. A. Chiocca, and T. S. Deisboeck, Emergence of a Subpopulation in a Computational Model of Tumor Growth, Journal of Theoretical Biology, 207, 431 (2000).
3. J. E. Schmitz, A. R. Kansal, and S. Torquato, A Cellular Automaton Model of Brain Tumor Treatment and Resistance, Journal of Theoretical Medicine, 4, 223 (2002).
Modeling Heterogeneous Tumor Growth in Confined Microenvironments

## Modeling Heterogeneous Tumor Growth in Confined Microenvironments

We generalize the cellular automaton model to study the effects of vasculature evolution on early brain tumor growth as well as proliferative tumor growth in confined microenvironments in Refs. [1] and [2], respectively. In Ref. [3], we formulate a cellular automaton model of tumor growth that accounts for several different inter-tumor processes and host-tumor interactions. In particular, the algorithm couples the remodeling of the microvasculature with the evolution of the tumor mass and considers the impact that organ-imposed physical confinement and environmental heterogeneity have on tumor size and shape. Furthermore, the algorithm is able to account for cell-level heterogeneity, allowing us to explore the likelihood that different advantageous and deleterious mutations survive in the tumor cell population. This computational tool we have built has a number of applications in its current form in both predicting tumor growth and predicting response to treatment. Moreover, the latent power of our algorithm is that it also suggests other tumor-related processes that need to be accounted for and calls for the conduction of new experiments to validate the model’s predictions.

Figure 1: Simulated growing tumors in homogeneous extracellular matrix. (a) Noninvasive tumor cells. (b) Invasive tumor with cellular motility. (c) Invasive tumor with greater cellular motility. (d) Invasive tumor with even greater cellular motility. Necrotic cells are black, quiescent cells are yellow, proliferative cells are red, invasive tumor cells are green, and the degraded extracellular matrix is blue. Figure taken from Ref. [4].

Understanding tumor invasion and metastasis is of crucial importance for both fundamental cancer research and clinical practice. in vitro experiments have established that the invasive growth of malignant tumors is characterized by the dendritic invasive branches composed of chains of tumor cells emanating from the primary tumor mass. The preponderance of previous tumor simulations focused on non-invasive (or proliferative) growth. The formation of the invasive cell chains and their interactions with the primary tumor mass and host microenvironment are not well understood. In Ref. [4], we present a novel cellular automaton model that enables one to efficiently simulate invasive tumor growth in a heterogeneous host microenvironment. By taking into account a variety of microscopic-scale tumor-host interactions, including the short-range mechanical interactions between tumor cells and tumor stroma, degradation of the extracellular matrix by the invasive cells and oxygen/nutrient gradient driven cell motions, our cellular automaton model predicts a rich spectrum of growth dynamics and emergent behaviors of invasive tumors. We also predict nontrivial coupling between the growth dynamics of the primary tumor mass and the invasive cells. The capability of our cellular automaton model suggests that sophisticated in silico tools could eventually be utilized in clinical situations to predict neoplastic progression and propose individualized optimal treatment strategies. In Ref. [5], we generalize our cellular automaton model to account for the mechanical pressure a confined solid tumor exerts on its surrounding microenvironment. For a comprehensive review of our work on computational modeling of tumors, see Ref. [6]. In Ref. [7], we further generalize our cellular automaton model to incorporate a variety of cell-level tumor-host interactions and different mechanisms for tumor dormancy, e.g., the effects of the immune system.

1. J. L. Gevertz and S. Torquato, Modeling the Effects of Vasculature Evolution on Early Brain Tumor Growth,Journal of Theoretical Biology 243, 517 (2006).
2. J. L. Gevertz, G. Gillies, and S. Torquato, Simulating Tumor Growth in Confined Heterogeneous Environments, Physical Biology, 5, 036010 (2008).
3. J. Gevertz and S. Torquato, Growing Heterogeneous Tumors in Silico, Physical Review E, 80, 051910 (2009).
4. Y. Jiao and S. Torquato, Emergent Behaviors From a Cellular Automaton Model for Invasive Tumor Growth in Heterogeneous Microenvironments, PLoS Computational Biology, 7, e1002314 (2011).
5. Y. Jiao and S. Torquato, Diversity of Dynamics and Morphologies of Invasive Solid Tumors, AIP Advances, 2, 011003 (2012).
6. S. Torquato, Toward an Ising Model of Cancer and Beyond, Physical Biology, 8, 015017 (2011).
7. D. Chen, Y. Jiao, and S. Torquato, A Cellular Automaton Model for Tumor Dormancy: Emergence of a Proliferative Switch, PLoS ONE, 9, e109934 (2014).
Structural Characterization and Modeling of Biological Systems

## A Novel Three-Phase Model of Brain Tissue Microstructure

In Ref. [1], we propose a novel biologically constrained three-phase model of the brain microstructure. Designing a realistic model is tantamount to a packing problem, and for this reason, a number of techniques from the theory of random heterogeneous materials can be brought to bear on this problem. Our analysis strongly suggests that previously developed two-phase models in which cells are packed in the extracellular space are insufficient representations of the brain microstructure. These models either do not preserve realistic geometric and topological features of brain tissue or preserve these properties while overestimating the brain’s effective diffusivity, an average measure of the underlying microstructure. In light of the highly connected nature of three-dimensional space, which limits the minimum diffusivity of biologically constrained two-phase models, we explore the previously proposed hypothesis that the extracellular matrix is an important factor that contributes to the diffusivity of brain tissue. Using accurate first-passage-time techniques, we support this hypothesis by showing that the incorporation of the extracellular matrix as the third phase of a biologically constrained model gives the reduction in the diffusion coefficient necessary for the three-phase model to be a valid representation of the brain microstructure.

## Spatial Organization and Correlations of Cell Nuclei in Brain Tumors

Accepting the hypothesis that cancers are self-organizing, opportunistic systems, it is crucial to understand the collective behavior of cancer cells in their tumorous heterogeneous environment. In Ref. [2], we ask the following basic question: Is this self-organization of tumor evolution reflected in the manner in which malignant cells are spatially distributed in their heterogeneous environment? We employ a variety of nontrivial statistical microstructural descriptors that arise in the theory of heterogeneous media to characterize the spatial distributions of the nuclei of both benign brain white matter cells and brain glioma cells as obtained from histological images. These descriptors, which include the pair correlation function, structure factor and various nearest neighbor functions, quantify how pairs of cell nuclei are correlated in space in various ways. We map the centroids of the cell nuclei into point distributions to show that while commonly used local spatial statistics (e.g., cell areas and number of neighboring cells) cannot clearly distinguish spatial correlations in distributions of normal and abnormal cell nuclei, their salient structural features are captured very well by the aforementioned microstructural descriptors. We show that the tumorous cells pack more densely than normal cells and exhibit stronger effective repulsions between any pair of cells. Moreover, we demonstrate that brain gliomas are organized in a collective way rather than randomly on intermediate and large length scales. The existence of nontrivial spatial correlations between the abnormal cells strongly supports the view that cancer is not an unorganized collection of malignant cells but rather a complex emergent integrated system.

## Structural Characterization and Statistical-Mechanical Model of Epidermal Patterns

In proliferating epithelia of mammalian skin, cells of irregular polygon-like shapes pack into complex, nearly flat two-dimensional structures that are pliable to deformations. In this work, we employ various sensitive correlation functions to quantitatively characterize structural features of evolving packings of epithelial cells across length scales in mouse skin. We find that the pair statistics in direct space (correlation function) and Fourier space (structure factor) of the cell centroids in the early stages of embryonic development show structural directional dependence (statistical anisotropy), which is a reflection of the fact that cells are stretched, which promotes uniaxial growth along the epithelial plane. In the late stages, the patterns tend toward statistically isotropic states, as cells attain global polarization and epidermal growth shifts to produce the skin’s outer stratified layers. In Ref. [3], we construct a minimalist four-component statistical-mechanical model involving effective isotropic pair interactions consisting of hard-core repulsion and extra short-range soft-core repulsion beyond the hard core, whose length scale is roughly the same as the hard core. The model parameters are optimized to match the sample pair statistics in both direct and Fourier spaces. By doing this, the parameters are biologically constrained. In contrast with many vertex-based models, our statistical-mechanical model does not explicitly incorporate information about the cell shapes and interfacial energy between cells; nonetheless, our model predicts essentially the same polygonal shape distribution and size disparity of cells found in experiments, as measured by Voronoi statistics. Moreover, our simulated equilibrium liquid-like configurations are able to match other nontrivial unconstrained statistics, which is a testament to the power and novelty of the model. The array of structural descriptors that we deploy enable us to distinguish between normal, mechanically deformed, and pathological skin tissues. Our statistical-mechanical model enables one to generate tissue microstructure at will for further analysis. We also discuss ways in which our model might be extended to better understand morphogenesis (in particular the emergence of planar cell polarity), wound healing, and disease-progression processes in skin, and how it could be applied to the design of synthetic tissues.

1. J. L. Gevertz and S. Torquato, A Novel Three-Phase Model of Brain Tissue Microstructure, PLoS Computational Biology, 4, e1000152 (2008)
2. Y. Jiao, H. Berman, T-R. Kiehl and S. Torquato, Spatial Organization and Correlations of Cell Nuclei in Brain Tumors, PLoS ONE, 6, e27323 (2011).
3. D. Chen, W. Y. Aw, D. Devenport, and S. Torquato, Structural Characterization and Statistical Mechanical Model of Epidermal Patterns, Biophysical Journal, 111, 2534-2545 (2016).

Questions concerning this work should be directed to Professor Torquato.