article: A major problem in biology is to understand how complex tissue shapes may arise through growth. In many cases this process involves preferential growth along particular orientations raising the question of how these orientations are specified. One view is that orientations are specified through stresses in the tissue (axiality-based system). Another possibility is that orientations can be specified independently of stresses through molecular signalling (polarity-based system). The axiality-based system has recently been explored through computational modelling. Here we develop and apply a polarity-based system which we call the Growing Polarised Tissue (GPT) framework. Tissue is treated as a continuous material within which regionally expressed factors under genetic control may interact and propagate. Polarity is established by signals that propagate through the tissue and is anchored in regions termed tissue polarity organisers that are also under genetic control. Rates of growth parallel or perpendicular to the local polarity may then be specified through a regulatory network. The resulting growth depends on how specified growth patterns interact within the constraints of mechanically connected tissue. This constraint leads to the emergence of features such as curvature that were not directly specified by the regulatory networks. Resultant growth feeds back to influence spatial arrangements and local orientations of tissue, allowing complex shapes to emerge from simple rules. Moreover, asymmetries may emerge through interactions between polarity fields. We illustrate the value of the GPT-framework for understanding morphogenesis by applying it to a growing Snapdragon flower and indicate how the underlying hypotheses may be tested by computational simulation. We propose that combinatorial intractions between orientations and rates of growth, which are a key feature of polarity-based systems, have been exploited during evolution to generate a range of observed biological shapes.
  Although there have been many experimental and theoretical studies on patterns of gene activities and their establishment in animals and plants [1]–[6] much less is known about how patterns of activity are linked to tissue growth and deformation. Addressing this problem represents a challenge because final form is usually not a direct readout of locally specified properties, but depends on mechanical constraints from neighbouring regions. For example, if the margin of a leaf has a higher specified growth rate than the centre, a wavy edge will emerge. The wavy edge is not directly specified but is a feature that emerges through the interaction between patterns of specified growth and the mechanical constraints of tissue continuity [7]. In such cases we may distinguish between specified growth, which is the growth that would be attained if each region grew independently of its neighbours (i. e. in mechanical isolation), and resultant growth, which is the growth observed when mechanical constraints of neighbours are taken into account (i. e. mechanically connected tissue). In animal systems a similar distinction is made between an imposed active deformation, and an elastic passive deformation [8]. Resultant growth can be measured experimentally by tracking tissue deformations over time [9]–[13]. However, to understand the mechanisms by which resultant growth arises we need to know how genes influence specified growth. Where specified growth is isotropic, genes need to control a single parameter, the local rate of growth. However, in many cases specified growth may be anisotropic requiring orientations as well as rates of growth to be under genetic control. Controlling orientations of growth requires a local axis to be defined (i. e. axiality, represented as a field of lines). In this respect growth is similar to stress which also has axiality. This similarity has led to the suggestion that stresses provide the primary cues for orienting growth. According to such a stress-based axiality mechanism, gene activity influences stresses in the tissue, the orientations of which are transduced to influence molecular properties of cells such as the cytoskeleton. These in turn modulate growth orientations which may further feed back to influence the pattern of stresses [14]–[18]. Recent support for such mechanisms in plants have come from studies of the effect of stresses on microtubule patterns [19]. A different way of specifying orientations of growth is through differential concentrations of signalling molecules. The varying concentrations define a local (cellular) polarity which includes both axiality and directional components (represented by a field of arrows). The axiality component is then used to orient growth. In this polarity-based axiality system, genes influence the distribution of signalling molecules which define a coordinated field of polarities. Incorporation of mechanical constraints then leads to resultant growth, which may feed back to influence, for example, tissue polarity orientations. In support of this system, there is considerable evidence that polarity is prevalent in biological tissues and may modulate growth [20]–[24]. For example, planar cell polarity (PCP) systems have been described in animals and implicated in processes such as growth of wings in Drosophila and convergent-extension in vertebrates [25]–[27]. Similarly, the polarised distribution of auxin transporters (PIN molecules) has been shown to be important for outgrowths of primordia in plants [28]. Elements of the stress-based axiality system have recently been modelled [29], [30]. Here we describe a framework and software implementation for the alternative polarity-based axiality approach, which we call GPT-framework. ur software, called GFtbox, is a MATLAB application available from http: //www. uea. ac. uk/cmp/research/cmpbio/GFtbox. This framework was developed with plant growth in mind, although it may also be useful for modelling animal systems where cell movement is limited. In accompanying papers we show how a biologically relevant model can be derived using the GPT-framework [31], [32]. This model of Snapdragon flower development is constrained by a range of experimental data including gene expression patterns, mutant phenotypes, clonal analysis, growth dynamics and changes in geometry. It provides a working hypothesis for how growth is specified and shows how reorientation of growth can account for key observations. In this paper we explore a series of simplified models which illustrate how growth and polarity may interact combinatorially during morphogenesis to generate a wide range of forms. The results highlight the value of being able to specify orientation independently of stresses in the generation of complex tissue shapes. In addition, we provide the theoretical foundations on which our modelling depends.
  Modelling the genetic control of tissue growth requires the incorporation of gene regulatory networks and signal propogation within a growing, mechanically connected, tissue. In the GPT-framework, tissue is treated as a continuous sheet of material with two surfaces and a thickness, here termed the canvas. Biologically, the canvas may correspond to a sheet of cells, single cells or subcellular components (e. g. walls). Regulatory factors are distributed over the canvas and may interact and propagate, allowing particular patterns and local polarities to be specified. Regulatory factors can be classified into two types. Identity factors do not propagate within the canvas, while signalling factors can. The regulatory factors specify a growth tensor field which describes the specified rates of growth parallel and perpendicular to the local polarity. Elasticity theory is used to compute the resultant deformation of the canvas. This deformation modifies the relationships within the canvas and thus feeds back to influence the regulatory factors. Our implementation (GFtbox) is specialised towards tissues that grow as sheets, such as petals or leaves, but the basic concepts are also applicable to bulk three-dimensional and flat two-dimensional tissues. In the results we study the interactions between tissue polarity and differential growth in the generation of shape through a series of models. For convenience each example has a setup phase during which the shape of the initial canvas and distribution of regional identities and signalling factors is established, and three components that form the model. (1) A Polariser Regulatory Network (PRN) controls the activity of various organisers from which tissue polarity information propagates. There are two types of organiser, termed organiser and organiser. As a convention, we show polarity pointing away from organisers, and towards organisers. Polarity propagation is implemented through a signalling factor called POLARISER (POL), the gradient of which defines local polarity. The PRN controls production and degradation of POL at organisers that anchor the polarity. POL may also be produced and degraded at a background rate throughout the canvas. (2) A gene regulatory network (GRN) controls the activity of identity or signalling factors encoded by genes. (3) A growth rate regulatory network (KRN) determines how identity or signalling factors influence specified growth rates parallel to, , and perpendicular to, , local polarity. The KRN also specifies the growth in thickness, . The specified growth rates for a region of the canvas are equivalent to the growth that would arise without the constraints of surrounding material (see Methods). In the first time step the specified growth field is applied to the initial canvas which may then distort through mechanical interactions in the continuous material (modelled according to elasticity theory, see Methods). The result is a slight deformation of the canvas (resultant growth field) that takes the regions of identity factors with it. Where a region containing an identity factor expands, that region inherits the properties of the parent region, so maintaining boundaries. In such new volumes, the concentrations of signalling factors are interpolated from the parent surrounding regions and then further adjusted according to their production, dilution, propagation and decay rates. The deformed canvas and expression pattern provides the starting point for the next time step and the sequence is reiterated. To verify the computational correctness of GFtbox, results were computed for several situations where analytical solutions are possible (see Text S1). In the following we explore combinatorial interactions between polarity and growth through a series of simple cases. We first consider deformations in 2D. A key feature of the polarity-based axiality system is that orientation and growth rates can be specified independently and then combined in various ways. This combinatorial aspect is unlike the stress-based axiality system where orientations can only be specified once stresses have been generated in the tissue. These stresses will depend on the pattern of specified growth rates and the geometry of the tissue. To illustrate the combinatorial interactions within a polarity-based axiality system we first model simple anisotropic growth (Case A) and differential isotropic growth (Case B) separately. We then combine them in different ways (Cases C-I). We use an initially square canvas marked with black discs (simulating cells that produce marked clones) and a grid to show the geometrical transformations [33]. In all Cases the total areal increase (accumulated growth) is the same. The state of the canvas before and after growth is illustrated in Figure 1 for each Case. Case A: Uniform polarity field with spatially uniform anisotropic specified growth rates. A gradient of POL is established during the setup phase through two organisers, (organiser) and (organiser) at the bottom and top boundaries respectively. The PRN involves promoting production of POL while promotes degradation of POL, forming a proximodistal gradient of POL (arrows). After the setup phase the POL gradient is frozen (fixed to the canvas so that the gradient deforms with the canvas). An identity factor is expressed uniformly throughout the canvas. The KRN is (the value of is indicated by the intensity of orange). The resultant growth transforms the square into a vertically stretched rectangle. The black discs become vertically oriented ellipses. The specified growth pattern underlying this transformation is straightforward to implement using the polarity-based axiality system. By contrast, an stress-based axiality system would require an additional step that generates vertically oriented stresses and thus an additional deformation. Moreover, the pattern of stresses would need to be maintained during growth unless there was a mechanism for fixing the axiality. Case B: Spatially varying isotropic specified growth rates. Differential growth is achieved by promotion of specified growth rates towards the right side of the square. This involves establishing an identity factor during the setup phase that is most strongly expressed along the right edge from where it declines gradually. The KRN involves promoting the specified growth rates equally in all directions (,). This leads to a gradient of locally isotropic specified growth that increases from left to right. The overall result is a curved fan. Curvature is not directly specified but arises through differential growth and mechanical constraints inherent in the canvas. Case C: A combination of Cases A and B: uniform polarity field with spatially varying anisotropic specified growth rates. The PRN and KRN are the same as in Case A while the pattern of is the same as in Case B. That is, specified growth rate is oriented parallel to the POL gradient and increases towards the right. The result is a convex fan with much stronger curvature than Case B. Thus anisotropic specified growth, which on its own produces no curvature (Case A), reinforces the curvature arising through differential growth. In principle this reinforcement may arise from two causes. 1) Because there is no, the gradient in is greater than in Case B. 2) Because polarity is local, the directions of specified growth rotates with the canvas, enhancing curvature. To separate the contributions of these two components, we fix the direction of specified growth by using an external (global) frame of reference, as shown in Case D. Case D: A combination of Cases A and B but using an external field to specify growth orientations. The gradient of POL is determined by an external frame of reference (y axis) instead of being embedded in the tissue. Biologically, external polarity information could be obtained from, for example, the effect of gravity. The result is a fan with reduced curvature compared to Case C. Note that ellipse orientations still deviate from the vertical because, even though growth is specified to be vertical, at each step mechanical constraints force the canvas to curve. The enhanced curvature of Case C over Case D reveals the contribution of orientations being specified internally (2) rather than externally (2). Another way of reducing curvature is by using a local polarity field that re-adjusts dynamically as the structure grows, as will be shown in Case E. Case E: The same model as Case C but allowing POL to continue diffusing rather than being frozen after the setup phase. As with Case D, the resulting curvature is less than Case C, particularly near the extremities. This is because growth orientations turn less near the extreme positions of the canvas. The previous Cases considered uniform polarity fields with differential growth. This raises the question of how non-uniform polarity fields may influence shape. We first consider these when combined with uniform growth rates (Cases F and G, Figure 2) and then with differential growth rates (Cases H and I, Figure 2). Case F: Similar to Case A but setting a spatially varying polarity field. A gradient of POL is established during the setup phase by, which is expressed along the horizontal midline, and, which is expressed in the top, bottom, and right edges, increasing toward the right corners. The resulting POL gradient is shown by the arrows. The polarity field is frozen (fixed to the canvas) after the setup phase. The distribution of is spatially uniform as in Case A. Growth at the top and bottom edges is oriented by the organisers producing a strongly concave right edge. Thus, curvature is generated as a result of non-uniform specified orientations of growth. The curvature is even stronger if the polarity field is not frozen after the setup period, as shown in Case G. Case G: Similar to Case F but allowing POL to continue diffusing rather than being frozen after the setup phase. The result is more concave than Case F. This is because of feedback between canvas geometry and the polarity field. We now look at the effect of introducing differential growth rates. Case H: Similar to Case F but with a gradient of specified growth rate. The PRN is the same as Case F leading to a polarity field pointing to the right corners. The KRN and distribution of are the same as Case C leading to increasing values of towards the right edge. The result is intermediate between Case C and Case F because the diagonally specified growth orientations counteract the curvature induced by differential growth. Thus, unlike Case C where local specification of orientation reinforces tissue curvature, here it antagonises curvature. This effect is still stronger when the POL gradient is not frozen as shown in Case I. Case I: Similar to Case H but allowing POL to continue diffusing rather than being frozen after the setup phase. The right edge grows to be almost vertical showing that an appropriate specified local polarity can antagonise curvature arising from differential growth (Case B). The main conclusion to emerge from Cases A to I is that the ability to combine specified growth rates with separately specified orientations provides an effective control mechanism for generating shape transformations. The shapes that emerge reflect interactions between specified orientations, differential growth and mechanical constraints. Depending on the spatial distribution of organisers and the dynamics of polarity propagation, tissue polarity can reinforce or antagonise curvatures resulting from differential growth or may generate curvature even in the context of uniform growth. So far we have only considered combinatorial interactions within the plane of the canvas. We next consider deformations out of the plane. We again consider a series of simplified cases (Figures 3 and 4) in which polarity and differential growth are treated separately (Cases J, K, O) and in combination (Cases L, M, N, P, Q). In each Case the up-down symmetry is broken by the centre of the initial canvas being slightly bowed upward. To simulate the presence of tissue beyond the boundaries of the initial canvas, the edges of tissue are prevented from moving or rotating out of the plane. Case J: Spatially varying specified orientation with a uniform areal growth rate. The PRN involves an organiser (), expressed in the middle of the canvas (blue). An outer region is defined by which keeps POL levels at zero. This leads to a divergent polarity field near the centre (arrows). POL continues to diffuse after the setup phase (i. e. the gradient is not frozen). is spatially uniform as in Case A. The KRN involves anisotropic growth in the polarised region (and). By default, growth is isotropic where the POL gradient is zero (). The result is a small spike. As with Case F tissue curvature has arisen through variations in specified growth direction even when areal growth rate is uniform. However, in Case J the curvature occurs out of the plane as well as in the plane. Case K: Spatially varying isotropic specified growth rates. During the setup phase is established in the centre of the canvas from where it declines in a graded fashion. As with Case B, the KRN has setting the specified growth rates, , . This leads to a gradient of locally isotropic specified growth rate that increases towards the centre. The result is a puffball-like central bulge exhibiting curvature both in and out of the plane of the canvas. The rounder shape compared to Case K illustrates the limitations of isotropic specified growth in creating elongated outgrowths. However, by combinging Cases J and K the outgrowth can be further exaggerated as shown in Case L. Case L: Spatially varying anisotropic specified growth rates (combining Case J and K). The PRN and KRN are the same as Case J leading to radially directed growth. The distribution of is the same as Case K leading to increased anisotropic specified growth towards the centre. The result is a tall central spike with a sharp tip showing how differential growth and anisotropy act in combination. In many biological structures, such as a growing plant apex, protrusions have rounded tips rather than sharp points. This can be achieved by reducing growth in the central region, as shown in Case M. Case M: Spatially varying anisotropic specified growth rates with a central region of no growth. This is similar to Case L, except that additional identity factor sets and to zero in a small central region. The final shape is a rounded projection similar to what might be observed in a plant apex. Such a model is also consistent with the observation that growth rates tend to be lower in the central region of plant apices [12]. We conclude that a range of outgrowths can be readily obtained by combining specified growth rates and orientations. As for the 2D cases, deformations lead to changes in orientations of the polarity field which feed back to influence further deformations. So far we have considered the effects of uniform and divergent polarity fields. A further elaboration is to combine these two as illustrated in Figure 4. Case N: A uniform polarity field with spatially varying anisotropic specified growth rates (combination of Cases A and K). The PRN and KRN are the same as in Case A leading to a left-right polarity field and anisotropic growth. The pattern of is the same as in Case K leading to enhanced growth in the centre. POL continues to diffuse after setup. The result is a thin bulge with grooves at each end. As with Case C, the polarity field is acting as a modulator rather than generator of curvature (no curvature is produced by the polarity field when combined with uniform anisotropic growth, Case A). We next look at the effect of combining the polarity fields in Cases A and J. Case O: Interacting polarity fields with spatially uniform anisotropic specified growth rates (combination of Cases A and J). The PRN and KRN are the same as Case A except that the additional organiser from Case J is included. The new organiser distorts the polarity field shown in Case N inducing a saddle point upstream. As result following growth, the canvas widens slightly in the centre and forms a central ripple. Thus, as with Cases F and J, some curvature arises even with uniform specified areal growth rates. Next we combine this polarity field with centrally increased specified growth rates. Case P: A combined polarity field with spatially varying anisotropic specified growth rates (combination of Cases K and O). The PRN and KRN are the same as Case O while the pattern of is the same as Case K. The result is an asymmetric spur reflecting the interactions between tissue polarity and growth. The asymmetry arises because the POL gradient generated by the central organiser flows in the same direction as the background POL gradient on one side but in the opposite direction on the other, creating a region of counterflow (arrowed). Disorganisation of growth in the counterflow region reduces growth along the main axis of the tissue. Asymmetry induced in this way is a feature of simple polarity-based axiality systems that would not occur in simple stress-based axiality systems. The orientation of the spur can be reversed by using a organiser instead of a organiser in the centre as shown in Case Q. Case Q: A combined polarity field with spatially varying anisotropic specified growth rates. The PRN and KRN and the pattern of are the same as Case P except that the central organiser is replaced by a organiser. This time the asymmetric spur points in the opposite direction to Case P because the counterflow region is on the other side. We conclude that combining polarity fields provides a further richness by generating asymmetries. The above Cases illustrate some basic combinatorial interactions between polarity and growth. To see how the same principles may apply to a biological example, we consider a simplified model of the Snapdragon corolla tube. To simplify the Snapdragon tube we assume the initial canvas comprises an initial cylindrical canvas closed at one end. As a first step we study locally isotropic specified growth (Case R) and then explore the effect of introducing specified anisotropic growth (Cases S and T). Case R: Spatially varying isotropic specified growth rates. An early step in the development of the Snapdragon flower is arching over of the tube through differential growth. We simplify this process by restricting growth rates in opposite regions of the cylinder and also at the base. This is achieved by having a general background level of which is inhibited in the base by and is also inhibited by a diffusing signal which is generated along opposite sides of the cylinder by. In this Case, specified growth is isotropic, . The result is a ballooned out bowl (Figure 5B, C) rather than an arched over tube. Some areas of the canvas near the base show anisotropic resultant growth, evident from elongated ellipses. This is shown more clearly Figure 5 C where the principal directions of resultant growth () are shown with short lines and the rate of anisotropic growth () is shown in magenta. As with curvature, resultant anisotropy is not specified directly but arises through the interaction between differential growth and mechanical constraints. However, the pattern and extent of resultant anisotropy is inconsistent with experimental observations of clones in the Snapdragon tube, which are highly elongated along the proximodistal axis [31]. To address this discrepency we introduce specified anisotropic growth through a polarity field as shown in Case S. Case S: A uniform polarity field with spatially varying anisotropic specified growth rates. A gradient of POL is established through two organisers, and - at the base and rim respectively (Figure 5A, arrows). The KRN is the same as Case R except that the specified growth rate is now anisotropic, and. Compared to the output from Case R, the sides of the cylinder curve towards each other rather than ballooning outwards (Figure 6 B). Thus, introducing specified anisotropy has a major effect, leading to a more closed shape. It also generates much more elongated clones matching experimental observations. However, continuation of growth leads to the two sides arching further (Figure 6 C) rather than creating the elongated shape that is observed experimentally. (In our implementation which does not currently include collision detection the two sides grow through each other. For clarity we therefore only show one side in Figure 6.) To address this discrepancy, we exploit the potential to reorient growth within the GPT-framework by modulating the polarity field as shown in Case T. Case T: Initially the same as Case S followed by reorientation of tissue polarity. There are two phases of growth, early and late. During the early phase the cylinder grows as in Case S. At the start of the late phase, the polarity field is modulated by restricting the spatial region of the organiser. This is achieved by activating an identity factor in the lateral regions of the cylinder which inhibits, restricting the distal organisers (cyan) to small regions at the apex of each arch (Figure 6D, E). The reoriention of polarity leads to vertical elongation of the arch rims, maintaining the closed shape, rather than the sides continuing to arch over. This captures an essential feature of Snapdragon corolla tube growth.
  We model growth through the accumulation of a series of small deformations of the tissue (canvas). Stresses are generated during the process as the canvas is mechanically interconnected. This may lead to anisotropic resultant growth even when growth is specified to be isotropic (e. g. Case R). In principle, such resultant stresses could be used, through stress-based axiality, to orient all forms of anisotropic growth. However, this would mean that specified orientations of growth would be dependent on differential rates of growth, precluding the possibility of independent control. By contrast, we show how a polarity-based axiality system allows diverse forms to be generated through combinatorial interactions between specified orientations and rates of growth. In this system, a key element in controlling growth orientations is the distribution of polarity organisers. These are of two types, or, allowing polarity fields to be anchored at both ends. Even when specified anisotropic growth is uniform over the canvas, a range of forms can be generated by varying the pattern of organisers. For example, starting from an initial square canvas it is possible to generate rectangles (Case A), concavities (Case F), small spikes (Case J) and ripples (Case O). In these Cases polarity was fixed after a setup period. Biologically, this would correspond to an initial period when polarity propagates across the tissue (when the tissue is small), followed by polarised cells maintaining their polarity and passing it on to their daughters. Another possibility is that polarity continues to propagate during growth leading to slight modifications of the resulting shape (compare Cases F and G). The range of shapes may be greatly extended by combining polarity fields with differential growth rates. For example, tissue polarity may reinforce or antagonise curvature arising through differential growth (Cases C and I). Both aspects are incorporated into the growing Snapdragon tube - reinforcement of curvature during the early phase leading to arching over (Case S), followed by antagonism of curvature leading to straightening (Cases T). It is also straightforward to generate extended outgrowths and apices by combining a single organiser with enhanced growth (Cases L and M). A further feature of polarity-based axiality systems is the emergence of asymmetries through interactions between polarity fields. For example, asymmetric spurs may arise because of counterflowing polarity on one side (Cases P and Q). The asymmetry of the outgrowths in these Cases results from the underlying polarity interactions and would not have arisen from a simple system with only stress-based axiality. In these examples only a few organisers are needed to achieve major shape transformations. To test whether the same simplicity might underly more complex biological transformations, we modelled growth of the Snapdragon flower [31]. This model is constrained by a range of experimental data. The expression pattern of the genes DIV, CYC, DICH and RAD are set according to experimental observations. The model has to not only account for the wild-type phenotype but also double (cyc, dich) and triple (cyc, dich, div) mutants. The model is also constrained by the observed changes in 3D shape determined by optical projection tomography at several developmental stages. In addition the pattern of growth rates and directions over each model petal need to be similar to those observed by clonal analysis. The model starts with an initial cylindrical canvas with five lobes and a proximodistal pattern of polarities established through two polarity organisers (and) (Figure 7 A). During the early phase of growth the ventral region of the tube arches over through differential anisotropic growth. To account for the observed pattern of clones a third organiser () is introduced (Figure 7 C). In the absence of this organiser the tube bulges out (Figure 7 F) similar to what happens in the simplified corolla with no reorientation of growth (Case S). However, with the introduction of the organiser the tube automatically straightens out during later stages, consistent with experimental observations. Thus, this biologically relevant case provides evidence for three organisers underlying major shape transformations and growth dynamics. In the Snapdragon model, the reorientation of growth is under the control of DIV, a gene that encodes a Myb-like transcription factor that affects flower shape and symmetry [31]. As well as its effect on organiser activity, DIV also influences growth rates. Thus, although rates and orientations of growth are specified separately in the model they can be regulated by a common gene. The polarity-based axiality system has the flexibility to account for global shape changes, observed growth patterns and clones without invoking large numbers of polarity organisers. This alone does not demonstrate the validity of invoking tissue polarity for the control of growth orientations. Nevertheless, tissue polarity is commonly observed in animals, for example, polarised cell movements [27] and in plants where the polar distribution of molecules within cells, such as PIN auxin transporters, suggests that cell polarity is also common [23], [24]. It has also been proposed that an auxin concentration maximum at a vascular boundary in the root tip establishes a distal polarity organiser in the root [20]. The GPT-framework allows hypotheses on polarity-based axiality growth to be established that can be subjected to further tests such as mechanical or genetic perturbations. The Snapdragon model, for example, was evaluated against predictions of shapes of multiple mutants not used to build the model [31], [32]. The results showed a good, quantitative, fit between predicted and observed shapes. The model also makes important predictions about the location of polarity organisers. Polarity markers are predicted to show reversals (i. e. arrows pointing away or towards each other) at these locations. In all our Cases we make the simplifying assumption that the tissue is linear over small deformations and has isotropic material properties. An elaboration of the GPT-framework would be to incorporate non-uniform properties, although this would also require these properties to be measured across the tissue during growth. The GPT-framework is consistent with current hypotheses regarding the mechanisms in which plant tissue grows under turgor pressure through the loosening and formation of bonds (Theorems 1 and 2, Methods). Loosening bonds in the cell wall allows the tissue to grow. If new material is inserted that restores the properties of the cell wall then the residual strain returns to zero (‘snip and fill’ [31]). Biologically this would require some form of feedback between resultant stresses (or strains) and cellular properties [34]. Feedback from stresses to microtubule patterns has been proposed [19], and this can be interpreted as reflecting the need to dissipate residual stresses rather than being the primary way of orienting specified growth. Cutting provides a convenient experimental way to evaluate the extent to which residual stresses accumulate or dissipate in a given biological system. Often they accumulate in certain regions in later developmental stages. For example, the dorsal and ventral petals of the adult Snapdragon flower press against each other holding the flower shut (not a part of the model in Green et al [31]). The observation that the accumulation of residuals varies systematically from region to region suggests that the process of dissipating or accumulating residuals is under genetic control. Stresses that are accumulated can be modelled with the GPT-framework and, to enable direct comparison with experimental results, the resulting shapes can be cut allowing the structure to spring into a new shape. The GPT-framework assumes that regions (e. g. cells) in a tissue do not slide or move past each other. This is valid for plants [35], making them particularly appropriate for this approach. The GPT-framework may also be applicable to some aspects of animal development. For example, finite element models have been used to capture deformations during Drosophila ventral furrow formation driven by apical constriction and apicobasal elongation of cells [8]. Comparable deformations can also be generated using GPT-framework by using a posterior-anterior polarity field [36] and incorporating negative growth (contraction) on one side of the canvas (Figure 8). Although this model does not incorporate all biologically relevant features such as constraints of the external vitelline membrane, it illustrates the flexibility of the approach. Clones generated in early wing development of Drosophila often stay as contiguous patches, indicating that connectivity is broadly maintained and extensive mixing of cells does not occur [37]. Greater cell mixing is observed for clones in developing mammalian tissues such as the heart or limb, although even in these cases cell movements are not sufficient to disrupt formation of clonal clusters or patches[38]. At the tissue scale it may therefore be reasonable to model many animal structures with the framework described here, particularly as orientated cell behaviours are thought to play a critical role [39], [40]. As well as multicellular tissues the canvas could represent a region of a plant or bacterial cell wall. By extension of the GPT-framework it may also be possible to capture the growth of compartments enclosed by a canvas (e. g. cells with their walls) or growth of a bulk solid. Thus, the GPT-framework provides a general approach that can be applied to growing tissues at many scales. The GPT-framework with its assumption of tissue polarity as a key component of growth specification provides an economical way of generating diverse shapes and forms. We hypothesise that this combinatorial richness is not only computationally attractive but has also been exploited during evolution to generate a range of observed biological shapes.
  The following theory covers the local specification of growth, how to compute the resulting growth given the mechanical properties of the canvas, how to handle residual growth, and how modelling using the GPT-framework relates to modelling growth in terms of turgor pressure and modifications to the mechanical properties of the cell walls. We distinguish two types of growth, specified and resultant. Resultant growth is the growth that can be directly observed by tracking or clonal analysis. Specified growth is the growth that would happen to an element of the canvas if it grew in isolation. Resultant growth emerges as result of specified growth in different regions interacting through connected tissue. This is illustrated in Figure 9. Panel (A) shows the initial state of a square tissue, divided into a number of small tiles. If we apply a radially increasing field of locally isotropic growth, then in (B) we have an exploded view of how this would affect each tile individually, if it were not attached to its neighbours. It is clear that without some further deformation, these tiles cannot fit together into a continuous tissue without gaps. This conflict between the specified growth field and the continuity of the tissue leads to an equilibrium compromise between the two shown at (C). It is mathematically determined by the partial differential equations of elasticity theory, and numerically computed by the finite element method, both of which we shall briefly summarise. Suppose that at a given time, each point at position in a tissue is moving with an instantaneous velocity. The resultant growth rate in the neighbourhood of is the gradient of the velocity field with respect to. This is the second rank tensor field (a two-dimensional matrix at each point of the tissue) whose components are, where and range from 1 to 3. This velocity gradient tensor represents both the change of shape and size and the rigid rotational motion of the material in the neighbourhood of the point. These are respectively its symmetric and skew-symmetric parts: , where and. is called the resultant strain rate tensor field, and the resultant vorticity. The vorticity field describes the angular velocity at each point. When the vorticity component of a tensor field is zero, the field is called irrotational. To avoid subscripts we abbreviate the definition of to, where is the differential operator defined by. The rate of resultant growth of the material in any particular direction is the sum. Because the resultant strain rate tensor at a point is symmetric, it can be diagonalised by suitably rotating the local frame of reference. The resulting three diagonal components are the principal rates of resultant growth, in three perpendicular directions. These are the eigenvalues of, and the principal growth directions are parallel to its eigenvectors. The growth directions and rates will in general vary over the tissue. To explain how resultant growth may be calculated from specified growth, it is convenient to think in terms of small displacements rather than velocities, by considering the effect over a small time. This is also how the computational implementation (to be discussed below) works, iterating through time in small steps. “Small” here means small enough that first-order approximations apply. In time a velocity field produces a small displacement field, and a growth rate or strain rate tensor field produces an amount of growth or strain, which we shall denote by the same symbols as before. At each point in the growing canvas, let be a specified strain tensor at that point, being the product of a strain rate tensor by a small time. This is the growth that would occur in a small region around in time if it were mechanically isolated from the rest of the tissue. Let be the displacement field that will result from this pattern of growth if the tissue remains in mechanical equilibrium, and the associated growth tensor field. Except in some special cases, such as uniform isotropic growth, will differ from. Even if the rotational component of is ignored, its strain component will still in general differ from: there may be no displacement field of which is the strain field. This is due to the constraint of physical continuity that we mentioned above. (For clarity, the amount of growth shown in Figure 9 has been made far greater than we would normally compute in a single time step.) Physical continuity is expressed mathematically by the St. Venant compatibility constraints [55]. If is a strain field of the form, then it necessarily satisfies the following partial differential equation: This can be verified by substituting for and (somewhat laboriously) finding that all of the terms in the resulting sum of third derivatives of components of cancel out. It is a deeper result that the St. Venant conditions are sufficient for such a velocity field to exist. If, on replacing by in the above equation, it fails to hold, then whatever deformation is applied, the material must remain in a state of frustration. There will be unrelieved residual strain given by. When the material is in mechanical equilibrium, the displacement field will be such as to minimise the energy contained in that residual strain. To calculate, we use the principle of virtual work: if the material is in equilibrium, and any additional infinitesimal displacement is applied, then it will do zero work against the stresses in the material ([56], ch. 2). These stresses are given by a tensor field calculated from the strain and the elasticity properties by the constitutive equation of the material: (1) The subscripts all range over the spatial dimensions 1–3. is the elasticity tensor or stiffness tensor, a 4th rank tensor field representing the elasticity properties of the substance [57]. The work done by any small strain against any stress is, and the total work done for strain and stress fields is found by integrating this over the whole tissue. This is the linear elastic constitutive model, which we are assuming to be valid for small strains. For some biological tissues this assumption may not be accurate, for example as noted in [58] for the mouse ventricle, which also notes that determining a more accurate constitutive model is experimentally challenging. To avoid writing explicit summations, we shall adopt the notations that if and are second rank tensors and and are fourth rank tensors, then: The work done by the strain against the residual stress is then: (2) where the integration is over the whole volume. For to be the equilibrium deformation we must have: (3) Except for degenerate situations (such as the initiation of buckling [59]), this determines up to a rigid translation or rotation of the whole object. We have omitted from equation 3 the possibility of external forces acting on the substance, since there are no such forces present in the applications used in this paper and the Snapdragon model [31]. Boundary conditions can also be applied which stipulate that some parts of the substance remain stationary. We describe how these are handled when we discuss numerical methods. Both the specified growth field and the resultant strain field are by definition irrotational. However, the resultant growth field in general does include rotations. Leaving aside rigid rotations of the whole tissue, the relative rotations between different parts of the tissue are entirely determined by the irrotational tensor. That is, relative rotations are caused solely by differential local growth and the continuity constraints, not by any explicit specification: rotations are always resultant, never specified. Since, the whole analysis carries back to the description in terms of velocities, strain rates, and growth rates. In plants, specified growth rates are always positive, but in animal tissue this is not always so. Both positive and negative growth rates in any direction can be handled computationally without difficulty. Figure 1 shows a simple model in which the shape changes with negligible change of volume. The residual strain is given by the tensor, which is the symmetric part of the residual growth tensor. Most of the examples in this paper discard the residual strain after each time-step of the simulation. In biological terms this is consistent with the observations of [19] that imply a feedback mechanism that acts to absorb stresses. To illustrate the effect of discarding or retaining residual strains we consider several cases in which we cut the canvas after growth or constrain the canvas during growth and then release the constraint. We contrast the effect of discarding residual strains (Cases U and V) with accumulating strains (Cases W and X). These are illustrated in Figure 10. Case U: Dissipation of residual strain with a non-uniform pattern of growth followed by cutting. This is identical to Case C in which residual strain is dissipated at each step. As expected, cutting the canvas induces no further changes of shape as there is no accumulated residual strain. We next show the result of constraining the canvas so that it cannot grow followed by releasing the constraint. Case V: Dissipation of residual strain with a non-uniform pattern of growth that is constrained, followed by release. The model is the same as Case C except that all the boundary points are fixed during growth (column 2). When these constraints are released the shape does not change (column 3) as there is no accumulated residual strain. We next consider the effect of accumulating the residual strain. Case W: Accumulating residual strain with non-uniform pattern of growth followed by cutting. The model is the same as Case C except that residual strain is accumulated at each step. The result is very similar to Case C but there is a small accumulated residual strain (column 2, blue). Cutting and allowing the canvas to relax releases some of this accumulated strain leading to a curve along the line of the cut compared to the straight line in Case U. Case X: Accumulating residual strain with non-uniform pattern of growth that is constrained, followed by release. The model is the same as Case V except that strain is accumulated during growth (blue shows accumulated strain). Releasing the constraints allows a shape to emerge similar to Case W uncut. We may also illustrate the effect of retaining residual strain with a 3D example. For this we use the simplified Snapdragon tube (Case T), but allow residual strain to accumulate on one side. This is illustrated in Figure 11. Case Y: This is the same as Case T, except that residual strain is retained on the right side. Figure 11 B shows how the resultant shape of the right side differs from the left. The residual strain is shown in blue (Figure 11 D). A further difference between the two sides is revealed by making vertical cuts and allowing the mechanical system to relax to a new geometry (Figure 11 C, E). As expected, cutting makes little difference to the left side as there is little residual strain. However, the right side springs apart revealing some of the stored residual strain. In the above examples the release of residual strain by cutting involves large displacements and rotations of the material. However, our computational methods are based on the linear elasticity theory of small displacements, and never directly solve large-displacement problems. The deformation resulting from the release of residual strain is, therefore, computed incrementally, by iteratively applying a small fraction of the residual strain, computing the resulting small deformation, transforming the remainder of the residual strain according to the new orientations of every part of the tissue, and repeating until an equilibrium is reached. The stiffness tensor is a fourth-rank tensor, which in three-dimensional space has components at each point. However, it satisfies certain symmetry properties which imply that it has at most 21 independent components. For isotropic materials, further symmetries imply that is determined by just two values: the bulk modulus and Poisson’s ratio. is the ratio of applied pressure to relative change in volume. We will see later (see TextS1, Equation (S1) et seq.) that its value is irrelevant for the calculations we require: it cancels out of the equations. When a block of material is compressed by external forces in one direction, Poisson’s ratio is the ratio of its transverse expansion to its longitudinal compression. In practice lies between 0 and 0. 5. As the value approaches 0. 5, while the bulk modulus is held constant, the resistance to unidirectional stretching and compression decreases towards zero. If the limit is instead approached by keeping the shear modulus constant, then the bulk modulus tends to infinity. In the former case, the material’s resistance to shears vanishes and it approaches the state of a liquid, while in the latter it approaches an incompressible solid with a finite elastic resistance to everywhere volume-preserving deformations. However, since there are no applied forces (such as gravity) in our models, but only growth described as a change in the resting shape of the material, the difference is more apparent than real. The elasticity tensor computed from (the shear modulus) and (Poisson’s ratio) is equal to the tensor computed from and multiplied by. As mentioned above, any such factor in the elasticity tensor cancels out (Text S1, Equation (S1) ), because all of the forces that we consider result from the material acting against itself. Both methods of computing the elasticity tensor for any value of less than 0. 5 give identical solutions to the equation, solutions which are independent of or. At exactly 0. 5 the equations become highly degenerate, and a different analysis is required to calculate the physical behaviour in the limit. Any value above 0. 5 is physically impossible for isotropic substances, as it would imply that the volume increased under compression, violating conservation of energy. Few experimental determinations of Poisson’s ratio for living plant tissues have been made. They range from 0. 18 to 0. 4 for onion epidermis [60]-[62]. We find that the growth behaviour of a model is insensitive to the precise value of (also see Case 6 in Text S1), and have generally set it equal to 0. 3 in our simulations. In our current models, for simplicity we have taken the elasticity properties to be uniform throughout the tissue and over the time of its development. However, elasticity that varies over the tissue and over time can also be described using the GPT-framework. The analysis so far has assumed that the deformations to be computed are always small. Growth by large amounts can be computed iteratively, by growing in a series of small time steps, in each of which the growth causes only a small deformation. The result is to produce a plastic flow of the material over large time intervals, computed by the theory of small deformations of purely elastic material. Plant growth is thought to occur from a transient reduction in the stiffness of cell walls allowing them to stretch under turgor pressure, new material being added to restore the stiffness [16], [17], [29], [63]–[67]. When the process is anisotropic, it may be because the cell wall fibres typically have coherent directionality, or because the weakening is distributed non-uniformly over the walls of a cell. Most studies have simplified the process by assuming that it is equivalent to increasing the amount of material in a region and then relaxing the shape. This is also the approach used in the GPT-framework: the specified-strain model. The simplification avoids the need to measure relative stiffness and consider turgor pressure. We show below how this approach is related to turgor-based systems. Suppose the tissue has stiffness tensor field, and turgor pressure field. As a result, it must be in some state of strain. The resulting stress in the tissue is. The condition for the tissue to be in equilibrium is that for any small displacement, . Now suppose we change the rest state. For any small piece of the tissue, its rest state is the shape it would be in if the turgor were removed (and the mechanical linkage to the rest of the tissue ignored). is the transformation from that state to the state that it takes up under turgor. Changing the rest state means applying a strain. (The minus sign is due to the fact that we want positive values of to model an increase in the resting size, but the effect of increasing the resting size is to put the tissue into a state of compression, which is described by negative values of strain.) When turgor is reapplied, the resulting strain is. (To validly add strains like this, we are assuming that all of the strains are small.) The applied strain will produce some equilibrium displacement field. The condition for the new equilibrium is. We can subtract from this the original equilibrium condition, leaving, our original equations (2) and (3). This means that the effect of a specified strain field is independent of the turgor, and we can ignore the turgor in our calculations. The following two theorems explore the relationship between the method of specifying the strain and the method of modulating the stiffness tensor. Theorem 1 Let a tissue have a stiffness tensor field, a turgor pressure field and a strain field, such that the tissue is in mechanical equilibrium. Suppose that the stiffness field is then changed to, where is small compared with. Let a second tissue of identical geometry have a stiffness tensor field and be in equilibrium under a strain field. Then there is a specified strain tensor field such that the deformation of the second tissue resulting from applying is the same as the deformation of the first tissue resulting from the change in stiffness. The strain field can be split into two parts: the strain due to turgor, which is (where is the compliance tensor, i. e. the inverse of), and a residual strain. The residual stress field in the first tissue is the residual strain multiplied by, which is. For the tissue to be in equilibrium in this state, the work done by any infinitesimal displacement field against the residual stress must be zero. This work is where the integration is over the whole tissue. Recall that is the differential operator that computes the strain tensor field of a displacement field. If the stiffness is changed to, a new equilibrium configuration will be established by a displacement field. The residual stress field is then, and in equilibrium we have for all. Subtracting the previous virtual work equation gives. is equal to. If we assume that is small in comparison with and is small in comparison with, then the last term can be omitted as being of second order, leaving an effective residual strain of. In the second tissue, the residual stress is initially, and equilibrium implies that is zero for all virtual displacement fields. When the strain is applied, it will produce a displacement field, and a residual stress. As for the first tissue, to determine the equilibrium value of we need only consider the effective residual stress. To prove the theorem it is sufficient to find a value for such that when is taken equal to, the residual strains and are identical at every point. Thus we require to satisfy: Since is invertible, its inverse being a compliance tensor, we can immediately calculate, where is the identity matrix. This proves the Theorem. If the first tissue of this theorem is a biologically accurate description of an increment of growth in terms of the tissue’s background stiffness, turgor, and change in stiffness, then the theorem tells us that we can find another description in terms of a specified strain which gives the same deformation. Furthermore, we have a free choice of the background stiffness. In particular, we can choose to be uniform and isotropic, and constant over time. However, the relationship between and is somewhat complex. When using specified strain to model the result of growth by stiffness modulation, we would like to obtain a closer connection between and, which we now proceed to do. Firstly, if we take, then the expression for simplifies to, and we need no longer calculate. Now suppose that is orthotropic. That is, at each point there are three orthogonal axes such that the change in stiffness is symmetric under a half-turn about each of them. These are called the principal axes of. Under certain extra conditions, we find that the principal axes of coincide with those of. Thus the same distribution of polarisation can be used for either description of growth. Theorem 2 Under the conditions of Theorem 1, suppose that the following conditions hold: 1. , , and are everywhere isotropic. 2. is orthotropic. Then by taking, the principal axes of the specified strain given by Theorem 1 coincide with those of. We have seen already that if, then. and are isotropic by the first condition. By the second condition, is the stress associated with the isotropic strain given the change in orthotropic material properties. Such a stress has the same principal axes as. Multiplying by the isotropic compliance leaves the axes unchanged. We now turn to how the specified growth tensor field is determined by concentration fields of growth factors, and how such concentration fields can be created by defining methods of production, consumption, diffusion, and interaction. A specified growth tensor has three principal axes at right angles to each other. When the tissue is a curved canvas of finite thickness, we assume that although the two sides of the canvas may grow at different rates, they have the same directions of principal growth axes, one of which is always perpendicular to the mid-plane of the canvas, the other two being parallel to it. These axes and the corresponding growth rates are determined by concentration fields of factors. We assume that factor concentrations do not vary through the thickness of the tissue, and therefore represent them computationally by their values on a two-dimensional mesh of triangles, being the midplanes of the mesh of pentahedra used for the elasticity calculation. In the GPT-framework factors can be classified into two types. Identity factors do not propagate within the canvas, while signalling factors can. The specified growth tensor at each point of the canvas is parameterised as follows. The specified principal directions of growth within the plane of the canvas are determined by the gradient of a signalling factor called POL. The specified rates of growth parallel to these directions on the two surfaces of the canvas are given by factors called and. Likewise, the specified growth rates perpendicular to the polarising gradient are given by factors and. The rate of growth of thickness of the canvas is specified by a factor. Propagation of signalling factors may occur through a variety of mechanisms, such as diffusion or active transport. Here we implement diffusion which biologically may be a proxy for a variety of underlying mechanisms. The evolution of a concentration field is modelled by the following equation: (5) The four terms on the right hand side represent respectively diffusion (with a diffusion constant), production at a rate, decay at a uniform rate, and the diluting effect of growth. Here is the volumetric rate of resultant growth. Our implementation handles dilution as a separate step and does not include it in the equation (see later). The values of, , and may vary in space. They may also vary in time, but we assume not rapidly, so that they can be assumed constant over a single timestep. To obtain the concentration field at time, we can make a linear approximation and write a forward Euler equation: (6) By methods similar to those for elasticity, we can discretise this relationship, representing the concentration distribution by its values at the vertices, and obtain an equation similar in form to equation (4): . and are calculated from the geometry of the mesh, the diffusion constant, the production and decay rates, and the current distribution. This set of equations can be solved to give the new concentration distribution. As for the calculation of displacements, the equations allow boundary conditions to be added stipulating that the concentration remains fixed at some nodes. Unlike the case of elasticity, here has full rank and the solution is uniquely determined. The sizes of and are and respectively, where is the number of vertices of the triangular mesh. This is of the value of for the elasticity computation, resulting in a much faster solution. A comparison of a computed diffusion pattern with its analytical solution is considered in Text S1. In the case of morphogens which do not diffuse, it is not necessary to solve the diffusion equation, and the effects of production and decay can be calculated directly, vertex by vertex. We compute diffusion separately from elastic deformation. In principle, the diffusion problem could be solved for a growing and deforming canvas, but over a short time interval only second-order effects arise from the interaction between growth and diffusion, except for the dilution effect mentioned above. When a material expands, the concentration of a physical substance spread through it must decrease in proportion. We make this correction as a separate step: after the diffusion and elasticity calculations, the concentration for each factor subject to dilution by growth is reduced at each point by the proportional expansion at that point, . It is unrealistic to assume that cells can detect the directions of arbitrarily shallow gradients; these also pose numerical problems. There are various options for dealing with very shallow gradients. (1) Generate new sources or sinks for signalling factors as space is created through growth. This would enable patterns to be continually elaborated as the shape expands. (2) Fix the pattern before it becomes too shallow. In the particular case of tissue polarity, for example, polarity may be frozen when the magnitude of the POL gradient falls below a certain threshold. This would be equivalent to a cell becoming polarised when the tissue is small enough for gradients to be measurable, and then retaining its polarity when the gradient falls below the threshold of detectability. Alternatively, (3) the polarity can disappear, resulting in isotropic growth. Case Z: Partitioning a canvas using a diffusing signal. Figure 13 shows a hypothetical example of how diffusion and thresholds can lead to the canvas being partitioned into regions to create a new central region from a peripheral one. An identity factor is expressed at the rim of the disc-shaped canvas (blue in Figure 13 A, C). The expression of an identity factor is represented by the value and non-expression by the value. (could represent a transcription factor expressed only at the rim.) Initially a signalling factor (blue in Figure 13 C) is present everywhere (in this model an initial value of) and its rate of production is promoted by. Diffusion and decay cause the level of to drop in the centre until it reaches a steady state, bowl shape, as shown in Figure 13 C. Wherever drops below a threshold the identity factor is expressed (is set to the value) so defining a new central region. In this example, can be considered as a regional organiser as it provides a source of the signalling factor, , that enables regions to be elaborated. This patterning process could also be used to control the timing of particular events. For example, there can be a further factor, with a high propagation rate, which is generated where falls below a second threshold. When this threshold is reached, will propagate rapidly to activate or inhibit factors bringing in a new patterning phase. Identity, signalling, and growth parameters may interact in many different ways. Rather than assume a fixed set of possible interactions, the software allows the user to write a general function to model interactions, the “interaction function”. The function is called on every iteration of the simulation, before the calculation of diffusion and growth. To simplify the task, a few standard functions are provided to model promotion and inhibition of one factor by another. These are: (7) (8) We use boldface for vectors of values, one value per mesh vertex, and italic for scalar values. Multiplication and division of vectors are to be understood elementwise. and both tend to vectors of 1’s as the components of tend to zero. If factor is to be assigned the value of factor promoted by factor by an amount, one writes the MATLAB equivalent of (i. e. y = z. *pro (k, x);). If is to be inhibited by, then. It is convenient to express inhibition and promotion in this way because the overall effects of different factors (say) that may be expressed in different regions can be obtained by multiplication (e. g.). The iterative loop of the simulation combines the regulatory and mechanical systems as follows.
summary: How do genes control the growth of cells into complex tissue shapes such as flowers, wings or hearts? A key requirement is that genes must be able to modulate growth along particular directions. Two mechanisms have been proposed for how this may work; one based on the directions of mechanical stresses in the tissue and the other on molecular signals that propagate and provide local polarities. Here we show how a polarity-based system has the advantage of being able to act in combination with growth rates to generate a wide range of shapes. By applying this system to the development of the Snapdragon flower, we show, by comparison of computational simulations with actual flower development, how a simple set of polarity controls may underlie the formation of complex biological structures.
keywords: growth control
  plant growth and development
  plant biology
  developmental biology
  plant science
  organism development
  morphogenesis
  pattern formation
  organogenesis
  biology
  systems biology
  biophysic al simulations
  computational biology
  genetics and genomics
title: Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth