Oriented clonal cell dynamics enables accurate growth and shaping of vertebrate cartilage
Abstract
Cartilaginous structures are at the core of embryo growth and shaping before the bone forms. Here we report a novel principle of vertebrate cartilage growth that is based on introducing transversally-oriented clones into pre-existing cartilage. This mechanism of growth uncouples the lateral expansion of curved cartilaginous sheets from the control of cartilage thickness, a process which might be the evolutionary mechanism underlying adaptations of facial shape. In rod-shaped cartilage structures (Meckel, ribs and skeletal elements in developing limbs), the transverse integration of clonal columns determines the well-defined diameter and resulting rod-like morphology. We were able to alter cartilage shape by experimentally manipulating clonal geometries. Using in silico modeling, we discovered that anisotropic proliferation might explain cartilage bending and groove formation at the macro-scale.
https://doi.org/10.7554/eLife.25902.001Introduction
Cartilage is an essential skeletal and supportive tissue in our body. The shape and size of each cartilage element results from complex developmental processes; mesenchymal cells initially condensate, differentiate into chondrocytes, and then an orchestrated growth of the entire structure occurs (Goldring et al., 2006). Often, cartilage plays an important role as a developmental intermediate, such as during the endochondral growth of the long-bones (Mackie et al., 2008). Cartilage elements vary widely in their shapes: they may be simple shapes like rods or bars (Meckel, cartilage templates of the future long bones and ribs) or sheet-like structures (in the head), but can be extremely complicated with a huge number of irregular shapes (for instance, in the inner ear or pelvis). The geometrical properties of cartilage elements must be fine-tuned during the growth because cartilage provides indispensable structural support to the body during development. Yet, how this is achieved despite drastic changes in size is unclear.
After early cartilage forms from mesenchymal condensations, growth typically occurs in all dimensions. However, the diversity of cell dynamics controlling precise early growth and shaping is not well studied. At the same time, the late growth of long rod-shaped cartilage elements in limbs is achieved through a mechanism of endochondral ossification that includes oriented cell dynamics in growth plate-like zones (Vortkamp et al., 1996). In the germinal zone of a growth plate, chondrocytes proliferate and produce progenies that form long streams oriented along the main axis of the forming skeletal element. Inside such streams, chondrocytes undergo flattening, oriented cell divisions and hypertrophy before dissipating and giving place to the forming bone (Nilsson et al., 2005), a process which is controlled by many signals (Kronenberg, 2003). This cell dynamic enables efficient extension of the skeletal element in a specific direction that coincides with the orientation of cell divisions in the proliferative zone (Abad et al., 2002). Growth plate disorders may result in dwarfism and other illnesses (De Luca, 2006).
Some parts of the cartilaginous skull (e.g. the basisphenoid of the chondrocranium) also undergo endochondrial ossification in synchondroses, and significant growth of the cranial base is achieved through a similar mechanism (Hari et al., 2012; Wealthall and Herring, 2006).
Synchondroses are mirror-image growth plates arising in the cranial base, which primarily facilitate growth in the anterio-posterior direction (Kettunen et al., 2006; Laurita et al., 2011; Nagayama et al., 2008; Young et al., 2006). Disorders in the development of synchondroses severely impact the elongation of the cranial base and often result in short-faced mutants and a general decrease of the cranial length (Ford-Hutchinson et al., 2007; Ma and Lozanoff, 1999). Insufficient or abnormal development of a cartilage element is one of the reasons for human craniofacial pathologies, providing a connection between the chondrocranium and facial bone geometry, size and placement (Wang et al., 1999).
The growth mechanism operating in growth plates and synchondroses involves the transformation of the cartilage into the bone. Since growth plates or synchondroses are oriented towards a specific direction, the expansion of a cartilage in other dimensions is not clear from the mechanistic point of view and requires further investigation. For example, although it is well known that the mouse chondrocranium develops as 14 independent pairs of cartilage elements that form one united structure, the logic behind further shaping and scaling remains unclear (Hari et al., 2012). How these initially separated large cartilaginous elements form, grow and fine-tune their geometry, thickness and smoothness during development is still not completely understood. We hypothesized that accurate cartilage growth might require alternative cell dynamics that do not involve hypertrophy, ossification or growth plates.
Such alternative cell dynamics may also contribute to the accuracy of scaling during cartilage growth. Scaling is a process of growth that maintains both the shape and the proportions of the overall structure. In nature, scaling often involves sophisticated principles of directional growth and a number of feedback mechanisms (Green et al., 2010). For instance, during bird development, the diversity in beak shape is constrained by the dynamics of proliferative zones in the anterior face (Fritz et al., 2014). Furthermore, scaling variations of beaks with the same basic shape result from signaling that controls the growth of the pre-nasal cartilage and the pre-maxillary bone (Mallarino et al., 2012). Indeed, in order to accurately scale a pre-shaped 3D-cartilaginous template both local isotropic and anisotropic cell dynamics may be required.
To assess changes in the complete 3D anatomy of the face following cellular-level mechanistic studies we used a variety of approaches including micro-computed tomography (μ-CT), genetic tracing with multicolor reporter mouse strains, multiple mutants and mathematical modelling.
Most importantly, we reveal here how oriented clonal behavior in the chondrogenic lineage controls the overall geometry of the cartilage elements, and show that this geometry can be manipulated with molecular tools at various levels.
Results
Cartilage elements form and grow in all parts of the vertebrate body. The developing face provides a remarkable variety of cartilage geometries and sizes and, therefore, may serve as a sophisticated model system to study the induction of complex cartilaginous structures.
The developing cartilaginous skull, the chondrocranium, displays a very complex geometry of mostly sheet-like cartilages that result from coordinated anisotropic growth in all dimensions. Such expansion of sheet-like cartilaginous tissue during embryonic development involves several mechanisms that were proposed in the past, including the formation and growth of cartilage at synchondroses, as well as at the apical growth zone.
To understand the changes in dimensions of chondrocranium growth at major developmental stages, we took advantage of 3D reconstructions using μ-CT enhanced with soft tissue contrasting (Figure 1). This approach enables the identification of various tissues and cell types in the embryo based on differential uptake of tungsten ions. We validated the μ-CT visualization of embryonic cartilage by directly aligning stained histological sections with the 3D models (Figure 1—figure supplements 1 and 2).
We analyzed the expansion of the chondrocranium due to synchondroses and found that despite a significant anterio-posterior elongation, synchondroses cannot entirely explain the growth dynamics in all directions: anterio-posterior, latero-medial and dorso-ventral vectors of growth (Figure 1A–C). Specifically, we found a complete absence of synchondroses and other endochondrial ossifications in the growing nasal capsule, even at the earliest postnatal stages, while membranous ossifications appeared well developed. The stereotypical clonal cell dynamics found in synchondroses (Figure 1D–I) did not appear during the development of the nasal capsule. Therefore, during the entire embryonic development, chondrocranium growth and shaping is largely aided by additional and unknown mechanisms of growth.
To investigate another possible mechanism of growth, we examined the apical growth zone of the nasal capsule. To understand growth dynamics there, we birth-dated different regions of the facial cartilage using genetic tracing in Col2a1-CreERT2/R26Confetti and Sox10-CreERT2/R26Confetti embryos (Figure 2—figure supplement 1). Both Col2a1-CreERT2 and Sox10-CreERT2 lines recombine in committed chondrocyte progenitors and in mature chondrocytes. 3D analysis following tamoxifen injections at different developmental stages allowed us to identify the parts of the cartilage that develop from pre-existing chondrocytes and the regions generated from other cellular sources (Figure 2, Figure 2—figure supplement 1). As an example, after genetic recombination induced at E12.5, locations with high amount of traced cells show structures that come from pre-existing cartilage, whereas areas comprising from non-traced cells present structures originating from de novo mesenchymal condensations. We discovered that important and relatively large geometrical features are produced from waves of fresh mesenchymal condensations induced directly adjacent to larger pre-laid cartilage elements between E13.5 and E17.5: this includes the frontal nasal cartilage, nasal concha, labyrinth of ethmoid and, consistent with previous suggestions, cribriform plate (Figure 2C and Figure 2—figure supplement 1A–K). These results cannot be safely inferred from 2D traditional histological atlases because of the complex geometry. Our results are complementary to the findings of McBratney-Owen and Morris-Key with coworkers, who demonstrated that the complete chondrocranium (including skull base) develops from 14 pairs of early independently induced large cartilaginous elements that fuse together during later development (McBratney-Owen et al., 2008). Here, we demonstrated how new adjacent mesenchymal condensation can increase the geometrical complexity of a single solid cartilaginous element.
To substantiate our results, we took advantage of Ebf2-CreERT2/R26Tomato transgenic mouse line that can genetically label only a few selected patches of early mesenchyme in the cranial region. We wanted to test if some of these labelled mesenchymal patches can undergo chondrogenesis independently and much later than most of the chondrocranium structure. If that would be the case, we could expect the formation of very sharp borders between the labelled and non-labelled cartilage due to the fusion of newly produced labelled cartilage with the old unlabeled one. If the local cartilage would form from labelled and unlabeled mesenchyme at the same time, the border would not form due to mesenchymal clone mixing that we observe when we label early neural crest. We injected Ebf2-CreERT2/R26Tomato animals with tamoxifen at E12.5 and analyzed the embryos at E17.5 (Figure 2—figure supplement 2). As a result, we discovered that the cartilage element connecting the inner ear with the basisphenoid was genetically traced, and demonstrated a very sharp border with non-traced cartilage (Figure 2—figure supplement 2C–D). μ-CT data confirmed that this element develops entirely after E14.5 from newly formed mesenchymal condensations adjoining the chondrocranium (Figure 2—figure supplement 2A–B), and this might be related to differential regulation at the neural crest-mesodermal border (McBratney-Owen et al., 2008; Thompson et al., 2012). At the same time, the main structure of the chondrocranium is expanded in a very precise and symmetrical way due to unknown cellular and molecular mechanisms that cannot be explained by the freshly induced condensations, the apical growth zone, or even cell dynamics in synchondroses. Our μ-CT results (Figure 2) show that various parts of the chondrocranium develop due to the growth of pre-existing cartilage not involving ossifications, while only additional features are induced in waves as de novo mesenchymal condensations that fuse with the main element during their maturation or expand in the process of ossification.
We further focused on the developing nasal capsule because its growth does not involve synchondroses while the apical growth zone and adjoining mesenchymal condensation only partly provide for the growth and shaping modifications.
The results obtained from comparisons of cartilaginous nasal capsules from different developmental stages showed that the shape of the structure is generally established by E14.5 (Figure 3, Figure 3—figure supplement 1, Figure 3—figure supplement 2, Video 1). Nevertheless, from E14.5 until E17.5 the cartilaginous nasal capsule is accurately scaled up with significant geometrical tuning (Figure 3A–B). Previous knowledge suggests that the underlying growth mechanism should be based on appositional growth of the cartilage during its transition to bone (Hayes et al., 2001; Li et al., 2017), however, numerous facial cartilages never ossify, but continue to grow.
Tomographic reconstructions of sheet-shaped cartilage elements in the nasal capsule revealed extensive expansion of the cartilage surface area and overall volume (Figure 3E–F). Surprisingly, the thickness of the cartilaginous sheets did not change as much as the other dimensions during nasal capsule growth (Figure 3C–F, Figure 3—figure supplement 1, Video 1). Thus, the sheet-shaped cartilage expands mostly laterally (within the plane) during directional growth. Therefore, we expected that clonal analysis of the neural crest progeny (with Plp1-CreERT2/R26Confetti) and of early chondrocytes (with Col2a1-CreERT2/R26Confetti or Sox10-CreERT2/R26Confetti) would reveal clonal units (so called clonal envelopes) oriented longitudinally along the axis of the lateral expansion of the cartilage. Surprisingly, and contrary to this, clonal color-coding and genetic tracing demonstrated transversely oriented clones represented by mostly perpendicular cell columns or clusters formed by traced chondrocytes (Figure 4, especially A-C, Figure 4—figure supplement 1).
To understand this process more in depth, we started with genetic tracing of the neural crest cells and their progeny in the facial cartilage with Plp1-CreERT2/R26Confetti (tamoxifen injected at E8.5). Clonal analysis and color-coding of neural crest-derived chondrogenic and non-chondrogenic ectomesenchyme showed intense mixing of neural crest-derived clones in any given location (Figure 4A–C, Figure 4—figure supplement 1D–G) (Kaucka et al., 2016). At the same time, chondrogenic ectomesenchyme demonstrated the presence of transversely oriented doublets of genetically traced and also EdU-labeled cells already at E13.5 (Figure 4—figure supplement 1G (inserts) and H). Next, analysis of neural crest progeny in established cartilage highlighted the presence of perpendicularly oriented clonal doublets and columns (Figure 4A–C). Further analysis of EdU incorporation and genetic tracing with chondrocyte-specific Col2a1-CreERT2/R26Confetti and Sox10-CreERT2/R26Confetti lines confirmed the existence of transversely oriented products of cell proliferation in the mature (E14.5-E17.5) cartilage (Figure 4D–F for EdU, Figure 4G,H–L and Figure 4—figure supplement 2 for lineage tracing). These results imply that cells in the sheet-shaped cartilage do not allocate daughter cells in lateral (longitudinal) dimensions as would be intuitively expected.
Thus, simple lateral or unidirectional proliferation cannot account for the accurate scaling of the sheet-shaped cartilage in the face. Instead, the cartilage development from chondrogenic condensations is achieved by a cellular mechanism that involves intercalation of columnar clonal units.
It was unclear to us why column-like structures, and no other shapes, are integrated into the sheet-shaped cartilage and how the fine surface is maintained during this mechanism of growth. To better understand possible mechanisms of accurate sheet-shaped cartilage surface development we modelled individual cell dynamics, in silico in 4D (3D + time) (Figure 5) (Hellander, 2015). We used this modelling to address two questions: firstly, under what conditions are clonal columns observed? Secondly, how is the sheet-like shape achieved by polarized or non-polarized cell divisions of single-cell thick layers and what are the controlling mechanisms? We tested a group of variables including: cell division speed, allocation of daughter cells in random- or defined directions, orientation cues in the tissue (equivalent to molecule gradients), as well as pushing/intercalating of the daughter cells during proliferation. We qualitatively compared the results from in silico simulations to our experimental clonal analysis from various genetic tracing experiments, in order to identify conditions in the model that were compatible with patterns observed in vivo.
The results of the mathematical modelling suggested that the clonal dynamics observed in natural conditions requires polarity cues in the system, specifically, a two-sided gradient of signals would be required to precisely fine- tune cartilage thickness (Figure 5A–J). At the same time, some yet to be identified mechanism controls the average number of cell divisions in a column, further controlling columnar height and undoubtedly regulating the local thickness of the cartilage. Combined with the observed introduction of the transverse clonal columns, oriented cell proliferation can provide fine surface generation and scaling (Figure 4—figure supplement 2). Moreover, the model highlighted the elegance of cartilage design involving transverse columnar clones in the sheet-shaped elements: this logic enables the uncoupling of thickness control (depends on cell numbers within a clone) and lateral expansion (depends on the number of initiated clones), which are likely two molecularly unrelated processes in vivo. The absence of a gradient during in silico simulations led to the generation of 3D asymmetrical clusters instead of straight columns (even in conditions of highly synchronized cell divisions, and starting from a laterally space-constrained initial configuration - suggesting the promotion of vertical growth due to space-exclusion in the lateral direction) (Figure 5I,J). This, in turn, led to the formation of surface irregularities in the cartilage with subsequent loss of local flatness (heat-map diagram in Figure 5I,J).
Importantly, lineage tracing also showed that for cartilaginous structures in the head with asymmetrical or complex irregular geometries, such as areas where several sheet-shaped cartilage elements were merged, clones were not constructed to perpendicular columns. In such locations, we identified irregular clonal clusters or randomly oriented clonal doublets, in accordance with the modelling results (Figure 6A–J). Thus, the shape and orientation of clones corresponds to the local geometry of the cartilage element.
Next, we attempted to target a molecular mechanism that controls the flatness and sheet-like shape of the facial cartilaginous sheets. We discovered that activation of ACVR1 (BMP type one receptor, ALK2) in developing cartilage leads to a phenotype with targeted clonal micro-geometries (Figure 6K–P,R). We utilized a constitutively activated caALK2 transgene (Fukuda et al., 2006) together with genetic tracing in a way that every GFP-expressing cell is carrying constitutively active ACVR1. This experiment revealed a dramatic change of the shape of clonal envelopes, changing from straight perpendicular columns to disorganized spherical clusters inside the sheet-shaped cartilages of transgenic Sox10-CreERT2/R26caALK2-IRES-GFP embryos (Figure 6K–N). The ectopically activated ACVR1 resulted in the presence of clonal spherical clusters that interfered with the cartilage borders and caused the formation of ectopic bumps, swellings and other abnormal local shapes - in accordance with the mathematical modelling predictions (substantially resembling the condition with no gradient, see Figure 5I) (Figure 6O–P). All recombined cells in this caALK2 experiment became Sox9+ chondrocytes. There were no other cell types found to be GFP+, including perichondrial cells. This result indicates that BMP family ligands either produce the gradient that directs the orientated behavior of chondrocytes inside of the cartilage or, alternatively, that an experimental increase of BMP signaling renders the cells insensitive to the gradient established by other molecules. In any case, ACVR1 mutation can be used as a tool to change columnar arrangements into clusters (Figure 6N,R). The activation of ACVR1 by Sox10-CreERT2 starting from E12.5 occurred both in perichondrial cells and in chondrocytes (based on our genetic tracing results using Sox10-CreERT2/R26Confetti). This later coincided with clonal bumps and bulging regions positioned mainly at the surface of sheet-shaped cartilaginous sheets (Figure 6L,N,P). These data also support the hypothesis that integration of clonal chondrocyte clusters into existing cartilaginous sheets likely depends on clonal shape and originates from the periphery of the cartilage. When this column-inserting process fails, the progeny of cells at the periphery of the cartilage forms ectopic bumps outside of the normal cartilage borders, and disrupts the flatness and straightness of cartilage surfaces.
Next, we attempted to block the planar cell polarity (PCP) pathway to challenge the system and disrupt the formation of perpendicular columns in the flat or curved cartilaginous sheets. To do this we performed µ-CT and EdU-incorporation analysis on Wnt/PCP mutants. Wnt/PCP pathway is well known for driving the cell and tissue polarity, and distinct facial phenotypes have appeared in Ror2, Vangl2 and Wnt5a homozygous mutants (Figure 7A). When EdU was administered 24 hr before embryo harvest, subsequent analysis showed no differences in the EdU-positive perpendicular clonal columns which formed within sheet-shaped facial cartilage of Wnt5a knockout mutants or wild type controls (Figure 7B–D). µ-CT analysis at early developmental stages showed that as early as at E12.5, Wnt5a mutants had abnormal shape and placement of the mesenchymal condensations that create a template for future cartilaginous structures (Figure 7E–F). Although µ-CT analysis of Wnt5a, Ror2 and Vangl2 homozygous mutants at later developmental stages confirmed that chondrocranium shape was heavily affected (with generally shortened nasal capsules as compared to both wild-type and heterozygous controls) (Figure 7A), we did not detect any defects in cartilage micro-geometry, including thickness or surface organization. Altogether, these results indicate that Wnt5a, Ror2 and Vangl2 do not control cartilage growth and shaping per se (via the insertion of perpendicular columns). Instead, they influence the position and shape of chondrogenic condensations, which define the future geometry of the facial chondrocranium (Figure 7E,F).
Following the prediction from our mathematical model that the thickness of the cartilage can be controlled by the number of cells in the inserted clonal column, we searched for the molecular mechanisms which control this. Knowing from our results that proliferation rate drops in the mature cartilage, we hypothesized that chondrocyte maturation speed may influence the number of cell divisions within a column. To test this suggestion, we analyzed G-protein stimulatory α-subunit (Gsα) knockout embryos (Figure 8). Inactivation of Gsα, encoded by Gnas, is known to lead to accelerated differentiation of columnar chondrocytes, without affecting other aspects of cartilage biology (Chagin et al., 2014). We analyzed three different locations in the developing chondrocrania, and observed a significant reduction of cartilage thickness in absolute metrics (Figure 8A,B,J), as well as in terms of the number of cells within each column (Figure 8B–H). Thus, the Gsα knockout is a perfect tool to test whether the modulation of differentiation speed can be used to create a variation of local cartilage thickness. The result of this experiment demonstrated that sheet-shaped cartilages in Gsα knockout embryos are thinner than that of littermate controls, while other parameters (including general size and shape of nasal capsule and other locations in the head together with the transverse orientation of chondrocyte columns) remain largely unaffected (Figure 8A,B,I). Thus, these data experimentally validated mathematical predictions and confirmed that the thickness of cartilage is determined by the number of cell divisions within a transverse clone, and that this is uncoupled from lateral expansion.
Next, we wanted to know how clonal cell dynamics accounts for the shape development in rod-shaped cartilages. For this we investigated the clonal dynamics in Meckel, rib and limb cartilages with the help of Confetti-based genetic tracing as well as EdU incorporation. The clonal arrangements appeared highly oriented and strongly resembled the clonal columns we observed in the facial cartilage. The columnar clones were oriented mostly transversally in the plane of a rod diameter and could not explain the early growth along the main axis of the skeletal element (Figure 9). These tracing results suggested that longitudinal extension is based on continuous development of chondrogenic mesenchymal condensations on the distal tip and is followed by the transverse proliferation of chondrocytes, which accounts for the proper diameter of a cartilaginous rod. The logic of oriented cell dynamics in sheet-shaped and rod-shaped cartilages is summarized in Figure 10.
Since the integration of clonal units is likely to be uneven in the cartilage, we questioned how the anisotropy of local proliferation can impact the shaping processes on a macro-scale. Starting from E14.5, the olfactory capsule is already formed of mature chondrocytes. Indeed, in this structure, proliferation was localized to specific regions, but remained generally low elsewhere (Figure 11A–B) according to the analysis of EdU incorporation. As we demonstrated above, proliferative regions expand due to the active integration of new clonal columns and clusters. We projected the low- and high proliferative zones onto the 3D structure of the nasal capsule at E13.5-E15.5 to understand not only the dynamics of lateral expansion, but also to see how the local expansion of cartilage may influence bending and geometrical changes on a large scale (Figure 11C–F). Since proliferative zones in nasal capsule are restricted and have defined edges, they inevitably induce tension and bending of the surrounding cartilage sheet.
In order to address the logic of distributed proliferative zones and its role in shape transitions between stages we took advantage of the mathematical model developed by the Enrico Coen and Andrew Bangham laboratories. This model has been efficiently validated and applied for advanced simulations of complex 4D plant organ development (Green et al., 2010; Kennaway et al., 2011). To simulate in silico nasal capsule shape transition from E13.5 to E14.5, we generated a basic E13.5-like shape by converting a sheet-shaped growing trapezoid into a corresponding 3D structure (Figure 12A, central part and Video 2). The result was considered as a simplified starting condition for further simulations. Next, two lateral zones with a low rate of proliferation were introduced according to their original position in E13.5 nasal capsule. Further simulations of the growth showed that these low proliferative zones impose a characteristic bending on the sides of the simulated structure. This bending corresponds to the lateral transformations observed in embryonic development of the nasal capsule between E13.5 and E14.5 (Figure 12B–C). This characteristic lateral bending did not depend on anterio-posterior polarity in the cartilage or formation of the groove at the midline (Figure 12D). According to the model, the polarity only affected the potential for the anterior elongation due to the anisotropic growth of the entire cartilaginous structure. Our results also suggested that the nasal septum functions as a slower proliferating anchoring point to the roof of the nasal capsule, which is necessary for the formation of the midline groove at E14.5. A simulated groove at the midline provided for the general bend and flattened shape of the in silico cartilage, similar to the native E14.5 nasal capsule and contrary to the model without the simulated midline groove (Figure 12C–D).
To validate the general rules of in silico transformations, we performed material modelling using plastic film to simulate anisotropic expansion and bending due to integration of local growing zones with attached borders. This simple material modeling demonstrated that growth zones/local expansions in the flat planes generate mechanical tensions which bend the structure (Figure 12E). We then performed another material modelling experiment using isotropic thermal expansion/constriction of a plastic film. For this purpose, we drew black regions (analogous to the lateral low proliferative zones in E13.5 nasal capsule) onto white plastic film that was cut in a shape of a trapezoid capable of transforming into a nasal capsule-like dome. Under the heating provided by a thermal infrared lamp, the black zones received more heat and isotropically shrunk. Shrinkage of the black zones created physical tensions that eventually bent the structure in a way similar to the original nasal capsule geometry at E14.5 (Figure 12F). The model with shrinking zones is comparable to the real growth conditions as the nasal capsule expands faster than spatially distributed slow proliferative regions (simulated as shrinking zones inside of the non-expanding plastic material). These results, combined with analysis of proliferation and 3D visualizations, strongly suggest that the distribution of uneven proliferative zones plays an important role in the shaping of the facial chondrocranium during embryonic development.
Taken together, we reveal a set of principles contributing to the accurate scaling and shaping of cartilage tissue during growth. The reverse engineering of this process highlights the involvement of highly specialized systems that control the directional growth at the levels of micro- (clonal shapes) and macro-geometries (proliferative regions in nasal capsule). Our results show that allometric growth of complex 3D cartilage elements is not achieved by simple, evenly distributed and/or unidirectional proliferation, but is sculpted by precisely localized proliferation.
Discussion
Here we report the discovery of how oriented cell behavior and molecular signals control cartilage growth and shaping. Previously, the use of chimeric avian embryos demonstrated the competence of facial mesenchyme in producing species-specific shapes and sizes of cartilage elements (Eames and Schneider, 2008), while facial epithelium and brain provided the instructive signals guiding generalized shaping of the face (Chong et al., 2012; Foppiano et al., 2007; Hu et al., 2015). Knowledge of how the facial cartilaginous elements are shaped has been rather restricted, and mainly concerned with the correct formation of chondrogenic mesenchymal condensations.
The accurate expansion of the chondrogenic condensation or cartilage during growth is no trivial matter. The general shape should be both preserved and modified at the same time. We show that anisotropic proliferation and oriented clonal cell dynamics are implemented to achieve the necessary outcome. The reverse engineering of this process highlighted the involvement of highly specialized systems that control the directional growth at the levels of micro- (clonal shapes) and macrogeometries (proliferative regions).
Allometric growth of complex 3D structures requires certain cellular logics and cannot efficiently proceed with equally distributed and/or unidirectional proliferation inside of the mesenchymal condensation or cartilage element. On the other hand, we did not observe the formation of growth plate-like zones in early sheet-shaped (nasal capsule) or rod-shaped cartilages, nor uniform expansion of cartilage in all directions. Thus, the underlying growth and shaping mechanisms required an explanation.
To test various strategies of cellular behavior during cartilage growth we devised a model simulating different aspects of multicellular dynamics in 3D together with lineage tracing of individual clones. Most of the currently existing models of cell dynamics and tracing operate in 2D space, which often limits the predictions (Jarjour et al., 2014). Our model suggested that a gradient-controlled orientation of clonal expansion can explain the biological observations (i.e. it is consistent with ordered columnar growth and its disruption results in spherical microdomains rather than columns), and showed the relation between the geometries of clonal domains (envelopes), the overall shape and the fineness of the surface. We confirmed the predictions from the model in a series of experiments involving tracing with multicolor reporters and manipulating the cartilage with mutations. Our results showed that the formation of oriented clones of chondrocytes with clonal envelope shape corresponds to the geometry of the analyzed locality. The sheet-shaped cartilage elements consisted of transversely oriented clonal columns, while asymmetric complex geometries revealed a variety of clonal shapes ranging from spherical to particularly oriented.
Genetic tracing initiated during transition of condensations into cartilage resulted in clonal columns within both sheet- and rod-shaped cartilage elements. This confirms that chondrogenic condensations undergo complex oriented cell dynamics during their development. Importantly, tracing of chondrocranium cartilage showed formation of transverse clonal columns as growth proceeded. Intercalation of newly born columns into pre-existing cartilage provided for the expansion potential in the sheet-shaped cartilage. This growth mechanism is very original and is not reported elsewhere so far.
A few studies have demonstrated how clonal envelopes form in accordance with the general shape of the structure. These were mainly conducted on Drosophila imaginal wing disc or growing flower petals. In all cases the authors highlighted that the shape of clonal geometries correlates with the major vector of expansion in the growing structure (Green et al., 2010; Repiso et al., 2013; Strutt, 2005). This implies the presence of polarized activity that directs the shaping of the tissue. Here, we provided the first experimental evidence of how the control of the directional clonal expansion influences the shape of a vertebrate tissue on a large scale. Moreover, in the sheet-shaped cartilage the orientation of clonal domains, i.e. the columns, does not correspond to the vectors of major expansion, but rather serves for uncoupling lateral expansion control and thickness tuning. In line with that, the number of chondrocytes comprising the clonal column or cluster depends on Gsα-mediated signals. Variations in this number do not significantly affect the lateral dimensions of the whole sheet-shaped cartilage structure: the thickness of the cartilage becomes less while the general geometry and size stay preserved. Additionally, the shape and orientation of clonal envelopes in cartilage is partially controlled by BMP signaling, since micro-geometries of clones depend on activation of ACVR1. Based on these results, we assume that BMP ligands (because of cAlk2/ACVR1 phenotype affected clonal orientation) expressed around the regularly shaped cartilages may play a role similar to the in silico predicted gradients. Indeed, the expression of INHBA, BMP5 and BMP3 fit this expression profile quite well (according to Allen Developing Mouse Brain Atlas (http://developingmouse.brain-map.org) and Eurexpress (http://www.eurexpress.org) in situ public databases). At least, BMP5 is clearly expressed at the cartilage periphery and has been shown to affect the cartilage shape by David Kingsley lab (Guenther et al., 2008).
Our experimental manipulations of planar cell polarity (PCP) pathway did not affect microgeometries and clonal domains, but strongly affected the chondrocranium shape on the macroscopical scale in several different ways. These phenotypes appeared to be rooted in pre-chondrogenic or early chondrogenic stages, and are based on distorted placement of mesenchymal condensations in the very early head. These experiments with Wnt/PCP mutants may potentially provide a better understanding of species-specific mechanisms of control and evolution of the facial shape on a macro scale.
Regular shapes require regular cellular arrangements and clonal cell dynamics. It is not only sheet-shaped cartilage in the head that demonstrate geometric regularity; rod-shaped cartilage (Meckel, embryonic ribs and long cartilages in limbs) also has a regular shape. Regular clonal patterns, conceptually similar to those found in sheet-shaped cartilage, explain conservative tissue dynamics during formation and growth of cartilaginous rods. Indeed, genetic tracing experiments suggested that formation of clonal columns is important for the diameter control, while chondrogenic condensations at the very tip of the rod-shaped growing structures enable elongation. Similar to the cell dynamics in the sheet-shaped cartilage, this mechanism may provide for uncoupling of length versus diameter control. Such uncoupling may generally enable developmental and evolutional plasticity of cartilage size and shape.
The mechanism controlling the thickness or diameter of sheet-shaped and rod-shaped cartilage elements not only includes spatially orientated behavior, but also involves the regulation of cell number within each chondrogenic clone. Immature chondrocytes are proliferatively active, while more mature chondrocytes show decreased proliferation. Therefore, differentiation speed emerges as a concept which could regulate the organ shape by impinging on clone size, thereby altering the thickness or diameter of the cartilage. This concept is known to operate in the brain and other tissues with classical stem cell/transiently amplifying cell arrangements (Díaz-Flores et al., 2006).
Clonal genetic tracing and EdU labeling experiments suggested that the origin of clonal columns and clusters might be represented by the cells located at the periphery of forming cartilage. The spherical clusters of chondrocytes forming at the periphery of the cartilage in cALK2 mutant mice may suggest that the cell source is also located at the periphery and might be a perichondrial cell. Clonal relationships between perichondrial cells and columns of chondrocytes also support the hypothesis of perichondrial cells acting as a stem population during cartilage expansion. In general, the heterogeneity and multipotency of perichondrial cells is still unclear, although there are multiple studies showing the perichondrium as a source of chondrocytes and osteoblasts (Kobayashi et al., 2011; Li et al., 2017; Maes et al., 2010).
In addition to this, the perichondrium might mediate non-autonomous effects in the cartilage in case of cAlk2 and GSα experiments. Genetic tracing shows that some perichondrial cells always recombine with Sox10-, Plp1- and Col2a1-CreERT2 lines, and, in case of functional experiments, may indirectly control some evens in more mature layers. Also, it is not clear how the fine border of the cartilage is set, and whether the perichondrial layer may play a key border-setting role during development and regeneration. This should be investigated further.
Next, our results show that tuning of macro-geometries on a large scale can be achieved through a stage-specific placement of proliferative hot zones where new clonal domains intercalate into the main cartilage structure. Anisotropic heterogeneous proliferation is a powerful tool, which, together with polarity in the tissue and local patterning, can drive the organ shape development (Ben Amar and Jia, 2013; Campinho and Heisenberg, 2013). The localized growth zones provide for the general expansion and also bend the cartilage by creating local tensions that require mechanical relaxation and influence further development of the overall shape (Schötz et al., 2013). For probing such transformations of the sheet-shaped facial cartilage we applied an in silico model that was already successfully validated in a number of growth, shaping and scaling tasks (Green et al., 2010; Kennaway et al., 2011). Such a model was necessary to understand why the high and low proliferation zones are positioned in such a specific way. Indeed, the discovered distribution of proliferative zones in the whole nasal capsule did not help us per se with intuitive explanations of geometrical changes on the macro-scale. Despite this counter-intuitive dataset, the mathematical model provided an insight into the logic of the high and low proliferation zones in relation to a transition between investigated cartilage shapes.
For example, it turned out that the position of lateral slow proliferation zones enables the generation of the symmetrical bends at the sides of the nasal capsule during transition from E13.5 to E14.5 developmental shapes. Furthermore, real material modelling confirmed the results predicted by the mathematical model, and generated lateral bends similarly to the native structure. The molecular mechanism controlling the dynamic distribution (patterning) of fast/slow proliferative zones in the cartilage is still unknown. It is likely linked to developmental signals from other tissues such as the olfactory epithelium or the mesenchyme surrounding the cartilage. Identification and validation of these signals will be essential in future studies and would involve a substantial combination of screening and functional approaches with transgenic animal models.
The anisotropic proliferation can be an important evolutionary mechanism that is directly responsible for the differences in snout geometry in a variety of phylogenetic groups. Additionally, it might be important for understanding the development of the facial shape variation in humans (Sheehan and Nachman, 2014) as well as numerous pathologies (Afsharpaiman et al., 2013).
One alternative way to fine-tune macro-geometry of a cartilage element is to continuously add on pre-shaped chondrogenic mesenchymal condensations from the pool of competent progenitors that are retained until late developmental stages. As we demonstrated, the formation of adjoining mesenchymal condensations occurs in sheet-shaped cranial cartilage. In the developing face, new chondrogenic condensations are responsible for introducing geometrically complicated fine details. Such mechanisms may also operate during amphibian metamorphosis, when most of the postmetamorphic cranial cartilage develops de novo and not from the pre-metamorphic cartilaginous elements (Kerney et al., 2012).
Taken together, we discovered important novel principles explaining the growth and shaping of cartilaginous structures. Further studies should focus, amongst other things, on the soluble signals emanating from other embryonic structures which influence the oriented behavior or proliferation of chondrogenic clones.
Materials and methods
Mouse strains and animal information
Request a detailed protocolAll animal (mouse) work has been approved and permitted by the Ethical Committee on Animal Experiments (Norra Djurförsöksetiska Nämd, ethical permit N226/15 and N5/14) and conducted according to The Swedish Animal Agency´s Provisions and Guidelines for Animal Experimentation recommendations. Genetic tracing mouse strains Plp1-CreERT2 (RRID:MGI:4837112) and Sox10-CreERT2 were previously described (Laranjeira et al., 2011; Leone et al., 2003; Yu et al., 2013). Plp1-creERT2, Sox10-creERT2 and Col2a1-CreERT2 (RRID:IMSR_JAX:006774) (Nakamura et al., 2006) (obtained from laboratory of S. Mackem, NIH) strains were coupled to R26Confetti (RRID:IMSR_JAX:017492) mice that were received from the laboratory of Professor H. Clevers (Snippert et al., 2010). The Stopflowed/floxedcaAlk2-IRES-GFP strain from the laboratory of Y. Mishina (Fukuda et al., 2006) was coupled to Sox10-CreERT2. The Ebf2-CreERT2 (RRID:MGI:4421811) strain was obtained from the laboratory of H. Qian, KI, and was coupled to R26Tomato. The Gsαfloxed/floxed strain was obtained from the laboratory of L. Weinstein (Sakamoto et al., 2005). Female mice which were homozygous for the reporter allele [Gt(ROSA)26Sortm4(ACTB-tdTomato,-EGFP)Luo/J; Jackson Laboratories] (Muzumdar et al., 2007) were coupled to homozygous Col2a1::creERT males [FVB-Tg(Col2a1-cre/ERT)KA3Smac/J; Jackson Laboratories] (Feil et al., 1997; Nakamura et al., 2006). To induce genetic recombination to adequate efficiency, pregnant females were injected intraperitoneally with tamoxifen (Sigma Aldrich, St.Louis, MO, T5648) dissolved in corn oil (Sigma Aldrich, C8267). Tamoxifen concentration ranged from 1.5 to 5.0 mg per animal in order to obtain a range of recombination efficiencies. Wnt5a, Vangl2 and Ror2 full knock-out embryos were obtained from heterozygous parents (Gao et al., 2011; Yamaguchi et al., 1999) at the expected Mendelian proportions.
Immunohistochemistry
Request a detailed protocolFor embryo analyses, heterozygous mice of the relevant genotype were mated overnight, and noon of the day of plug detection was considered E0.5. Mice were sacrificed with isoflurane (Baxter, Deerfield, IL, KDG9623) overdose, and embryos were dissected out and collected into ice-cold PBS. Subsequently, the samples were placed into freshly prepared 4% paraformaldehyde (PFA) and depending on the developmental stage they were fixed for 3–6 hr at +4°C on a roller. Embryos were subsequently cryopreserved in 30% sucrose (VWR, Radnor, PA, C27480) overnight at +4°C, embedded in OCT media (HistoLab, Serbia, 45830) and sections cut of between 14 µm to 200 µm on a cryostat (Microm International, Germany), depending on the following application. If needed, sections were stored at −20°C after drying for 1 hr at room temperature, or processed immediately after sectioning. Primary antibodies used were: goat anti-GFP (FITC) (Abcam, UK, 1:500, RRID:AB_305635), rabbit anti-Sox9 (Sigma Aldrich, 1:1000, RRID:AB_1080067), rabbit anti-Sox5 (Abcam, 1:500, RRID:AB_10859923), sheep anti-ErbB3 (RnD Systems, Minneapolis, MN, 1:500, RRID:AB_2099728). For detection of the above-mentioned primary antibodies we utilized 405, 488, 555 or 647-conjugated Alexa-fluor secondary antibodies produced in donkey (Invitrogen, Carlsbad, CA, 1:1000, RRID:AB_162543, RRID:AB_141788, RRID:AB_141708, RRID:AB_142672, RRID:AB_2536183, RRID:AB_141844,). Sections were mounted with 87% glycerol mounting media (Merck, Germany) or in Vectashield Antifade Mounting Medium with DAPI (Vector Laboratories, Burlingame, CA, RRID:AB_2336790).
EdU incorporation analysis
Request a detailed protocolEdU (Life Technologies, Carlsbad, CA) was injected intraperitoneally into the pregnant females (65 µg per gram of body mass) either 6- or 24 hr before the embryos were harvested. Cells with incorporated EdU were visualized using Click-iT EdU Alexa Fluor 647 Imaging Kit (Life Technologies) according to the manufacturer’s instructions.
Microscopy, volume rendering, image analysis and quantifications
Request a detailed protocolConfocal microscopy was performed using Zeiss LSM710 CLSM, Zeiss LSM780 CLSM and Zeiss LSM880Airyscan CLSM instruments. The settings for the imaging of Confetti fluorescent proteins were previously described (Snippert et al., 2010). The imaging of the confocal stack was done with a Zeiss LSM780 CLSM, Plan-Apochromat 3 10x/0.45 M27 Zeiss air objective.
Histological staining
Request a detailed protocolSlides were stained for mineral deposition using von Kossa calcium staining: 5% silver nitrate solution was added to the sections at a room temperature and exposed to strong light for 30 min. After that the silver nitrate solution was removed, and slides were washed with distilled water for three times during 2 min. 2.5% sodium thiosulphate solution (w/v) was added to the sections and incubated for 5 min. Slides were again rinsed for three times during 2 min in distilled water. The sections were then counterstained using Alcian blue. Alcian blue solution (0.1% alcian blue 8GX (w/v) in 0.1 M HCl) was added to the tissue for 3 min at room temperature and then rinsed for three times during 2 min in distilled water. Slides were then transferred rapidly into incrementally increasing ethanol concentrations (20%, 40%, 80%, 100%) and incubated in 100% ethanol for 2 min. Finally, the slides were incubated in two xylene baths (for 2 min and then for 5 min) before mounting and analysis.
Statistics
Statistical data are represented as mean ± s.e.m. Unpaired version of Student’s t-test was used to calculate the statistics (P value). All results were replicated at least in three different animals. Statistical analysis and graphs were produced in GraphPad Prism (La Jolla, CA, RRID:SCR_002798) or Oriana Software (Kovach Computing Services, UK). Spearman coefficient was used for correlation assessment of microgeometries corresponding to different locations in the cartilage.
In Figure 8 the difference between control (mean = 5.9, sem = ±0.23, n = 4) and mutant (mean = 4.3, sem = ±0.25, n = 3) olfactory cartilage thickness is significant (p=0.0053). The difference between control (mean = 10.6, sem = ±0.83, n = 3) and mutant (mean = 5.7, sem = ±0.61, n = 3) basisphenoid cartilage thickness is significant (p=0.0087).
Tissue contrasting for µ-CT scanning
Request a detailed protocolOur staining protocol has been modified from the original protocol developed by Brian Metscher laboratory (University of Vienna, Austria). After dissection, the embryos were fixed with 4% aqueous solution of formaldehyde in PBS for 24 hr at +4°C, with slow rotation. Samples were then dehydrated by incubation in incrementally increasing concentrations of ethanol in PBS (30%, 50%, 70%); samples were incubated at +4°C for two days in each concentration to minimize the tissue shrinkage.
We found that the best signal to noise ratio on scans results from contrasting the samples with 0.5–1.0% PTA (Phosphotungstic acid, Sigma Aldrich) in 90% methanol. After sample dehydration, the tissue-contrasting PTA solution was added to the samples and then changed every day with the fresh solution. E12.5 embryos were contrasted with 0.5% PTA for four days while E15.5 embryos were stained in 0.7% PTA for six days. E16.5 and E17.5 embryos were decapitated, and the contrasting procedure was extended to 9–15 days in 1% PTA to ensure the best penetration of the contrasting agent. Subsequently, tissues were rehydrated through a methanol gradient (90%, 80%, 70%, 50% and 30%), to sterile distilled water. After that, rehydrated embryos were embedded in 0.5% agarose gel (A5304, Sigma-Aldrich) and placed in polypropylene conical tubes (0.5, 1.5 or 15 ml depending on the sample size) to minimize the amount of surrounding agarose gel, and to avoid movement artifacts during X-ray computed tomography scanning.
µ-CT analysis (micro computed tomography analysis)
Request a detailed protocolThe µ-CT analysis of the embryos was performed using laboratory system GE phoenix v|tome|x L 240, equipped with a 180 kV/15W maximum power nanofocus X-ray tube and high contrast flat panel detector DXR250 with 2048 × 2048 pixel, 200 × 200 µm pixel size. The exposure time was 900 ms in all 2000 positions. The µ-CT scan was carried out at 60 kV acceleration voltage and with 200 μA X-ray tube current. The voxel size of obtained volumes appeared in the range of 4 µm - 6 µm depending on a size of an embryo. The tomographic reconstructions were performed using GE phoenix datos|x 2.0 3D computed tomography software.
The cartilage of the olfactory system was segmented manually using Avizo - 3D image data processing software (FEI, Hillsboro, OR). The volumetric data of a segmented region were transformed to a polygonal mesh that describes the outer boundary of the region. The polygonal mesh consisting of triangles is a digital geometrical representation of the real object. The polygonal mesh of the olfactory system was imported to VG Studio MAX 2.2 software (Volume Graphics, Germany) for surface smoothening. The analysis of wall thickness at different embryonic stages was performed in order to show the differences or similarities in the thickness of the cartilage structures (Tesařová et al., 2016). The results are shown on the polygonal mesh by a colour map. The growth zones in facial chondrocranium at different stages were outlined on top of the 3D polygonal mesh based on the EdU analysis and confocal microscopy results.
Computer simulations of shape transitions of nasal capsule structure
Models were developed using the growing polarised tissue (GPT) framework and implemented in the MATLAB application GFtbox (Kennaway et al., 2011; Kuchen et al., 2012) (RRID:SCR_001622). In this method, an initial finite element mesh, also termed the canvas, is deformed during growth. The pattern of deformation depends on growth-modulating factors, whose initial distribution was established during setup. Factors have one value for each vertex and values between vertices are linearly interpolated across each finite element. In the models described here, the initial canvas is oriented with regard to the external xy-coordinate system such that the canvas base is parallel to the x-axis and the midline is parallel to the y-axis. The initial nasal capsule-line canvas consists of 1800 elements. Elements were not subdivided during the simulations.
Each model has two interconnected networks: the Polarity Regulatory Network (PRN) specifies tissue polarity and hence specified orientations of growth, and the Growth-rate Regulatory Network (KRN) determines how factors influence specified growth rates. In total, growth interactions are specified by three equations, one for the PRN and two for the KRN. These networks determine the specified polarity and growth fields across the canvas. Growth rates are influenced by factors distributed across the canvas. Growth can be promoted in a region by the pro function or inhibited by the inh function as follows:
Due to the connectedness of the canvas, this specified growth differs from the resultant growth by which the system is deformed.
Fixed midline models
Request a detailed protocolThese models question how the structure can transform given that the septum actively anchors the midline and the central groove.
Setup
Request a detailed protocolThe initial set-up phase runs from 0 to 12 time steps and during this phase the canvas deforms from a square sheet into the starting shape for the nasal capsule-like structure. Factor MID is expressed along the proximal-distal midline and used to anchor the midline vertices in the z-plane. Factor CHEEKS is expressed either side of the midline.
PRN
Request a detailed protocolA proximo-distal polarity field is set up and used to define the orientations of growth. This field is specified as being oriented parallel to the midline throughout growth by the gradient of a polarity factor, POLARISER (POL). POL has a linear gradient across the canvas with the highest level of one at the proximal base and zero at the distal tip.
KRN
Request a detailed protocolThe growth phase occurs after the initial setup phase at time step 13. During this phase there are options for specifying either isotropic growth or anisotropic growth.
During isotropic growth, the growth rate K was set to:
K = 0.05 · inh(100, iCHEEKS)
During anisotropic growth the specified growth rate parallel to the polarity field, Kpar, was defined as:
Kpar = 0.05 · inh(100, iCHEEKS)
while growth perpendicular to the polarity field, Kper, was set to zero.
Non - Fixed midline models
Request a detailed protocolThese models are aiming to simulate what happens to the shape transition when the midline and corresponding central groove are not fixed in space (and can bend or change in any other way) following tensions in the whole simulated structure. We used this approach to question how much the roof of the nasal capsule is anchored by the nasal septum.
Setup
Request a detailed protocolAs with the Fixed-midline model, an initial setup phase runs for 0–12 time-steps in which a square sheet is deformed into an alternative starting shape for the nasal capsule-like structure. In this model the proximo-distal midline was allowed to deform in the z-plane. Factor CHEEKS is expressed either side of the midline and offset slightly distally.
PRN
Request a detailed protocolA proximo-distal polarity field is set up as in the Fixed mid-ridge model.
KRN
Request a detailed protocolThe growth regulatory network is defined as in the Fixed mid-ridge model.
Mathematical model
Request a detailed protocolFor detailed description, please see the Appendix, Appendix 1—figure 1 and 2 and also (Kaucka et al., 2016)
Appendix 1
Individual based model (IBM) for cartilage dynamics
In order to model and illustrate the growth of the cartilage on the cellular level, we developed an individual-based model incorporating cell proliferation (including displacement of surrounding cells via pushing) (Hellander, 2015). The model is stochastic because we want to be able to capture effects of e.g. synchronicity in cell division and the degree of determinism needed to achieve and ordered columnar growth of the structure. We do this by letting the cell division times, direction of allocation of daughter cells after cell division, etc., be random variables. This document describes the details of the model and its implementation.
There are several popular modeling frameworks for simulating interacting cells. In the cellular potts model (Graner and Glazier, 1992) a single biological cell can be composed of multiple lattice sites making it possible to use a more detailed description of cell shape and include more detailed description of more mechanical properties. In off-lattice center based models cells are often modeled as spheres with pair-wise interactions and a force-based description to evolve the system dynamics, for an overview see (Van Liedekerke et al., 2015). Vertex models can offer even more realistic models of cell mechanics (Fletcher et al., 2013) but become expensive and complicated in three space dimensions.
Rather than these more comprehensive mechanical models, we use a simplistic rule-based, on-lattice stochastic cellular automaton (CA). In the language of a recent review (Van Liedekerke et al., 2015) our model falls into the category of a Type B CA. These types of models are widely used in e.g. cancer tumor modeling and for simulation of monolayers and spheroids (Radszuweit et al., 2009). The simulation code is written in Python, relies on the PyURDME package for spatial stochastic simulations (www.pyurdme.org) and it is freely available for download from www.github.com/ahellander/multicell (Hellander, 2015) under the GPLv3 license. A copy is archived at https://github.com/elifesciences-publications/multicell.
The basic entities in our simulation are Agents and Events. An Agent is a model (implemented as a Python class) of an individual (cell). Events simulate discrete state changes involving one or several agents. They occur at a certain time (assuming no other event involving the same agents occurs first). They rely on rules that specify how and under what conditions the event is to be executed. A simulation is initialized by creating the initial population of agents and events, and then creating a priority queue (in our case implemented using a heap data structure using the Python module’ heapq’). In each iteration of the algorithm, the event with the shortest time is popped from the queue and executed (assuming that all of its rules and conditions can be satisfied), the system time updated, new events derived from any newly created agents are created and inserted into the queue, and all existing events affected by changes in the agents or the system state are updated.
Each agent occupies one voxel of a tessellation of 2D or 3D space, and each lattice site can only accommodate one agent. Following the recommendation in (Radszuweit et al., 2009) a Delaunay triangulation is used. The mesh resolution is chosen such that the average voxel size is close to the desired cell size (∼7 µm radius) taken from the experiments. Being a lattice model, the shape and volume of the cells are a lattice property and are given by the dual grid (Voronoi cells in the case of a Delaunay triangulation). This is illustrated in Appendix 1—figure 1A. The interpretation is that individual agents occupy the dual elements (dashed lines) of an unstructured triangular (2D) or tetrahedral primal mesh (3D) (solid lines).
The individual agents – colored coded cells
Each individual cell is modeled as an individual agent with the following properties:
Color (a label used to track the lineage).
Mean cell division time, .
Variance in cell division time,
When visualized in 2D, we draw cells as polygons (the actual dual cells) and in 3D for practical reasons we visualize them as spheres centered on the vertices of the primal mesh, with radius chosen such that the volume corresponds to the volume of the corresponding dual element. On the unstructured mesh, there will be a size distribution for the mesh elements, i.e. there will be a small variation in the size associated with each lattice site, see Appendix 1— figure 1B.
Cell proliferation
Cell division time
The time until a cell, or individual, divides, is assumed to be a random variable. Although a multi-stage model of cell division can give rise to an an Erlang distribution (Radszuweit et al., 2009) which is found to match experimental data for another system, we are not calibrating our model to experimental data on cell division time distributions. The dividing cell (referred to as the mother cell) create a daughter cell after a normally distributed waiting time .
where the mean division time and the variance are parameters of the model to be supplied as input to the simulation. As a measure of the degree of variability in cell division times, we use the standard deviation over the mean,
The smaller value of , the more deterministic and synchronized the cell cycles of individual cells are. The way our model is set up, there are no events that leads to the recalculation of a cell’s division time. A cell gets assigned a time to division at creation and each cell then divide according to its internal clock irrespective of if it gets pushed etc.
After division, the daughter cell needs to be deposited on the grid. The division direction, or receiving voxel, is sampled according to a discrete distribution. In the simplest case, all directions of division are equally probable and the direction distribution is uniform. In the general case, weights are assigned according to an external, deterministic gradient.
Cell division direction
The division direction is also a random variable and the number of possible directions are given by the connections to the neighbors on the grid. Each individual has a property that sets its polarization, represented by a normalized vector pointing in the preferred direction of division. With no polarization, each possible division direction is equally probable. With polarization, the probability to divide in a certain direction is biased by the gradient. The weights for sampling the division direction are taken to be
where is the position of the vertex in the grid for which the agent resides, and is the length of the edge connecting grid points and and is a given concentration profile. The parameter dictates how perfectly the cells become polarized by the concentration profile . A value leads to equal probabilities for all directions, and very large value of means that the division direction will always be in the direction of the maximal value of the gradient (the division direction becomes deterministic in the direction of the maximal gradient in the limit . Values in between the extremes describes an increasing precision in polarization axis alignment with the gradient field.
Cell pushing
If the receiving lattice for the daughter cell site is empty (i.e. occupied by matrix), it is simply deposited there. When the daughter cell cannot be placed on a free lattice site, there is an attempt to reorganize the structure by pushing neighboring cells to make room for it. The procedure is illustrated in Appendix 1—figure 2. The probability for the displaced cell to move to a given neighboring grid point depends on the direction of pushing. Let be the vector along the edge connecting the mother cell and the daughter cell , pointing towards the daughter cell. Let be the unit vector along the edge connecting the daughter cell and one of its neighbors, . The weight for moving the displaced cell to the neighbor with index is given by
where
That is, to account for the directionality implied by the originally sampled division direction , the probability of the pushed cell to move to an adjacent lattice site is proportional to the scalar projection of the vector connecting the pushed cell and the vertex at index k onto the direction vector of the push. is a penalty parameter in the interval [0, 1] that dictates the degree of resistance in pushing the cell into an already occupied site. If, it is not possible to push a cell into an occupied site, and if , there is no resistance what so ever and only the direction of the push affects the displacement direction. Any value in between is a tradeoff between the two extremes. In Appendix 1—figure 2, for example, there is a high probability to push to , due to being almost parallel to (high ), but depending on the value of , the site may be sampled instead since that site is empty, even if the directionality contributions is smaller . Once a neighbor has been displaced, the pushed cell moves into that lattice site and the daughter cell gets deposited on the now free lattice site of the displaced cell. The displaced cell, may in turn then go on to displace additional cells and this procedure is repeated until a cell gets pushed into an unoccupied site.
Measure of order in the cartilage model
We are interested in assessing what factors are the main determinants to the degree of ordered columnar growth in the cartilage sheet patches. To that end, we postulate that a perfectly ordered structure consists of clonal columns growing straight and directed along the axis perpendicular to the initial condition starting plane (Figure 5G). We then use the following metric to quantify the degree of order in the structure
with
where is a unit vector perpendicular to the initial condition plane, and are normalized vectors joining two consecutive points (sorted by y-coordinate) in clones with the same color. is the number of cells of a given color minus one (the number of vectors), and is the number of unique clones tracked. This is illustrated graphically in Figure 5A. With this metric, a score of would mean that all columns are perfectly aligned to the main growth axis and a score of would mean that they are all perpendicular to it. We use this metric to score realizations of the process either in the absence of a gradient, or when the gradient is uniform in planes parallel to the center plane, so that in the case of a perfect polarization, cells should all deposit their daughter cells perpendicular to the center plane.
References
-
The role of the resting zone in growth plate chondrogenesisEndocrinology 143:1851–1857.https://doi.org/10.1210/endo.143.5.8776
-
Respiratory difficulties and breathing disorders in achondroplasiaPaediatric Respiratory Reviews 14:250–255.https://doi.org/10.1016/j.prrv.2013.02.009
-
The force and effect of cell proliferationThe EMBO Journal 32:2783–2784.https://doi.org/10.1038/emboj.2013.225
-
Signaling by SHH rescues facial defects following blockade in the brainDevelopmental Dynamics 241:247–256.https://doi.org/10.1002/dvdy.23726
-
Adult stem and transit-amplifying cell locationHistology and Histopathology 21:995–1027.https://doi.org/10.14670/HH-21.995
-
Regulation of cre recombinase activity by mutated estrogen receptor ligand-binding domainsBiochemical and Biophysical Research Communications 237:752–757.https://doi.org/10.1006/bbrc.1997.7124
-
Implementing vertex dynamics models of cell populations in biology within a consistent computational frameworkProgress in Biophysics and Molecular Biology 113:299–326.https://doi.org/10.1016/j.pbiomolbio.2013.09.003
-
Inactivation of Pten in osteo-chondroprogenitor cells leads to epiphyseal growth plate abnormalities and skeletal overgrowthJournal of Bone and Mineral Research 22:1245–1259.https://doi.org/10.1359/jbmr.070420
-
The control of chondrogenesisJournal of Cellular Biochemistry 97:33–44.https://doi.org/10.1002/jcb.20652
-
Simulation of biological cell sorting using a two-dimensional extended Potts modelPhysical Review Letters 69:2013–2016.https://doi.org/10.1103/PhysRevLett.69.2013
-
The development of articular cartilage: evidence for an appositional growth mechanismAnatomy and Embryology 203:469–479.https://doi.org/10.1007/s004290100178
-
Signals from the brain induce variation in avian facial shapeDevelopmental Dynamics : An Official Publication of the American Association of Anatomists.https://doi.org/10.1002/dvdy.24284
-
Fate mapping reveals origin and dynamics of lymph node follicular dendritic cellsThe Journal of Experimental Medicine 211:1109–1122.https://doi.org/10.1084/jem.20132409
-
Cartilage on the move: cartilage lineage tracing during tadpole metamorphosisDevelopment, Growth & Differentiation 54:739–752.https://doi.org/10.1111/dgd.12002
-
Histological development and dynamic expression of Bmp2-6 mRNAs in the embryonic and postnatal mouse cranial baseThe Anatomical Record. Part A, Discoveries in Molecular, Cellular, and Evolutionary Biology 288:1250–1258.https://doi.org/10.1002/ar.a.20402
-
Glial cells in the mouse enteric nervous system can undergo neurogenesis in response to injuryJournal of Clinical Investigation 121:3412–3424.https://doi.org/10.1172/JCI58200
-
Tamoxifen-inducible glia-specific cre mice for somatic mutagenesis in oligodendrocytes and schwann cellsMolecular and Cellular Neuroscience 22:430–440.https://doi.org/10.1016/S1044-7431(03)00029-0
-
Endochondral ossification: how cartilage is converted into bone in the developing skeletonThe International Journal of Biochemistry & Cell Biology 40:46–62.https://doi.org/10.1016/j.biocel.2007.06.009
-
Development and tissue origins of the mammalian cranial baseDevelopmental Biology 322:121–132.https://doi.org/10.1016/j.ydbio.2008.07.016
-
Wnt/beta-catenin signaling regulates cranial base development and growthJournal of Dental Research 87:244–249.https://doi.org/10.1177/154405910808700309
-
Endocrine regulation of the growth plateHormone Research in Paediatrics 64:157–165.https://doi.org/10.1159/000088791
-
Chondrocyte-specific knockout of the G protein G(s)alpha leads to epiphyseal and growth plate abnormalities and ectopic chondrocyte formationJournal of Bone and Mineral Research 20:663–671.https://doi.org/10.1359/JBMR.041210
-
Glassy dynamics in three-dimensional embryonic tissuesJournal of the Royal Society Interface 10:20130726.https://doi.org/10.1098/rsif.2013.0726
-
Organ shape: controlling oriented cell divisionCurrent Biology 15:R758–R759.https://doi.org/10.1016/j.cub.2005.08.053
-
The origin of the stapes and relationship to the otic capsule and oval windowDevelopmental Dynamics 241:1396–1404.https://doi.org/10.1002/dvdy.23831
-
Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel resultsComputational Particle Mechanics 2:401–444.https://doi.org/10.1007/s40571-015-0082-3
-
Endochondral ossification of the mouse nasal septumThe Anatomical Record Part A: Discoveries in Molecular, Cellular, and Evolutionary Biology 288:1163–1172.https://doi.org/10.1002/ar.a.20385
-
A Wnt5a pathway underlies outgrowth of multiple structures in the vertebrate embryoDevelopment 126:1211–1223.
Article and author information
Author details
Funding
European Molecular Biology Organization (ALTF 216-2013)
- Marketa Kaucka
Svenska Sällskapet för Medicinsk Forskning
- Marketa Kaucka
National Institutes of Health
- Andreas Hellander
Svenska Forskningsrådet Formas
- Phillip T Newton
- Andrei S Chagin
- Kaj Fried
- Igor Adameyko
Karolinska Institutet
- Phillip T Newton
- Andrei S Chagin
- Kaj Fried
- Igor Adameyko
Bertil Hållstens Forskningsstiftelse
- Igor Adameyko
Åke Wiberg Stiftelse
- Igor Adameyko
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to thank to Olga Kharchenko for illustrations. We are grateful to Prof. Susan Mackem from NIH for the Col2a1-CreERT2 strain and to Prof. Hans Clevers for the R26Confetti strain. M K was supported by EMBO Long-term Fellowship and SSMF (Svenska Sällskapet för Medicinsk Forskning) fellowship. IA and KF were supported by grants from the Swedish Research Council and Karolinska Institutet. IA received support from The Hållsten Research Foundation and Åke Wiberg Foundation. T Z and J K were supported by the Ministry of Education, Youth and Sports of the Czech Republic under the project CEITEC 2020 (LQ1601). J J and A Ha. were supported by the Ministry of Education, Youth and Sports of the Czech Republic (project LQ1605). A He. was supported by the Swedish strategic research programme eSSENCE, the Swedish Research Council (grant #2015–03964) and National Institute of Health under Award Number 1R01EB014877-01. ASC, and PN were supported by the Swedish Research Council (grant # 2016–02835) and Karolinska Institutet. A Deo lumen, ab amicis auxilium. The contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Institute of Health. All authors declare no conflict of interest.
Ethics
Animal experimentation: All animal (mouse) work has been approved and permitted by the Ethical Committee on Animal Experiments (Norra Djurförsöksetiska Nämd, ethical permit N226/15 and N5/14) and conducted according to The Swedish Animal Agency´s Provisions and Guidelines for Animal Experimentation recommendations.
Copyright
This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Metrics
-
- 3,537
- views
-
- 693
- downloads
-
- 48
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Developmental Biology
The morphogen FGF8 establishes graded positional cues imparting regional cellular responses via modulation of early target genes. The roles of FGF signaling and its effector genes remain poorly characterized in human experimental models mimicking early fetal telencephalic development. We used hiPSC-derived cerebral organoids as an in vitro platform to investigate the effect of FGF8 signaling on neural identity and differentiation. We found that FGF8 treatment increases cellular heterogeneity, leading to distinct telencephalic and mesencephalic-like domains that co-develop in multi-regional organoids. Within telencephalic regions, FGF8 affects the anteroposterior and dorsoventral identity of neural progenitors and the balance between GABAergic and glutamatergic neurons, thus impacting spontaneous neuronal network activity. Moreover, FGF8 efficiently modulates key regulators responsible for several human neurodevelopmental disorders. Overall, our results show that FGF8 signaling is directly involved in both regional patterning and cellular diversity in human cerebral organoids and in modulating genes associated with normal and pathological neural development.
-
- Developmental Biology
Wnt signaling plays crucial roles in embryonic patterning including the regulation of convergent extension (CE) during gastrulation, the establishment of the dorsal axis, and later, craniofacial morphogenesis. Further, Wnt signaling is a crucial regulator of craniofacial morphogenesis. The adapter proteins Dact1 and Dact2 modulate the Wnt signaling pathway through binding to Disheveled. However, the distinct relative functions of Dact1 and Dact2 during embryogenesis remain unclear. We found that dact1 and dact2 genes have dynamic spatiotemporal expression domains that are reciprocal to one another suggesting distinct functions during zebrafish embryogenesis. Both dact1 and dact2 contribute to axis extension, with compound mutants exhibiting a similar CE defect and craniofacial phenotype to the wnt11f2 mutant. Utilizing single-cell RNAseq and an established noncanonical Wnt pathway mutant with a shortened axis (gpc4), we identified dact1/2-specific roles during early development. Comparative whole transcriptome analysis between wildtype and gpc4 and wildtype and dact1/2 compound mutants revealed a novel role for dact1/2 in regulating the mRNA expression of the classical calpain capn8. Overexpression of capn8 phenocopies dact1/2 craniofacial dysmorphology. These results identify a previously unappreciated role of capn8 and calcium-dependent proteolysis during embryogenesis. Taken together, our findings highlight the distinct and overlapping roles of dact1 and dact2 in embryonic craniofacial development, providing new insights into the multifaceted regulation of Wnt signaling.