University of Southern California Dissertations and Theses

Patterning and growth in biological systems: a computational exploration of biological robustness, epithelial growth dynamics and development

(USC Thesis Other)

Patterning and growth in biological systems: a computational exploration of biological robustness, epithelial growth dynamics and development

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content PATTERNING AND GROWTH IN BIOLOGICAL SYSTEMS: A COMPUTATIONAL EXPLORATION OF BIOLOGICAL ROBUSTNESS, EPITHELIAL GROWTH DYNAMICS AND DEVELOPMENT By George Courcoubetis A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) May 2021 Copyright 2021 George Courcoubetis ii To my family iii Acknowledgements I would love to thank, first and foremost, my advisor Prof. Stephan Haas for his mentorship and support throughout all of my doctoral study. Prof. Haas’s mentorship was invaluable when building intricate computational tools and approaching complex research questions. He was supportive during trying times and provided me with a solid foundation to build upon. The lessons learned are invaluable and will always be cherished. In addition, Prof. Haas’s collaborative nature allowed me to branch multiple collaborations and study a diverse number of systems. Prof. Sergey Nuzhdin played an integral role in introducing me to biology and challenging me with hypothesis and new ideas. Furthermore, I want to thank Prof. James Boedicker and Leonardo Morsut for introducing me to new projects and helping me develop further as scientist. My family has supported me throughout my academic career. I want to thank my mother Dora and my two brothers, Aristotelis and Leonardo, for always being there for me. I want to especially thank my father who nurtured my curiosity about how the world works from a young age and supported me throughout my doctoral study. Lastly, I want to thank my friends in LA and in Greece. Especially, I want to thank my colleagues Georgios and Vagelis for being great roommates and my support during the toughest times. Georgios and I shared many fruitful discussions iv often leading to helpful insights. Another thank you goes to my surfing and skateboarding group for the amazing adventures that kept me inspired. Finally, I want to thank my friends back home that always helped me remember the important things in life and to keep going. v Abstract This thesis is focused on computational modeling and quantitative analysis of patterning and growth in biological systems. Patterning is studied in the context of the Drosophila eye disc, an epithelial tissue that contains ~800 eye units arranged crystalline order, with an emphasis on photoreceptor differentiation. Coupled ODE equations are parametrically perturbed to determine and quantify the robustness of the developmental process. Generalization of the translational order parameter, utilizing the notion of probability distance, is used to quantify the transition of the pattern from order to disorder. The transition is found to exhibit a ubiquitous threshold response to noise, with important implications on cryptic variation and robustness of gene regulatory networks. Growth is studied in the context of the vertex model, in which a uniform growth of an epithelial monolayer and the Drosophila eye disc are simulated. During the growth process, novel critical behavior is identified, where the tissue grows in bursts. The magnitude of growth bursts, or avalanches, are shown to be distributed according to power law distributions and are found to match the critical behavior of stress avalanches in sheer induced amorphous materials. The avalanches are also found to be accompanied with collective cellular motion and vortex formation. The eye disc growth model dynamics are used to investigate correlations between chemical and signaling components. The onset of avalanches, signaling properties and global vi cellular packing are found to depend on the coordinated apical constriction present in the developmental process. vii List of Publications The work presented in this thesis contains material from the following publications and preprints: Courcoubetis G, Ali S, Nuzhdin SV, Marjoram P, Haas S. Threshold response to stochasticity in morphogenesis. PLOS ONE. 2019 Jan 30;14(1):e0210088. Courcoubetis G, Xu C, Nuzhdin SV, Haas S. Avalanches During Epithelial Tissue Growth; Uniform Growth and a Drosophila Eye Disc Model. bioRxiv. 2021 Mar 1;2021.03.01.433318. Other publications and preprints completed during the course of the Ph. D but not included in this thesis are: Courcoubetis G,. Gangan M., Lim S., Haas S., Boedicker J. Formation, collective motion, and merging of macroscopic bacterial aggregates (in preparation) Lam C, Morsut L. A Modular Computational Framework for the Design of Multicellular Genetic Circuits of Morphogenesis. bioRxiv. 2020 Jul 11;784496. (authorship added during rebuttal process in ACS Synthetic Biology) viii Table of Contents Acknowledgements ............................................................................................ iii Abstract ................................................................................................................ v List of Publications ............................................................................................ vii Introduction ......................................................................................................... 1 Chapter 1. Threshold response to stochasticity in morphogenesis ...................... 5 1.1 Introduction ..............................................................................................6 1.2 Methods..................................................................................................16 1.3 Quantitative measures of structural order ..............................................28 1.4 Results and interpretation ......................................................................38 1.5 Conclusions ............................................................................................54 1.6 Supporting information ..........................................................................57 Chapter 2. Avalanches During Epithelial Tissue Growth; Uniform Growth and a Drosophila Eye Disc Model .................................................................................. 62 2.1 Introduction ............................................................................................62 ix 2.2 Generic uniform growth model ..............................................................70 2.3 Drosophila eye disc growth model ........................................................77 2.4 Model Predictions ..................................................................................85 2.5 Discussion ..............................................................................................86 2.6 Computational Methods .........................................................................88 Conclusion ....................................................................................................... 104 References ....................................................................................................... 108 1 Introduction Biological organisms constitute the epitome of complexity. From the scale of proteins, lipids and DNA to the ensemble of cells to tissues and ultimately organisms, living systems exhibit remarkable emergent behavior. One valid question is, given the stochastic behavior of all the constituents, how are organisms able to reliably control emergence and take form during morphogenesis? In the first chapter, tools are developed to quantify the order of emerged patterns, which are tailored for general spatial configurations. The inherent robustness in a developmental model is quantified and tested up to the regime of disorder. This quantification uncovers a low-stochasticity regime where the final structure is unaltered but beyond a threshold, the structure becomes rapidly disordered. This sigmoid functional form, in response to stochasticity, points to cryptic generic variation and a prerequisite for morphogenetic process stability. The developmental system studied in Chapter 1 concerns the differentiation of R8 photoreceptor cells on the Drosophila eye disc. The Drosophila eye disc constitutes a unique biological tissue, as the ommatidia or eye units are positioned with crystal-like precision. In the eye disc, the R8 cells differentiate in the wake of a signaling and mechanical wave that propagates linearly though the tissue. To generate the crystal-like pattern, differentiated R8 cells excrete transcripts that 2 locally inhibit R8 differentiation in nearby cells while at the same time excreting long range signals that induce R8 cell differentiation. The gene regulatory network responsible for the formation of these highly ordered structures has been the focus of multiple studies. This has led to the refinement of the patterning mechanisms and the intricacies of the gene regulatory network at play. In addition to the components that generate the pattern, morphogens responsible for pattern stabilization have been introduced. Thus, it constitutes an ideal model to study the properties of robustness and disorder. To describe the order of the emergent patterns, novel quantitative measures are developed. Due to the small system size and the spatially correlated disorder present in the disordered patterns, traditional measures of order do not directly apply to the spatial patterns. Probability distance is a measure that quantifies the difference between two probability distributions. By fixing one of the two distributions on a reference distribution, and varying the target distribution, one can use probability distance as a measure of similarity. In this scenario, the probability distributions of nearest neighbor distances of the resulting point pattern are used. To quantify disorder, the reference distribution is taken from the ordered point pattern. Then, when stochasticity is introduced, the nearest neighbor distributions from the resulting point patterns are used as the target distributions. Thus, probability distance is used to describe the degree of disorder in the patterns. Interestingly, if 3 one chooses the uniform distribution as a reference, the translational order parameter is approached. Thus, the quantitative tool developed serves as a generalization of the translational order parameter, applicable to an any distribution of interest. Morphogenesis requires proper growth of cellular tissues, often concurrently with spatial patterning. Epithelial tissues, in tandem with mesenchymal tissues, are responsible for the formation of the majority of organs during embryogenesis. In the second chapter, an epithelial tissue model is developed to simulate growth and development. For large tissue size, the epithelium exhibits area discontinuities, or growth avalanches, which are a result of increased cell packing density. This is a novel emergent phenomenon, which has a dominant effect in tissue growth dynamics. The avalanches are accompanied with collective cell movement, with cells in the core moving the least and cell on the outer edges the most. The net radial motion of cells leads to the large-scale deformations and changes in cellular packing density. In addition, the cell trajectories are also found to form a vortex in the center of the tissue, a phenomenon that has been observed in vitro. Systems that display avalanche instabilities often exhibit self-organized criticality (SOC). SOC is a property present in dynamical systems, where slowly driven processes lead to a critical point without the need of tuning parameters. At criticality, the systems exhibit scale-invariance, with quantities distributed 4 according to power law distributions. Distinct processes, such as earthquakes and magnetic spin avalanches, that share the same exponent in their distributed quantities, belong in the same universality class. This correspondence is useful, as it allows for analysis and results from one system to translate to the other. In the context of epithelial growth, the avalanche magnitude probability density distributions obey a power law, a clear signature of self-organized criticality. In this case, the individual cell divisions slowly drive the system towards the critical point of avalanche formation. Furthermore, the exponent of the power law is shown match the exponent of the probability density distribution of stress avalanche magnitudes in sheered amorphous materials. The growth model was later enhanced with signaling, mechanical reaction of cells, and compartmental growth to simulate development during the third instar Drosophila eye disc. The Drosophila eye disc is an ideal system to study the effects of collective mechanical deformation on the dynamics of the growth process, as a signaling wave propagates through the tissue, accompanied with stern apical constriction. The coordinated mechanical deformations during this process are found to couple with the avalanche onset and to affect signaling propagation providing testable experimental predictions. The model was calibrated to match the cell number, growth rate and other signaling properties of the experimental system. 5 Chapter 1. Threshold response to stochasticity in morphogenesis 6 1.1 Introduction Deterministic outcomes from inherently stochastic components Biological systems are intrinsically noisy but nonetheless produce deterministic outcomes. During development, organisms utilize signaling molecules, i.e. morphogens, to generate a body plan and differentiate cells. With modern experimental techniques, it is possible to measure temporal concentrations of selected morphogens in each cell. The expression of a gene depends on the probabilistic outcomes of several factors, such as molecular binding affinities, processivity, and regulatory sequence interactions. Specifically, genetically identical cells produce morphogen transcripts in asynchronous bursts and in varying quantities (1). It is therefore essential to investigate how ordered deterministic structures are formed from underlying stochastic components, such as gene expression. In this study, as an approximation to the stochastic nature of morphogenesis, parametric variation is introduced in the partial differential equations of the mathematical model. Earlier studies have addressed robustness of developmental processes, also termed as canalization (2), which remain a subject of great interest today (3–6). Developmental robustness has been highlighted as a necessary condition that narrows down dramatically the search for plausible models 7 and one that can even “predict key mechanistic and molecular properties of the associated biochemical circuits” (3). For example, the Turing mechanism has been deemed inapplicable in many instances due to its sensitivity to noise common in developmental processes (6). In this chapter, we analyze the response of a developmental system to noise quantitatively from a statistical physics perspective. Recent studies in developmental biology have addressed the issue of robustness (7,8). For a model to be plausible, there are requirements that the suggested mechanism is robust when parametric noise is introduced. However, in previous works the definition of robustness has mostly been used in a qualitative manner, and the detailed properties of the disordering mechanisms have been left unexplored. In this contribution, novel quantitative measures of robustness are developed and, the value of their application is demonstrated in the analysis of pattern formation in finite-system simulations. These new tools enable the determination of the functional form of the model’s response to stochasticity, leading to a rigorous description of a noise threshold. Specifically, a low-noise regime of robustness is identified in which, even in the presence of stochastic input, there is a deterministic outcome, i.e. structural order. In addition, these new quantitative measures of structural order allow for a precise analysis of the behavior of the system after this threshold is crossed. This disordered regime is relevant to genetic mutations, environmental perturbations and cryptic genetic variation. 8 Pattern formation in the Drosophila eye disc is an ideal model system to study the effects of stochasticity on deterministic development The compound eye, found primarily in insects, consists of an array of repeating visual units. The Drosophila compound eye is made up of approximately 800 unit eyes, known as ommatidia. In wild-type Drosophila, the structure of the ommatidia resembles a near-perfect hexagonal lattice. This highly structured pattern is developed via the delicate coordination of cell signaling, proliferation, movement and apoptosis (9). Some of these cellular processes are guided through the communication of a few conserved molecules, known as morphogens, resulting in tissue morphogenesis. While the resulting functional eye emerges in the adult, the role of each cell within repeating ommatidial arrays, and other head structures, are specified during larval development (9). The larval eye-antenna imaginal disc, hereafter referred to as the eye disc, contains numerous cells that produce various morphogens. The Drosophila eye disc is one of the simplest systems for studying morphogenesis since this tissue develops in an approximately two-dimensional fashion. This is because the eye-antenna disc is derived from an epithelial monolayer (10). Furthermore, this is an ideal system to model morphogenesis, since cell lineage has 9 a minimal effect on pattern formation, and differentiation relies primarily on cell- to-cell signaling (11). The underlying developmental mechanisms that guide the formation of the eye disc from larva to adult have been studied extensively, allowing for the generalization and analysis of biologically realistic models. To investigate the stability of the developmental pathways in the Drosophila eye disc, numerical simulations are used, and the results are analyzed from a physics perspective. The spatial organization of the ommatidia begins with the specification of the R8 cells, the future centers of the ommatidia, on the undifferentiated cells of the eye disc, as shown in Fig 1. Mathematical modeling of the Drosophila eye disc pattern-formation mechanism has been the focus of previous investigations, (7,8) which proposed a model that reproduces the hexagonal lattice pattern of differentiated R8 cells, which is robust enough to be biologically plausible. Using this latest model as a basis, we examine and quantitatively test the robustness of the emerging spatially ordered patterns of differentiated R8 cells when transcriptional noise is introduced. Until today, only qualitative measures of spatial order of the pattern have been applied in this mathematical framework (7). To this end, we implement several appropriate measures of spatial order, testing the functional relationship of R8 cell pattern order with increasing stochasticity. 10 Fig 1. Experimental image of the positions of the R8 cells during morphogenesis. In this fluorescent microscopy image, one can see DE-Cadherin linked with GFP. This labels the top part of each cell shown as a ring. When multiple cells are in close proximity, they appear as one bright point, which in this case point the positions of the R8 cells. During the development of the eye disc, morphogens are produced and emitted by individual cells in quantities that are determined by morphogen inputs from other cells in their vicinity (9). In addition, the production and diffusion of morphogens is non-uniform, i.e. a given morphogen concentration in a cell does not result in a single production or diffusion rate. Therefore, the highly ordered structure of the ommatidia in the Drosophila eye disc, emerging from such stochastic constituents, requires robustness in the underlying gene regulatory networks. Our study aims to 11 quantify eye disc developmental robustness. We focus on one essential step out of many involved in eye disc development, in which additional robustness is introduced (12). However, the regulation of R8 cell distances is a crucial developmental step in precisely placing each repeating eye unit, thereby strongly influencing the organism’s overall visual acuity. Biological framework of R8 cell specification Each ommatidium is made of 20 cells, 8 of which are photoreceptors. The specification of the first photoreceptor, R8, guides the specification and orientation of the remainder of the cells within the ommatidium (9). Thus, the spacing between neighboring R8 cells is pivotal for positioning and refining geometric relationships in resulting ommatidia. Cellular defects that misplace a single ommatidium will influence the position of neighboring ommatidia, thus propagating further flaws in the lattice. If even a single gene involved in ommatidia formation is perturbed, all unit eyes can be affected, since the regulatory logic that generates each ommatidium is repeated using the same set of morphogens and gene networks (9). Prior to R8 cell specification, the eye disc is composed of tightly-packed undifferentiated cells dividing asynchronously. Differentiation starts with the initiation of the morphogenetic furrow (MF), a physical indentation in the eye disc 12 that sweeps through the tissue to dictate the pace at which ommatidia are specified and positioned. The MF advances anteriorly via the communication of several morphogens, primarily Hedgehog (Hh) and Decapentapegic (Dpp). The former morphogens have been identified as the long range diffusible activators for the production of atonal (ato), the transcription factor that initiates R8 cell differentiation (13,14). The MF is initiated in the posterior end of the eye disc and sweeps through the tissue towards the anterior end. As the MF sweeps through the eye disc, the cells posterior to the MF differentiate and commit to their respective role within the ommatidium. Anterior to the MF, the cells are asynchronously proliferating but at the anterior interface of the MF, cells are arrested in G1 to allow synchronous divisions (15). In the anterior part of the MF, the transcription factor atonal is expressed in intermediate levels in a stripe by all cells. As the MF progresses, the stripe of ato- expressing cells converges to clusters of 12 cells and finally, to a single ato- expressing cell, the R8 cell (9,16). The inhibition signals originating from the refining ato-expressing clusters define the template for the refinement of ato- expresing stripes to the subsequent column (14,16). Based on the results of (7), which deem the role of Scabrous as a short range diffusible activator, the diffusible inhibitor of stripe refinement has not yet been identified but is related to the 13 epidermal growth factor (14). Finally, the Notch-Delta pathway is known to be associated with both stripe and cluster refinement (14,17,18). Process leading to pattern formation in the Drosophila eye disc In this section, we provide a simplified description of the rationale behind the current mathematical model, describing the basic mechanism of pattern formation in the eye disc. The goal is to provide the reader with an intuitive understanding of the foundations of the mathematical model by introducing the minimum amount of elements and interactions. The details of the full model are elaborated in the methods section below. The cells of the eye disc are assumed to lie on a hexagonal grid, as appropriate for a hard disk tight packing. After the MF passes, a portion of these cells differentiate to become activated R8 cells, and the remainder will be left undifferentiated. At the end of the process, ato-expressing cells are positioned on a hexagonal lattice. This process is illustrated in Fig 2. 14 Fig 2. Visualization of the foundations of the R8 cell specification mechanism in the Drosophila eye disc. Simplified illustration of pattern formation mechanism in the Drosophila eye disc as a result of the competition between short-range inhibitor and long-range activator morphogens. Posterior region of the eye disc with an initial row of differentiated precursor cells, denoted by red dots. The morphogenetic furrow (MF), modeled by a plane wave front, moves to the right towards the anterior region. The gray circles represent the boundaries of regions affected by the short-range inhibitor, where the morphogenetic furrow (MF) will not initiate production of atonal. Therefore, differentiation can only occur in those locations which are not affected by the short-range inhibitor, leading to the hexagonal supper-lattice of differentiated R8 cells, which is experimentally observed. This diagram does not incorporate cluster formation. Some morphogens are diffusible, while others are cell-specific and stay within the boundaries of the cell. In this study, ato is the only non-diffusing, cell-specific morphogen, while all others can diffuse between cells. The various morphogens either have inhibitory or activating properties. The inhibitors are morphogens that, 15 when present in a cell, decrease the rate of production of another morphogen, and vice-versa for the activators. The following paragraph defines elements and their interactions in an approximate fashion in order to provide an understanding of the patterning process with minimum complexity. Initially, we will deal with the production of single R8 cells rather than ato expressing clusters, which are defined as cells that express ato after the MF has passed. There are two main ingredients that make the pattern of R8 cells form. First, there is an inductor (or activator) signal that represents the MF, which causes undifferentiated cells to produce atonal. It can be approximated as a rectangular wave form, with one edge expanding with a constant velocity that eventually moves through the entire hexagonal lattice. Second, there is an inhibition signal that blocks atonal when present above a threshold. When a cell expresses ato above a threshold it rapidly, compared to the speed of the MF, produces an inhibitory signal that represses the production rate of ato of all cells within its vicinity, depicted by the circular regions in Fig 2. Ultimately, the first cells that do not lie within the inhibitory regions, extending from the most anterior column of ato-expressing cells, will become R8 cells. As an initial condition, a periodic array of differentiated cells are placed as in the posterior region. The spacing of these cell clusters is chosen such that it will allow 16 for the pattern to propagate. In this simplified setting, the initial configuration is defined by differentiated cells separated evenly in a single row in the posterior- most region. The MF then starts to propagate from posterior to anterior and causes the pattern to propagate as shown in Fig 2. The actual model that leads to pattern formation is more elaborate. Now that there is a basic understanding of how the pattern propagates, some more intricate aspects can be introduced. In order for the model to be biologically plausible, cell clusters are formed instead of single cells. Even in this illustrative context, as it can be inferred schematically from Fig 2, the propagation of the pattern is very sensitive to the shape of the inhibitory regions. A recent refinement of the model eliminates this flaw by including a short ranged diffusible activator that is emitted rapidly, compared to the speed of the MF, by ato expressing cells. In simple terms, this forces the first uninhibited cell to receive the linearly propagating activator signal to create a circular activator region. This activator region defines the cluster size and shape and makes the model robust, preventing the creation of greatly elongated cell clusters, as explained in Ref. (7). 1.2 Methods Mathematical formulation 17 We start with a model of coupled differential equations that describe morphogenetic pattern formation for R8 cell specification on the Drosophila eye disc obtained from (7), an article expanding on the work of (8). The model relies on the simplifying assumption that the number of cells in the eye disc is fixed (7,8). The cells in the eye disc are arranged in a hexagonal grid, treating it as an underlying two- dimensional structure. The lattice constants for the hexagonal lattice are set to 1, and the points of the lattice are defined by the relation , where (i, j) are integers and are the perpendicular unit vectors in the x and y directions respectively. Every cell, due to the hexagonal geometry of the lattice, has six adjacent neighbors. The resulting coupled partial differential equations are of the form: 18 These represent the spatio-temporal evolution of the morphogen concentrations a, u, s, and h, responsible for pattern formation in the drosophila eye disc. The morphogen interactions are summarized in Fig 3(a). Here, the upper-case indexes label each cell. The pro-neural transcription factor marking the center of the future ommatidia is a (for atonal), u represents all diffusible inhibitors, s denotes short-range diffusible activator (scabrous), and h (for hedgehog) describes the morphogenetic furrow (MF) which activates production of atonal along a propagating wave front. MF propagation is mathematically derived to have this functional form from the underlying morphogen differential equations (7). The constant v is the velocity of the 19 MF, Dh is the diffusivity, Ph the production rate, and τh is the reaction time scale of hedgehog Fig 3. Interactions between morphogens and the simulation of the R8 cell ordered pattern. a) In this diagram we show the interactions defining the regulatory network between the four morphogens in the mathematical model. h acts solely as an activator to a, a increases production of itself and leads to production of u and s, s activates a and lastly, u inhibits the effect of h and s. b) Ordered R8 cell pattern obtained by the mathematical model. The pattern was obtained numerically using the parameters used in (7) (shown in S1 Table). The mathematical model contains multiple parameters, operators and functions that play important roles in the dynamics of the system. First, the characteristic reaction time scale of each morphogen is given by corresponding τ’s on the left hand side of the equations. In addition, each morphogen’s differential equation contains a term with a −λ coefficient that incorporates the spontaneous decay of the morphogens. This term, combined with τ, determines the mean lifetime of the morphogen via the ratio . In addition, for the morphogens that diffuse, there are Laplacian operators associated with diffusion, Δ which are discretized at the cell 20 level (see S1 Appendix for exact expression). For cells on the boundary, reflective boundary conditions are used. Note that although only the equations for u and s have this diffusion term, the differential equations that lead to the functional form of h are also based on diffusion. Each Laplacian operator is multiplied by a term D, the value of which (when divided by the respective τ) sets the scale for how fast the morphogen diffuses. The Heaviside functions are essential, as they set the morphogen concentration thresholds that start, accelerate and stop production. They are a simplified version of the more biologically plausible sigmoid functions, and are defined by The Heaviside function θ(h i − h1) in Eq 1 is responsible for initiating production of atonal, a, in the target cell with a rate 𝐺 /𝜏 𝑎 , once the concentration of h crosses the threshold h1. This term is multiplied by (1 − θ(u i − u1)), a function that goes to zero when u i > u 1 . This term introduces the functionality of u as a direct inhibitor of atonal in the model. The last term in Eq 1 introduces sas an activator of atonal. Furthermore, the first term in Eq 1, Pa θ(a i − aa), sets a threshold for atonal at which its production becomes refractive to inhibitory signals. Finally, the two remaining differential equations for the morphogens u and s in Eqs 2 and 3 have similar 21 functional forms. Here, the Heaviside functions specify the threshold value for atonal in the target cell that is needed for initialization of production. In this model parameter sets have to be tuned for the experimentally observed ordered pattern to emerge. Furthermore, as the placement of the ato expressing cells can vary in cluster size and separation, we use a biologically plausible and stable parameter set, determined by comparison of numeric simulations to experimental imaging data of developing drosophila ommatidia previous stability analyses adopted from (7,8), shown in S1 Table. More specifically, the decay length scales and time scales for morphogen concentration levels used in the numerics agree with the known experimental data (8). In Fig 3(b), we show simulation results for this ideal parameter set after the MF has passed the entire cell grid. The mathematical model is simulated numerically on a hexagonal lattice of 44 vs 120 cells. The size and shape of the simulation grid does not reflect the actual size of the eye disc. The mathematical model evaluated on the 44 vs 120 simulation grid, while using the parameter set shown in S1 Table, results in the differentiation of 85 ommatidia in contrast with the ∼800 ommatidia present in the wild-type eye disc. Moreover, the simulation output contains 20 columns of cells with alternating four to five cells on each column. In contrast, the real eye disc contains 32-34 columns and the number of ommatidia per column varies along the eye disc to 22 provide its elliptical shape (9). Therefore, the custom size used in the simulations is acceptable to probe how noise propagates between columns while keeping the system size small enough for computational efficiency. Placement of initial clusters, edge effects and the ordered pattern. In order for the pattern to propagate, one must place preexisting clusters of cells with non-zero atonal concentration as an initial condition. Without proper placement, the pattern does not propagate periodically. In this model, a cell is considered to be part of an ato expressing cluster when it expresses a concentration of a above a threshold value aa. Prior work has not specified a method for determining the initial differentiated cell placement other than trying all possible configurations. To determine the spacing of the differentiated cells, an isolated differentiated cluster was simulated. After the cell cluster reached a steady state, the shape of the resulting inhibition region was recorded. Then, the spacing between clusters or single cells was determined by requiring an inhibition region that resembles Fig 2. Note that the size of the cluster is determined by a combination of model parameters and was obtained by trials. 23 In this study, the emerging ordered pattern of R8 cells is not a perfect hexagonal lattice. It is defined as a pattern that retains its structure over long length scales; see S1 Fig for the quantitative definition used. Due to the discreteness of the lattice, for a given radius of inhibition, it is not necessary that a placement of initial clusters exists such that the pattern propagates without slight displacements. In addition, the edges of the simulation grid limit the flawless specification of cell clusters. Specifically, as can be seen in Fig 2, the placement of each ato expressing cell cluster is determined by the diffusion of the inhibitor by two ato expressing clusters from the previous column. At the edge, only one cluster determines the shape of the inhibitory region leading to a different separation length scale. Finally, even in the biological system, the pattern is not periodic down to single cell positions and a mechanism that requires such fine tuning for a perfect placement of clusters would be unrealistic. Introducing noise. Starting from validated simulations of the noiseless case (7,8), we are now well positioned to study the effects of parametric disorder on pattern formation. The motivation for this is that even in quasi-identical cells, it is experimentally observed that genes are expressed at appreciably different levels due to stochasticity. Here we investigate how such variations affect the outcome of pattern formation in the 24 Drosophila eye disc. To introduce effects of disorder to the system, one or two parameters of interest are chosen at a time. The reason for this choice is to probe the system response when one of the many components start to fail as often encountered in genetic mutations and other stresses. The analysis for adding equal amount of noise to many parameters and adding variable amounts of noise for a pair of parameters was included for completeness in S2 Fig. For the chosen model parameters, their value is picked from a normal distribution centered at the mean value that produces the ordered pattern of differentiated R8 clusters in Fig 4(a). For each cell of the system a unique value for the parameter is drawn. They are then kept constant during the time evolution. The width and center of each Gaussian distribution is chosen to be the same for all cells. When more noise is introduced, the width of the distribution is increased accordingly. More specifically, in the numerical data discussed below the widths are tuned from 0% to 60% of the mean. As an example, introducing noise in the production rate of u corresponds to picking values from the probability distribution: 25 where μ = 〈Pu〉. Introducing noise in the diffusivity of the morphogens is a bit more intricate and is explained in S1 Appendix. The effect of introducing noise in the diffusivity of u and s is presented in Fig 4. The interpretation of adding noise to the diffusivities corresponds to cell size variation. In this chapter, σ refers to the standard deviation of the normal distribution, as shown for Pu in Eq 6, and therefore defines the degree of noise. In the unlikely instance that a negative value is drawn from the distribution, the absolute value is considered. Fig 4. Pattern formation simulation results in the Drosophila eye disc, the morphogenetic furrow moves from left to right. Cluster positions and sizes of ato expressing are shown for 26 increasing noise σ in the diffusivity of the morphogens D u and D s. (a)-(c) are the final patterns for σ/ μ = 0, σ/ μ = 30%, and σ/ μ = 40% respectively. Here, μ is the mean of the corresponding normal distribution and the value that is used to generate the perfect pattern. Cluster refinement. Until now, in the parameter regime deemed biologically plausible, the mathematical model in Eqs 1-4 produces patterns containing clusters. However, the actual developmental process includes an extra step that reduces the clusters to single activated R8 cells, which then become the center of the future ommatidia. This process is associated with the Notch-Delta pathway, whereby the cell that produces the most Delta inhibits its nearest neighbors (18). In order to reproduce this refinement process, a single cell for each cluster is kept as the R8 cell. The cell is chosen to be the most central cell in the most populated row in each cluster, an approximation to the non-trivial pathways and models identified in the literature (7). This cluster refinement process is illustrated graphically in Fig 5. This extra step, of cluster refinement, has been analytically shown to increase robustness of the pattern (19). 27 Fig 5. Visualization of the cluster refinement step. The code takes the ato expressing cell clusters, as shown in (a), as an input and determines the single R8 cell that will become the center of the future ommatidium, as shown in (b). This example is shown for a noise level of σ/ μ = 40%, applied to the diffusivity of the morphogens u and s. Numerical evaluation Numerical forward integration on a 120 × 44 hexagonal lattice with temporal step size 10 −2 was used to evaluate the spatio-temporal evolution of the coupled differential Eqs 1-4. To obtain sufficient data for statistical analysis, the code was parallelized and run on USC’s high-performance supercomputer cluster center. 28 Each computation was assigned to multiple processors at once, running for around 10 hours on each node. To obtain reliable error bars, at least 50 realizations of each noise level were simulated. The total eye disc size width was chosen such to eliminate edge effects, allowing sufficient space for clusters to form on the edges. Also, a relatively long eye disk length of 120 lattice sites was used is to allow investigation of error propagation of clusters as a function of position from posterior (P) to anterior (A). The simulations were terminated once the morphogenetic furrow had passed throughout all of the eye disc. 1.3 Quantitative measures of structural order While previous quantitative studies have concentrated on characterizing the disorder of the R8 and avian point patterns of experiments and creating minimal underlying models (20,21), in this work we focus on the robustness of the mathematical model and its structural characteristics in the presence of stochasticity. More specifically, in this study the order of the pattern is quantified when varying the noise of parameters that represent regulatory and other processes, distinguishing from order quantification when changing morphological parameters directly (20,21). Curiously, while varying morphological parameters directly, a 29 threshold response is also apparent (21). Furthermore, the order parameter that is developed and applied is distinct from qualitative measures applied before (7,8), and is more closely related with the measures introduced in quantitative studies (20,21). This analysis will provide insight into pathways with which developmental systems either cope with or succumb to stochasticity beyond a threshold, ultimately leading to malformation. We wish to analyze the emerging point pattern using appropriate measures of structural order, including translational, bond orientation and variance of nearest neighbor distributions. Details of how these measures are designed and applied are discussed in this section. Traditional measures for determining structural order have been developed in solid state physics (22–25). However, in the context of the spatio-temporal formation of patterns in biological developmental models there are additional issues to address. Foremost, the degree of structural order in the final pattern is generally not homogeneous in the presence of stochasticity. While, per initial conditions, the posterior point pattern starts off close to an ordered state, it may result in a strongly disordered anterior region at later stages of development, depending on the degree of stochasticity. This type of disorder is correlated between rows because of the manifestly Markovian mechanism that produces the patterns, i.e. the points on a 30 given row are specified based on the geometry of the points in the previous row (26). This observation renders traditional solid state measures used to evaluate the structure of homogeneous crystals insufficient (23). Furthermore, the type of correlated disorder exhibited in the point patterns of the eye disc is unlike random thermal displacements found in atoms or infrequent impurities. Instead, such correlated disorder causes Bragg’s law to become inapplicable, since it is based on the assumption that the underlying unperturbed lattice is periodic (22). Specifically, the model system is far away from the two instances where a corrected Bragg’s law for imperfect lattices could still be applied. The Debye-Waller approximation (27,28) applies only for uncorrelated deviations from a perfectly periodic lattice, and the available corrections for correlated deviations (29) apply only in the limit of small deviations and short-ranged Gaussian correlations. In order to study the structure of amorphous solids and general highly disordered materials, as in the case of developmental models, it is therefore necessary to resort to the pair distribution function as a local measure of structural order (23). The definition of the pair distribution function is given by 31 In the context of developmental models, this function describes the probability of finding an activated cell at position r, given that another activated cell is located at position rj. The pair distribution when treated as function of the scalar distance, as in the case of isotropic patterns (23), is denoted by g2(r) and is called the radial distribution function. In the following analysis, we will use a variation of the radial distribution function. Practically, in numerical simulations of eye disk development one can calculate the radial pair distribution function only for small distances. Each simulation of the pattern only contains around 40 activated cells. Moreover, the calculation of g2(r) is limited by the open boundary conditions, because the application of periodic boundary conditions to a single realization and the combination of multiple realizations is unphysical. Therefore, here we focus on a variation of the radial distribution function, that includes only nearest neighbors. The nearest neighbors were identified via the Voronoi tessellation method (30). Furthermore, we analyze a scalar measure of the spatial order of the patterns, the translational order parameter defined as 32 where ηc is a cutoff limited by the simulation size (24). This order parameter can be interpreted as the Kolmogorov probability distance (31) between the radial distribution function of the target pattern and the Poisson random point process. (The radial distribution function of a Poisson random point process is the uniform distribution, which in the normalization used in Eq 8 is equal to 1). The translational order parameter, T, is a general order metric used to describe systems independently of the underlying crystal structure. In this setting, we examine how fast the pattern deteriorates compared to the perfect pattern, so we replace the uniform probability distribution with that of the perfect pattern. In addition, since we only considered the radial distribution of nearest neighbors, the cutoff c in the integration in Eq 8 was replaced with the furthest nearest neighbor distance. Finally, we study the bond angle order parameter. Contrary to the translational measure, this order parameter evaluates the spatial orientation of vectors connecting the nearest neighbors of all points. It is defined by 33 (25), taking a value of 1 for a perfect hexagonal point pattern and 0 for a completely random pattern. Here Nn is the total number of nearest neighbors in the point pattern, the j’s sum over all activated cells, and the k’s sum over all nearest neighbors of a given reference point. Lastly, θjk is the angle of the vector connecting each point with its nearest neighbor with respect to a fixed axis. The discussion in this section connects and motivates the order measures used in this article with the tools available in the literature. It is important to note that equations Eqs 7 and 8 apply to point patterns in continuous space. In this article, the point patterns are defined by positions of the activated R8 cells which are constrained to lie on a hexagonal lattice. This means that the argument of the pair distribution function becomes discrete. Thus, we are going to consider discrete probability functions. This leads to the scalar pair distribution function to be defined as: where Nbins is the number of bins used, nj the frequency of activated cells at the jth bin and N is the total number of R8 cells. Similarly, the translational order parameter is defined as: 34 Probability distribution functions of nearest-neighbor distances and angles. To quantitatively determine the degree of disorder, a post-processing code is used to record the activated R8 cell positions in the final patterns and generate the probability distributions of nearest-neighbor angles and distances. Firstly, the R8 clusters are refined to single R8 cells by the cluster refinement method explained above. For every R8 cell, the nearest neighbors are identified via the Voronoi tessellation (30). Next, nearest-neighbor distances and angles of each R8 cell are calculated, and a filter is applied to correct for boundary effects for the nearest neighbor angles. Combining this information for all random realizations, the probability distribution functions of nearest-neighbor distances and angles are determined. Since the cells in the simulation are positioned on a hexagonal lattice, the nearest-neighbor distances and angles take discrete values. As a consequence, the probability distributions are also discrete. They are shown in Fig 6 for the ordered case. These types of probability distribution functions are used to quantify the order of the R8 35 cell point patterns by calculating the variance and the distance between probability distributions. The same approach can be taken when interpreting experimental data obtained from imaging. Fig 6. Nearest neighbor (a) distance and (b) angle probability distributions of R8 cell point pattern in the Drosophila eye disc in the absence of stochasticity. Distances are naturally binned by the available sites on the underlying hexagonal lattice, whereas angular positions are collected in bins of 0.05 radians. Even in the zero-noise case, the histograms contain multiple peaks in Fig 6. The histogram peaks for the nearest neighbor distances are a consequence of three effects. First, there are two length scales associated with the nearest neighbor distances in the pattern, one for the x-axis and one for the y-axis. Even for a perfect hexagonal lattice there are two lattice constants, one for each dimension. Second, cells in the edges of the simulation grid are placed in a different separation than the cells in the bulk. The two former points result to four peaks. Third, by the definition of nearest neighbor dictated by the Voronoi tessellation, the edge cells along the y- 36 direction boundaries are nearest neighbors of each other creating a fifth peak. The very short peaks that remain in the rightmost side of Fig 6(a) are a consequence of edge effects and the Voronoi tessellation in the x-direction boundary. Finally, slight blurring of each peak is a consequence of clusters forming one or two cells away than their counterparts. Variance as a measure of order. The nearest-neighbor angles and distances can be thought of as samples of an underlying distribution. If this distribution has non-vanishing higher-order moments, it is not trivial to produce a measure of disorder. However, in the simple case where the distribution can be approximated by a Gaussian, its variance, can be used as a reliable measure of the spread of the distribution. The variance used in this context is the sample variance of the nearest neighbor distances and angles. It is calculated from N values, xi, using where is the sample mean. 37 Probability distribution distance as a measure of order: Fidelity and Kolmogorov distance. Here we introduce two new and useful measures of order which best address the needs of this study: fidelity and Kolmogorov distance. Both of these rely on the concept of distance between probability distributions, and are inspired by the translational order parameter T defined in Eq 8. The concept of probability distance applied as an order measure is, to the best of our knowledge, novel. The most novel adjustment is to apply the measure of probability distance not only on the nearest neighbor distances distributions but also the nearest neighbor angles distributions. Note that given any set of probability distributions that include an ordered case, of any dimension and type of variable, the order measures introduced here can be applied to quantify deviation from the ordered state. To define these measures of order, we use the zero noise case probability distribution as a reference and calculate its distance from each of the noisy probability distributions. This choice is intuitive, since we are addressing the question of how disordered the pattern is relative to the perfect pattern. The results did not depend strongly on the reference distribution, see S3 Fig for the results when using Eq 11 where the reference distribution is the uniform distribution. The fidelity of two discrete probability distributions p(xi) and q(xi) is defined as 38 where the sum runs over all the bins of the discrete distributions (31). The fidelity between two distribution falls into the range 0 ≤ F(p(xi), q(xi)) ≤ 1, where 1 is attained only when p(xi) = q(xi) ∀ xi. The Kolmogorov distance, in this context, is defined as where the sum again extends over all the bins of the discrete distributions (31). Akin to the fidelity, the Kolmogorov distance is bounded by zero and one. This measure represents the maximum deviation of the two probability distribution functions given that a collection of events occur. In order for the Kolmogorov distance to match the functional form of the fidelity, we will consider 1 − D(p(xi), q(xi)) instead. This quantity is 1 when the two probabilities are the same and 0 when they do not have a single common non-zero element. 1.4 Results and interpretation Threshold response 39 To quantify disorder in the emerging activated R8 cell patterns, the nearest neighbor distance and angle probably distribution functions were computed for various noise levels in the diffusion coefficients Du, Ds. (Later we will discuss the effect of stochasticity on other model parameters). Their histograms are shown in Fig 7. The binning used for the distance histograms reflects the discreteness of the underlying hexagonal lattice on which the cell centers are placed. For the distance diagrams, all possible nearest neighbor distances on the hexagonal lattice are used as bins. For the angle histograms, since there are many more possibilities, a constant bin size of 0.05 radians was chosen in order to appropriately resolve the distribution. 40 Fig 7. Nearest neighbor distance (left panels) and angle (right panels) probability distributions generated from the point-pattern in the Drosophila eye disc R8 cell specification for increasing noise levels. Stochasticity was introduced in the differential equations by drawing the value of the diffusion coefficients D u and D s from a normal distribution with mean μ and standard deviation σ. The nearest neighbor distance histograms in (a) to (d) and nearest neighbor angle histograms in (e) to (h) correspond Gaussian stochasticity with standard deviations σ/ μ = 0%, 30%, 40%, 60%. While 41 the first two histograms are almost identical, as noise in the diffusivity of the morphogens increases, the peaks get smeared out and the histograms appear wider. The angle histograms, unlike the distance histograms, are skewed towards larger angles. This asymmetry makes the analysis of the distributions more complicated, as in this case the variance, paints an incomplete picture. As discussed below, the probability distance order measure eliminates this issue. The angle histograms are skewed as a direct consequence of the cluster refinement process. As seen in the plots of atonal patterns (Fig 4), the onset of disorder is signaled by elongated clusters of activated R8 cells. Since these elongated clusters are brought down to a single cell in the cluster refinement step, and they cast large inhibition radii, the number of nearest neighbors found for the elongated clusters is less than the six found in the perfect pattern. A lesser number of neighbors trivially leads to higher angles in this case. Both histograms, distances and angles, start out sparse for the case without any stochasticity in Du, with empty bins between peaks, indicative of the highly ordered repeating lattice structure produced by the simulations in this case. As stochasticity increases, the resulting point patterns become aperiodic, leading to denser histograms. Furthermore, the histograms become wider with increasing stochasticity. This is captured in the variance, shown in Fig 8. 42 Fig 8. Disorder in the R8 cell point pattern of the developing Drosophila eye disc as a function of stochasticity. Here, the variance of the distributions of (a) nearest neighbor distances and (b) nearest neighbor angles is used to quantify the amount of disorder. As in the previous figure, stochasticity is introduced in the diffusivity of the morphogens D u and D s, and the x-axis represents the noise level in terms of the standard deviation of the underlying Gaussian distribution from which this parameter is drawn. As the noise level in the simulation is increased beyond a threshold, the variance grows. The variance plots in Fig 8 illustrate the generic nature of robust pattern formation in stochastic developmental models. They show that for small levels of stochasticity (here σ/ μ < σc/ μ = 20%) the ordered pattern remains basically unaffected. However, beyond this threshold the variance starts increasing linearly, both in the angle and distance histograms. This threshold response is a central finding of this study. As discussed below, this is a universal phenomenon, largely independent of which model parameters experience stochasticity and which measures of order are used to assess the system response. It implies that these stochastic models capture 43 biological resilience against stochastic variation up to a certain point, resulting in a regime of deterministic outcomes in spite of stochastic input. Let us now discuss how this type of threshold behavior is also picked up by the more general measures of order given by Eqs 13 and 14. The results of applying these probability distance measures is shown in Fig 9. Note that saturation observed in the nearest angle distribution variance plot is apparent in both cases, and the responses of both observables (distances and angles) to stochasticity exhibit a universal sigmoidal functional form, implying threshold behavior. These measures reveal the same interesting quality as the variance: the pattern order exhibits an initial resistance to weak input stochasticity that eventually gives way to structural malformation at larger noise levels. The sigmoidal functional form confirms that the mechanism does exhibit robustness for σ < σc. The fact that this is observed in the Kolmogorov measure, the fidelity, the bond orientation and variance measures supports the notion that such sigmoidal, threshold response to stochasticity is a universal feature (see supporting information, S4 Fig, for bond orientation sigmoidal response). 44 Fig 9. Probability distance measures applied to analyze response to stochasticity in R8 cell point pattern formation. (a)-(b) are generated using the fidelity, F, and (c)-(d) using the Kolmogorov distance, K, for nearest neighbor distances and angles respectively. Stochasticity was introduced in the model through drawing the parameter D u and D s in the differential equations Eqs 1-4 from a normal distribution with varying standard deviation, σ, and mean value, μ. The simulation results were compiled, and nearest R8 cell neighbor angle and distance distributions were obtained. Kolmogorov and fidelity of the resulting distributions where computed with reference to the perfect pattern case. Here, a value of 1 corresponds to a perfectly ordered pattern, and a value of 0 corresponds to a completely irregular pattern. This figure illustrates that the functional form is independent of the measure used. Next, we turn to the question of universality with respect to how stochasticity is introduced. We observe that stochasticity in other model parameters also produces 45 threshold behavior, indicating that the observed threshold response is a generic feature of this class of evolutionary models. Specifically, in Fig 10 we show the response to increasing stochasticity levels in the production rate of atonal Pa, the production rate due to the activator S, and the production rate caused by the morphogenetic furrow G. In all these cases, there is also a regime σ < σc where the ordered pattern remains intact. For the threshold response in the case of noise in all parameters, except the ones determining h, the two remaining production rates Pu and Ps and for more cases see S2 Fig. Fig 10. Universal threshold response to stochasticity in different model parameters. In each sub-plot, stochasticity is introduced in the following parameters of the underlying differential equations: (a) D u and D s, (b) P a, (c) S, and (d) G. Introduction of parametric noise in all of the former parameters, produces a sigmoid response. The order measure used is the 46 fidelity of the nearest neighbor distance probability distributions and sigmoid curves were used to fit the data. Next let us address the origin of the observed resilience of this model towards weak noise with σ < σc, signified by the initial plateaus observed in Fig 10. The main reason for these wide regimes, exhibiting robustness with respect to parameter stochasticity, is the presence of Scabrous as a diffusible activator (7). Introducing Scabrous decreases the sensitivity of the patterning mechanism to the shape of the inhibitory regions (7) visualized in Fig 2, and thus to slightly misplaced clusters, making the model more robust with respect to parametric perturbations. This effect can be identified in the ordered state where displaced cluster positions do not lead to dissipation of the order of the pattern, see S1 Fig. While this ingredient is not needed for the pattern to propagate, without this element even low stochasticity levels of Du cause it to deteriorate, making the model of pattern formation biologically implausible (7). This effect is quantitatively demonstrated in Fig 11, which compares the response to stochasticity with and without Scabrous as a diffusible activator. In the absence of Scabrous, with the same parameter set and with optimized initial conditions, there is an immediate deterioration of the patterned order upon introduction of infinitesimal noise. In contrast, the presence of Scabrous changes this functional relationship to the observed sigmoid response, 47 thus introducing a resilience scale to the pattern formation mechanism set by the Scabrous production rate Ps and the activation threshold s1. Fig 11. Robustness of Drosophila eye disc R8 cell pattern order versus stochasticity, with and without the diffusible activator Scabrous in the model. The analysis is applied in the same setting as Fig 9, i.e. stochasticity is introduced in the diffusivity of the morphogens D u and D s. In the absence of the diffusible activator s, there is a significantly increased sensitivity to stochasticity for infinitesimal non-vanishing σ. The order measure used is the fidelity of the nearest neighbor distance probability distributions. In addition to the robustness introduced by Scabrous, the threshold response of the order of the pattern to stochasticity in the various parameters shown in Fig 10 is a consequence of three further elements: (i) the spatial discreteness of the underlying lattice, (ii) the threshold activated response of cells to morphogen levels, and (iii) propagation of avalanches from anterior to posterior. As we now discuss, the 48 robustness for low stochasticity levels is a consequence of the lattice discreteness and the threshold response of the cells to morphogen concentrations, whereas the sharp decay is a consequence of propagation of avalanches. To understand how the discreteness of the underlying lattice contributes to robustness for low noise levels, consider the effect of altering the diffusivity of a morphogen. Since the local morphogen concentration drops off at distances on the order of a few cells, and the cells lie on a discrete hexagonal lattice, a small change in the diffusivity radius will not result in any difference unless it leads to the inclusion or exclusion of additional cells. A continuous change in the radius of inhibition will exclude or include new cells when it increases by a factor comparable to the lattice spacing, thus creating a threshold response. Alternatively, consider the effect of the morphogenetic furrow, mathematically denoted by Eq 4, which propagates with a constant speed through the lattice of cells. Two consecutive cells on along the direction of propagation of the morphogenetic furrow will receive the activation signal with a delay lag, δt. Only if the production rates, denoted by G in Eq 1, picked from the normal distributions differ above a finite threshold one will result to a misplaced cell. This effect is shown in the absence of s as a short range activator, when noise in introduced in G, in S5 Fig. An analytical approach to this effect and its contribution in error is carried out in (19). 49 Second, the threshold response of pattern formation to parameter stochasticity is linked to the threshold activation of cells in response to morphogen levels. The turn- on and turn-off of morphogens in the differential equations is governed by step functions with threshold values, aa, au, u1, h1, s1. This has the consequence that the effect of morphogen concentration is binary, i.e. the parameter variation must be sufficient such that the concentration threshold is crossed for activation or inhibition to occur. Nevertheless, even for more continuous cell responses to morphogen levels, there is always a threshold value and the same argument can be made. To rule out whether the threshold response to stochasticity is an artifact of the binary functional form of the Heaviside functions appearing in Eqs 1-4, further simulations were performed, see S2 Appendix. Finally, to explain the fast decay of the pattern beyond threshold stochasticity, the propagation of avalanches comes into play, an effect also identified by (7,8). Beyond σc, R8 clusters start to get significantly elongated and misplaced. The first such elongated R8 cell cluster causes the next row of R8 cell clusters to be misplaced, leading to a cascade onset of pattern disorder. This avalanche effect is quantitatively discussed next, where we analyze the order of the pattern as a function of position. Increasing disorder of the pattern from posterior to anterior 50 Here we use the fidelity measure to evaluate error propagation as a function of position on the spatial posterior-to-anterior axis. The mechanism of pattern formation relies delicately on carefully spaced precursor R8 cells that are placed in the posterior region before the growth process described by the differential equations is initiated. Their positions define the subsequent cluster spacing of activated R8 cells for the entire eye disc. In the presence of stochasticity, this information dissipates during the propagation of the morphogenetic furrow from posterior to anterior, since the ommatidial pattern on a column is defined by the shape of the inhibition signals emitted by the ommatidial row before it. It is thus expected that as the morphogenetic furrow travels from posterior to anterior, errors in the pattern order grow in an avalanche fashion. This effect is observed in the simulations, as shown in Fig 12, where the disorder of the pattern clearly increases from posterior to anterior. 51 Fig 12. Decay of R8 pattern order as a function of position on the eye disc (same parameter set as in Fig 5). The local order measure is calculated as a function of position, from posterior to anterior. The x-axis refers to regions of the simulated eye disc 0-40, 40-80, 80-120 respectively, whereas the y-axis refers to the fidelity probability distance measure applied to nearest neighbor distance distributions. The three colored lines correspond to various stochasticity levels in D u and D s, from σ/ μ = 15% to σ/ μ = 25% in steps of 5%. Local order generally decreases from posterior to anterior. Once the threshold of σ c/ μ ≈ 20% is crossed, there is a very sharp spatial decay of the order of the pattern, signaling the underlying avalanche effect. Within the noise regime studied, the order parameter F decreases by two effects. First, the pattern becomes increasingly disordered as it propagates. Second, a larger part of the pattern becomes disordered. Interestingly, in the high noise regime, the pattern saturates to a value for F that does not increase further from posterior to anterior, see S6 Fig. The higher the noise, the faster it reaches to the saturated value of F. Thus, in the regime of high noise the disorder is largely due to a larger part of 52 the pattern being in the disordered state. Conversely, for low noise regimes, the dominant effect is the intensification of disorder from posterior to anterior. These two effects are also universal, observed for all parameters, since they are the precursors of a threshold response to stochasticity. Previous studies have not, so far, quantitatively analyzed the consequences of finite system sizes on the order of the pattern. To illustrate these finite size effects on the threshold response quantitatively, in Fig 13 we plot the order parameter obtained from numerical simulations for different eye disc sizes, illustrating the effect of size on the robustness of the system. Here one clearly observes that the larger eye-disc, the sharper the threshold response to stochasticity. Since our numerical simulations were performed on lattices much smaller than the actual eye disc, the observed threshold effect can be extrapolated to be much more pronounced for realistic eye disc sizes. 53 Fig 13. Pattern order as a function of eye disc size (same parameter sets as in Fig 8) The threshold response is more pronounced as the length of the eye disc in the direction of the morphogenetic furrow is increased. The order parameter is calculated for cell number 40, 80 and 120 while keeping the size perpendicular to the propagation direction of the morphogenetic furrow at 44 cells. The order measure used is the fidelity of the nearest neighbor distance probability distributions and noise is added on the diffusivity of the morphogens. This variation in the degree of dependency of the pattern order measure upon stochasticity levels across different eye disk sizes is analogous to the finite size scaling of order parameters commonly observed in interacting physical many body systems, such as Ising models on lattices. In analogy to phase transitions in the thermodynamic limit of such models, i.e. for infinite system sizes, here we observe precursor sigmoidal response to increasing stochasticity levels that becomes more pronounced for larger system (eye disk) sizes, ultimately culminating into a true 54 phase transition as Lx → ∞. The structure of the underlying system of coupled dynamical equations suggests that this transition belongs to the universality class of directed percolation (32). While a complete scaling analysis is beyond the scope of this study, it will be provided in a future publication. 1.5 Conclusions The model of coupled differential equations in Eqs 1-4, describing eye disc R8 cell differentiation, provides robust developmental mechanisms that lead to resilience against stochastic perturbations. This robustness can be attributed to the discrete nature of the underlying lattice and the threshold activation of cells to morphogens. In addition, the introduction of Scabrous, as rigorously illustrated in (7), leads to increased resilience against stochasticity, especially in stochasticity in the diffusivity of the inhibitor. Beyond a critical noise threshold σc, there is an acute loss of order with increasing stochasticity due to avalanches of misplaced differentiated cells. This threshold value can be obtained by fitting the response with a sigmoid and is found to occur when the relative noise variation of a parameter reaches a certain value, distinguishing the work from more qualitative observations of robustness (7,8). Furthermore, the preliminary results of introducing noise to multiple parameters, see S2 Fig, suggest an additive effect of the respective variables on the system’s response to noise. This suggests that the 55 number of parameters subject to fluctuations in the regulatory network change the quantitative threshold but the sigmoid response is conserved. Finally, when the R8 pattern formed in the absence of noise, this and previous studies (7) observe small variations from column to column in the pattern. In this analysis, we verified that these variations do not lead to loss of of order even for large system sizes and thus, strengthened the plausibility of this mathematical formulation. This threshold response to noise is an essential characteristic of developmental systems. A lack of resilience would deem biological systems unable to produce complex, nearly perfect structures in an inherently noisy background. This is especially true in eye disc formation with large numbers of ∼ 800 ommatidia. The introduction of Scabrous creates a redundancy that supports the threshold response. The model system studied here exhibits biologically plausible robustness, identifying the sigmoid functional form as one of the emergent responses to noise of developmental models. The expression of a gene depends on a probabilistic outcome determined by the upstream regulatory sequences that act in cis, and the binding and processivity of other molecules that act in trans. Thus, mutational changes in cis-regulatory elements and trans-regulatory factors can affect transcription levels and transcriptional noise (33). These changes in transcriptional noise reveal some of the 56 evolutionary constraints on cis-regulatory elements (34). The sigmoid functional response buffers the system against perturbations. This functional form indicates the effect of Cryptic Genetic Variation (CGV). CGV is largely neutral until it is exposed in certain genetic backgrounds. Robustness in development allows for the accumulation of CGV without any observable changes in phenotype. Once the genetic background changes, for example through mutational perturbation, this previously-neutral bottled-up CGV can be released to produce strikingly different phenotypes (35). Our analysis further illustrates that as the R8 clusters are specified, errors propagate, leading to increased irregularity from posterior to anterior. This rate is important as it suggests that there is an interplay between how large an eye disc is with how perfectly ordered it can be, with larger eye discs being more likely to accumulate some positional errors. This has interesting implications on limiting eye disc size, a hypothesis that can be investigated experimentally. The structural order measures outlined in this chapter can be used to quantify order of patterns in developmental systems containing local point patterns. The application of the radial pair correlation function for nearest neighbors, as specifically performed in this study, is useful for data sets that have nontrivial periodic structure and are small in size. Specifically, the measure is useful for 57 analysis of experimental data of the eye-disc, since the effect of curvature of the eye-disc is eliminated and restrictive boundary effects are resolved. The point- pattern order measure used here, based on probability distance, is independent of the functional form of the underlying nearest neighbor distributions and boundary conditions, thus making it applicable independent of the form of the point pattern. 1.6 Supporting information S1 Appendix. Exact expression of the Laplace operator and noise. https://doi.org/10.1371/journal.pone.0210088.s001 (PDF) S2 Appendix. Threshold response is not an artifact of the Heaviside functions. https://doi.org/10.1371/journal.pone.0210088.s002 (PDF) S1 Table. Numerical parameters of complete mathematical model. This parameter set is taken from previous studies, plugging it in the mathematical equations produces the desired pattern in a biologically plausible regime (7). https://doi.org/10.1371/journal.pone.0210088.s003 (PNG) S1 Fig. Quantification of order in the R8 cell pattern in the absence of noise as a function of anterior to posterior. 58 In this plot, the fidelity order measure for nearest neighbor distances is applied on an R8 cell pattern generated by the mathematical model. The mathematical model is evaluated on an elongated cell grid of 200 vs 44 cells in the case of no noise with the parameters in S1 Table. The order measure used is the fidelity of the nearest neighbor distance distributions. In this case, the fidelity order measure uses as a reference the pattern of the first 0-40 cells in the direction of the morphogenetic furrow (chosen as the x-direction in this chapter) and uses as a target the pattern of the 40-80, 80-120 and 160-200 cells as shown in the x-axis of the plot. The first data point starts at 1 and then saturates to a value of 9.85 for the rest of the slices. The plot supports the claim that, in the zero noise case, the resulting pattern retains its degree of order as it propagates. The reason that the first data point is higher than the other is that the first column of cells is placed manually as an initial condition and it is placed slightly to the right to avoid the effects of reflective boundary conditions. Details on the quantitative measure of order can be found in the Methods section of the main text. https://doi.org/10.1371/journal.pone.0210088.s004 (PDF) S2 Fig. Parametric variation leads to a threshold response when it is applied to a variety of parameters and in different combinations. (a) Parametric variation is applied in equal degree for all parameters except the ones that determine h propagation. (b) Parametric variation is applied in variable 59 amplitude for the pair Pa and Du, Ds. (c) Parametric variation on Pu. (d) Parametric variation on Ps. (e) Parametric variation in equal degree for all the production rates and Du, Ds as in (7). The order measure used is the fidelity of the nearest neighbor distance distributions. Parametric variation was not applied on the elements that define h as it is approximated analytically and the pattern exhibited great sensitivity to such variations. https://doi.org/10.1371/journal.pone.0210088.s005 (PDF) S3 Fig. Probability distance order measures lead to threshold response when uniform distribution is used as a reference. (a)-(b) are generated using the fidelity, F, and (c)-(d) using the Kolmogorov distance, K, for nearest neighbor distances and angles respectively. The probability distance measures are applied on the R8 point pattern with noise added in the model for the parameters Du and Ds. The probability distance order measures need a reference distribution to quantify order. In contrast to the rest of the plots in this article, the reference distributions were chosen to be the uniform distributions. The threshold response to stochasticity is apparent in the plots. As the pattern deteriorates with increasing noise, it approaches the uniform distribution. Thus, the fidelity between the distributions increases and the Kolomogorov distance decreases. https://doi.org/10.1371/journal.pone.0210088.s006 60 (PDF) S4 Fig. Sigmoid response of bond orientation order parameter. The bond orientation order parameter also exhibits the sigmoid response to noise. This further supports the fact that the response of the system to noise is physical. To calculate the bond orientation order, the Voronoi diagram method was used to precisely determine the nearest neighbors in the point pattern. https://doi.org/10.1371/journal.pone.0210088.s007 (PDF) S5 Fig. The threshold response is apparent, although at a lower noise level, even in the absence of s as a local inhibitor. https://doi.org/10.1371/journal.pone.0210088.s008 (PDF) S6 Fig. Fidelity order measure as a function of position on the simulated eye- disc. In this plot, noise was introduced in Du and Ds. Similarly to Fig 12, the plot shows the order of the pattern as a function of position from anterior to posterior. The x- axis refers to regions of the simulated eye disc 0-40, 40-80, 80-120 respectively, whereas the y-axis refers to the fidelity probability distance measure applied to nearest neighbor distance distributions. The conclusion is that the pattern saturates to a value of F. As the noise is increased, this saturation happens earlier in the eye disc. 61 https://doi.org/10.1371/journal.pone.0210088.s009 (PDF) Acknowledgments We thank Avishai Gavish, first author of (7), for answering multiple questions regarding the mathematical model adopted in this chapter. In addition, Alexander Samsonov, Simon Restreppo and George Styliaris for valuable discussions. Lastly, we would like to acknowledge the high performance computing center at USC (HPC) for their support and computing resources. Computation for the work described in this chapter was supported by the University of Southern California’s Center for High-Performance Computing (hpc.usc.edu). 62 Chapter 2. Avalanches During Epithelial Tissue Growth; Uniform Growth and a Drosophila Eye Disc Model 2.1 Introduction Ensembles of epithelial cells give rise to solid and liquid structural properties. Epithelial tissues are often considered viscoelastic materials, which exhibit elastic material properties on short time scales and liquid properties on long time scales (36). A prevailing notion is that the viscoelastic structure of tissues is a consequence of being in a liquid phase, but in proximity to a transition to a glass phase (37,38). Notably, growing and space adaptive tissues do not exhibit glassy behavior, since this would imply solid material properties on long time scales (37). Computational models that use Self Propelled Particle (SPP) models to describe epithelial dynamics capture a density dependent glass transition (37,39– 41). However, such a transition can also be achieved by tuning cell structural parameters and fluctuations (42). In this computational study, we find that the viscoelastic nature of tissue growth can be accompanied by sequences of cell flow avalanches. The tissue remains in an elastic solid structural state during growth. As the cells proliferate, the cell 63 density increases until the internal pressure is sufficient to cause a large scale rearrangement event, which we define as an avalanche. During this rearrangement, the cells collectively flow radially, leading to a sharp increase in size. The properties of the identified avalanches are investigated in the context of a uniformly growing epithelial tissue and a model of Drosophila eye disc growth during the third instar. In nature, avalanches are ubiquitous. From earthquakes to ferromagnets, large scale structural rearrangements arise, whereby a small perturbation leads to a pronounced collective response (43,44). Systems characterized by such avalanche instabilities exhibit self organized criticality (SOC) (45). A main result is that SOC systems give rise to power law scaling signatures in observables, and by fitting to power laws, universal exponents can be obtained. These inferred exponents can be used to identify the strength of correlations in the system and the universality class of the process. Universality classes are groups of processes that share the same scaling properties. As an example, amorphous materials, which include glassy systems, exhibit avalanche behavior when exposed to strain (46,47). A successful characterization of such avalanches utilizes the physics of self organized criticality and phase transitions (46,48). Recently, sheer stress induced avalanches were identified in the context of a vertex model, and both their universality class and material properties were mapped to amorphous 64 materials (49). Furthermore, signatures of criticality were identified in vitro for the Drosophila wing disc epithelium, suggesting that the nonlinear properties of collective processes in amorphous materials and the vertex model are at play in biological tissue growth (49). Epithelial tissue growth and migration has been the subject of several recent in vitro studies. These generally branch into two categories, either focused on confined tissues (50–53), or on freely expanding tissues (54–57). In the former category, processes that involve density induced motility transitions, front propagation, and spatial sensing are probed. In the latter category, processes involving mesenchymal-to-epithelial transition, contact inhibition, size dependence of growth and boundary formation are investigated. In both categories, cell motility and collective phenomena are a central theme. In the context of growth, a ubiquitous collective phenomenon, identified in vitro for both freely expanding and constrained growing cell tissues, is the formation of vortices by cell trajectories (53,55). Similar to our discovery of avalanches during growth, which is reported here, the vortices are a density driven collective phenomenon (53). Interestingly, the cell trajectories during avalanches also exhibit vortex signatures. However, the avalanches also involve a coordinated radial motion of cells not reported in (53,55). In our study, the underlying model treats the tissue in a confluent state, as in epithelial tissues in vivo, where growth 65 occurs in compartments separated by boundaries. Thus, it provides an appropriate framework linking the phenomena probed with what takes place in vivo. Biophysical properties, such as mechanical forces, impact cellular processes during growth phases in development (58). Epithelial cells have been shown to regulate their shape and division rate by responding to mechanical force stimuli (50,59–61). Furthermore, control of the cell division plane via cell-shape and mechanical considerations has been experimentally detected in animal and plant epithelial cells (59,62). Numerical implementations of epithelial tissue growth, which include mechanosensing and division plane control, have shown to generate improved epithelial mosaics that represent in vitro structures more closely (63,64). Single cell resolution computational models are a valuable tool to model growth, since each cell can respond to its distinct state and incorporate phenomena at the single cell level (65,66). We use a vertex model (67) to simulate the spatial dynamics of epithelial tissues. Vertex models have provided successful mechanical simulations of confluent epithelial tissues with single cell resolution (68,69). This type of epithelial tissue computational modeling treats cells as polygons and uses energy considerations to track the spatial arrangement of cells, vertices, and edges. It has been successful in simulating apoptosis, mitosis, and cell shape changes (68,69). In the vertex models, 66 each edge, perimeter and area of the cells contribute to the energy and affect the structure of the tissue. There are two main categories of dynamics implemented in vertex models. The first kind is the active vertex model, which explicitly updates the structure using the equations of motion for each cell or vertex in the overdamped regime (70). The other kind is the quasistatic vertex model, which assumes that relaxation occurs in very short time scales and drives the system to the first accessible local energy minimum using gradient descent (67). The local energy minima correspond to mechanically stable configurations, as the negative of the gradient of the energy corresponds to the forces on the vertices. These tissue models have been verified in Drosophila development with experimental data and laser aberration experiments, and have been recently successfully implemented in describing complex phenomena such as epithelial tissue folding and growth (67,70,71). Epithelial tissues, across metazoa, have characteristic cellular packing geometries and distributions captured by the vertex model. One characteristic of epithelial geometry is formulated in Lewis' law, which states that the average cell size for cells with n neighbors, 𝐴 𝑛 , is proportional to the number of neighbors (72). This characteristic distribution has been found in a variety of metazoan tissues (67,73– 75). Another characteristic signature prevalent across growing epithelial tissues is the distribution of cell sides (63,64,67,76). The cells are predominantly 67 hexagonal, with pentagons and heptagons following suit. These experimentally observed cellular distributions have been reproduced by models that combine single cell geometrical configurations and proliferation (63,65,67) . To investigate epithelial tissue growth in the context of an in vitro system, we simulated epithelial growth in the Drosophila eye disc during the third instar. The Drosophila eye disc is part of the eye-antennal disc, which begins with approximately 70 founder cells and develops to a final size of approximately 44,000 cells (77). During the first half of the second instar, the eye-antennal disc grows via rapid cell proliferation, and at the late second instar, it segregates into the retinal progenitor and the antennal/head epidermal fields (78). Photoreceptor development initiates at the beginning of the third instar, on an elliptical epithelial sheet, the retinal progenitor, composed of undifferentiated cells. The size of the retinal progenitor starts at approximately 10,000 𝜇 𝑚 2 (79). At that stage, a visible indentation of the tissue, termed the morphogenetic furrow (MF), begins propagating linearly from the posterior to the anterior of the eye disc, leaving behind differentiating photoreceptors (9,15,78,80). Cells within the MF region are constricted, with their apical surface reduced by a factor of about 15, from 9𝜇 𝑚 2 ahead of the furrow to 0.6𝜇 𝑚 2 (15). There are two mitotic waves, the first mitotic wave occurs anterior to the MF, where cells differentiate asynchronously, and the second mitotic wave affects undifferentiated cells between photoreceptor 68 cells behind the MF (15,81). The first mitotic wave produces the majority of cells of the final eye disc and takes place ahead of the MF, whereas the second mitotic wave ensures that there are sufficient cells for proper ommatidial formation (77). The second mitotic wave is more controlled, takes place posterior to the MF, and occurs in a synchronized and organized manner (82). The contribution to growth of the second mitotic is negligible compared to the first mitotic wave and will not be included in the model (79). The growth dynamics of the third instar Drosophila eye disc have been described quantitatively. The growth rate over time has been determined by analyzing experimental tissue size data over different developmental time points (79). Both exponential decay and inverse area functional forms fitted the data exceptionally well (79). Furthermore, the area rule has been the focus of a subsequent study, which identifies a diluting and decaying transcript, cytokine unpaired, as a possible candidate for cell growth regulation (83). Thus, functional forms and coefficients that describe the time dependent eye imaginal disc growth rate quantitatively are available and were used to simulate the process. In order to simulate the propagation of the MF, cellular signaling was incorporated on top of the vertex model. In vitro, the MF is a result of signal transduction of multiple morphogens with multiple roles. The primary signal that is associated as 69 the initiator and driver of the MF is hedgehog (Hh) (78,80,84–86), a morphogen that is sourced from differentiating neural cells posterior to the MF. Hh diffuses in the epithelial tissue up to and including the constricted cells in the morphogenetic furrow (80). The propagation of the MF is mediated by Decapentaplegic and regulated by Hemothorax, Hairy and extra macrochaetae (80,87–89). In the model herein, a simplified gene regulatory network was used, see (90) for a model with a detailed treatment of the gene regulatory network. The propagation of the MF was simulated by introducing an effective morphogen that induces its own expression while diffusing from posterior to anterior (see methods for details). Few studies have so far been published that combine signal transduction processes, transcriptional responses and vertex models (66,91–94). In the model by Wartlick et al, the authors were able to reconstruct the observed growth of the wing disc via an integrated signal transduction and the vertex model by Farhadifar et al (94). Schilling et al. developed a model for the production, diffusion, and sensing of Hh in the Drosophila wing primordium and explained how anterior/posterior cell sorting can be explained by a homotypic boundary model (92) . We will use a similar approach, i.e. a finite-volume method defined by a tessellation of the centroids of the cells. Notably, an alternative approach has been introduced utilizing the Lagrangian–Eulerian (ALE) formulation and finite element method to describe morphogen transport within a tissue coupled to a vertex model (93). Furthermore, 70 an earlier study was able to model wing disk growth by considering a complex mechanical feedback, with cell constriction affecting the cell cycle and a biased angle of mitosis (65). Finally, a report on the zebrafish retina was able to model highly ordered packing via coupling planar cell polarity proteins with edge tension (66,68). The eye disc model developed here belongs in the class of models that couple signaling and vertex models. In the context of the eye disc, the model is novel in its treatment of the epithelial structure and a natural extension of the viscous fluid based treatment of drosophila eye disc growth seen in (90). Results 2.2 Generic uniform growth model In order to study cellular avalanches during epithelial tissue growth, a vertex model is used, which is based on an energy functional that optimizes cell surface area, and minimizes edge tension as well as contractility of the cell perimeter. In our initial implementation, the tissue is initialized to a small circular size and is allowed to proliferate into empty space. The cell doubling rate is set to be the same for each cell and constant over time. In addition, the doubling rate is also set to weakly depend on cell area (64) , approximating cellular mechanosensing , 71 with bigger than average cells having a higher division rate compared to smaller cells (see methods for details). Contrary to the periodic boundary condition case, we discovered the addition of size dependent doubling rate was necessary for open boundary growth (SI Fig. 1). Without mechanosensing, the cells become unrealistically packed, there is significant apoptosis and the increase of tissue size is suppressed. Importantly, avalanches were detected with or without the size dependent doubling rate. Furthermore, the tension at the tissue periphery is set to an increased value in order to enforce a smooth boundary. The approximately circular symmetry of this system allows for probing the radial profile of structural quantities. Three frames at different times of a simulation realization are shown in Figs. 1A-C. In Fig. 1C, the tissue increases significantly in size, as a constant cell division rate leads to exponential growth. A complete video of a simulation of this growth process can be found in SI Vid 1. Fig. 1. Uniform growth model. A-C: Successive frames of the tissue growth process at different time points, T=20.0h, 30.0h and 40.0h respectively. The tissue starts at a small size composed of 50 cells and is simulated until it reaches 350 times its original size, allowing the understanding of different growth regimes. Cell 72 division takes place at a constant doubling rate of ~6hrs. The cells are color coded with respect to their area. The tissue retains a smooth boundary during growth and an approximately circular shape. In early times, the tissue grows continuously and exhibits discontinuities, or avalanches, at later times of ~40 hrs. Next, we probe aspects of the collective cell motility during a characteristic avalanche. In Fig. 2A, the cell trajectories during an avalanche event are visualized. During the avalanche, the cells move collectively, with displacement increasing in a radial fashion. Interestingly, the trajectories of the cells also have a significant tangential component. The tangential component is visualized in Fig. 2B by color coding the trajectories with their respective angles. Trajectories of cells in the middle of the tissue have a significant tangential component, contrary to the radial direction of trajectories on the periphery. This signature indicates the presence of a cellular vortex. Interestingly, tissue scale vortex formation has been experimentally observed in growing epithelial monolayers (55). In addition, the radial profile of the vortex, formed during an avalanche, is consistent with the profile reported experimentally (55). For the cell trajectories without the onset of avalanches, see SI Fig. 2. 73 Fig. 2. Single cell trajectory characteristics during avalanches. A: Cell trajectories during an avalanche, color coded with their respective length. The cells in the outer region displace the most, whereas the trajectories are shorter at the center of the tissue. The tissue boundary before and after the avalanche is outlined, visualizing the net expansion of the tissue. B: Cell trajectories during an avalanche, color coded with their respective angle relative to the x-axis, as defined by the inverse tangent function. On the outer region of the tissue, the trajectories are approximately radial. The vorticity is evident towards the center of the tissue, where the trajectory angles start to change. The time dependence of the evolving total area and average cell area in this realization is depicted in Figs. 3A&B. Initially, the tissue grows at a predominantly continuous rate, as cells can easily rearrange within the tissue. As it becomes larger, at approximately 30h, more pronounced discontinuities in the area growth are detected. These avalanches correspond to large scale cell rearrangements. Due to the initial small size of the tissue, the increased boundary tension dominates the energy term, and the tissue begins its growth at a high cell 74 packing state. As the tissue grows further, the average cell area increases, as the boundary tension energy becomes negligible compared to the bulk energy. Fig. 3. Uniform growth model observables. A: Total tissue area versus time. The tissue is initiated as a small size circle, and the cells undergo mitosis with a constant doubling rate. At early times, the tissue follows predominately continuous growth. After ~30 h, the tissue displays growth predominately with discontinuities. B: Average cell area versus time. The discontinuities in average cell area become more pronounced and regular at ~30 h. C-D: Avalanches in the uniform growth model. C: Average cell area versus radial distance from the tissue center at consecutive time points. The spatial signature of the average cell area before an avalanche occurs is depicted. As cells divide, the average cell area decreases globally, with higher packing from the center to the tissue boundaries. D: Average cell area versus radial distance from the tissue center at consecutive time points. The spatial signature of the average cell area is depicted shortly before and after an avalanche occurs. Before the avalanche, there is a positive gradient from the tissue center to the tissue boundaries for cell size. After the avalanche, the average area increases with a discontinuity. The magnitude of the gradient diminishes. 75 The radial profile of the average cell area is probed to understand the onset of avalanches. Before the avalanche, the cell density increases as cell division takes place (Fig. 3C). The spatial signature suggests increased cell packing in the center of the tissue, which diminishes outwards. Right after the avalanche, the radial cell density discontinuously increases (Fig. 3D), and the gradient of cell packing is less pronounced after the avalanche. This is consistent with the radial increase of trajectory length shown in Fig. 2A. Based on these measured quantities, we can address what causes the avalanches and what controls their magnitude. To this end, we repeat 150 simulations of the process and isolate the avalanches. To do so, we skip the first 30 h of development and record an avalanche if a cell density discontinuity is equal to 0.1 𝜇 𝑚 2 or greater. To quantify the spatial inhomogeneity of cell packing, we use the standard deviation of the average area as a function of radial distance from the tissue center. Interestingly, Fig. 4A shows that there is a strong correlation between the average cell area before the avalanche and the jump magnitude, with a spearman’s rank correlation coefficient of ρ=-0.83 and a p-value of 0. Furthermore, SI Fig. 3 illustrates that there is weak correlation between spatial density inhomogeneity and jump magnitude. Thus, avalanche onset and strength 76 are more sensitive to the global density of cells rather than to the details of their spatial distribution. Fig. 4. Avalanche Characteristics. A: Average cell area versus avalanche magnitude. The avalanche magnitude is quantified as the percent discontinuity in the total area. The more tightly packed the epithelium, right before the onset of the avalanche, the larger the total area discontinuity. B: Probability density of avalanche magnitudes. The functional form is a power law distribution, with a fit shown by the red dotted line. The exponent found is τ=1.27±0.04. Such power law distributions are a universal characteristic of avalanches found in nature. To understand whether the system exhibits self organized criticality (SOC), the probability distribution of avalanche magnitudes is shown in Fig. 4B. To provide more avalanches and solidify the fit, the growth model was run 150 times for an extended runtime, until the tissues reached 300,00 𝜇 𝑚 2 . The magnitudes are distributed according to a power law distribution, with an exponential cutoff, which is a ubiquitous property of SOC. The scaling exponent is found to be 77 τ=1.27, in agreement with analogous exponents observed in sheer induced, stress avalanches in overdamped elastoplastic models (46,95) . More information about the fitting procedure and comparison with exponential fit is shown in SI Fig. 4. This result reinforces the notion that the vertex model and epithelial tissues share common material properties and universality classes with plastic amorphous systems (49). Varying the boundary tension and the mechanosensing amplitude did not cause the exponent τ to change ( SI Fig. 5) . 2.3 Drosophila eye disc growth model In order to study the cellular avalanches in the context of a particular developmental process that is experimentally accessible, a more detailed model was developed. In this version of the model, the cell doubling rate is set by experimental data, the tissue is initialized in an elliptical shape, and there is signaling of a single effective morphogen which drives the propagation of the morphogenetic furrow (MF) see methods for details. The boundary tension at the boundary is again set to the same increased value as the original uniform growth model to ensure a smooth boundary. A complete video of a simulation of this growth process can be found in SI Vid 2. 78 A visualization of the simulation of the growth process at three different time steps, initiation, middle, and end of growth are shown in Fig. 5. Initially, the transcription factor is induced in the posterior boundary cells (Fig. 5A). This will eventually cause the auto-activation of transcript production in neighboring cells, and will lead to the initiation of the MF. At this time point, all the cells are dividing in the tissue. At the midpoint, the MF has traversed approximately halfway through the tissue (Fig. 5B). Cells entering the MF immediately stop dividing and become immobile 2.5 h after the MF has passed, to simulate differentiation. At the midpoint, the growth rate has decreased by approximately 80%, and the majority of cell divisions has taken place. Finally, at the end point (Fig. 5C), the growth process is complete, as the MF has propagated to the anterior tissue boundary, and all cells have arrested their cell cycle. 79 Fig. 5. Drosophila eye disk tissue at different times of growth. A-C: Vertex model simulation of tissue growth, with cells in the MF highlighted in yellow. A: Time T=3.4 h; the MF is initialized at the posterior tissue boundary. B: Time T=47.0 h; the majority of cell divisions have taken place, while the MF has linearly propagated to the bulk of the tissue. C: The MF reaches the anterior boundary of the tissue, signaling the cell cycle arrest of all cells and termination of growth. During the growth process, the structural details of the cellular mosaic of the tissue were measured. To quantitatively describe the mosaic of the epithelial cells, average areas versus number of sides and the distribution of polygon classes were plotted (Fig. 6). The measured quantities match experimentally established, system independent, epithelial structure signatures. The emergent relationship between average cell size of a given polygon class, versus number of edges is linear (Fig. 6A), in accordance with Lewis’s law. Furthermore, the distribution of cell polygon classes matches the corresponding distribution found in vitro (Fig. 6B) with P=0.98. Thus, the model gives rise to realistic epithelial packings. 80 Fig. 6. Quantifying epithelial packing during the simulated growth process. A: Average cell area per number of cell sides, divided by the mean area, is plotted versus cell number of sides. The values are obtained by combining data from all frames in the simulation. A linear relationship is observed, in accordance with Lewis' law. B: Distribution of number of cell sides. The simulation values are obtained by combining data from all frames in the eye disc simulation. The experimental values were obtained from (76) and correspond to Drosophila melanogaster third instar larval wing disc. Aspects of the tissue growth process can be captured by quantifying the temporal evolution of key observables. The anterior area depicted in Fig. 7A initially grows, reaches a plateau and then decreases until it vanishes, as needed for growth termination. There are three effects contributing to the anterior area profile. Firstly, new cells are introduced due to cell mitosis. Second, since the MF separates the anterior from the posterior compartment, as the MF moves across the tissue, it turns anterior cells to posterior cells. Third, the rate of mitosis in the anterior decreases over time. The posterior area depicted in Fig. 7C grows solely via the cells that exit the anterior area due to MF propagation. Initially, the eye disc is small, and the MF passes through a small number of cells. At the midpoint, the posterior area grows the fastest, as the MF is at the center of the tissue and passes through a large number of cells. Towards the end of the growth process, it tapers down, as it reaches the anterior edge. The total area depicted in Fig. 7B is the sum of the anterior and posterior compartments. Finally, the posterior length depicted in Fig. 7D, is shown to increase linearly over time. It is defined as the 81 length from the posterior edge midpoint to the midpoint of the MF. Initially, there is a plateau, which corresponds to the initialization time needed for the MF to begin propagating. Casting aside the initialization time, the linear nature of the posterior length (𝑅 2 = 0.99 ) indicates that the posterior length indicates that the MF is propagating at a constant speed throughout the growth process. The value of the slope of the posterior length corresponds to the MF speed, which is 3.4 𝜇𝑚 /ℎ𝑟 , matching the experimental value. Fig. 7. Drosophila eye disc growth dynamics. A: Time dependence of anterior area, defined as the portion of the tissue ahead of the MF. Cell division takes place only on the anterior portion of the eye disc. B: Time dependence of total area. The final value is in the same order of magnitude found in wild-type Drosophila strains. C: Time dependence of posterior area, defined as the portion of the tissue before the MF. Cell division is arrested in the posterior portion of the tissue. The posterior grows solely by cells exiting the propagating MF. D: Time dependence of the posterior length, defined as the distance of the posterior edge to 82 the midpoint of the MF. Aside from initialization time, the posterior length increases linearly, with slope of 3.4μm/hr, in accordance with MF propagation experimental data. Similarly to the simpler model, the quantification of the growth process reveals discontinuities, i.e., avalanches, in the temporal evolution of the tissue area. This can be seen in the anterior area in Fig. 7A, and consequently also in the total area (Fig. 7B). To understand better what is taking place at the single cell level and across tissue compartments, the average cell area during the growth process is plotted in Fig. 8A. In Fig. 8B we visualize the tissue area avalanche magnitude for each tissue compartment. Comparing Fig. 8A to Fig. 8B, one can see that each avalanche is accompanied by a discontinuity in the average cell area. Furthermore, the discontinuities take place after the average cell area decreases and always lead to an increase of the average cell area. Therefore, the avalanches are a consequence of increased cell packing in the eye disc that leads to a global rearrangement of cells. Also, the avalanches mainly take place in the anterior of the eye disc (Fig. 8B). 83 Fig. 8. Avalanches in eye disc growth dynamics. A: Time dependence of the average cell area. The avalanches are captured by discontinuous transitions from high cell packing to low. B: Time dependence of tissue area jumps for total, anterior and posterior area, marked by green, red and blue respectively. The avalanches lead to tissue area discontinuities, with the majority occurring in the anterior. Changing apical constriction affects the onset of the avalanches. In the model, we can tune the strength of apical constriction in the MF. In Fig. 9, the MF cell constriction is incrementally increased from no constriction to high constriction, by changing the preferred area of cells in the MF. We define the apical constriction factor as the percent decrease in the preferred cell area of cells in the MF. Fig. 9A shows the effect of the constriction strength on the cell number, with higher constriction leading to significantly larger cell numbers. On the other hand, Fig. 9B shows total area, which reflects that eye discs are larger when the apical constriction in the MF is increased. However their relationship is not a scalar multiple, as average cell size is variable between different apical constriction strengths. This can be seen by inspection of two simulation examples (Figs. 9C & 9D) of the average cell size, one for no apical constriction and one for strong 84 apical constriction in the MF. When the apical constriction is increased, the avalanches become less frequent and the average cell area before they take place is decreased. Thus, apical constriction makes it harder, in the sense that cells reach a higher packing state, for global cell rearrangements take place. Finally, since the reaction-diffusion equation is discretized on the lattice of the cells, independent of packing, a decreased cell size leads to a decrease in the MF speed (SI Fig. 6). Since the growth rate decays based on the posterior length and thus the MF speed, apical constriction makes the growth rate decay slower. Fig. 9. Avalanches depend on the magnitude of apical constriction in the MF. A: Time dependence of the cell number. The larger the apical constriction, the higher the cell number. B: Time dependence of the tissue area for different constriction strengths. The larger the apical constriction, the bigger the tissue area. C-D: Average area for no constriction and moderate constriction. The case of no apical constriction corresponds to 85 more frequent avalanches that take place in higher densities. This highlights the observation that apical constriction in the MF makes avalanches harder to achieve, and leads to higher cell packing. 2.4 Model Predictions The model predicts that proliferating epithelial tissues, with short relaxation time scale dynamics, exhibit discontinuous spurs of growth. The onset and strength of the discontinuities, which we term avalanches, are positively correlated with increased cellular packing. Thus, the model predicts that the average cellular area depends strongly on whether the tissue growth is before or after an avalanche. Furthermore, the model predicts the formation of tissue scale vortices during growth, in agreement with experiments (55). In the context of the Drosophila eye disc development, with the assumption that the dynamics fall into a regime of short relaxation time scales, we expect that the average cell density will exhibit considerable variation during growth. Since avalanches take place less frequently at later times during tissue growth, the variations should be less intense at later developmental times, consistent with a more densely packed eye disc. In addition, the stochastic onset and strength of avalanches suggests that cellular density should vary between eye discs at the 86 completion of eye disc growth. Increase of the apical constriction in the MF makes avalanches harder to attain, with cells needing to be more packed. Thus, mutant eye discs with lower apical constriction should, on average, exhibit a larger average anterior cell size compared to the wild-type. 2.5 Discussion To date, most insights regarding the material properties of epithelial tissue have been based on a viscoelastic description. Epithelial tissues have been shown to give rise to both solid and liquid material properties. In addition, numerous cellular properties have been identified that coordinate and regulate cell proliferation. Our vertex model implementation predicts that, as tissues grow, the growth process is dominated by avalanches of cells. However, do epithelial tissues exhibit avalanches in vivo under proliferation stress? And which biological and biophysical processes control the onset and magnitude of the avalanches? Addressing these issues is imperative for developing realistic integrated biophysical computational models that accurately describe developmental processes. Especially, the interplay between growth and cellular signaling requires a valid description of the underlying cellular mosaic and its dynamics. 87 Similar to the signatures identified in our model simulation of avalanche dynamics, time resolved evidence of cell area oscillations have been reported in expanding epithelia (51). However, so far, experiments that involve freely expanding cell cultures have not yet reported this phenomenon. One of the reasons may be the boundary conditions: in freely expanding epithelia, the boundary expands by forming protrusions with cells, often losing confluence. This is in contrast with the simulated system, where the tissue retains its confluency and smooth boundary. The computational model treats the tissue dynamics within a quasistatic paradigm. That is, it is assumed that the energy relaxation time scale is very short. Changing the solver from gradient descent to overdamped dynamics, will affect the prevalence of avalanches. The magnitude of viscosity, cell properties and the nature of cellular fluctuations control the strength and prevalence of avalanches. When the relaxation time scale increases, the avalanches are smoothed out until they cannot be resolved. In vitro, the relaxation time depends on details of the cells, such as cell cycle progression, expression profiles of cadherins and acto- myosin activity. Therefore, the choice of dynamical update rules have profound and fundamental effects on collective phenomena and growth processes, stressing the importance of experimentally quantifying fluctuations and relaxation 88 timescales (96,97). Our computational study illustrates the intricacy of growth processes, their dependence on the dynamics and emergent phenomena when coupling signaling and mechanical responses. 2.6 Computational Methods Spatial Structure There are many manifestations of vertex models. Here, we adapt a model introduced by Farhadifar et al, which has been successful in emulating the growth of drosophila wing disk epithelia. In this model, the tissue is assumed to evolve quasistatically, i.e., the tissue re-arranges instantaneously to the nearest local minimum via gradient descent (67). This approach originally uses a 2D description of the tissue, but it can be generalized and applied in 3D (71). The energy function considers cell elasticity, actin myosin bundles, and adhesion molecules. Mathematically, it contains an elastic term for area, tensions for each edge and contractility for the perimeter, 89 , where α is the cell index, 𝐾 𝛼 is the elastic coefficient for compressibility, 𝐴 𝛼 is the area and 𝐴 𝛼 0 the preferred area , <i,j> is the summation over all edges,𝛬 𝑖𝑗 is the edge tension and similarly 𝛤 𝛼 is the contractility of the perimeter 𝐿 𝛼 . Determining the epithelial tissue structure in a given conformation of cells is more intricate than sole gradient descent energy minimization. Cellular structures undergo topological transitions, if they further minimize the energy of the configuration. The possible transitions have been addressed in the context of foam physics (98). For instance, cells can undergo neighbor exchange via a T1 transition (67,69). Cells do not have a fixed number of edges, and if it is energetically favorable, new edges can form via T3 transitions (67,71). Contrary, cells that shrink undergo inverse T1 transitions, which reduce their number of edges. When a cell obtains three edges or reaches a critical apoptotic size, it can be removed from the structure via a T2 transition (67,98). These transitions conserve the topology of the tissue while implementing energy minimization. The coding environment for the spatial structure was taken from the tyssue repository on github. The tyssue library (DOI 10.5281/zenodo.4475147) is a product of code refactoring originally from (71). 90 Neglecting the attachment of the eye imaginal disc to the antenna disc and the peripodial membrane (78), we approximate the eye disc to be suspended and with open boundary conditions. Due to the asymmetry of the growth process, as the MF linearly propagates, one cannot implement periodic boundary conditions, as it is often done in other systems (94). To approximate the effect of boundary cells, increased tension was implemented at the cell edges along the boundary. The boundary tension parameters were chosen to enforce circularity of the final structure SI Fig. 7A. The increased boundary tension did not lead to significant changes in the growth process SI Fig. 7B. This is motivated by the established differential adhesion on tissue boundary layer, often thought as creating an effective surface tension (99,100). To simulate the structural component of the morphogenetic furrow, a simplistic approach was taken. After determining the cells that make up the MF, their preferred area was instantly set to a decreased value. We term apical constriction factor the percent decrease in preferred cell area of cells in the MF. Furthermore, the preferred area for cells that exit the MF was set back to the default value. In the case of the eye disc model, 2.5 hrs after the cells exhibit the MF, they are artificially made immotile by no longer considering their contribution to the energy term, to simulate the differentiating region. 91 Growth In the vertex model, cell division can be achieved by directly introducing a new edge in a target cell. Some implementations include area growth, but it has been shown that area does not correlate with volume (64), so we use an alternate implementation. In this implementation, the new edge is chosen by drawing a line through the cell centroid with random orientation (65). The two daughter cells inherit the concentration of the mother cell. To simulate the cell cycle, a variable V, which corresponds to volume, is assigned to each cell (64). This variable increases over time, and when it crosses a constant threshold, arbitrarily chosen to be 1, the cell divides. To model the mechanical feedback of the cells, also known as mechanosensing, the magnitude of increase over time includes an area dependent term (64). This formulation has been shown to recreate epithelial cell topologies similar to the ones seen in experiments. Finally, to match the experimental cell doubling, the step increase of this is normalized by 𝑡 𝑛𝑜𝑟𝑚 (𝑘 ), which depends on the growth rate k. The variable and its step increase are shown in the following equation: 92 In this equation 𝛥 𝑉 𝛼 is the volume contribution in a single iteration for cell α, 𝑉 𝑑 is the threshold volume for division, 𝛽 is a constant growth contribution rate, 𝜇 determines the strength of mechanosensing contribution to growth and 𝑡 𝑛𝑜𝑟𝑚 is a factor that matches the experimental growth rate. Finally 𝐴 𝛼 , 𝐴 , 𝐴 (0) represent the area of cell α, the average cell area in the tissue and the preferred area of the cell respectively. Note that the preferred area of a cell is a structure variable in the vertex model implementation. The functional form of 𝑡 𝑛𝑜𝑟𝑚 depends on the growth profile we want to use. For the exponential growth profile used to model the eye disc, 𝑘 = 𝑘 0 𝑒 −𝛿 𝑃𝐿 𝐿 𝑝 (79,90), where 𝛿 𝑃𝐿 is a constant and 𝐿 𝑝 is the posterior length (used interchangeably with developmental time since 𝐿 𝑝 = 𝑣 𝑀𝐹 𝑡 ) . As an approximation, we assume that the target cell has an average area simplifying the equation to: Note that since cells start at 𝑉 𝑑 /2 as a starting volume, they require 𝑉 𝑑 /2 to divide. Now, 𝑡 𝑛𝑜𝑟𝑚 can be derived from the equation: Where in the above, we assume that β and 𝑘 0 are in units 1/𝑠 . 93 Eye Disc Model: Signaling In order to describe the progression of the MF, a simplified approach was taken. The complete complex interactions of multiple morphogens and transcripts are not needed for this study. To that end, signaling was modeled by considering only a single effective transcription factor obeying an auto-activation regulatory loop. This transcription factor obeys the reaction diffusion equation shown below. The differential equation was discretized on the cell center lattice generated by the vertex model (94). The lattice was approximated by a hexagonal packing and the appropriate discretization method was implemented (101). In this equation, c is the concentration of the transcription factor 𝐷 𝑐 is the diffusion coefficient, 𝑃 𝑐 is the production rate and 𝜆 is the decay rate. The θ depicts a step function which becomes 1 when the concentration crosses 𝑐 0 . The emergent dynamics of this equation, after initializing a concentration or source, give rise to a non-dissipating linearly propagating concentration gradient. This is ideal for the eye disc, as the MF is known to propagate at a constant speed of 3.4𝜇𝑚 /ℎ𝑟 , as measured by the dorsal-ventral point of the MF (79). Finally, cells 94 that constitute the MF were defined as cells that have a concentration between a range, chosen to be 0.4 and 1. Given the autoregulatory loop for MF propagation, an initial flux or concentration of the effective morphogen is necessary. In vitro, Hh diffuses through the posterior tissue boundary and initiates the signaling cascade necessary for MF propagation. To initialize the propagation of the MF in the eye disc simulation, we follow a similar approach with (90). The posterior boundary cells obeying 𝑥 (𝑡 ) ≤ 0.2 𝐿 𝐴𝑃 , where 𝐿 𝐴𝑃 is the total anterior to posterior length, received a boundary flux of the morphogen equal to 40% of the autoproduction rate 𝑃 𝑐 . Varying this parameter simply changes the time scale of the initialization phase. Furthermore, the boundary flux was tapered down by a time dependent term 𝜏 (𝑡 ) = 1 − 𝜃 (𝑡 − 10ℎ𝑟 ), such that the boundary flux ceases after 10 hrs. Model calibration and initialization The specifications to initialize the eye disc size and shape were taken from (90). An elliptical tissue was extracted from a previously grown epithelial tissue under periodic boundary conditions as in (64,67). For the uniform growth model, a circular tissue with a radius of 11μm was extracted from the same previously grown epithelial tissue. At the start of both uniform growth and eye disc simulations, cells are set to be at a random point of their cell cycle. That is done 95 by setting the variable V, defined above, at a random value from 0.5 to 1. Note when V=1, a cell divides and V is halved, making the minimum V be 0.5. Next, we discuss how we set the spatial and temporal scales in the simulation. Firstly, the time and space scales used in the generic uniform growth model were taken from the eye disc. Given that the average cell size in the posterior portion of the eye disc is known to be 9𝜇 𝑚 2 (15), obtaining the average cell area from a sample simulation and equating it to the experimental value provides the length scale. For an intermediate value of constriction, the MF propagation speed was measured. Then, it was equated to the experimentally derived value of 3.4 𝜇 𝑚 /ℎ𝑟 (79) to obtain the time scale. This allowed for the conversion of simulation time to developmental time and thus, the calibration of the growth rate. The apical constriction is known to cause the 9𝜇 𝑚 2 cells ahead of the furrow to 0.6𝜇 𝑚 2 (15). The two dimensional vertex model cannot support such an intense constriction. To this end the apical constriction is varied and the corresponding emergent behaviors are investigated. Coupling signaling and tissue structure iterations A custom script was developed for the eye disc model, which combined reaction diffusion based signaling and the vertex model. Firstly, the tyssue library was used for the vertex model and its dynamical updates. For signaling, which 96 translates to ordinary differential equations under discretization, python’s scipy.integrate.solve_ivp() function was used. The two aspects were updated sequentially. The differential equations were time-evolved for a constant interval of time. Then, the properties of the cells, a subset of which depend on the morphogen concentration, were updated. Afterwards, the tyssue solver was used to relax the system to the nearest energy minimum. Then, the configuration of the resulting structure was plugged in the differential equations and they were time propagated. The time resolution of this transition between signaling evolution and mechanical update was increased until the macroscopic system observables were invariant. The mechanical structure and cell state were updated every 4 minutes of developmental time. 2.7 Supplemental information Supplemental Movies 97 SI Vid 1 winglike_ growth.avi 98 SI Vid 2 eye_growth.avi 99 Supplemental Figures SI Fig. 1. Mechanosensing analysis in the uniform growth model. A: Cell death during growth for different mechanosensing amplitudes. When mechanosensing is toggled off, at zero mechanosensing magnitude (μ), there is significant cell death due to proliferation stresses. When non-zero values of mechanosensing are used, cell death becomes negligible. B: Area versus time for different mechanosensing amplitudes. The mechanosensing amplitude significantly affects the growth process. Since this is exponential growth, small variations intensify exponentially. Close to μ=0.04, which is the parameter used and is taken from (64), the area does not deviate when changing mechanosensing. Growth discontinuities or avalanches are present for all mechanosensing parameters. C) Division number versus time for variable mechanosensing magnitudes. When mechanosensing is turned off, the cells divide periodically since all daughter cells grow instantaneously. As the mechanosensing amplitude increases, the division number increases. That is the case because ΔV is bounded by 0, meaning cells can only grow or arrest their cell cycle, making the ΔV contribution due to mechanosensing biased to be positive. D) Average cell area versus time for different 100 mechanosensing amplitudes. When mechanosensing is off, the tissue does not reach the preferred area of 9𝜇 𝑚 2 . SI Fig. 2. Cell trajectories between two successive avalanches. The movement of cells recorded after an avalanche at T=44.1h and before the next avalanche at T=44.9h. Cell movement does not exhibit the collective behavior seen during an avalanche. There are disparate regions of the tissue where cells exhibit significant displacement and the boundary does not significantly expand. 101 SI Fig. 3. Standard deviation of the radial average cell area right before an avalanche versus the tissue jump magnitude after the onset of the given avalanche. The y-axis provides a measure of cell packing inhomogeneity. The correlation was found to be weak but non-negligible, with a spearman’s rank correlation coefficient of ρ=0.17 and a p-value of 0. SI Fig. 4. Fitting considerations. A: Log-log plot of the probability density versus tissue jump (or avalanche) magnitude. The coefficient was obtained to be τ=1.27 by a linear fit. The power law has an exponential cutoff due to finite system size. B: Semi-log plot for probability density versus tissue jump magnitude. The linear fit is poor as the distribution does not exhibit an exponential signature. 102 SI Fig. 5. Exponent dependence on parameter variation. A: Log-log plot of the probability density versus tissue jump magnitude when mechanosensing amplitude, μ, is one quarter of the original value. The slope coefficient was obtained to be invariant with τ=1.27. B: Log-log plot of the probability density versus tissue jump magnitude when boundary tension is changed to the bulk edge tension. The slope coefficient was obtained to be τ=1.25 and considering the error bar, within the 1.27 value. SI Fig. 6. Effect of apical constriction on MF speed and cell packing. A: The MF speed decreases as apical constriction increases. B: The cellular mosaic becomes increasingly packed with increasing apical constriction. This speed dependence on apical constriction is a consequence of the signaling differential equations, which are discretized on the cellular lattice, and are independent of the physical distance between cells. The blue dots involve data taken from simulations were the MF reached the anterior portion and the growth process was concluded. The red dots indicate data from simulations that were not completed since the tissue reaches a very large number of cells, making them very computationally costly. 103 SI Fig. 7. Boundary tension analysis. A: Circularity index, defined as 𝑐 0 = 4𝜋𝐴 𝑃 −2 , where A is the tissue area and P the perimeter . The tissue became increasingly circular as the boundary tension was increased. The effect of increasing boundary tension is more prominent for low values and saturates ~1 for larger values. Boundary tension value used throughout this chapter was set to 0.7. B: Cell number versus time. The changes in boundary tension do not significantly affect the growth rate of the tissue. Acknowledgements We would like to thank Guillaume Gay, Yehuda Ben-Zion, Konstantin Kozlov and Svetlana Surkova for useful discussions. Research reported in this publication was supported by the National Institutes of Health under award number NIH GM128193. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. Computation for the work described in this chapter was supported by the University of Southern California’s Center for High-Performance Computing (carc.usc.edu). 104 Conclusion During development of biological organisms, multiple complex structures are formed. In many instances, these structures need to exhibit a high degree of order to be functional, although many of their constituents are intrinsically stochastic. Hence, it has been suggested that biological robustness ultimately must rely on complex gene regulatory networks and clean-up mechanisms. Here we explored developmental processes that have evolved inherent robustness against stochasticity. In the context of the Drosophila eye disc, multiple optical units, ommatidia, develop into crystal-like patterns. During the larva-to-pupa stage of metamorphosis, the centers of the ommatidia are specified initially through the diffusion of morphogens, followed by the specification of R8 cells. Establishing the R8 cell is crucial in setting up the geometric, and functional, relationships of cells within an ommatidium and among neighboring ommatidia. Here we studied an PDE mathematical model of these spatio-temporal processes in the presence of parametric stochasticity, defining and applying measures that quantify order within the resulting spatial patterns. We observed a universal sigmoidal response to increasing transcriptional noise. Ordered patterns persist up to a threshold noise level in the model parameters. In accordance with prior qualitative observations, as the noise is further increased past a threshold point of no return, these ordered patterns rapidly become disordered. Such robustness in development allows for 105 the accumulation of genetic variation without any observable changes in phenotype. We argue that the observed sigmoidal dependence introduces robustness allowing for sizable amounts of genetic variation and transcriptional noise to be tolerated in natural populations without resulting in phenotype variation. In the physicists' perspective, epithelial tissues constitute an exotic type of active matter with non-linear properties reminiscent of amorphous materials. In the context of a circular proliferating epithelium, modeled by the quasistatic vertex model, we identified novel discrete tissue scale rearrangements, i.e. cellular flow avalanches, which are a form of collective cell movement. During the avalanches, the cellular trajectories are radial in the periphery and form a vortex in the core. After the onset of these avalanches, the epithelial area grows discontinuously. The avalanches were found to be stochastic, and their strength is determined by the density of cells in the tissue. Overall, avalanches regularize the spatial tension distribution along tissue. Furthermore, the avalanche distribution was found to obey a power law, with an exponent consistent with sheer induced avalanches in amorphous materials. To decipher the role of avalanches in organ development, we simulated epithelial growth of the Drosophila eye disc during the third instar using a computational model, which includes both signaling and mechanistic signaling. During the third instar, the morphogenetic furrow (MF), a ~10 cell 106 wide wave of apical area constriction propagates through the epithelium, making it a system with interesting mechanical properties. These simulations were used to understand the details of the growth process, the effect of the MF on the growth dynamics on the tissue scale, and to make predictions. The avalanches were found to depend on the strength of the apical constriction of cells in the MF, with stronger apical constriction leading to less frequent and more pronounced avalanches. The results herein highlight the dependence of simulated tissue growth dynamics on relaxation timescales, and serve as a guide for in vitro experiments. This thesis identifies novel phenomena, develops novel tools and models to study development and epithelial dynamics but also provides the foreground for further explorations. Firstly, the question arises: How does robustness in development depend on the nature of the noise that is introduced? For instance, one could extend the robustness study of photoreceptor differentiation process to include dynamical noise, as a way to model transcriptional bursting present in developmental systems. Furthermore, a logical next step would be to study the robustness of a collection of simulated developmental processes and to provide a cross-system analysis to assess the universality of the threshold response to noise. Second, the avalanches present during epithelial growth do depend on the nature of the implemented solver. That is, if the dynamics are changed from quasistatic 107 to overdamped viscosity dynamics, the cells will be allowed to displace up to a fixed distance. Since the avalanches involve large cellular displacements, the details of the viscous solver would dampen their onset. Thus, the question arises, where does this transition takes place? Furthermore, one can investigate the effect of fluctuations on the dynamics. Specifically, one can introduce dynamical random noise on the cellular vertices and study how the noise amplitude affects the dynamical evolution of the system. Lastly, a logical extension of the Drosophila eye disc growth model would be to include the gene regulatory network responsible for the patterning of photoreceptors described in Chapter 1. This would provide a novel model system, perfect to investigate the effects of coupling growth and patterning. 108 References 1. Chubb JR, Trcek T, Shenoy SM, Singer RH. Transcriptional Pulsing of a Developmental Gene. Curr Biol. 2006 May 23;16(10):1018–25. 2. Waddington CH. Canalization of Development and Genetic Assimilation of Acquired Characters. Nature. 1959 Jun;183(4676):1654–5. 3. Averbukh I, Gavish A, Shilo B-Z, Barkai N. Dealing with noise: The challenge of buffering biological variability. Curr Opin Syst Biol. 2017 Feb 1;1:69–74. 4. Barkai N, Shilo B-Z. Robust Generation and Decoding of Morphogen Gradients. Cold Spring Harb Perspect Biol. 2009 Nov 1;1(5):a001990. 5. Gursky VV, Surkova SYu, Samsonova MG. Mechanisms of developmental robustness. Biosystems. 2012 Sep 1;109(3):329–35. 6. Keller EF. Developmental Robustness. Ann N Y Acad Sci. 2002;981(1):189– 201. 7. Gavish A, Shwartz A, Weizman A, Schejter E, Shilo B-Z, Barkai N. Periodic patterning of the Drosophila eye is stabilized by the diffusible activator Scabrous. Nat Commun. 2016 Feb 15;7(1):1–10. 8. Lubensky DK, Pennington MW, Shraiman BI, Baker NE. A dynamical model of ommatidial crystal formation. Proc Natl Acad Sci. 2011 Jul 5;108(27):11145–50. 9. Kumar JP. Building an ommatidium one cell at a time. Dev Dyn. 2012;241(1):136–49. 10. Haynie JL, Bryant PJ. Development of the eye-antenna imaginal disc and morphogenesis of the adult head in Drosophila melanogaster. J Exp Zool. 1986;237(3):293–308. 11. Lawrence PA, Green SM. Cell lineage in the developing retina of Drosophila. Dev Biol. 1979 Jul 1;71(1):142–52. 12. Wolff T, Ready DF. Cell death in normal and rough eye mutants of Drosophila. Dev Camb Engl. 1991 Nov;113(3):825–39. 109 13. Greenwood S, Struhl G. Progression of the morphogenetic furrow in the Drosophila eye: the roles of Hedgehog, Decapentaplegic and the Raf pathway. Dev Camb Engl. 1999 Dec;126(24):5795–808. 14. Baonza A, Freeman M. Notch signalling and the initiation of neural development in the Drosophila eye. Dev Camb Engl. 2001 Oct;128(20):3889– 98. 15. Wolff T, Ready DF. The beginning of pattern formation in the Drosophila compound eye: the morphogenetic furrow and the second mitotic wave. Development. 1991 Nov 1;113(3):841–50. 16. Voas MG, Rebay I. Signal integration during development: Insights from the Drosophila eye. Dev Dyn. 2004;229(1):162–75. 17. Cagan RL, Ready DF. Notch is required for successive cell divisions in the developing Drosophila retina. Trends Genet. 1989 Jan 1;5:364. 18. Baker NE, Yu S, Han D. Evolution of proneural atonal expression during distinct regulatory phases in the developing Drosophila eye. Curr Biol. 1996 Oct 1;6(10):1290–302. 19. Gavish A, Barkai N. A two-step patterning process increases the robustness of periodic patterning in the fly eye. J Biol Phys. 2016 Jun;42(3):317. 20. Jiao Y, Lau T, Hatzikirou H, Meyer-Hermann M, Joseph C. Corbo, Torquato S. Avian photoreceptor patterns represent a disordered hyperuniform solution to a multiscale packing problem. Phys Rev E. 2014 Feb 24;89(2):022721. 21. Kim S, Cassidy JJ, Yang B, Carthew RW, Hilgenfeldt S. Hexagonal Patterning of the Insect Compound Eye: Facet Area Variation, Defects, and Disorder. Biophys J. 2016 Dec 20;111(12):2735–46. 22. Kittel C. Introduction to Solid State Physics. Wiley; 2004. 704 p. 23. Egami T, Billinge SJL. Underneath the Bragg Peaks: Structural Analysis of Complex Materials. Elsevier; 2003. 424 p. 24. Truskett TM, Torquato S, Debenedetti PG. Towards a quantification of disorder in materials: Distinguishing equilibrium and glassy sphere packings. Phys Rev E. 2000 Jul 1;62(1):993–1001. 110 25. Halperin BI, Nelson DR. Theory of Two-Dimensional Melting. Phys Rev Lett. 1978 Jul 10;41(2):121–4. 26. Markov AA, Nagorny NM. The Theory of Algorithms [Internet]. Springer Netherlands; 1988 [cited 2021 Mar 2]. (Mathematics and its Applications). Available from: https://www.springer.com/gp/book/9789027727732 27. Debye P. Interferenz von Röntgenstrahlen und Wärmebewegung. Ann Phys. 1913;348(1):49–92. 28. Waller I. Zur Frage der Einwirkung der Wärmebewegung auf die Interferenz von Röntgenstrahlen. Z Für Phys. 1923 Dec 1;17(1):398–408. 29. Lindenmeyer PH, Hosemann R. Application of the Theory of Paracrystals to the Crystal Structure Analysis of Polyacrylonitrile. J Appl Phys. 1963 Jan 1;34(1):42–5. 30. Okabe A, Boots B, Sugihara K, Chiu SN. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. John Wiley & Sons; 2009. 697 p. 31. Rachev ST, Klebanov L, Stoyanov SV, Fabozzi F. The Methods of Distances in the Theory of Probability and Statistics [Internet]. New York: Springer- Verlag; 2013 [cited 2021 Mar 2]. Available from: https://www.springer.com/gp/book/9781461448686 32. Statistical Mechanics of Interacting Systems: Theory of Phase Transitions. In: Principles of Equilibrium Statistical Mechanics [Internet]. John Wiley & Sons, Ltd; 2000 [cited 2021 Mar 2]. p. 253–5. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/3527603158.part3 33. Metzger BPH, Duveau F, Yuan DC, Tryban S, Yang B, Wittkopp PJ. Contrasting Frequencies and Effects of cis- and trans-Regulatory Mutations Affecting Gene Expression. Mol Biol Evol. 2016 May;33(5):1131–46. 34. Metzger BPH, Yuan DC, Gruber JD, Duveau F, Wittkopp PJ. Selection on noise constrains variation in a eukaryotic promoter. Nature. 2015 May;521(7552):344–7. 35. Gibson G, Dworkin I. Uncovering cryptic genetic variation. Nat Rev Genet. 2004 Sep;5(9):681–90. 111 36. Forgacs G, Foty RA, Shafrir Y, Steinberg MS. Viscoelastic Properties of Living Embryonic Tissues: a Quantitative Study. Biophys J. 1998 May 1;74(5):2227–34. 37. Schötz E-M, Lanio M, Talbot JA, Manning ML. Glassy dynamics in three- dimensional embryonic tissues. J R Soc Interface. 2013 Dec 6;10(89):20130726. 38. Angelini TE, Hannezo E, Trepat X, Marquez M, Fredberg JJ, Weitz DA. Glass-like dynamics of collective cell migration. Proc Natl Acad Sci. 2011 Mar 22;108(12):4714–9. 39. Henkes S, Fily Y, Marchetti MC. Active jamming: Self-propelled soft particles at high density. Phys Rev E. 2011 Oct 12;84(4):040301. 40. Berthier L. Nonequilibrium Glassy Dynamics of Self-Propelled Hard Disks. Phys Rev Lett. 2014 Jun 4;112(22):220602. 41. Berthier L, Kurchan J. Non-equilibrium glass transitions in driven and active matter. Nat Phys. 2013 May;9(5):310–4. 42. Bi D, Lopez JH, Schwarz JM, Manning ML. A density-independent glass transition in biological tissues. Nat Phys. 2015 Dec;11(12):1074–9. 43. Ruina A. Slip instability and state variable friction laws. J Geophys Res Solid Earth. 1983;88(B12):10359–70. 44. Hubert A, Schäfer R. Magnetic Domains: The Analysis of Magnetic Microstructures [Internet]. Berlin Heidelberg: Springer-Verlag; 1998 [cited 2021 Jan 8]. Available from: https://www.springer.com/gp/book/9783540641087 45. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: An explanation of the 1/ f noise. Phys Rev Lett. 1987 Jul 27;59(4):381–4. 46. Karimi K, Ferrero EE, Barrat J-L. Inertia and universality of avalanche statistics: The case of slowly deformed amorphous solids. Phys Rev E. 2017 Jan 12;95(1):013003. 47. Salerno KM, Maloney CE, Robbins MO. Avalanches in Strained Amorphous Solids: Does Inertia Destroy Critical Behavior? Phys Rev Lett. 2012 Sep 6;109(10):105703. 112 48. Bohn F, Durin G, Correa MA, Machado NR, Pace RDD, Chesman C, et al. Playing with universality classes of Barkhausen avalanches. Sci Rep. 2018 Jul 26;8(1):1–12. 49. Popovic M, Druelle V, Dye N, Jülicher F, Wyart M. Inferring the flow properties of epithelial tissues from their geometry. New J Phys [Internet]. 2020 [cited 2021 Feb 17]; Available from: http://iopscience.iop.org/article/10.1088/1367-2630/abcbc7 50. Streichan SJ, Hoerner CR, Schneidt T, Holzer D, Hufnagel L. Spatial constraints control cell proliferation in tissues. Proc Natl Acad Sci. 2014 Apr 15;111(15):5586–91. 51. Serra-Picamal X, Conte V, Vincent R, Anon E, Tambe DT, Bazellieres E, et al. Mechanical waves during tissue expansion. Nat Phys. 2012 Aug;8(8):628– 34. 52. Vedula SRK, Leong MC, Lai TL, Hersen P, Kabla AJ, Lim CT, et al. Emerging modes of collective cell migration induced by geometrical constraints. Proc Natl Acad Sci. 2012 Aug 7;109(32):12974–9. 53. Doxzen K, Vedula SRK, Leong MC, Hirata H, Gov NS, Kabla AJ, et al. Guidance of collective cell migration by substrate geometry. Integr Biol. 2013 Aug 22;5(8):1026–35. 54. Puliafito A, Hufnagel L, Neveu P, Streichan S, Sigal A, Fygenson DK, et al. Collective and single cell behavior in epithelial contact inhibition. Proc Natl Acad Sci. 2012 Jan 17;109(3):739–44. 55. Heinrich MA, Alert R, LaChance JM, Zajdel TJ, Košmrlj A, Cohen DJ. Size- dependent patterns of cell proliferation and migration in freely-expanding epithelia [Internet]. eLife. eLife Sciences Publications Limited; 2020 [cited 2021 Jan 15]. Available from: https://elifesciences.org/articles/58945 56. Huergo MAC, Pasquale MA, González PH, Bolzán AE, Arvia AJ. Dynamics and morphology characteristics of cell colonies with radially spreading growth fronts. Phys Rev E. 2011 Aug 11;84(2):021917. 57. Lee RM, Kelley DH, Nordstrom KN, Ouellette NT, Losert W. Quantifying stretching and rearrangement in epithelial sheet migration. New J Phys. 2013 Feb;15(2):025036. 113 58. Eder D, Aegerter C, Basler K. Forces controlling organ growth and size. Mech Dev. 2017 Apr 1;144:53–61. 59. Pan Y, Heemskerk I, Ibar C, Shraiman BI, Irvine KD. Differential growth triggers mechanical feedback that elevates Hippo signaling. Proc Natl Acad Sci. 2016 Nov 8;113(45):E6974–83. 60. Irvine KD, Shraiman BI. Mechanical control of growth: ideas, facts and challenges. Development. 2017 Dec 1;144(23):4238–48. 61. Ehrig S, Schamberger B, Bidan CM, West A, Jacobi C, Lam K, et al. Surface tension determines tissue shape and growth kinetics. Sci Adv. 2019 Sep 1;5(9):eaav9394. 62. Mao Y, Tournier AL, Hoppe A, Kester L, Thompson BJ, Tapon N. Differential proliferation rates generate patterns of mechanical tension that orient tissue growth. EMBO J. 2013 Oct 30;32(21):2790–803. 63. Gibson WT, Veldhuis JH, Rubinstein B, Cartwright HN, Perrimon N, Brodland GW, et al. Control of the Mitotic Cleavage Plane by Local Epithelial Topology. Cell. 2011 Feb 4;144(3):427–38. 64. Aegerter-Wilmsen T, Smith AC, Christen AJ, Aegerter CM, Hafen E, Basler K. Exploring the effects of mechanical feedback on epithelial topology. Development. 2010 Feb 1;137(3):499–506. 65. Aegerter-Wilmsen T, Heimlicher MB, Smith AC, de Reuille PB, Smith RS, Aegerter CM, et al. Integrating force-sensing and signaling pathways in a model for the regulation of wing imaginal disc size. Development. 2012 Sep 1;139(17):3221–31. 66. Salbreux G, Barthel LK, Raymond PA, Lubensky DK. Coupling Mechanical Deformations and Planar Cell Polarity to Create Regular Patterns in the Zebrafish Retina. PLOS Comput Biol. 2012 Aug 23;8(8):e1002618. 67. Farhadifar R, Röper J-C, Aigouy B, Eaton S, Jülicher F. The Influence of Cell Mechanics, Cell-Cell Interactions, and Proliferation on Epithelial Packing. Curr Biol. 2007 Dec;17(24):2095–104. 68. Alt Silvanus, Ganguly Poulami, Salbreux Guillaume. Vertex models: from cell mechanics to tissue morphogenesis. Philos Trans R Soc B Biol Sci. 2017 May 19;372(1720):20150520. 114 69. Fletcher AG, Osterfield M, Baker RE, Shvartsman SY. Vertex Models of Epithelial Morphogenesis. Biophys J. 2014 Jun;106(11):2291–304. 70. Barton DL, Henkes S, Weijer CJ, Sknepnek R. Active Vertex Model for cell- resolution description of epithelial tissue mechanics. PLOS Comput Biol. 2017 Jun 30;13(6):e1005569. 71. Monier B, Gettings M, Gay G, Mangeat T, Schott S, Guarner A, et al. Apico- basal forces exerted by apoptotic cells drive epithelium folding. Nature. 2015 Feb;518(7538):245–8. 72. Lewis FT. The correlation between cell division and the shapes and sizes of prismatic cells in the epidermis of cucumis. Anat Rec. 1928;38(3):341–76. 73. Heller D, Hoppe A, Restrepo S, Gatti L, Tournier AL, Tapon N, et al. EpiTools: An Open-Source Image Analysis Toolkit for Quantifying Epithelial Growth Dynamics. Dev Cell. 2016 Jan 11;36(1):103–16. 74. Sánchez‐Gutiérrez D, Tozluoglu M, Barry JD, Pascual A, Mao Y, Escudero LM. Fundamental physical cellular constraints drive self‐organization of tissues. EMBO J. 2016 Jan 4;35(1):77–88. 75. Escudero LM, Costa L da F, Kicheva A, Briscoe J, Freeman M, Babu MM. Epithelial organisation revealed by a network of cellular contacts. Nat Commun. 2011 Nov 8;2(1):1–7. 76. Patel AB, Gibson WT, Gibson MC, Nagpal R. Modeling and Inferring Cleavage Patterns in Proliferating Epithelia. PLOS Comput Biol. 2009 Jun 12;5(6):e1000412. 77. Kumar JP. My what big eyes you have: How the Drosophila retina grows. Dev Neurobiol. 2011 Dec;71(12):1133–52. 78. Kumar JP. The fly eye: Through the looking glass. Dev Dyn. 2018;247(1):111–23. 79. Vollmer J, Fried P, Sánchez-Aragón M, Lopes CS, Casares F, Iber D. A quantitative analysis of growth control in the Drosophila eye disc. Development. 2016 May 1;143(9):1482–90. 80. Treisman JE. Retinal differentiation in Drosophila. Wiley Interdiscip Rev Dev Biol. 2013;2(4):545–57. 115 81. Baker NE. Cell proliferation, survival, and death in the Drosophila eye. Semin Cell Dev Biol. 2001 Dec;12(6):499–507. 82. Crossman SH, Streichan SJ, Vincent J-P. EGFR signaling coordinates patterning with cell survival during Drosophila epidermal development. PLOS Biol. 2018 Oct 31;16(10):e3000027. 83. Vollmer J, Fried P, Aguilar-Hidalgo D, Sánchez-Aragón M, Iannini A, Casares F, et al. Growth control in the Drosophila eye disc by the cytokine Unpaired. Development. 2017 Mar 1;144(5):837–43. 84. Curtiss J, Mlodzik M. Morphogenetic furrow initiation and progression during eye development in Drosophila: the roles of decapentaplegic, hedgehog and eyes absent. Development. 2000;127(6):1325–36. 85. Domínguez M, Hafen E. Hedgehog directly controls initiation and propagation of retinal differentiation in the Drosophila eye. Genes Dev. 1997 Dec 1;11(23):3254–64. 86. Heberlein U, Moses K. Mechanisms of Drosophila retinal morphogenesis: the virtues of being progressive. Cell. 1995 Jun 30;81(7):987–90. 87. Pichaud F, Casares F. homothorax and iroquois-C genes are required for the establishment of territories within the developing eye disc. Mech Dev. 2000 Aug;96(1):15–25. 88. Roignant J-Y, Treisman JE. Pattern formation in the Drosophila eye disc. Int J Dev Biol. 2009 May 22;53(5–6):795–804. 89. Brown NL, Sattler CA, Paddock SW, Carroll SB. Hairy and Emc negatively regulate morphogenetic furrow progression in the drosophila eye. Cell. 1995 Mar 24;80(6):879–87. 90. Fried P, Sánchez-Aragón M, Aguilar-Hidalgo D, Lehtinen B, Casares F, Iber D. A Model of the Spatio-temporal Dynamics of Drosophila Eye Disc Development. Thieffry D, editor. PLOS Comput Biol. 2016 Sep 14;12(9):e1005052. 91. Okuda S, Miura T, Inoue Y, Adachi T, Eiraku M. Combining Turing and 3D vertex models reproduces autonomous multicellular morphogenesis with undulation, tubulation, and branching. Sci Rep. 2018 Feb 5;8(1):2386. 116 92. Schilling S, Willecke M, Aegerter-Wilmsen T, Cirpka OA, Basler K, Mering C von. Cell-Sorting at the A/P Boundary in the Drosophila Wing Primordium: A Computational Model to Consolidate Observed Non-Local Effects of Hh Signaling. PLOS Comput Biol. 2011 Apr 7;7(4):e1002025. 93. Smith AM, Baker RE, Kay D, Maini PK. Incorporating chemical signalling factors into cell-based models of growing epithelial tissues. J Math Biol Heidelb. 2012 Sep;65(3):441–63. 94. Wartlick O, Mumcu P, Kicheva A, Bittig T, Seum C, Jülicher F, et al. Dynamics of Dpp Signaling and Proliferation Control. Science. 2011 Mar 4;331(6021):1154–9. 95. Talamali M, Petäjä V, Vandembroucq D, Roux S. Avalanches, precursors, and finite-size fluctuations in a mesoscopic model of amorphous plasticity. Phys Rev E. 2011 Jul 29;84(1):016115. 96. Fodor É, Mehandia V, Comelles J, Thiagarajan R, Gov NS, Visco P, et al. Spatial Fluctuations at Vertices of Epithelial Layers: Quantification of Regulation by Rho Pathway. Biophys J. 2018 Feb 27;114(4):939–46. 97. Curran S, Strandkvist C, Bathmann J, de Gennes M, Kabla A, Salbreux G, et al. Myosin II Controls Junction Fluctuations to Guide Epithelial Tissue Ordering. Dev Cell. 2017 Nov 20;43(4):480-492.e6. 98. Stavans J. The evolution of cellular structures. Rep Prog Phys. 1993 Jun 1;56(6):733–89. 99. Sussman DM, Schwarz JM, Marchetti MC, Manning ML. Soft yet Sharp Interfaces in a Vertex Model of Confluent Tissue. Phys Rev Lett. 2018 Jan 29;120(5):058001. 100. Dahmann C, Oates AC, Brand M. Boundary formation and maintenance in tissue development. Nat Rev Genet. 2011 Jan;12(1):43–55. 101. Fabero JC, Bautista A, Casasús L. An explicit finite differences scheme over hexagonal tessellation. Appl Math Lett. 2001 Jan 1;14(5):593–8.

Abstract (if available)

Abstract This thesis is focused on computational modeling and quantitative analysis of patterning and growth in biological systems. Patterning is studied in the context of the Drosophila eye disc, an epithelial tissue that contains ~800 eye units arranged crystalline order, with an emphasis on photoreceptor differentiation. Coupled ODE equations are parametrically perturbed to determine and quantify the robustness of the developmental process. Generalization of the translational order parameter, utilizing the notion of probability distance, is used to quantify the transition of the pattern from order to disorder. The transition is found to exhibit a ubiquitous threshold response to noise, with important implications on cryptic variation and robustness of gene regulatory networks. Growth is studied in the context of the vertex model, in which a uniform growth of an epithelial monolayer and the Drosophila eye disc are simulated. During the growth process, novel critical behavior is identified, where the tissue grows in bursts. The magnitude of growth bursts, or avalanches, are shown to be distributed according to power law distributions and are found to match the critical behavior of stress avalanches in sheer induced amorphous materials. The avalanches are also found to be accompanied with collective cellular motion and vortex formation. The eye disc growth model dynamics are used to investigate correlations between chemical and signaling components. The onset of avalanches, signaling properties and global cellular packing are found to depend on the coordinated apical constriction present in the developmental process.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Avalanches during epithelial tissue growth; uniform growth and a drosophila eye disc model

PDF

Robustness and stochasticity in Drosophila development

PDF

Microbial interaction networks: from single cells to collective behavior

PDF

Metabolomic and proteomic approaches to understanding senescence and aging in mammary epithelial cells

PDF

Adaptive patterning during skin formation: assembly of vasculature and adipose tissue

PDF

Exploring the molecular and cellular underpinnings of organ polarization using feather as the model system

PDF

Ancestral inference and cancer stem cell dynamics in colorectal tumors

PDF

Physical principles of membrane mechanics, membrane domain formation, and cellular signal transduction

PDF

Structure – dynamics – function analysis of class A GPCRs and exploration of chemical space using integrative computational approaches

PDF

Developmental trajectories of sensory patterns in young children with and without autism spectrum disorder: a longitudinal population-based study from infancy to school age

PDF

Signaling mechanisms governing intestinal regeneration and gut-glia cross-talk in Drosophila

Asset Metadata

Creator Courcoubetis, George (author)

Core Title Patterning and growth in biological systems: a computational exploration of biological robustness, epithelial growth dynamics and development

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 04/14/2021

Defense Date 03/16/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag avalanche,biophysics,compound eye,Computational Biology,developmental biology,Drosophila eye disc,epithelial growth,morphogenesis,OAI-PMH Harvest,patterning,self-organized criticality,spatial order metrics,vortex formation

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haas, Stephan (committee chair), Boedicker, James (committee member), Haselwandter, Christoph (committee member), Morsut, Leonardo (committee member), Nuzdhin, Sergey (committee member)

Creator Email courcoub@usc.edu,gcourcou@hotmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-443180

Unique identifier UC11666653

Identifier etd-Courcoubet-9457.pdf (filename),usctheses-c89-443180 (legacy record id)

Legacy Identifier etd-Courcoubet-9457.pdf

Dmrecord 443180

Document Type Dissertation

Rights Courcoubetis, George

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA