A multi-scale model for hair follicles reveals heterogeneous domains driving rapid spatiotemporal hair growth patterning

Abstract

The control principles behind robust cyclic regeneration of hair follicles (HFs) remain unclear. Using multi-scale modeling, we show that coupling inhibitors and activators with physical growth of HFs is sufficient to drive periodicity and excitability of hair regeneration. Model simulations and experimental data reveal that mouse skin behaves as a heterogeneous regenerative field, composed of anatomical domains where HFs have distinct cycling dynamics. Interactions between fast-cycling chin and ventral HFs and slow-cycling dorsal HFs produce bilaterally symmetric patterns. Ear skin behaves as a hyper-refractory domain with HFs in extended rest phase. Such hyper-refractivity relates to high levels of BMP ligands and WNT antagonists, in part expressed by ear-specific cartilage and muscle. Hair growth stops at the boundaries with hyper-refractory ears and anatomically discontinuous eyelids, generating wave-breaking effects. We posit that similar mechanisms for coupled regeneration with dominant activator, hyper-refractory, and wave-breaker regions can operate in other actively renewing organs.

Keywords: computational biology; developmental biology; hair follicle; mouse; pattern formation; skin; stem cells; systems biology.

Figure 1.. Model recapitulates hair cycling and its associated activator and inhibitor signaling dynamics.

(A) Schematic depiction of HF growth dynamics during telogen and anagen. Telogen and anagen HFs are shown on the left and in the center, respectively. In both hair cycle phases, Region I (purple) represents bulge and Region II (orange) represents DP with HG during telogen phase, and DP with matrix during anagen phase. On the right, schematic drawing of diffusive activator (Act. L. in green) and inhibitor (Inh. L. in red) interactions with their corresponding receptors (Act. R and Inh. R, not depicted) that form ligand-bond-receptors (Act. LR and Inh. LR) and their coupling with physical growth of the HF (blue) is shown. (B) Typical noise-free dynamics of the activator (green) and inhibitor (red) and cyclic HF growth (blue) are shown. X-axis is time in simulated days. Y-axis for activator and inhibitor shows simulated signaling levels, and for HF growth – simulated length of the HF. Grey area demarcates one modeled hair growth cycle. (C) The duration of ~anagen and ~telogen phases as the function of inhibitor signaling strengths. X-axis shows modeled inhibitor levels with ‘0’ being an arbitrary baseline levels. Y-axis shows time in simulated days. Upon stronger inhibitory signaling (high Inh. L level) ~anagen shortens (yellow) and ~telogen lengthens (purple). The entire cycle (blue) becomes longer either with stronger or weaker inhibitory signaling. When inhibitory signaling becomes either very strong or very weak, the excitability of the system breaks down and HFs equilibrate in one state (grey regions). Also see Appendix 2—tables 1, 2 and 4. (D–E’’) A total of 236 putative activator genes (green) and 122 putative inhibitor genes (red) available from a whole skin microarray dataset were identified to recapitulate temporal dynamics of the simulated activator (D) and inhibitor (E), respectively. Multiple WNT pathway members are in the putative activator gene set (D’, D’’), while BMP pathway members are among the putative inhibitor genes (E’, E’’). See gene list in Dataset 1. For all genes log-transformed, zero-mean expression profile values were calculated using colorimetric ratio-scale algorithm as reported in (Lin et al., 2009). DOI: http://dx.doi.org/10.7554/eLife.22772.003

Figure 2.. Spatiotemporal patterning of early hair cycles.

(A–D) Analysis of the whole mount dorsal skin samples from P0 (n = 3) (A) and P1 WT mice (n = 3) (A’) reveals subtle head-to-tail and lateral-to-medial hair cycle asynchronies. Asynchronies were inferred from examining the size of pigmented HFs. Larger HFs result from earlier anagen onset. (B) Heatmaps of skin samples from A and A’ built based on black pixel density (reflecting pigmented anagen HFs). (C) Quantification of anagen HFs at different phases confirms head-to-tail pattern asynchrony. Morphological definition of anagen phases used for this analysis is provided at the bottom on the panel. (D) Analysis of the whole mount dorsal skin samples from P0 220bpMsx2-hsplacZ mice, where lacZ reporter activates in anagen HFs starting from phase IIIb, confirms head-to-tail asynchrony. (E, F) Modeling rapid hair growth pattern evolution in the context of two heterogeneous domains. (E) Schematic depiction of the modeling conditions with central Dorsal domain flanked by two lateral Ventral sub-domains with coupling between Dorsal and Ventral HFs. (F) Compared to Dorsal domain HFs, Ventral domain HFs were assigned with higher levels of total available activator and inhibitor receptors, allowing shorter ~anagen and ~telogen duration. Furthermore, hair cycle asynchrony was introduced into Dorsal domain to model the initial head-to-tail asynchrony. In simulations, interactions between HFs across domain boundaries result in bilateral symmetry during the second cycle (simulated time t68-78.5). Also, initial asynchrony breaks down in the cycle 3 (t130), and partial bilateral symmetry maintains into the late cycles (see Appendix 2—video 4). Scale bars: A, A’, D – 5 mm. Images on A, A’ and D are composites. DOI: http://dx.doi.org/10.7554/eLife.22772.004

Figure 3.. Dorsal-ventral HF interactions produce bilateral symmetry.

(A) Time-lapse bioluminescence in dorsal and ventral skin of the representative Flash mouse between days P5-P48. Bioluminescent signal is color-coded according to the colorimetric scale on the right. (B) Combined temporal dynamics (from 6 Flash mice) of the bioluminescent signal-based anagen measurements over four hair cycles (days P5-P119). Dorsal skin dynamics are in brown and ventral skin dynamics are in blue. Prominent temporal advancement of ventral over dorsal anagen initiation can be seen during second, third and fourth cycles (light blue areas). This advancement is accompanied by dominant ventral-to-dorsal anagen wave spreading. (C, C’) Mapping of Flash-based anagen reveals ventral-to-dorsal hair growth wave propagation and bilateral pattern symmetry. New anagen areas for each time point are color-coded. Second anagen initiation is shown on panel C, and third anagen initiation on panel C’. Also see Appendix 1—figure 6. (D–G) Hair growth distribution patterns on P17 (D), P21 (E), P39 (F) and P55 (G). Three mice were analyzed at each time point. Inverted whole mount skin samples from representative mice are shown. Schematic pattern maps are provided with color-coded anagen (green), catagen (yellow) and telogen (red) regions. Also see Appendix 1—figures 7–12. (H, I) HF cycling dynamics in chin skin grafts remain faster compared to dorsal skin grafts. After transplantation, first anagen initiated similarly in both chin and dorsal skin grafts, however, second anagen initiated significantly faster in chin grafts. Representative chin and dorsal skin grafts are shown on (H). Combined temporal dynamics of skin grafts in anagen and telogen are shown on (I). Dorsal graft dynamics are in brown and chin graft dynamics are in blue. Temporal advancement of chin over dorsal second anagen initiation is highlighted with light blue color. Also see Appendix 1—figure 13. DOI: http://dx.doi.org/10.7554/eLife.22772.005

Figure 4.. BMP and WNT signaling differences underlie regionally specific telogen phase duration.

(A–B’’) PCA analysis reveals largely non-overlapping transcriptomic trajectories across six hair cycle stages in chin (B), ventral (B’) and dorsal domains (B’’). Combined, deconstructed PCA plots are shown with all data points marked as grey dots and domain-specific data points outlined and color-coded. Color-coding is based on the hair cycle timeline from Appendix 1—figure 14A; transcriptomic trajectories are drawn with dark lines. (C) Deconstructed PCA plot for refractory (early second) telogen is shown with domain-specific data points highlighted and color-coded based on the schematic drawing on A. (D–F) Analysis of refractory telogen data identified 1407 differentially expressed genes across the three domains (D), with each domain showing enrichment for distinct gene ontologies (E). Multiple putative hair cycle activator and inhibitor genes show domain-specific differential expression (F). Putative activators are in green and putative inhibitors are in red. For each gene, relative fold changes for ventral over chin and dorsal over chin expression levels are indicated. Genes that show cyclic expression patterns are highlighted with blue. See additional expression data analysis on Appendix 1—figures 14–19 and in Dataset 2. Asterisk marks non-canonical WNT ligand. (G, H) Analysis of Axin2-lacZ skin during second telogen reveals faster activation of WNT signaling in chin and ventral HFs over dorsal HFs. At P36 majority of HFs in chin and ventral skin have WNT-active DPs. In dorsal skin, the number of HFs with WNT-active DPs is low at P36, but increases by P42. (I) Analysis of P42 BRE-gal skin shows that many more dorsal HFs have BMP-active bulges as compared to chin and ventral HFs. Also see Appendix 1—figure 22. (J) Overexpression of Wnt7a results in disruption of ventral-to-dorsal hair growth wave dominance and spontaneous anagen appears in the dorsal domain at P60. (K) Overexpression of Bmp4 results in stalled ventral-to-dorsal hair growth wave spreading and patchy, asymmetric hair growth at P57. Scale bars: G-I – 200 um, J-K – 500 um. DOI: http://dx.doi.org/10.7554/eLife.22772.006

Figure 5.. Ear skin shows hyper-refractory properties with telogen arrested HFs.

(A–C) Morphogenesis and physiological cycling of ear HFs. (A) Analysis of ear tissue histology shows that developing HFs first appear on day P4, and progress toward mature anagen by P7. They enter catagen around P15 and first telogen by P17 (based on three mice for each time point). (B, C) Whole mount ear skin analyses show that ear HFs fail to enter second coordinated anagen and, instead, remain in an extended telogen. Seldom, isolated anagen HFs can be found (see P95 sample on B). Data are based on three mice for each time point. (D) HFs along medial side of the ear re-enter anagen after plucking (also see Appendix 1—figure 25). Experiment is based on five mice for each time point analyzed. Representative ear skin image and accompanying heatmap is shown. Heatmap criteria are shown at the bottom. (E) Unlike in dorsal skin (see Appendix 1—figure 27A), ear HFs poorly respond to topical SAG treatment. Anagen induction is limited to the medial edge of the ear. (F) Unlike in dorsal skin (see Appendix 1—figure 27B), ear HFs fail to re-enter anagen in response to topical cyclosporin A treatment. Experiments on E and F are based on three mice for each time point analyzed. Representative ear skin images and accompanying heatmaps are shown. Scale bars: A – 100 um. DOI: http://dx.doi.org/10.7554/eLife.22772.007

Figure 6.. WNT and BMP signalings modulate ear HF hyper-refractory state.

(A–C) Transcriptomes of first telogen ear skin, first telogen dorsal skin and ear cartilage/muscle complex are distinct, as revealed by PCA analysis (A). They contain 1334 differentially expressed genes (B), spanning distinct gene ontologies (C). (D, E) These tissues show differential expression of multiple ligands and antagonists for several major signaling pathways, prominently WNT and BMP. Putative activators are in green and putative inhibitors are in red. For each gene, relative fold changes for ear skin over dorsal skin and cartilage/muscle complex over dorsal skin expression levels are indicated. Select genes are highlighted. (F, F’) BRE-gal reporter reveals high BMP activity in telogen ear HFs and in the adjacent cartilage/muscle complex (n = 8). (G, G’) Axin2-lacZ reporter reveals near absence of WNT activity in ear HFs and cartilage/muscle complex. Seldom, sites of dermal reporter activity can be found (n = 8). (H–H’’) Compared to wild type mice (n = 4) (H), ears of Krt14-Noggin (n = 4) (H’) and Krt14-Wnt7a mice (n = 4) (H’’) show prominent increases in spontaneous anagen frequency. Cumulative heatmaps from four individual ear samples are shown. Also see Appendix 1—figure 28. Scale bars: F, G – 500 um; F’, G’ – 100 um. DOI: http://dx.doi.org/10.7554/eLife.22772.008

Figure 7.. Hair growth waves distort at the non-propagating boundaries.

(A, B) Introduction of a non-propagating barrier (A) or an aperture (B) into the model produces simulations with distorted anagen spreading wave front (green). (C, E) Distortions in the geometry of hair growth waves are commonly seen in the head region at the boundaries with the hyper-refractory ears and eyelids, the physical breaks in the skin. Seldom, similar distorted patterns can be seen around limb skin (D). Hair growth patterns on C-E are accompanied by color-coded schematic drawings. Colors are defined at the bottom. Hair growth distortion patterns shown were documented in ten mice each. DOI: http://dx.doi.org/10.7554/eLife.22772.009

Appendix 1—figure 1.. BMP changes affect hair length.

Grossly, Krt14-Noggin mice display longer than normal pelage (A). In contrast, Krt14-Bmp4 mice display generalized pelage shortening. Also see Appendix 1—figure 2. (B) Compared to control, hairs are longer in Krt14-Noggin and shorter in Krt14-Bmp4 mutants. Krt14-Noggin hairs lengthen by up to ~12%, while Krt14-Bmp4 hairs shorten by 35–46%. Arrows mark hair ends. Scale bars: B – 1 mm. Images on B are composites. DOI: http://dx.doi.org/10.7554/eLife.22772.013

Appendix 1—figure 2.. Changes in BMP signaling result in morphological pelage length defects.

Compared to control, Krt14-Noggin mice display generalized increase in pelage length, prominently on the tail and paws. Note fused paws, short digits, lack of claws and polydactyly, phenotypes that were reported previously. Krt14-Bmp4 mice display visibly short pelage, patches of prominent thinning on the trunk, as well as short hairs on the tail and paws. Also note other paw defects, longer digits and long, curved claws. DOI: http://dx.doi.org/10.7554/eLife.22772.014

Appendix 1—figure 3.. BMP modulates anagen phase duration.

Repression and activation of inhibitory BMP leads to increased and decreased anagen duration, respectively. Histologically, day P19 Krt14-Noggin mice display delayed (A), while day P15 Krt14-Bmp4 mice show premature anagen termination (B). Scale bars: A, B – 1 mm (whole mount) and 200 um (histology). DOI: http://dx.doi.org/10.7554/eLife.22772.015

Appendix 1—figure 4.. Downregulation of WNT leads to hair length shortening.

Krt14-Cre;Wnt7b^fl/fl mice display prominent pelage shortening during the first cycle (A’), which then becomes largely restored during the second cycle (A). Indeed, second cycle hairs in the same mutant HFs are ~15–23% longer than first cycle hairs (B’). Note that hairs on B) were dyed red. In Krt5-rtTA;tetO-Dkk1 mice, induction of Dkk1 results in prominent pelage shortening (A’’) and individual hair shorten by 22–25% as compared to non-induced control (B’’). Scale bars: B-B’’ – 1 mm. Images on B-B’’ are composites. DOI: http://dx.doi.org/10.7554/eLife.22772.016

Appendix 1—figure 5.. Downregulation of WNT leads to decreased anagen phase duration.

Histologically, day P14 Krt14-Cre;Wnt7b^fl/fl mice (A) and induced day P44 Krt5-rtTA;tetO-Dkk1 mice (B) display premature anagen termination. Scale bars: A – 1 mm (whole mount), A, B – 200 um (histology). DOI: http://dx.doi.org/10.7554/eLife.22772.017

Appendix 1—figure 6.. Spatial luminescence signal mapping across hair cycles.

Spatial maps of the luminescence signal in individual Flash mouse across three consecutive hair cycles are shown. Luminescence signal is represented as a ‘heatmap’ ranging from low (blue) to strong (red). For each time point, both ventral and dorsal views of the animal are provided. Images are grouped into two types of sets: (i) ‘propagation waves’, which show telogen-to-anagen wave spreading and (ii) ‘regression waves’, which show anagen-to-catagen-to-telogen wave spreading. During each hair cycle, both propagation and regression waves demonstrate notable temporal phase advancement in the ventral as compared to dorsal skin. Also see main Figure 3. DOI: http://dx.doi.org/10.7554/eLife.22772.018

Appendix 1—figure 7.. Distribution of the hair cycle stages in P0 mouse skin.

Day P0 skin displays prominent spatial distribution of anagen HFs (based on macroscopic evaluation of HF pigmentation). Prominent patterning features at this stage include: (a) head-to-tail anagen phase advancement in the dorsal skin (compare insert 1 with insert 5) and (b) dorsal-to-ventral anagen phase advancement reflective of the temporary delayed HF morphogenesis in the ventral skin (compare insert 2 and 4 with inserts 1 and 3). See main Figure 2A–2D for further details. DOI: http://dx.doi.org/10.7554/eLife.22772.019

Appendix 1—figure 8.. Distribution of the hair cycle stages in P2 mouse skin.

On day P2, initial head-to-tail hair cycle asynchrony in the dorsal skin becomes less prominent, while dorsal-to-ventral hair cycle asynchrony is maintained (inserts 2 and 4). DOI: http://dx.doi.org/10.7554/eLife.22772.020

Appendix 1—figure 9.. Distribution of the hair cycle stages in P7 mouse skin.

By day P7, HFs in both dorsal and ventral skin are in mature anagen. DOI: http://dx.doi.org/10.7554/eLife.22772.021

Appendix 1—figure 10.. Distribution of hair cycle stages in P15 mouse skin.

By day P15, HFs in the chin domain and most of the cranial domain undergo catagen-to-telogen transition, whereas HFs in the dorsal and ventral domains are in anagen. At this time, mouse-to-mouse pattern variability starts to become prominent, with more or less (as compared to the example pattern shown here) of the chin and cranial domains having transitioned to telogen. DOI: http://dx.doi.org/10.7554/eLife.22772.022

Appendix 1—figure 11.. Distribution of the hair cycle stages in P20 mouse skin.

By day P20, all but the most posterior portion of the dorsal domain HFs (insert 5) progressed into the first telogen. DOI: http://dx.doi.org/10.7554/eLife.22772.023

Appendix 1—figure 12.. Distribution of the hair cycle stages in P29 mouse skin.

By day P29, all ventral and most of the dorsal skin is in second anagen. One prominent exception is the cranial domain, situated in the space between eyes and ears, which at this time is still in the first telogen (insert 1). In some cases (as in the skin sample shown here), small regions of skin in the most posterior portion of the dorsal domain are also in the first telogen (insert 6). DOI: http://dx.doi.org/10.7554/eLife.22772.024

Appendix 1—figure 13.. Hair growth cycle dynamics in skin grafts.

Upon grafting onto the dorsum of pigmented SCID host mice, WT chin skin grafts (A–C) cycle faster compared to WT dorsal skin grafts (D–F). While first post-transplantation anagen starts similarly in both graft types, second anagen starts substantially faster in chin grafts. Grafts in telogen are outlined. DOI: http://dx.doi.org/10.7554/eLife.22772.025

Appendix 1—figure 14.. Domain-specific transcriptomic dynamics during hair cycle.

(A) PCA analysis reveals largely non-overlapping transcriptomic trajectories across six hair cycle stages in chin (purple), ventral (blue) and dorsal domains (green). Deconstructed PCA plots are shown with all data points marked as grey dots and domain-specific data points highlighted for each of the six hair cycle stages. (B) Multiple putative hair cycle activator and inhibitor genes show domain-specific differential expression at each hair cycle stage. Putative activators are in green and putative inhibitors are in red. For each gene, relative fold changes for ventral over chin and dorsal over chin expression levels are indicated. Genes that show cyclic expression patterns are highlighted in blue. DOI: http://dx.doi.org/10.7554/eLife.22772.026

Appendix 1—figure 15.. Double-normalized fold change calculation algorithm for RNA-seq data.

Double-normalized fold change on Appendix 1—figures 17–19 was calculated as follows: if [y/EA] <2, then Y = y/EA if [y/EA] >2, then Y = [y/EA]/2, where EA is gene expression value for early anagen. DOI: http://dx.doi.org/10.7554/eLife.22772.027

Appendix 1—figure 16.. Putative hair cycle activator and inhibitor genes.

Overlap between model-predicted activators (A) and inhibitors (B) from the microarray (see main Figure 1D–1E) and RNA-seq datasets. DOI: http://dx.doi.org/10.7554/eLife.22772.028

Appendix 1—figure 17.. Model-predicted hair cycle activators.

Multiple putative activator show cyclic temporal expression patterns. Some activator genes show model-predicted cyclic expression patterns in all three domains (A), while other activator genes show cyclic patterns only in one domain (B–D). Select genes are listed. For A), activator genes with distinctly faster expression dynamics in chin/ventral and dorsal domains are listed. (*) Relative fold change was calculated as shown on Appendix 1—figure 15. DOI: http://dx.doi.org/10.7554/eLife.22772.029

Appendix 1—figure 18.. Model-predicted hair cycle inhibitors.

Multiple putative inhibitors show cyclic temporal expression patterns. Some inhibitor genes show model-predicted cyclic expression patterns in all three domains (A), while other inhibitor genes show cyclic patterns only in one domain (B–D). Select genes are listed. For A), inhibitor genes with distinctly faster expression dynamics in chin/ventral and dorsal domains are listed. (*) Relative fold change was calculated as shown on Appendix 1—figure 15. DOI: http://dx.doi.org/10.7554/eLife.22772.030

Appendix 1—figure 19.. Additional activators and inhibitors.

Additional domain-specific activators (A–A’’) and inhibitors (B–B’’) identified on RNA-seq. (*) Relative fold change was calculated as shown on Appendix 1—figure 15. DOI: http://dx.doi.org/10.7554/eLife.22772.031

Appendix 1—figure 20.. Domain-specific expression of select markers.

Immunostaining on early telogen skin reveals thicker perilipin-expressing dermal fat layer in chin and ventral domains. Despite having a prominent subcutaneous fat layer, dorsal skin has thin dermal fat layer (B). Phalloidin staining shows distinctly thicker panniculus carnosus muscle in dorsal and ventral domains as compared to chin domain (C). αSMA staining shows more contractile cells around telogen HFs in chin domain (D). Scale bars: A-D – 50 um. DOI: http://dx.doi.org/10.7554/eLife.22772.032

Appendix 1—figure 21.. Domain-specific differences in macrophage abundance.

FACS analysis on early telogen skin reveals that chin and ventral domains contain significantly more CD11b⁺;F4/80⁺ macrophages as compared to dorsal domain. Example FACS plots are shown on A–C. Statistical analysis is shown on D and domain-specific data is color-coded as follows: chin (purple), ventral (blue) and domains (green). Average percentage of macrophages is indicated within each bar (** p<0.01). DOI: http://dx.doi.org/10.7554/eLife.22772.033

Appendix 1—figure 22.. Domain-specific expression of WNT and BMP reporters.

(A) In P36 Axin2-lacZ WNT reporter mice, the majority of chin and ventral telogen HFs contain WNT-active DPs, while many dorsal telogen HFs have WNT-negative DPs. (B) By P42, the number of telogen HFs with WNT-active DPs in dorsal skin raises to the levels compatible to ventral skin. By P42, chin HFs in some animals having already entered early anagen. (C) In P42 BRE-gal BMP reporter mice, dorsal skin contains many more HFs with BMP-active bulges as compared to chin skin. Ventral skin shows intermediate percentage of BMP-active HFs. Also see main Figure 4. Scale bars: A-D – 100 um. DOI: http://dx.doi.org/10.7554/eLife.22772.034

Appendix 1—figure 23.. Distribution of the hair cycle stages in P23 Krt14-Wnt7a mouse skin.

On day P23, Krt14-Wnt7a mouse shows disruption of ventral-to-dorsal dominance and spontaneous anagen initiation occurs in several locations: chin domain (insert 3) and several sites along the midline in the dorsal domain (inserts 2, 4 and 6). In addition, anagen is prominent in the cranial domain (insert 1). DOI: http://dx.doi.org/10.7554/eLife.22772.035

Appendix 1—figure 24.. Differences in club hair density between the domains.

(A, B) In four-month-old mice, a large portion of chin (54.6%) and ventral telogen HFs (41.7%) contain three club hairs. However, in the same mice, only 0.4% of dorsal HFs contain three club hairs and the majority of HFs only have two clubs (C). Quantification of club hair density is shown on the right. Pie charts are color-coded according to the top label: single club HFs (purple), double club HFs (blue), triple clubs HFs (green) and multiple club HFs (more than 3, orange). Scale bars: A-C – 200 um. DOI: http://dx.doi.org/10.7554/eLife.22772.036

Appendix 1—figure 25.. Anagen induction in response to hair plucking.

Club hair plucking activates anagen in dorsal skin (D) as well as in the ear skin along the medial side of the caudal pinna (A, B). Also see main Figure 5D. Hair plucking across the whole ear (C) reveals that HFs along the lateral side and in the center of the pinna either do not re-enter anagen (center) or re-enter it with a significant delay (lateral side). Hair cycle activation on D) was accessed based on the re-appearance of pigmentation. Ear hair plucking experiments on A–C) are based on five mice for each time point analyzed. Back hair plucking experiment on D) is based on five mice. Representative images and accompanying heatmaps are shown. DOI: http://dx.doi.org/10.7554/eLife.22772.037

Appendix 1—figure 26.. Lack of hair growth propagation across the ear pinna.

Anagen activation in response to plucking in the ear skin remains restricted to the plucked site (caudal pinna skin medially and rostral pinna skin laterally) and does not propagate across the ear onto the opposite side. Caudal pinna skin is the same skin sample as on main Figure 5D. DOI: http://dx.doi.org/10.7554/eLife.22772.038

Appendix 1—figure 27.. Hair growth activation in response to SAG and cyclosporin A in dorsal skin.

Topical SAG (A) and Cyclosporin A treatment (B) induce rapid hair cycle activation in dorsal telogen skin. Hair cycle activation was assessed based on the re-appearance of pigmentation, a feature of early anagen. Each experiment shown is based on three mice. Representative mice are shown. DOI: http://dx.doi.org/10.7554/eLife.22772.039

Appendix 1—figure 28.. Partial rescue of the hair cycle in the ear skin of mutant mice.

Individual ear skin samples from WT control (A), Krt14-Noggin (B) and Krt14-Wnt7a mice (C) are shown. Each image is accompanied by a heatmap color-coded according to the label in the top-right corner. (D) Statistical analysis is shown on the number of active HFs (these with early and late anagen and catagen morphologies) per ear skin sample. Average number of active HFs is indicated within each bar. DOI: http://dx.doi.org/10.7554/eLife.22772.040

Appendix 1—figure 29.. Landscape model of hair growth pattern formation.

Mouse skin is represented as a landscape with peaks (ventral and chin skin), valleys (dorsal and cranial skin) and obstacles (eye openings, ear pinnae and limbs). In analogy to water streaming downward, hair growth waves (green) preferentially spread from ventral to dorsal skin, producing bilateral symmetry. In analogy to water waves diffracting around physical obstacles, hair growth waves (red) diffract around hyper-refractory ear pinnae, limb skin, and eye openings. DOI: http://dx.doi.org/10.7554/eLife.22772.041

Appendix 2—figure 1.. Geometry of a single hair follicle.

DOI: http://dx.doi.org/10.7554/eLife.22772.043

Appendix 2—figure 2.. Illustration of the hair growth cycle phases.

DOI: http://dx.doi.org/10.7554/eLife.22772.044

Appendix 2—figure 3.. Temporal patterns of activator and ligand production.

(A) The temporal patterns of activator (green line) and inhibitor ligands (red line) produced from Appendix 2—table 4 and calculated by Appendix 2-Equations (8). (B) The temporal patterns of activator and inhibitor ligands produced from Appendix 2-Equations (9, 10). Comparing (A) to (B), we see they follow the similar patterns, indicating that our modeling of activator and ligand productions (Appendix 2-Equations (9, 10)) is reasonable. DOI: http://dx.doi.org/10.7554/eLife.22772.045

Appendix 2—figure 4.. Illustration of 2D and 3D HF model.

DOI: http://dx.doi.org/10.7554/eLife.22772.046

Appendix 2—figure 5.. When activator/inhibitor signals are either too strong or too weak (Figure 1C, grey regions), the excitability breaks down and instead of cycling, the HF reaches a stable equilibrium state.

(A) When inhibitor is very strong, signal difference per se cannot reach the threshold ${\mathrm{\Delta }}_{+}$ required for spontaneous excitation, and the HF equilibrates in a prolonged telogen-like state (Figure 1C, Equilibrium II). (B) When inhibitor is very low, it fails to catch up with the activator, and the HF equilibrates in a prolonged anagen-like state instead (Figure 1C, Equilibrium I). DOI: http://dx.doi.org/10.7554/eLife.22772.047

Appendix 2—figure 6.. Effects of inhibitor ligand level and noise on ∼anagen and ∼telogen lengths.

(A, B) Mean and variances of $\sim$anagen and $\sim$telogen durations as functions of different inhibitor ligand levels and different levels of multiplicative noise. (A) Here, multiplicative noise of medium strength was added to both activator and inhibitor equations (Appendix 2-Equations (3, 4)). Similar changes in the mean of $\sim$anagen and $\sim$telogen durations as those in the noise-free simulation on Figure 1C. (B) With a medium level of inhibitory signaling, increasing the multiplicative noise strength shortens both $\sim$anagen and $\sim$telogen durations and increases their variances. (C, D) Mean and variances of $\sim$anagen and $\sim$telogen durations as functions of different inhibitor ligand levels and different levels of additive noise, the results are qualitatively the same with (A, B) obtained from mutiplicative noise. All simulations are based on 50 samples. DOI: http://dx.doi.org/10.7554/eLife.22772.048

Appendix 2—figure 7.. Dorsal-ventral interaction simulated with additive noise with noise strength λA=λI=0.05.

Similar to the results obtained with multiplicative noise (Figure 2F) and biological observations, dorsal-ventral interaction lasts for mostly 2 or 3 cycles, after which it breaks into spontaneous waves showing some bilateral symmetries. DOI: http://dx.doi.org/10.7554/eLife.22772.049

Appendix 2—figure 8.. Dependence of ∼anagen and ∼telogen length on receptor levels.

In noise-free simulations we adjusted the levels of total available activator and inhibitor receptors ($R_{\text{tot}}^A$ and $R_{\text{tot}}^I$). This allowed differential HF sensitivity to the same ligand level changes. (A, B) Comparative changes in $\sim$anagen and $\sim$telogen duration under the conditions of low (solid bars) and high levels (checkered bars) of $R_{\text{tot}}^A=R_{\text{tot}}^A=R$ under three different levels of inhibitor signaling. (C-F) Temporal dynamics of noise-free simulations performed using varying levels of inhibitor signaling under low $R$ (C, D) vs. high $R$ regime (E, F). DOI: http://dx.doi.org/10.7554/eLife.22772.050

Appendix 2—figure 9.. Separate roles of RtotA and RtotI on ∼anagen and ∼telogen lengths.

DOI: http://dx.doi.org/10.7554/eLife.22772.051

Appendix 2—figure 10.. With inhibitory pathway from inhibitor to activator incorporated to the model, the dynamics are still stable.

(A, B) The inhibitor pathway is incorporated either as inhibitor LR inhibiting the activator ligand-receptor (A, Appendix 2-Equation (14)) or inhibitor ligands inhibiting the production of activator ligands (B, Appendix 2-Equation (15)). (C, D) In either type of inhibitory pathway, $\sim$anagen (green line) is shortened and $\sim$telogen (red line) is lengthened; moreover, as the inhibition becomes stronger (larger $a$), $\sim$anagen length quickly drops to a minimum while $\sim$telogen gradually reaches a maximum. (E) With very strong inhibition, the activator LR is pushed almost flat while the inhibitor LR still shows the hill-valley pattern, however, the HF growth (blue line) still cycles. (F) With the inhibitory pathway incorporated into our model (in the activator ligand-receptor binding with strength $a=10$), changing the inhibitor ligand levels affects the $\sim$anagen, $\sim$telogen and cycle lengths in a qualitative similar way as no pathway incorporated (Figure 1C), although the excitable regime is being switched from ${\beta }_I\in [-9,0.5]$ down to ${\beta }_I\in [-13,-4]$. DOI: http://dx.doi.org/10.7554/eLife.22772.052

Appendix 2—figure 11.. 2D simulation of dorsal HF waves, when all HFs have equal probability of anagen initiation, allowing spontaneous anagen initiations.

DOI: http://dx.doi.org/10.7554/eLife.22772.053

Appendix 2—figure 12.. 3D simulation of dorsal HF waves, when all HFs have equal probability of anagen initiation, allowing spontaneous anagen initiations.

Here we see spontaneous anagen initiations over the dorsal domain from the very first cycle, against the experimental observations of early head-to-tail wave degenerating into later spontaneous initiations showing some bilateral symmetries. The full time course is shown in Appendix 2—video 1. See Appendix 2—tables 1 - 3 and 5 for parameter values, except ${\lambda }_A={\lambda }_I=0.035$. DOI: http://dx.doi.org/10.7554/eLife.22772.054

Appendix 2—figure 13.. 2D (A) and 3D (B) simulations of dorsal HF waves, when HFs near the head have higher build-in levels of activator signaling.

Here we see highly successive steady waves characterized as the ‘target-like wave’ pattern, also against experimental observations. The full time course of the 3D simulation is shown in Appendix 2—video 2. See Appendix 2—tables 1 - 3 and 5 for parameter values, except ${\lambda }_A={\lambda }_I=0.02$ and ${\mathrm{\Delta }}_{+}=1.8$ in the 3D simulation 13B. DOI: http://dx.doi.org/10.7554/eLife.22772.055

Appendix 2—figure 14.. 2D simulation of dorsal HF waves, when HFs near the head have build-in temporal phase advancement.

DOI: http://dx.doi.org/10.7554/eLife.22772.056

Appendix 2—figure 15.. 3D modeling of HF growth pattern evolution with head-to-tail hair cycle asynchrony.

Initial asynchrony was introduced into the field by assigning a group of HFs on one end with a one-time increase in activator ligands. Simulations faithfully reproduce head-to-tail asynchrony of the first few cycles. However, this asynchrony fails to degrade quickly, and persists during multiple consecutive cycles. In addition, no clear bilateral symmetry is observed. This differs from rapid HF growth pattern evolution observed in mice. Also see Appendix 2—video 3 and see Appendix 2—tables 1 - 3 and 5 for parameter values, except ${\lambda }_A={\lambda }_I=0.035$. DOI: http://dx.doi.org/10.7554/eLife.22772.057

Appendix 2—figure 16.. The high level of inhibitor in ear skin results in hair cycling stops and HFs equilibrate in an extended telogen.

Inhibitor level is high in this whole hyper-refractory domain. Since initially the activator level is slightly higher than later, noise is able to randomly initiate the first HF growth cycle. After that, HFs stay in an extended telogen. See Appendix 2—tables 1 - 3 and 5 for parameter values, except ${\mathrm{\Delta }}_{+}=3$, ${\lambda }_A={\lambda }_I=0.05$, ${\alpha }_I^1=22$ and ${\beta }_A(t)$ given by Appendix 2-Equation (16). DOI: http://dx.doi.org/10.7554/eLife.22772.058

Appendix 2—figure 17.. Anagen waves spreading from dorsal skin could not propagate into the ear skin due to the hyper-refractory state.

We simulate a domain of dorsal skin, with a small patch (the region inside the box) having high inhibitor level in simulating the hyper-refractory ear skin. Steady HF waves propagate through the dorsal region for several cycles together with a few stochastic initiations, however, no HF wave is able to propagate into the ear skin due to the hyper-refractory state. DOI: http://dx.doi.org/10.7554/eLife.22772.059

Appendix 2—figure 18.. Simulation results of a wave-breaker region where HF growth is disabled while molecules are still allowed to diffuse in this region.

The growth-impaired region creates a distortion-like effect, which breaks the HF growth wave propagating. DOI: http://dx.doi.org/10.7554/eLife.22772.060

Appendix 2—figure 19.. Increasing Δ+ greatly increases ∼telogen length, slightly increases ∼anagen length.

DOI: http://dx.doi.org/10.7554/eLife.22772.061

Appendix 2—figure 20.. For different values of Δ+, increasing inhibitor level will always result in shorter ∼anagen and longer ∼telogen, although to different extents.

DOI: http://dx.doi.org/10.7554/eLife.22772.062

Appendix 2—figure 21.. At the single HF level, increasing Δ- greatly decreases ∼anagen length while slightly decreasing ∼telogen length.

DOI: http://dx.doi.org/10.7554/eLife.22772.063

Appendix 2—figure 22.. For different values of Δ-, increasing inhibitor level will always result in shorter ∼anagen and longer ∼telogen, although to different extents.

DOI: http://dx.doi.org/10.7554/eLife.22772.064

Appendix 2—figure 23.. With fixed activator and inhibitor levels, changing either of the four parameters (H+,H-,m+,m-) has little effect on ∼anagen/∼telogen length.

DOI: http://dx.doi.org/10.7554/eLife.22772.065

Appendix 2—figure 24.. Increasing inhibitor level will always result in shorter ∼anagen and longer ∼telogen, and the pattern does not change much under different values of H+,H-,m+,m-.

DOI: http://dx.doi.org/10.7554/eLife.22772.066

Appendix 2—figure 25.. At the single HF level, increasing hmax will shorten ∼telogen length, and will also shorten ∼anagen length within a certain region.

DOI: http://dx.doi.org/10.7554/eLife.22772.067

Appendix 2—figure 26.. For different values of hmax, increasing inhibitor level will always result in shorter ∼anagen and longer ∼telogen, although to different extents.

DOI: http://dx.doi.org/10.7554/eLife.22772.068

Appendix 2—figure 27.. Larger |hmax| (i.e., longer HF) results in a longer wave interval (B) and a wider wave (C), but the wave speed is hardly affected (A).

Simulations are noise-free, on a $60\times 100$ uniform domain, with two values of ${\mathrm{\Delta }}_{+}=2.6$ or $2.7$. DOI: http://dx.doi.org/10.7554/eLife.22772.069

Appendix 2—figure 28.. Dorsal-ventral interaction with RtotA=RtotI=15 for the dorsal while RtotA=RtotI=20 for the ventral.

First few HF waves behave similar to Figure 2F where $R_{\text{tot}}^A=R_{\text{tot}}^I=15$ for the dorsal while $R_{\text{tot}}^A=R_{\text{tot}}^I=50$. However, with the 15/20 profile, there is no clear delay between dorsal and ventral HF wave development during later cycles. DOI: http://dx.doi.org/10.7554/eLife.22772.070

Appendix 2—figure 29.. Dorsal-ventral interaction with RtotA=RtotI=15 for the dorsal while RtotA=RtotI=30 for the ventral.

First few HF waves behave similar to Figure 2F where $R_{\text{tot}}^A=R_{\text{tot}}^I=15$ for the dorsal while $R_{\text{tot}}^A=R_{\text{tot}}^I=50$. With the 15/30 profile, there is a little delay between dorsal and ventral HF wave development during later cycles. DOI: http://dx.doi.org/10.7554/eLife.22772.071

Appendix 2—figure 30.. Dorsal-ventral interaction with RtotA=RtotI=15 for the dorsal while RtotA=RtotI=40 for the ventral.

First few HF waves behave similar to Figure 2F where $R_{\text{tot}}^A=R_{\text{tot}}^I=15$ for the dorsal while $R_{\text{tot}}^A=R_{\text{tot}}^I=50$. With the 15/40 profile, there is clearly a delay between dorsal and ventral HF wave development during later cycles. DOI: http://dx.doi.org/10.7554/eLife.22772.072

Appendix 2—figure 31.. Periodic and excitable regimes in the FHN model.

(A) Nullcline $\dot{w}=0$ (red dash-dot line) intersects with $\dot{v}=0$ (green dash line) in its first piece ${\displaystyle (v_{ss}<v_0)}$. The system is excitable in this case. (B) With nullclines shown in (A), it requires strong enough noise to pass over the threshold so end the competent telogen. When noise is not strong enough, the system stays in competent telogen for a long time. Stochastic trajectory (black line) presented. (C) Nullcline $\dot{w}=0$ (red dash-dot line) intersects with $\dot{v}=0$ (green dash line) in its second piece ${\displaystyle (v_0<v_{ss}<v_1)}$. The system is periodic in this case. (D) With nullclines shown in (C), the system undergoes cyclic dynamics automatically without help from noise, as there is no threshold to pass over. Deterministic trajectory (black line) presented. DOI: http://dx.doi.org/10.7554/eLife.22772.073

Appendix 2—figure 32.. In the FHN model, changes in anagen length (solid line, TP⁢A), refractory telogen length (dashed line, TR), competent telogen length (dash-dot line, TC) to the activator (I) or inhibitor (J) production rate, under different values of positive feedback in the activator pathway (α1).

DOI: http://dx.doi.org/10.7554/eLife.22772.074

Appendix 2—figure 33.. In the FHN model, changes in anagen length (solid line, TP⁢A), refractory telogen length (dashed line, TR), competent telogen length (dash-dot line, TC) to the positive feedback in the activator pathway (α1) under different values of activator (I) and inhibitor (J) production rates.

DOI: http://dx.doi.org/10.7554/eLife.22772.075

Appendix 2—figure 34.. The nulleclines (w˙=0 red dash-dot lines, and v˙=0 green dash lines) and stochastic trajectories (black lines) under different values of I, J and α1 in FHN model.

(A, E) $I=0.015$, $J=0.017$, as in Appendix 2—figure 33D. (B, F) $I=0.0158$, $J=0.017$, as in Appendix 2—figure 33E. (C, G) $I=0.06$, $J=0.017$, as in Appendix 2—figure 33F. (D, H) $I=0.017$, $J=0.017$, as in Appendix 2—figure 33B. DOI: http://dx.doi.org/10.7554/eLife.22772.076

Appendix 2—figure 35.. Dorsal-ventral interaction simulation by the FHN model.

In comparison to Figure 2F. Dorsal temporal advancement is built-in, the HFs have shorter telogen in the ventral domain than in the dorsal. Figure shows the activator activities (white-low, black-high). Interaction between the domains with the dorsal domain showing head-to-tail asynchrony caused by the built-in temporal advancement last for several cycles, until it breaks into waves propagating from ventral to dorsal showing bilateral symmetry, which lasts during the following few cycles. Using the FHN model, we cannot obtain the quick degradation of the early dorsal head-to-tail asynchrony as early as the third or fourth cycle, as observed in experiments or simulated by our model. DOI: http://dx.doi.org/10.7554/eLife.22772.077

Appendix 2—figure 36.. Dorsal wave propagation with temporal phase advancement built-in in HFs near the head, simulated by the FHN model.

In the homogeneous dorsal domain, the early head-to-tail asynchrony quickly degrades into spontaneous initiations, and no bilateral symmetry is observed. DOI: http://dx.doi.org/10.7554/eLife.22772.078