PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of elifeeLifeRecent contentAbout eLifeFor authorsSign up for alerts
 
eLife. 2017; 6: e22772.
PMCID: PMC5610035

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

Valerie Horsley, Reviewing editor
Valerie Horsley, Yale University, United States;

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.

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

Research Organism: Mouse

eLife digest

Skin includes hundreds of thousands of hair follicles that cycle through different stages of activity. Each follicle grows hair, sometimes (in the case of long hairs like human head hair and horse tail hairs) for several years, before losing it. The follicle then goes through a resting stage before starting to grow another hair. To achieve high hair density, the follicles need to coordinate their hair-making activities. If they all worked independently from one another, bald patches would inevitably form that would compromise how effectively the skin works.

Groups of cells can communicate using a variety of chemical signals. It was not known whether cells in hair follicles from different regions of the skin rely on the same signals to communicate, and whether follicles in neighboring regions are able to ‘understand’ one another.

Through a combination of mathematical modeling and experimental results from mice, Wang, Oh et al. now show that hair follicles across the body use a common signaling system. This system consists of a pair of signals: ‘activators’ that stimulate hair growth, and ‘inhibitors’ that prevent it. The balance between these two signals affects the pattern of hair growth. For example, higher levels of activators allow fur to grow thickly on the belly of the mouse, likely to protect against heat loss and injuries from the ground. By contract, higher levels of inhibitors make the hairs on the ear sparse, which may prevent them from interfering with hearing.

There is little evidence that hair follicles on the scalp communicate in adult humans. Learning to activate and control communication between these follicles could provide a way to treat male pattern baldness and similar conditions. Understanding how hair follicles communicate may also help researchers to develop ways of regenerating other fast-renewing organs, such as the gut and bone marrow.

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

Introduction

Featuring prominent growth cycles, the hair follicle (HF) is a model system of choice for studying tissue regeneration. At the level of cellular activities, the hair growth cycle consists of three consecutive phases: anagen, phase of active proliferation; catagen, apoptotic involution phase; and telogen, relative quiescence phase (Al-Nuaimi et al., 2010; Paus and Foitzik, 2004; Schneider et al., 2009; Stenn and Paus, 2001). Cyclic regeneration is sustained by the bulge stem cells, located at the base of the permanent HF portion (Cotsarelis et al., 1990). During anagen initiation, signals from the niche, including the dermal papilla (DP), stimulate bulge stem cells and adjacent hair germ (HG) progenitors to proliferate (Enshell-Seijffers et al., 2010; Greco et al., 2009; Legrand et al., 2016). Activated progenitors generate all lower HF structures, including the outer root sheath (ORS) and hair matrix. During catagen, a widespread apoptotic program remodels the HF back toward a telogen state (Botchkarev et al., 2001b; Fessing et al., 2006; Foitzik et al., 2000; Lindner et al., 1997; Mesa et al., 2015). Conceptually, since the bulge produces downward migrating progeny (Hsu et al., 2011), it effectively serves as a progenitor source, while the matrix functions as a sink, and the ORS as a channel for progenitors transiting between them.

The signaling mechanisms that time these coordinated cellular activities during hair regeneration remain incompletely understood (Al-Nuaimi et al., 2014; Bernard, 2012; Lin et al., 2009; Paus et al., 1999). The putative ‘hair cycle clock’ is thought to be composed of one or several activator/inhibitor pairs acting to time key cycle phase transitions at set thresholds of their activities. Accordingly, cycle pace will depend on the speed at which activators and inhibitors reach their respective thresholds (Chen et al., 2015). Importantly, HFs exist as large populations and at least in the dorsal skin they interact to coordinate growth cycles (Hodgson et al., 2014; Plikus et al., 2011, Plikus and Chuong, 2008a, Plikus et al., 2008b). Such coordination implies that at least some of the activators and inhibitors should be present between HFs, in the so-called skin macro-environment. Previous work on dorsal skin indicates that BMP and WNT pathways constitute important components of the hair cycle clock. Indeed, defects in either of these pathways can dramatically change hair cycle progression (Botchkarev et al., 2001a; Botchkarev and Sharov, 2004; Choi et al., 2013; Enshell-Seijffers et al., 2010; Kandyba and Kobielak, 2014; Kandyba et al., 2013; Kobielak et al., 2003, 2007; Sharov et al., 2005, 2006), and ligands and antagonists for both pathways mediate macro-environmental coordination between HFs (Chen et al., 2014; Plikus et al., 2011; Plikus and Chuong, 2014). Additionally, FGF, PDGF, TGFβ, TNFα and other pathways can modulate hair cycle (Chen et al., 2015; Festa et al., 2011; Higgins et al., 2014; Ito et al., 2003; Kimura-Ueki et al., 2012; Leishman et al., 2013; Oshimori and Fuchs, 2012; Plikus, 2012; Rivera-Gonzalez et al., 2016). Importantly, the combined signaling activities for the above pathways partition the hair cycle in the dorsal skin into four functional phases, each with its distinct activator/inhibitor profile: propagating anagen (P), autonomous anagen (A), refractory telogen (R) and competent telogen (C) (Hodgson et al., 2014; Plikus and Chuong, 2014; Plikus et al., 2008b). Interactions between HFs enable hair regeneration across dorsal skin to self-organize into dynamic patterns. Critical for this self-organization are the following HF-to-HF interactions: P-phase HFs can induce neighboring C-phase HFs to enter anagen via diffusible activators, leading to hair growth coupling, while A-phase anagen or R-phase telogen HFs cannot couple due to high levels of inhibitors (Murray et al., 2012; Plikus et al., 2011; Plikus and Chuong, 2014; Plikus et al., 2008b). It remains unknown, however, whether this self-organization mechanism and its underlying WNT/BMP signaling activities is a general feature of all body skin or a special case for dorsal skin only.

Integrative understanding of large-scale hair regeneration requires a systems biology approach. Previous modeling on HFs include cellular automaton models (Halloy et al., 2000; Plikus et al., 2011), feedback-control model (Al-Nuaimi et al., 2012) and the FitzHugh-Nagumo (FHN) excitable medium model (Murray et al., 2012). Here, we present a unified three-dimensional and stochastic modeling framework for the HF that captures: (i) activator/inhibitor signaling dynamics in a single HF, (ii) cyclic growth of a single HF, and (iii) coupling between multiple HFs through diffusive signals. Using this model, we reveal that skin as a whole behaves as a heterogeneous regenerative field, where: (a) dominant hair cycle waves start in the ventrum, (b) propagate dorsally in a bilateral pattern, (c) stop at the boundary with hyper-refractory ear skin, and (d) break at non-propagating anatomical landmarks, such as eyelids and ears. We also show that WNT and BMP serve as a universal activator/inhibitor signaling pair, whose varying activities underlie distinct hair regeneration dynamics in all anatomical locations studied. These results provide new understanding of how the entire skin of the animal manages all of its hair regeneration.

Results

A multi-scale model recapitulates a single growing HF, as well as HF-to-HF communication

First, in modeling the geometry of a single HF, we considered four key expression sites for activator/inhibitor ligands, antagonists, and receptors along the HF axis: bulge, HG, matrix, and DP (Figure 1A). During the cycle, the bulge (assigned as Region I) remains relatively static, whereas the DP moves up and down along the HF axis. Also dynamic are the HG and matrix. The former only exists during telogen, while the latter only exists during anagen. The HG grows down to make matrix during anagen onset, whereas during catagen, the matrix collapses, and a new HG reforms. Simplistically, cyclic HG→matrix→HG dynamics are coordinated with the DP; thus, in the model we identify them jointly as Region II. Next, we considered that both regions produce signaling factors. Although a biological simplification, we assumed that Region I does so at a rather constant rate, while Region II shows distinct temporal dynamics (Appendix 2—table 4). We also assumed that Region II is essential for sending hair cycle-promoting signal(s), while Region I is the primary signal target. In short, we hypothesized that the essential temporal molecular dynamics in the HF operate as follows: Region II generates a signaling ligand (L) gradient; Region I detects it and transmits it into ligand-bound receptors (LR) that then, through a series of intermediate signaling steps not captured in the model directly (such as activities of the downstream signaling pathways and involvement of additional cell populations), regulates cyclic HF growth (Figure 1A). The molecular signaling events, either activating or inhibitory, can be summarized as:

(1)t[L] = Diffusion+Production+Reaction of L and R(2)t[LR] = Reaction of L and R +Degradation +Extra Source
Figure 1.
Model recapitulates hair cycling and its associated activator and inhibitor signaling dynamics.

where L, R and LR stand for ligands, receptors, and ligand-bound receptors, respectively. In the dynamics of LR (Equation 2), the ‘Extra Source’ describes stochastic signaling effects due to noise, and potential signaling contributions from Region I (Appendix 2-Governing equations for activators and inhibitors). As Equations 1 and 2 show, ligand-receptor interactions in the model take place only for the same signaling pathway, and no direct pathway cross-talk is set to occur. This, again, is a biological simplification. Recently, evidence for pathway interactions have emerged (Kandyba et al., 2013), and its effect is explored in Appendix 2-Possible interactions between the activator and inhibitor pathways do not qualitatively alter the HF dynamics.

Our model integrates key signaling features of the hair growth cycle: strong activator signals enhance HF growth, while strong inhibitor signals prevent it. We modeled HF growth through the spatial average of LR concentration differences between the levels of activator and inhibitor in Region I (Equation 7 in Appendix 2-Modeling HF phases by concentration difference). We assumed the hair cycle has two critical ‘checkpoints’: (i) the event in late competent telogen, when production of activator starts to increase (Chen et al., 2014; Greco et al., 2009; Oshimori and Fuchs, 2012; Plikus et al., 2008b), and (ii) the event of anagen termination, when the HF starts to involute. Thus, our model recognizes two phases determined by these checkpoints: ~anagen, starting from the moment of activator amplification until anagen termination, and ~telogen, lasting until the next activator amplification event. In the context of the conventional hair growth cycle, ~anagen incorporates the late portion of competent telogen and the entire anagen, while ~telogen includes catagen, refractory telogen and the remainder of competent telogen (Plikus et al., 2011; Plikus and Chuong, 2014; Plikus et al., 2008b) (Appendix 2-Modeling HF phases by concentration difference; Appendix 2—figure 2).

Model simulations produce several emergent behaviors. The cycle becomes autonomous – that is, it displays stable periodicity and excitability emerges naturally without a built-in ‘clock’ (Figure 1B). Cycling is maintained within a range of parameter values, allowing testing for various intrinsic and extrinsic signaling scenarios (Figure 1C). Associated with these dynamics are periodic changes in the system’s geometry – the signaling source in Region II moves cyclically. Simulations indicate that the moving HF geometry in the model is critical, greatly contributing to the regulation of the cycle. In a single HF model, activator/inhibitor diffusion occurs only along the HF axis. When a HF population is modeled, hair-to-hair communication emerges naturally as ligand diffusion from neighbors supplements intrinsic HF ligand levels. As such, hair cycle pace depends on interactive signaling between neighboring HFs – a feature that we explore below.

HF cycling emerges from the growth-mediated coupling of activator and inhibitor

Our model predicts that HF cycling occurs only within a certain range of signal strengths, that is the excitable regime (Figure 1C, white region). Within this regime, activator and inhibitor are predicted to inversely modulate duration of both ~telogen and ~anagen phases. At certain, either too high or too low signal strengths, the excitability is predicted to break down and the HF is expected to enter a non-cycling state of equilibrium (Figure 1C, grey regions). For example, when inhibitor levels are very high, the HF is predicted to equilibrate in an extended telogen (Appendix 2—figure 5A), while extended anagen is predicted for the opposite signaling condition (Appendix 2—figure 5B).

Next, we used bioinformatic and experimental approaches to validate the model’s key prediction that the same activator or inhibitor pathway can inversely modulate telogen and anagen phase duration. Considering the established roles for BMP and WNT as respective inhibitor and activator pathways regulating telogen duration in the dorsal skin, we explored if they can also regulate anagen duration in the same skin region in a model-predicted fashion. First, we found that model-predicted temporal dynamics for inhibitor and activator during ~anagen (Figure 1D and 1E) match the actual anagen expression dynamics for multiple BMP and WNT pathway members established on a highly temporally resolved whole-tissue dorsal skin microarray dataset (Lin et al., 2009) (Figure 1D–E’’; Appendix 1-Identifying model predicted hair cycle activators and inhibitors). We also show that perturbing BMP (for details see Appendix 1-Validating model-predicted roles for BMP signaling in hair cycle control) and WNT in transgenic mice (for details see Appendix 1-Validating model-predicted roles for WNT signaling in hair cycle control) alters dorsal anagen phase duration and leads to hair length defects in a way that is consistent with the model’s predictions. Overall, this data shows that our model generates biologically meaningful outcomes and that its predictive power is robust.

Model reveals skin is a heterogeneous regenerative field

Next, we set out to explore novel aspects of hair regeneration at the population level. For this purpose, we modeled a linear array of HFs (i.e. two-dimensional organization; Appendix 2—figure 4A) and a grid of HFs (i.e. three-dimensional organization; Appendix 2—figure 4B). In both cases, the diffusion of activators and inhibitors accompanying each HF during growth naturally led to HF coupling (Appendix 1-Validating model-predicted roles for BMP signaling in hair cycle control) and emergence of several known features of collective hair growth behavior, including spontaneous anagen initiation and anagen wave spreading (Appendix 2—figures 11, 12). We then focused on the phenomenon of bilaterally symmetric hair growth that is prominent in young mice (Plikus et al., 2008b) yet remains unexplained. Conventionally, first anagen in the dorsal skin of newborn mice is considered synchronous. On the other hand, adult mice display fully asynchronous and asymmetric dorsal hair growth patterns (Chen et al., 2014; Plikus and Chuong, 2008a; Plikus et al., 2009). This, however, is preceded by prominent bilateral symmetry, which often persists into the fourth hair cycle (Plikus and Chuong, 2008a). We now show that in the three-dimensional model where all HFs are assumed to be identical, full asynchrony evolves within just one cycle, and bilateral symmetry cannot be achieved (Appendix 2-Dorsal and ventral HF patterns; Appendix 2—figures 11, 12; Appendix 2—video 1). Therefore, we hypothesized that first anagen is inherently asynchronous as a result of spatially patterned HF development. Indeed, spatial distribution of early anagen HFs in the dorsal skin of newborn mice (Figure 2A–D) reveals head-to-tail and subtle lateral-to-medial asynchronies. We modeled the impact of these asynchronies on hair growth pattern evolution. Simulations reproduced head-to-tail asynchrony (Appendix 2-Dorsal and ventral HF patterns; Appendix 2—figures 14, ,15;15; Appendix 2—video 3); however, it persisted for at least 10 cycles, which is far more than the 3–4 cycles observable in mice. Moreover, prominent bilateral symmetry failed to form.

Figure 2.
Spatiotemporal patterning of early hair cycles.

We note that the above and previous simulations (Murray et al., 2012; Plikus et al., 2011) were performed on homogenous HF populations, where all HFs are assumed to be identical. We then considered that novel patterns might develop upon interaction of two or more HF populations, whose activator/inhibitor signaling levels are inherently different. In principle, dorsal skin HFs can interact with HFs from other body regions, such as ventral skin, where hair cycle dynamics are potentially distinct. Because all skin is continuous and forms an approximation of a cylinder, we modeled it as an unrolled sheet, where two Ventral sub-domains flank a rectangular Dorsal domain (Figure 2E). For initial modeling conditions (Appendix 2-Dorsal and ventral HF patterns), we considered that: (i) the first cycle on the dorsal skin has built-in head-to-tail asynchrony, and that (ii) ventral HFs develop with a 3- to 4-day delay relative to dorsal HFs (Appendix 1—figures 79). Because ventral HFs are known to produce distinctly shorter hairs (Candille et al., 2004), in the model we assumed that they have faster cycle dynamics compared to dorsal HFs (Appendix 2-Changes in the total amount of activator and inhibitor receptors results in different sensitivity of ~anagen and ~telogen lengths to signaling changes, Appendix 2-Dorsal and ventral HF patterns; Appendix 2—figures 8, ,9).9). Indeed, in this configuration, our model readily reproduced patterns with aspects of bilateral symmetry already in the second cycle as the result of dominant waves spreading from the Ventral to the Dorsal domain (Figure 2F, t68-78.5). Importantly, after the second cycle, the effect of the initial built-in head-to-tail asynchrony started to disappear. Instead, the interaction between Ventral and Dorsal HFs continued to produce prominent bilateral symmetry in the third (Figures 2F, t130-145) and later cycles (Appendix 2—video 4). Taken together, the model predicts that rapid hair growth pattern evolution requires interaction of two or more skin domains with distinct hair cycle parameters.

Ventral-dorsal interactions produce bilaterally symmetric hair growth patterns

Next, we imaged Flash mice, whose luciferase reporter produces skin-specific WNT activity signal and allows to sensitively and non-invasively determine the location and percentage of anagen HFs across the entire body (Hodgson et al., 2014) (Figure 3A–C). Luminescence levels were measured both dorsally and ventrally and mice were followed up until day P119, encompassing up to five hair growth cycles. Combined analysis from multiple mice reveals prominent phase advancement in ventral over dorsal anagen, specifically during the second, third, and fourth hair cycles (Figure 3B, blue area). Additionally, the spatial luminescence signal mapping reveals distinct ventral-to-dorsal anagen propagation with features of bilateral symmetry during second (Figure 3C; Appendix 1—figure 6) and third cycles (Figure 3C’), supporting the patterning mechanism predicted by the model (Figure 2F). We also mapped body-wide hair growth patterns on the basis of anagen HF pigmentation between days P0-P55 (Figure 3D–G; Appendix 1—figures 712). This analysis confirms ventral over dorsal anagen phase advancement starting from the second cycle and provides the following additional insights:

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

(i) Ventrally, anagen phase is the shortest in the ‘chin domain’, ending around P10. It is longer in the ‘ventral domain’ proper, ending in the genital area around P14 and in the chest area around P17 (Figure 3D).

(ii) Dorsally, anagen is the shortest in the ‘cranial domain’, ending around P14. In the ‘dorsal domain’ proper it ends as a head-to-tail wave between P15-P20 (Figure 3D; Appendix 1—figures 10, ,1111).

(iii) First ventral telogen is shorter than dorsal telogen. Second anagen initiates in the chin and ventral domains already between P21-24 and then spreads toward ventral-dorsal boundaries in form of two bilaterally symmetric waves (Figure 3E; Appendix 1—figure 12). Second anagen also ends faster in the ventral skin, maintaining ventral-dorsal asynchrony and bilateral symmetry (Figure 3F).

(iv) Third anagen initiates the fastest in the chin domain, as early as P42 (Figures 3G and and4H4H).

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

When transplanted onto the back of pigmented SCID mice, chin skin grafts (n = 8) showed faster cycling compared to dorsal skin grafts (n = 8). While first post-transplantation anagen started with similar timing in both chin and dorsal grafts, consecutive anagen started significantly faster in chin grafts (Figure 3H and I; Appendix 1—figure 13). Furthermore, in many instances, grafts induced anagen in the surrounding dorsal host skin. Taken together, these data support that dominant ventral-to-dorsal hair wave spreading drives rapid hair growth pattern evolution and bilateral symmetry. Underlying this behavior are faster hair growth cycle dynamics in chin and ventral HFs, a property that is partially maintained upon skin grafting.

Next, we asked if faster hair cycle dynamics in chin and ventral domains correlate with distinct molecular dynamics in putative activators and inhibitors. We performed RNA-seq profiling of whole skin from chin, ventral and dorsal domains at six hair cycle time points: first (aka competent) telogen, early anagen, mid-anagen, late anagen, catagen and early second (aka refractory) telogen. Analysis revealed non-overlapping transcriptomic trajectories of the hair cycle between the three domains (Figure 4A–B’’) and domain-specific expression patterns for multiple putative activator and inhibitor genes at all hair cycle time points (Appendix 1—figures 1419; Dataset 2). We then asked if refractory properties of early telogen differ between the domains. Differential gene expression analysis (Figure 4C–D) revealed enrichment in chin and ventral domains for gene ontologies related to macrophage function and lipid storage, and enrichment in chin domain for muscle-related genes (Figure 4E). Consistently, chin skin shows more contractile cells around HFs, and chin and ventral skin have thicker dermal adipose tissue and substantially more CD11b+;F4/80+ macrophages as compared to dorsal skin (Appendix 1—figures 20, ,21).21). Furthermore, dorsal early telogen skin shows gene expression changes consistent with higher refractivity – it is enriched for several BMP ligands, and depleted for BMP antagonists and WNT ligands (Figure 4F). Consistently, in Axin2-lacZ WNT reporter mice, many more HFs with WNT-active DPs are seen in chin and ventral as compared to dorsal skin at P36 (Figure 4G; Appendix 1—figure 22A). WNT activity increases in dorsal skin to the levels of ventral skin only by P42 (Figure 4H; Appendix 1—figure 22B). Furthermore, in P42 BRE-gal BMP reporter mice, many more HFs with BMP-active bulges are seen in dorsal as compared to chin and ventral skin (Figure 4I; Appendix 1—figure 22C). In Krt14-Wnt7a mice, spontaneous anagen initiation sites in dorsal skin overrun ventral-to-dorsal wave dominance (Figure 4J; Appendix 1—figure 23). In contrast, in Krt14-Bmp4 mice, ventral-dorsal hair growth waves stall and asymmetric anagen patches form instead (Figure 4K). Together, this data confirms that lower refractivity and the underlying differences in BMP and WNT activities form the bases for ventral-dorsal hair growth dominance.

Ear pinna behaves as a hyper-refractory skin domain

Our model also predicts conditions when hair cycling stops and HFs equilibrate in an extended telogen, such as due to high levels of inhibitors (Appendix 2-Hyper-refractory domain; Appendix 2—figure 5A, ,16).16). We profiled mouse skin for the existence of such behavior and found ears to match such prediction. In the ear skin, HF morphogenesis begins between days P2-P4, and HFs remain in anagen until about P15 (Figure 5A). After first anagen, and for at least three months, they remain in an extended telogen, while at the same time dorsal HFs have already reached their third cycle (Figure 5B and C). Seldom, solitary anagen HFs can be found, but no coordinated hair growth waves, characteristic to other skin regions, are observed (Figure 5B, day P95). Moreover, anagen waves spreading from cranial skin could not propagate into ear skin (Figure 7E). These observations are consistent with the possibility that ear skin is hyper-refractory. Next, we examined ear HFs’ responses to several potent anagen inducers: cyclosporin A (Maurer et al., 1997; Paus et al., 1989), smoothened agonist (SAG) (Paladini et al., 2005) and hair plucking (Chen et al., 2015). We show that while dorsal telogen HFs readily respond to cyclosporin A (Appendix 1—figure 27B), ear HFs remain quiescent even 3 weeks after treatment (Figure 5F). Anagen can be induced in response to SAG; however, this response occurs late, after 3 weeks, and remains restricted to the medial side of the ear (Figure 5E). This is contrasted by rapid SAG-induced anagen in dorsal skin (Appendix 1—figure 27A). Plucking induces anagen along the medial side of the ear; however, there is no anagen wave spreading into the unplucked region, a feature common in dorsal skin (Chen et al., 2015) (Figure 5D; Appendix 1—figures 25A-B, ,26).26). Furthermore, whole ear plucking experiments reveal very sparse anagen activation along the lateral side (Appendix 1—figure 25C). These data demonstrate that physiologically, adult ear HFs equilibrate in a hyper-refractory telogen state, yet in principle remain capable of regenerative cycling in response to selective external stimuli.

Figure 5.
Ear skin shows hyper-refractory properties with telogen arrested HFs.
Figure 7.
Hair growth waves distort at the non-propagating boundaries.

To understand how ear HF hyper-refractivity relates to activator and inhibitor signaling levels, we compared on RNA-seq refractory telogen dorsal skin with telogen ear skin and, additionally, cartilage/muscle complex, a structure unique to ears. We show that, transcriptionally, these three tissue types are distinct (Figure 6A), containing large number of differentially expressed genes (Figure 6B; Dataset 3) enriched for distinct gene ontologies (Figure 6C). Analysis of the signaling pathways implicated in the hair cycle control revealed a number of differentially expressed WNT and BMP pathway ligands and antagonists (Figure 6D). Compared to dorsal skin, ear skin is enriched for transcripts for several WNT antagonists, including Dkkl1, Dkk2 and Sfrp2, as well as collagen Col17a1, implicated in HF stem cell maintenance (Matsumura et al., 2016). Cartilage/muscle complex is prominently enriched for Bmp5, and multiple WNT antagonists, including Frzb, Sfrp2, Sfrp5 and Wif1. Additionally, it showed upregulated expression of other known hair cycle inhibitors Fgf18 (Kimura-Ueki et al., 2012; Leishman et al., 2013) and Ctgf (Liu and Leask, 2013).

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

We validated WNT and BMP changes from RNA-seq by studying relevant pathway reporters and measuring changes in ear hair cycling in mutant mouse models. Using Axin2-lacZ reporter mice, we show isolated sites of WNT activity in ear skin dermis, and a lack of activity in telogen HFs as well as in the cartilage and muscle (Figure 6G and G’). Using BRE-gal reporter mice, we show high levels of BMP activity in telogen ear HFs (in the bulge), as well as in the cartilage and muscle (Figure 6F and F’). Overexpression of the BMP antagonist Noggin in Krt14-Noggin mice partially rescued the hyper-refractory state – substantially more spontaneous anagen HFs can be found in Krt14-Noggin ears as compared to wild-type control (Figure 6H and H’; Appendix 1—figure 28B, D). Wnt7a overexpression in Krt14-Wnt7a mice also reactivated anagen in ear skin, albeit to a lesser extent compared to Noggin overexpression (Figure 6H’’; Appendix 1—figure 28C, D). Together, these results support that hyper-refractivity of ear HFs depends on higher levels of BMP ligands and WNT antagonists, in part produced by the cartilage/muscle complex (Figure 6E).

Hair growth waves distort around hyper-refractory and hairless skin regions

Lastly, our model predicts that hair growth waves can form distorted patterns around non-propagating skin regions, such as hyper-refractory hair-bearing skin or hairless skin (Figure 7A and B; Appendix 2-Hyper-refractory domain and Wave breaker; Appendix 2—figure 17; Appendix 2—video 5). We considered that pattern distortion could occur in the cranial skin at the boundaries with hyper-refractory ears and eyelids – naturally occurring physical breaks in the skin. Indeed, we observe that hair growth waves prominently break around the eyelids and ears – anagen waves propagate faster through the hair-bearing skin around the eyelids and ears, and then distort into the spaces in front of these anatomical structures (Figure 7C and E). Similar patterns are also observed for the ventral-to-dorsal hair growth wave around the limbs (Figure 7D). We conclude that distortions of hair growth waves around anatomical structures with temporary or permanent non-propagating properties contribute to rapid body-wide hair growth pattern evolution.

Discussion

Growth-regulated parallel signaling makes the HF an excitable medium

Previous mathematical models have recapitulated cycling of a single HF (Al-Nuaimi et al., 2012; Halloy et al., 2000) or in HF populations in two dimensions (Murray et al., 2012; Plikus et al., 2011). Here, we developed a multiscale model where coupling of activator and inhibitor signals with the movements of a HF in a three-dimensional space simulates cyclic growth and communication between neighboring HFs. In a single HF regime, our model faithfully predicts the effects that changes in WNT and BMP signaling can exert on the length of the anagen phase of the hair cycle.

Similar to the FHN generic excitable media model (Murray et al., 2012), our model also recapitulates several known population-level features of the HF system such as spontaneous hair growth initiation and hair wave spreading. Importantly, however, only our model allows incorporation of differential HF growth in space, a feature required for simulating heterogeneous skin properties such as interactions between skin domains with different hair cycle frequencies or the hair wave distortion effect. Thus, while the multiscale nature and non-linearity make our model more difficult to derive analytical results, its heterogeneous domain feature allows studying complex skin-wide hair growth dynamics (see Appendix 2-Comparison with FitzHugh-Nagumo (FHN) model).

HF morphogenesis across mouse skin is spatially asynchronous

Hair growth in newborn mice is commonly thought to occur simultaneously across the entire skin. In fact, we show that the first cycle is already distinctly patterned: at birth, anagen HFs in dorsal skin have head-to-tail and lateral-to-medial asynchronies, while first anagen entry by ventral HFs is delayed by approximately 3 days and proceeds as a concentric lateral-to-midline wave. Similarly delayed by 6 days are ear HFs. First anagen naturally follows the process of HF morphogenesis, which is known to be temporarily asynchronous, and to occur, at least in the dorsal skin, in three successive waves (reviewed in Clavel et al., 2012). Pattern-wise, development of HFs relies on reaction-diffusion (Sick et al., 2006) and on space-filling expansion-induction mechanisms (Cheng et al., 2014). Importantly, models for both mechanisms assume spatially synchronous HF morphogenesis. Our findings of spatial asynchrony of the first anagen indicate spatial asynchrony of HF morphogenesis. Future studies will be required to understand the modeling and signaling aspects of such phenomenon.

Hair cycle patterns evolve from the interaction of heterogeneous skin domains

Our data reveal prominent regional differences in hair cycle dynamics and show that interaction between HFs across domain boundaries drives rapid evolution of complex hair growth patterns. Specifically, we show that during early postnatal cycles, chin and ventral domains become the dominant sources of skin-wide anagen waves. Such dominant behavior of chin and ventral domains is accompanied by distinct activity dynamics for WNT and BMP, putative hair cycle activators and inhibitors, respectively. Transgenic mouse studies further confirm the functional importance of differential WNT and BMP activities in setting distinct hair growth pace across discrete anatomical skin regions. Admittedly, an in-depth follow-up study will be necessary to identify and verify the major site-specific cellular sources for WNT and BMP ligands and antagonists.

We also show that ear skin behaves as a hyper-refractory domain, where telogen HFs are resistant to anagen-inducing stimuli and cannot participate in hair growth wave propagation. We reveal that such hyper-refractivity relates to high levels of BMP ligands and WNT antagonists, in part produced by the cartilage/muscle complex, a structure unique to the ear skin. Thus, novel behaviors can be produced by the cooption of signals from new tissue modules, rather than by the modification of preexisting ones. This finding parallels the modulatory effects of non-HF cell types on the dorsal skin hair cycle, including adipose progenitors (Festa et al., 2011; Rivera-Gonzalez et al., 2016), mature adipocytes (Plikus et al., 2008b), and resident macrophages (Castellana et al., 2014; Chen et al., 2015). Finally, we show that anatomically defined structures that cannot propagate hair growth waves, namely ears and eyelids, can generate a ‘wave-breaker’ effect. Similar distortion effects are likely to occur around other anatomical structures, such as the tail and genitals, and around skin defects, such as scars, and can jointly contribute to rapid hair growth pattern evolution.

Taken together, our study reveals that the skin as a whole functions as a complex regenerative landscape with regions of fast, slow, and very slow hair renewal (Appendix 1—figure 29). We show that this behavior produces a fur coat with variable hair density, which likely serves an adaptive role, such as in thermoregulation. Mechanistically, we show that the WNT/BMP activator/inhibitor signaling pair modulates hair regeneration in all skin regions studied. This suggests that the WNT/BMP ‘molecular language’ for hair growth is general, rather than a special case for a specific body site. Its generality allows for hair-to-hair communications to arise across anatomic domain boundaries, which, in turn, enables novel hair growth dynamics not obvious from prior work – fast cycling skin regions (such as chin skin) function as a kind of hair growth pacemaker. Furthermore, our findings on ear hair cycle expand the repertoire of tissues with signaling macro-environment function to include any closely-positioned anatomic structures with signaling properties, such as cartilage.

We posit that some of the newly found hair regeneration features can have analogs in other organs. For instance, dominant anatomically defined pacemakers are common in the electrically coupled muscle-based tissues, including heart and stomach, where they generate directional contractile rhythmicity. Other actively regenerating organs, such as the intestines and bone marrow, can likely contain anatomic regions of faster and slower regeneration and, conceivably, they can be coupled to work in coordination. Knowledge learned from the skin system in the current study can guide the search for regenerative landscapes in these and other organs. Because coordination principles observed in the skin may be universal, the likelihood of them operating in other organs is substantial despite prominent anatomical differences.

Materials and methods

Computational modeling

The modeling framework is based on a hybrid approach, with individual HFs modeled as an expanding or contracting one-dimensional line and with the diffusive molecules described in reaction-diffusion equations (Appendix 2, Equations 1–4). The latter are solved using a finite difference scheme with the standard central difference approximation on the diffusion (see Appendix 2-1-dimensional (1D) HF model to Numerical methods in Appendix 2).

Experimental mouse models

Krt14-Noggin (Plikus et al., 2004, 2005), Krt14-Bmp4 (Guha et al., 2004), Krt14-Cre;Wnt7bfl/fl (Kandyba and Kobielak, 2014), Krt14-Wnt7a (Plikus et al., 2011), Krt5-rtTA;tetO-Dkk1 (Chu et al., 2004), 220bpMsx2-hsplacZ reporter (Brugger et al., 2004), BRE(Bmp response element)-gal BMP reporter (Javier et al., 2012), Axin2-lacZ (Lustig et al., 2002), and Flash WNT reporter mice (Hodgson et al., 2014) were used. For Dkk1 induction, P30 Krt5-rtTA;tetO-Dkk1 mice were placed on 2 mg/ml Doxycycline-containing water ab libitum, and skin was collected at P44 for histology and at P50 for hair length measurements.

Skin grafting

5 × 5 mm skin grafts from chin and dorsal domains of P21 C57BL/6J male mice were transplanted onto the dorsum of gender-matched pigmented P50 SCID recipients. At the time of grafting, donor skin was in first telogen and recipient skin was in second telogen.

Hair plucking

In dorsal skin, club hairs were plucked from 5 × 5 mm areas. In the ear pinna, plucking was done on the caudal skin. For quantitative plucking, approximately 500 club hairs we plucked along the medial ear side.

Topical drug treatment

Cyclosporin A: for the dorsal skin, 100 ul of Cyclosporin A solution (1, 5, and 10 mg/ml) was applied topically once a day for 7 days. For the ear pinna, caudal skin was treated with 100 ul of 10 mg/ml of Cyclosporin A once a day for 7 days. Smoothened agonist (SAG): for the dorsal skin, 120 uM of SAG in DMSO/acetone was applied topically once a day for 4 days as described (Paladini et al., 2005). For the ear pinna, caudal skin was treated with 25 ul of SAG solution once a day for 4 days.

Hair length measurements and club hair counting

Guard, awl, auchene and zigzag club hair types were photographed, traced and calibrated using Adobe Illustrator software. See Appendix 1—table 1. Club hair density was evaluated on whole-mount telogen skin samples that were pre-treated with 1 mg/mL Collagenase/Dispase and counterstained with hematoxylin.

Histology and immunohistochemistry

Histology was performed on 4% PFA-fixed sections. For BRE-gal and Axin2-lacZ specimens, whole mount lacZ staining was performed first followed by histology. The primary antibodies used were rabbit anti-keratin Krt5 (1:250, Abcam, UK), rabbit anti-perilipin (1:750; Cell Signaling), rabbit anti-αSMA (1:200; Abcam). Actin was detected with phalloidin (Alexa Fluor 488; Molecular Probes).

Whole mount in vivo bioluminescence imaging

Whole body imaging of Flash mice was performed as previously described (Hodgson et al., 2014). Briefly, mice were injected with 150 mg/kg of firefly D-luciferin substrate and imaged with the Xenogen IVIS Spectrum system.

FACS and analysis

Second telogen skin from C57BL/6J male mice was treated with Dispase to separate epidermis from dermis. Epidermis was digested with Accutase and dermis with Collagenase. Epidermal and dermal cell suspensions were combined and stained with anti-CD11b (eBioscience) and anti-F4/80 antibodies (eBioscience). Due to small tissue size, chin skin cells from three mice were combined for each experiment. FACS data were analyzed using FlowJo.

RNA-sequencing and analyses

Total RNA was isolated using the RNeasy Mini Kit (Qiagen). RNA samples with RIN >8.0 were considered for cDNA library preparation. Full-length cDNA library amplification and tagmentation was performed as previously described (Picelli et al., 2014). Libraries were multiplexed and sequenced as paired-end on an Illumina Next-Seq500 platform. Paired-end reads were aligned to the mouse genome (mm10/gencode.vM8) and quantified using the RNA-seq by Expectation-Maximization algorithm (RSEM) with standard parameters (version 1.2.25) (Li and Dewey, 2011). Samples were batch-effect corrected. EdgeR (version 3.14.0) was employed to identify differentially expressed genes (DEGs) across samples of interest. FPKM values were taken as inputs for PCA analysis and DEG analyses. Data is available at GEO: GSE85039.

Acknowledgements

MVP is supported by the NIH NIAMS grants R01-AR067273, R01-AR069653, Edward Mallinckrodt Jr. Foundation grant and Pew Charitable Trust grant. QN is supported by NSF grants: DMS 1161621, and DMS 1562176 and NIH grants: P50-GM076516, R01-GM107264, and R01-NS095355. JWO is supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (2016R1C1B1015211 and 2014R1A5A2010008). KK (Kobielak) is supported by NIH grants R01-AR061552 and National Science Centre, Poland, OPUS Grant #2015/19/B/NZ3/02948. KK (Khosrotehrani) is supported by the NHMRC Career Development Fellowship 1023371. HLL is supported by NIH NCI T32 training grant (T32-CA009054). CFGJ is supported by NSF-GRFP (DGE-1321846) and MBRS-IMSD training grant (GM055246). ZY is supported by the programs of the Major Project for Cultivation Technology of New Varieties of Genetically Modified Organisms of the Ministry of Agriculture (2014ZX08008001, 2013ZX08008-001); Beijing Nature Foundation Grant (5162018); State Key Laboratory Open Project Grant (2015SKLB6-16). BA is supported by the NIH NIAMS grants R01-AR056439. Authors thank Dr. John A. Kessler for Krt14-Bmp4 mice, Dr. Sarah E. Millar for Krt5-rtTA;tetO-Dkk1 mice, Hoang Ha, Buu Le Dao, Andy Lau, Kathleen Nguyen, Hyeon Sung Lee and Manda Nguyen for technical assistance.

Appendix 1

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

Identifying model predicted hair cycle activators and inhibitors

To find out which signaling pathway activities fit the temporal dynamics of activator and inhibitor predicted by the model (Figure 1D and E), we examined a highly temporally resolved whole-skin microarray dataset, which includes nine consecutive time points: five for anagen, three for catagen, and one for telogen (Lin et al., 2009). We identified two sets of 236 and 122 genes whose temporal dynamics recapitulate those of the simulated activator (Figure 1D’ and D’’) and inhibitor signals (Figure 1E’ and E’’; Dataset S1), respectively. Focusing on major signaling pathways, we show that multiple members of WNT and BMP are represented in each gene set. The putative activator set includes WNT ligands (Wnt2, Wnt9a, Wnt7b), soluble WNT activator r-Spondin1, multiple WNT-specific Tcf transcription factors, and soluble BMP antagonists, Follistatin and Sostdc1. The putative inhibitor set includes multiple BMP ligands (Bmp2, Bmp4, Bmp6, Bmp7), BMP receptors (Bmpr1a, Bmpr1b and Bmpr2), BMP-specific transcription factors (Smad2, Smad5 and Smad7), and soluble (Wif1, Sfrp1), transmembrane (Kremen1), and intracellular WNT antagonists (Nkd2, Prickle1).

Validating model-predicted roles for BMP signaling in hair cycle control

For BMP signaling, we examined Krt14-Bmp4 (Guha et al., 2004) and Krt14-Noggin mice (Plikus et al., 2004, 2005) overexpressing BMP ligand and soluble antagonist, respectively. Consistent with the notion of BMP acting as anagen inhibitor, Krt14-Bmp4 mice show shortened pelage, prominently on their ears, tail, and paws, and partial baldness on the trunk (Appendix 1—figure 1A, ,2).2). Fully grown dorsal hairs in Krt14-Bmp4 are significantly shorter compared to control across all hair types by 35–47% (p<0.05) (Appendix 1—figure 1B). This is accompanied by the shortened anagen duration as established by histology on day P15 (Appendix 1—figure 3B). In contrast, Krt14-Noggin mice show a general increase in hairiness, prominent on their ears, tail, and paws (Appendix 1—figure 1A, ,2),2), and significantly longer than WT dorsal hairs: 12% longer for guard (p<0.01) and 7% longer for auchene type (p<0.01) (Appendix 1—figure 1B). In parallel, we observe longer anagen phase duration as revealed by histology on day P19 (Appendix 1—figure 3A). Hair length changes in Krt14-Noggin mice were not statistically significant for the zigzag and awl types.

Validating model-predicted roles for WNT signaling in hair cycle control

For WNT signaling, we examined Krt14-Cre;Wnt7bfl/fl (Kandyba and Kobielak, 2014) and Krt5-rtTA;tetO-Dkk1 mice (Choi et al., 2013), with constitutive skin-specific deletion of the WNT ligand Wnt7b and conditional overexpression of the soluble WNT antagonist Dkk1, respectively. Wnt7b follows model-predicted activator expression pattern (Figure 1D’’) and consistent with the previous report (Kandyba and Kobielak, 2014), Wnt7b-deficient mice display short first anagen (Appendix 1—figure 5A) and short pelage (Appendix 1—figure 4A’, B’). We now show that, intriguingly, during the following hair growth cycle, Krt14-Cre;Wnt7bfl/fl mice compensate for the initial defect (Appendix 1—figure 4A), and hair length increases between cycles by an average of 22% (p<0.05) (Appendix 1—figure 4B, B’). In Krt5-rtTA;tetO-Dkk1 mice, induction of Dkk1 overexpression during telogen prominently blocks HFs in an extended telogen (Choi et al., 2013) – a scenario equivalent to ‘Equilibrium II’ in our model (Figure 1C). Induction of Dkk1 overexpression during anagen leads to early anagen termination (Appendix 1—figure 5B) and shortened hairs (Appendix 1—figure 4A’’, B’’). Compared to control, club hairs in Krt5-rtTA;tetO-Dkk1 mice across all four types shorten by 22–25% (p<0.05). Taken together, our modeling and experiments confirm that WNT and BMP serve as a core hair cycle activator/inhibitor pair, affecting both telogen and anagen phase timing.

Appendix 1—figure 1.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig1.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig2.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig3.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig4.jpg
Downregulation of WNT leads to hair length shortening.

Krt14-Cre;Wnt7bfl/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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig5.jpg
Downregulation of WNT leads to decreased anagen phase duration.

Histologically, day P14 Krt14-Cre;Wnt7bfl/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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig6.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig7.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig8.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig9.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig10.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig11.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig12.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig13.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig14.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig15.jpg
Double-normalized fold change calculation algorithm for RNA-seq data.

Double-normalized fold change on Appendix 1—figures 1719 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig16.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig17.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig18.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig19.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig20.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig21.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig22.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig23.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig24.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig25.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig26.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig27.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig28.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app1-fig29.jpg
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 1—table 1.

Numbers of club hairs used for analysis.

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

Hair typeControlKrt14-NogginKrt14-Bmp4ControlKrt14-Cre;
Wnt7b-/-
Krt5-rtTA; tetO-Dkk1 controlKrt5-rtTA; tetO-Dkk1 induced
 Auchene52432350253043
 Awl70554050253045
 Guard73484052253045
 Zigzag54644018281933

Appendix 2

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

1-dimensional (1D) HF model

Geometry of a hair follicle

The 1D HF model is on a computation domain along the z-axis: z[Z,0] (Appendix 2—figure 1). HF growth, characterized by elongation of the bottom part of the HF, ranges in a region [hmin,hmax][Z,0]. Region I is the bulge region, which does not move associatively with HF growth, and is located at [hI,hI+dI]. Region II includes HG and DP during telogen, or matrix and DP during anagen. Region II is the bottom region of the HF, hence it constantly moves in association with HF growth. Computationally, Region II is located at [h(t),h(t)+dII], where h(t) marks the bottom tip of the HP. We have the following relation: hmin=hI-dII.

Equations for activators and inhibitors

The dynamical system consists of the diffusion of activator (ActL) and inhibitor ligands (InhL), the reactions between ligands and receptors of the same species, where ActR/InhR denotes the free activator/inhibitor receptors; and ActLR/InhLR, which denotes the activator/inhibitor ligand-bound receptors. We assume that the total amount of free and bound receptors for either activator or inhibitor are preserved at any z-level, i.e.,

[ActR](z;t)+[ActLR](z;t)RtotA,[InhR](z;t)+[InhLR](z;t)RtotI

The dynamical system includes the following equations:

1. Two stochastic PDEs on the whole computational domain [Z,0], depicting the production (in region II) and the diffusion of the ligands, and their binding to corresponding receptors (in region I):

t[ActL](z,t)=DA2z2[ ActL]diffusion+IA(h(t))δII(z)production(1)+[konA[ActL](RtotA[ActLR])+koffA[ActLR]]δI(z)reactiont[InhL](z,t)=DI2z2[InhL]+II(h(t))δII(z)(2)+[konI[InhL](RtotI[InhLR])+koffI[InhLR]]δI(z)

where δI and δII indicates region I and II:

δI(z)={1if z[hI,hI+dI]0otherwiseδII(z;h(t))={1if z[h(t),h(t)+dII]0otherwise

Functions IA(h) and II(h) give the ligand production rates, which will be discussed in Appendix 2-Modeling the production of ligands.

2. Two stochastic ODEs in region I, evaluated at each z[hI,hI+dI][hmax,0], depicting the binding reaction between ligands and corresponding receptors, and the degradation of (ligand-bound) receptors:

ddt[ActLR](z;t)= konA[ActL](RtotA[ActLR])(3)(koffA+kdegA)[ActLR]+βA+σAddt[InhLR](z;t)= konI[InhL](RtotI[InhLR])(4)(koffI+kdegI)[InhLR]+βI+σI

where βA and βI represent the constant contributions due to the reactions between the ligands produced from and receptors located at region I, or any extra sources that may be added to the system; σA,σI are noise terms. In our simulations, unless specified, we use multiplicative noise, that is

σA=λA[ActLR]dωdt,σI=λI[InhLR]dωdt,

where λA,λI are the noise strength, and dt is the time step.

3. Boundary conditions for PDEs:

(5)At z=Z:[ActL]0,[InhL]0(6)At z=0:z[ActL]=0,z[InhL]=0

The parameter values related to the above equations are given in Appendix 2—table 1, temporal and geometric parameter values used in simulations are given in Appendix 2—table 3.

Modeling HF phases by concentration difference

Let Δ be the difference of the average of [ActLR] and [InhLR] in region I:

Δ(t)=[ActLR]dz[IhnLR]dz=1dIhIhI+dI[ActLR](z;t) dz1dIhIhI+dI[InhLR](z;t) dz
(7)

Modeled HF growth relies on Δ(t), subject to the following rules: during telogen, HF rests at its minimum length, meanwhile Δ progressively increases; once Δ exceeds a certain threshold Δ+ at very late telogen, the activator senses it and amplifies its own production (i.e. positive feedback). The latter results in a quick increase in Δ, which in turn induces anagen initiation and HF growth toward its maximal length. At that point in time, inhibitor starts to accumulate while activator degrades, leading to a decrease of Δ. Once Δ deceases back to Δ-=0 (i.e., activator and inhibitor levels balance out), HF starts to involute toward its original minimal length. In other words, HF completes its cycle and returns to telogen (Appendix 2—figure 2). There are two critical time points associated with these dynamics: t+, the moment when Δ reaches Δ+ and activator production amplifies; and t-, the moment of anagen termination, when Δ decreases back to Δ-. On these bases, we divided modeled hair growth cycle into two phases: the period from a t+ to its following t- and referred to as anagen, which includes very late competent telogen (C) and the full of anagen — both propagating (P) and autonomous anagen (A); and the period from a t- to the next t+ is referred to as telogen, which includes catagen (Cat), refractory (R) telogen and most of competent telogen (Appendix 2—figure 2).

Modeling the production of ligands

The production rates of ActL and InhL are modeled in relative to the follicle growth h(t) based on experimental observation presented in Appendix 2—table 4, where we qualitatively evaluate the activity strengths of activator/inhibitor ligands (L), antagonists (A), receptors (R), and ligand-bound receptors (LR), in region I and region II during different phases within a full HF cycle.

We simplify our model by eliminating the antagonist (A in Appendix 2—table 4), and estimate the net ligand production rate using the following formula:

Net Production Rate = LA+LR2
(8)

where the values of L, A and LR are from Appendix 2—table 4, depending on the phase of hair growth. Appendix 2—figure 3A shows the temporal pattern estimated from the above equation, produced by data from Appendix 2—table 4.

The production rates of activator (IA) and inhibitor (II) are given by the following equations:

(9)IA(h)= {αA+,1hmax+hhmax+hmin+αA+,0t+t<tαAtt<next t+(10)II(h)= αI0+αI1|h|

Appendix 2—figure 3B is produced from Appendix 2-Equations (9, 10). Comparing Appendix 2—figure 3A with Appendix 2—figure 3B, we see that they follow a similar pattern, suggesting that our modeling design of Appendix 2-Equations (9, 10) is reasonable.

The parameter values are given in Appendix 2—table 5.

Modeling of HF growth

h(t) gives the location of the bottom tip of the follicle, thus -h(t)=|h(t)| gives the length of the follicle. h(t) is modeled as follows, with parameter values given in Appendix 2—table 5:

1. A follicle has a minimum length -hmin=-hI+dII (i.e., the length of the follicle during telogen): h(t)hmin.

2. Upon t+, the follicle starts growing:

ddth(t)=H+(tt+)m+
(11)

3. When the follicle grows to the maximum length -hmax, the follicle stops growing.

4. Upon t-, the follicle dies:

ddth(t)=H(tt)m
(12)

until it returns to the minimum length -hmin.

For illustration of modeling of h, please see Appendix 2—figure 1 and and22.

2-dimensional (2D) and 3-dimensional (3D) HF models

The dynamics of 2D and 3D HF models are the same as the 1D HF model; however, in the 1D model we are considering a single HF in a computation domain of [Z,0], whereas in the 2D and 3D models we consider an array of HFs in a computation domain of [0,X]×[Z,0] (Appendix 2—figure 4A) and of [0,X]×[0,Y]×[Z,0] (Appendix 2—figure 4B), respectively.

The diffusion terms in Appendix 2-Equations (1, 2) will switch to:

2z2Δ={2x2+2z2in 2D2x2+2y2+2z2in 3D

The boundary conditions at z=0 and z=Z are the same as Appendix 2—Equations (5,6). In 2D model, we apply no leak boundary conditions on x=0,X:

At x=0,X:x[ActL]=x[InhL]=0

In 3D model, we apply periodic or no leak boundary conditions on x=0,X, and no leak boundary conditions on y=0,Y, depending on the detailed model. For example, when we apply periodic and no leak boundary conditions on x=0,X and y=0,Y, respectively, we have:

[ActL](x=0,y,z;t)=[ActL](x=X,y,z;t)[InhL](x=0,y,z;t)=[InhL](x=X,y,z;t)At y=0,Y:y[ActL]=y[InhL]=0

Parameter values used in 2D and 3D simulations are given in Appendix 2—table 1 - 3, ,55.

Numerical methods

In 1D and 2D model, we use implicit finite difference to discretize the diffusion term and forward Euler methods to discretize the production and reaction terms. For example, the discretized equation of Appendix 2-Equation (1) is in the form of

unm+1unmΔtDAun+1m+12unm+1+un1m+1Δz2=fnm
(13)

where unm stands for [ActL]nm, and f represents the production and reaction terms in Appendix 2-Equation (1). It is unconditionally stable, and the truncation error is O(Δt)+O(Δx2) in 1D, or δO(Δt)+O(Δx2)+O(Δy2) in 2D.

In 3D model, implicit finite difference method is not practical, since it requires solution of a system of Nx×Ny×Nz equations. Instead we use explicit finite difference method. The truncation error is δO(Δt)+O(Δx2)+O(Δy2)+O(Δz2). To make this method stable, we require that

ν=DΔt(1Δx2+1Δy2+1Δz2)12

For the discretization parameter values we used in simulations (Appendix 2—table 3), we have ν0.13, which satisfies the stability requirement.

Effects of noise on HF dynamics

As we discussed in the main text, cycling occurs only within a range of signal strengths (Figure 1C, white region), and when activator/inhibitor signals are either too strong or too weak (Figure 1C, grey regions), HF reaches a stable equilibrium state and fails to cycle. When inhibitor is very strong or activator is very weak, the signal difference threshold required for spontaneous excitation cannot be reached (Appendix 2—figure 5A), and HF equilibrates in a prolonged telogen-like state. On the other hand, when inhibitor is very low, it fails to catch up with the activator, and HF equilibrates in a prolonged anagen-like state instead (Appendix 2—figure 5B).

Within the excitable region on which the HF can cycle, we used noise-free model to show the effects of signal levels on anagen and telogen lengths: higher inhibitor level results in shorter anagen and longer telogen (Figure 1C, yellow and magenta lines). Qualitatively similar results were also obtained in the stochastic regime with multiplicative noise in both activator and inhibitor signaling (Appendix 2—figure 6A).

We investigate how noise in signaling might affect the HF growth by adding multiplicative noise with different strengths in both activator and inhibitor signaling. Simulations show that strong signaling noise shortens the average length of both anagen and telogen, and at the same time, increases their variability (Appendix 2—figure 6B). This is interesting, although not surprising, because both anagen and telogen checkpoints are largely determined by the differential level between the activator and inhibitor signals. Stronger noise is likely to increase the probability of the signal reaching the critical differential level, advancing HFs to the next phase - an irreversible process.

Next we repeat the above simulations with additive noise instead. With additive noise added to both activator and inhibitor LR equations, anagen is shortened while telogen is extended, the same as Appendix 2—figure 6A Appendix 2—figure 6B. However, the change of Std of additive noise is not as much as when it is from multiplicative noise (Appendix 2—figure 6C,D). On the other hand, for medium inhibitor level but different noise strength, stronger additive noise also leads to shorter anagen and telogen, but the change is not as large as it is for multiplicative noise - mostly due to the reason that we choose same noise strength λA and λI for both multiplicative and additive noise to make a comparison, in which case multiplicative noise would cause bigger effects on the activator and inhibitor levels.

Finally we would like to point out that we use multiplicative noise in all 2D and 3D simulations. However, if additive noise is adopted instead, the results will not change essentially. For example, we repeat the simulation of interaction between dorsal and ventral HF waves with additive noise having strength λA=λI=0.05. Similar to that with multiplicative noise (Figure 2F), interaction between dorsal and ventral HF waves last for two cycles, then breaks into more stochastic waves with bilateral symmetry (Appendix 2—figure 7).

Changes in the total amount of activator and inhibitor receptors results in different sensitivity of anagen and telogen lengths to signaling changes

While higher inhibitor levels or lower activator levels always result in shorter anagen and telogen, we find that by altering certain background parameter values, HF may achieve different sensitivity to such signal changes. One such pair of parameters we find is the total amount of activator and inhibitor receptors (RtotA and RtotI).

We investigate anagen and telogen lengths with either both low amounts of activator and inhibitor receptors (RtotA=RtotI=6), or both high amounts (RtotA=RtotI=60). The result is shown in Appendix 2—figure 8. First, under the same signaling levels, high RtotA and RtotI always result in both short anagen and telogen. This difference between receptor levels in ventral and dorsal domains was used in the 3D modeling of HF wave interactions across dorsal and ventral domains (see Appendix 2-Dorsal and ventral HF patterns); as shown by experimental observations, ventral HFs have both shorter anagen and telogen than dorsal HFs. Next, these two modeling HF types also show different sensitivity to inhibitor signaling changes. When the receptor levels are both low, both anagen and telogen phases react more dramatically to changes in signaling (Appendix 2—figure 8A,B, solid bars). For instance, after the same net increase in the inhibitor ligand value, simulated anagen phase decreased by 16 days and telogen phase increased by 23.8 days under low RtotA and RtotI regime. However, under high RtotA and RtotI regime, anagen phase decreased by 6.4 days and telogen phase increased by 14.8 days. Appendix 2—figure 8C,D exemplifies changes in activator/inhibitor LR profiles and accompanying changes in HF growth dynamics under low RtotA and RtotI regime, while Appendix 2—figure 8E,F under high RtotA and RtotI regime. Based on this, we speculate that higher RtotA and RtotI abundance confers HFs with better ability to read out signal changes, and increases the robustness of anagen and telogen durations.

We also explore all four possible combinations of low vs. high activator and inhibitor profiles. From the results in Appendix 2—figure 9 we see that each combination gives different anagen and telogen lengths, and that each combination results in different sensitivity to inhibitor level changes.

Possible interactions between the activator and inhibitor pathways do not qualitatively alter the HF dynamics

In our model the temporal patterns of activator and/or inhibitor ligand activities are dependent on HF growth, which is based on experimental observation. This is a phenomenological description, part of which may come out from any interaction between the pathways.

Recently, Kandyba et al. (2013) discovered the inhibition from BMP (inhibitor) to WNT (activator); other than this, so far there is no solid support that other interactions between activator and inhibitor exist. We incorporate this inhibitory interaction into our model in two ways:

1. the inhibitor LR inhibits the ligand-receptor binding of activator (Appendix 2—figure 10A):

ddt[ActLR](t)=konA[ActL](RtotA[ActLR])1+(a[InhLR])m(koffA+kdegA)[ActLR]+βA+σA
(14)

or

2. the inhibitor ligands inhibits the production of activator ligands (Appendix 2—figure 10B):

IAIA1+(a[InhL])m
(15)

However, in either way, please be aware that the original model which gives the phenomenological description already involves possibly existing interactions, and adding such an inhibitory feedback only enhances the effect.

In either Appendix 2-Equation (14) or (15), we run simulations with m=2,βI=-4 for a wide range of a, where small/large a corresponds to weak/strong inhibition of the inhibitor (BMP) to the activator (WNT), with a=0 indicating that the inhibition does not exist, and the dynamics return to Appendix 2-Equation (1-4). Simulation results show that in either mode of inhibition, anagen length will decrease and telogen length will increase. (Appendix 2—figure 10C,D). Moreover, as the inhibition becomes stronger, anagen length quickly reaches a minimum (Appendix 2—figure 10C,D, green lines); however, telogen length increases gradually (Appendix 2—figure 10C,D, red lines). Observing the dynamical patterns of HF growth and activator and inhibitor LR, we find that when the inhibition is extremely strong, the inhibition of the inhibitor to the activator push the activator LR almost flat (Appendix 2—figure 10E, green line), while the inhibitor LR still shows the hill-valley pattern (Appendix 2—figure 10E, red line).

When the activator is inhibited, a single HF is still able to present three dynamic regimes upon changes of signal levels: an excitable one and two equilibrium ones, the same as in Figure 1C. Appendix 2—figure 10F gives the profile of the excitable regime, from which we see that the dynamics react to inhibitor levels in a same way as before, with the only difference that the excitable regime requires a weaker inhibitor level now: βI[-13,-4] with the inhibitory feedback in L-R binding with a strength a=10, compared to βI[-9,0.5] without any feedback (Figure 1C).

2D and 3D modeling details

Dorsal and ventral HF patterns

It is observed on mouse dorsal skin that during the first 2-to-3 cycles, synchronized anagen waves propagate from head to tail, followed by more spontaneous, asynchronous waves from the fourth cycle onwards, while still showing some bilateral symmetry. To study the cause of such patterns, we propose four scenarios. The first three of them consider a uniform HF domain simulating the dorsal domain only, and they are simulated by both the 2D and 3D models. The fourth scenario involves interactions between different HF domains: that is simulating the interactions between dorsal and ventral domains; hence, it is only simulated in the 3D model. Below we list certain modeling details of the 3D simulations, with the 2D simulations designed in similar ways:

1. Consider only the dorsal domain, where all HFs have equal probability of anagen initiation. To simulate this scenario, we consider on a 60×100 domain, we add multiplicative noise terms σA and σI (λA=λI=0.035) to the activator and inhibitor LR Appendix 2-Equations (3, 4) (Appendix 2—figure 11 for 2D, Appendix 2—figure 12 and Appendix 2—video 1 for 3D simulation). In this case, we see spontaneous HF waves from the very first cycle without any bilateral symmetry.

2. Consider only the dorsal domain, where HFs near the head have higher built-in levels of activator signaling. Again we simulate this scenario on a 60×100 domain, adding multiplicative noise terms σA and σI (λA=λI=0.02) to the activator and inhibitor LR Appendix 2-Equations (3, 4). To simulate the built-in high level of activator near head, for all HFs on the top row (y=100), we give them constant increase in their activator production (IAIA+50). In addition, we set a small threshold Δ+ (Δ+=1.8) for anagen activation (Appendix 2—figure 13, Appendix 2—video 2). In this case, wee see highly successive wave patterns steadily propagate from head to tail, which may last for many cycles. This scenario is referred to as 'target pattern' (Murray et al., 2012).

3. Consider only the dorsal domain, where HFs near the head have built-in temporal phase advancement. To simulate this scenario, for all HFs on the top row (y=100), we give a one time increase (at t=0.1) in the activator production (IAIA+3000). The noise terms are the same as in the first scenarios (λA=λI=0.035). In both 2D and 3D simulations, the first few waves steadily propagate from head to tail (Appendix 2—figure 14 for 2D, Appendix 2—figure 15 and Appendix 2—video 3 for 3D simulations). Spontaneous initiations gradually show up in later waves. However, the initial head-to-tail wave pattern persists for at least 10 consecutive cycles despite constant disruption from spontaneous initiations, which violates biological observation where the steady head-to-tail pattern breaks down as early as the third wave. Moreover, 3D simulation shows no bilateral symmetry.

4. Finally we consider the scenario where dorsal HFs interact with ventral HFs, where dorsal HFs near the head have built-in temporal phase advancement than other dorsal HFs. Moreover, the first anagen entrance of ventral HFs occurs later than in dorsal HFs, and ventral HFs have shorter anagen and telogen phases. We consider a 100×100 domain, where the 60×100 area in the middle is to simulate the dorsal domain, with its HFs having less available receptors (RtotA=RtotI=15) resulting in longer anagen and telogen; the two side-domains, each with a size of 20×100 simulate the ventral domain, which flank the dorsal domain and are connected via periodic boundary condition, and HFs in them have more available receptors (RtotA=RtotI=50). In the beginning, for all HFs on the top row in the dorsal domain (20x80,y=100), we give a one time increase in the activator production (IAIA+1000), and we block the growth of ventral HFs by setting all HF length h=hmin in these areas during t14. Multiplicative noise terms σA and σI (λA=λI=0.02) are added to the activator and inhibitor LR Appendix 2-Equations (3, 4). In this case, we see wave patterns similar to those in our biological observations, shown in Figure 2F and Appendix 2—video 4 and discussed in the main text.

Hyper-refractory domain

In studying the ear pinna skin behaving as hyper-refractory HF domain, we model two phenomena: 1) the high level of inhibitor in ear skin results in hair cycling termination, and HFs equilibrate in an extended telogen; 2) anagen waves spreading from dorsal skin could not propagate into the ear skin due to its hyper-refractory state.

We model the first phenomenon on a 100×100 domain, with parameter values given in Appendix 2—table 1 (3D dorsal column), 1 and 5, except Δ+=3, λA=λI=0.05, αI1=22 and

βA(t)=2+1.8t+0.2
(16)

Noise is able to randomly initiate the first HF growth cycle, but the high inhibitor levels prevent HF wave propagating, and no other HF cycles appear, HFs stay in an extended telogen (Appendix 2—figure 16).

For the second phenomenon, that is no HF wave propagating from dorsal skin into the ear skin, we model on a 100×100 domain. Parameter values are given in Appendix 2—table 1 3D Dorsal column, with λA=λI=0.02, and the region [30,70]×[30,60] and with elevated inhibitor level (αI0=2 in this region) to simulate the ear skin, while the other part the dorsal domain. We stimulate a top-to-bottom HF wave by giving a one-time increase (at t=0.1) in the activator ligand production in the top row (IAIA+1000 at y=100). We see steady HF wave propagate through the dorsal region for several cycles along with stochastic initiations; however, no HF wave is able to propagate into the ear skin (Appendix 2—figure 17).

Wave breaker

Hairless skin domains and physical breaks in the skin can break homogenous wave spreading. To model this wave-breaker phenomenon, we consider a 100×100 domain to simulate the Cranial domain, with parameter values given in Appendix 2—table 1 (the 3D dorsal column), 3 and 5, with λA=λI=0.02. The [30,72]×[56,60] region is 'cut-off' (Figure 7A, black region), that is no HF growth permitted and no molecules are allowed to diffuse into this region. Additionally, we apply no-slip boundary conditions to the four boundaries of this strip. We stimulate the top row of HF growth at the beginning of the simulation by giving it a short-time increase in activator production (for 0t0.1, set IAIA+400 for y=100), which creates a steady HF wave propagating from top to bottom. The ‘eyelid' region breaks the wave propagation, creating the distortion-like effect which is particularly clear during the first few cycles.

In comparison, we also model another wave-breaker scenario: while still no HF growth is allowed in the wave-breaker region, molecules are allowed to diffuse into this region. This might be a similar scenario to a skin where a part of it has the HF growth impaired. We block the HF growth in this region by setting h(t)hmin at all times. Other parameter values are the same as in the previously described scenario. Simulation results show that such a growth-impaired region also creates a distortion-like effect which is quite similar to the 'eyelid' simulation results (Appendix 2—figure 18).

Sensitivity test of parameters

In this part, we test the sensitivity of several crucial parameters in our model.

Threshold Δ+

At single HF level, Δ+ has great effect on telogen length, but little effect on anagen length. For instance, whereas increasing Δ+ greatly increases telogen length, it slightly increases anagen length (Appendix 2—figure 19). For different values of Δ+, increasing inhibitor level will always result in shorter anagen and longer telogen, although to different extents (Appendix 2—figure 20).

Threshold Δ-

On the other hand, at the single HF level, increasing Δ- greatly decreases anagen length while slightly decreasing telogen length (Appendix 2—figure 21). For different values of Δ-, increasing inhibitor level will always result in shorter anagen and longer telogen, although to different extents (Appendix 2—figure 22).

Growth parameters H+,H-,m+,m-

Parameters H+,m+ determine the growth rate of a HF during phase P, while H-,m- determine its degeneration rate during catagen. The HF dynamics are not sensitive to these parameters. With fixed activator and inhibitor levels, changing either of these four parameters have little effect on anagen/telogen length (Appendix 2—figure 23). Increasing inhibitor level will always result in shorter anagen and longer telogen, and the pattern does not change much under different values of these parameters (Appendix 2—figure 24).

Maximum HF growth length hmax

At the single HF level, increasing hmax will decrease telogen length, with a relatively slight effect on anagen length (Appendix 2—figure 25). For different values of hmax, increasing inhibitor level will always result in shorter anagen and longer telogen, although to different extents (Appendix 2—figure 26).

At HF wave level, we test the effect of hmax in noise-free 3D model on a 60×100 domain, with parameter values given in Appendix 2—table 1 (3D dorsal column), 3 and 5, except λA=λI=0, and Δ+=2.6 or 2.7 as shown above. At t=0.1, we give a one time increase in the activator production (IAIA+3000) for all HFs on the top row (y=100) to generate steady head-to-tail HF waves. We tracked the time when the second HF wave reaches y=60 and y=40 (denoted as T602 and T402); when it leaves y=40 (denoted as T¯402); and when the third HF wave reaches y=40 (denoted as T403). We measured the wave speed, wave interval and wave length in the following way:

WaveSpeed=20T402T602WaveInterval=T403T402WaveLength=WaveSpeed(T¯402T402)

for different values of hmax. Simulation results show that longer HF (i.e., larger |hmax|) results in a longer wave interval and a wider wave, but the wave speed is hardly affected (Appendix 2—figure 27).

Total available receptors RtotA and RtotI

As we discussed in section Section 5, the total amount of activator and inhibitor receptors results in different anagen and telogen lengths, and we use this property to model the difference between dorsal and ventral HFs. In Figure 2F, the dorsal-ventral interaction is simulated with RtotA=RtotI=15 for the dorsal and RtotA=RtotI=50 for the ventral domains, which leads to an interaction lasting for 2–3 cycles. To test the effects of RtotA and RtotI on the dorsal-ventral interaction behavior, we keep RtotA=RtotI=15 for the dorsal while using RtotA=RtotI=20,30,40 for the ventral domain.

Simulations show that with RtotA=RtotI=20,30,40 for the ventral domain, the HF wave interactions during the first 3 cycles are similar to Figure 2F where RtotA=RtotI=50 for the same domain. However, during later cycles, with large ventral RtotA and RtotI, there is clearly a delay between ventral and dorsal HF wave development (Appendix 2—figure 28, ,29,29, ,3030).

Comparison with FitzHugh-Nagumo (FHN) model

The FitzHugh-Nagumo model is a 2D model that describes an excitable medium, and it has also been used in the study of HF behavior (Murray et al., 2012). Below we provide a detailed comparison between the FHN model and our model. While there are many similarities between the two models in predicting experimental observations, we do find there are several important differences between them at both single and population HF levels. The codes and parameters of the FHN model are based on the original ones from Murray et al. (2012).

Single HF

Main differences between our model and the FHN model at the single HF level:

1. FHN model cannot predict reasonable anagen and telogen times as observed in experiment. Plikus et al. (2008b) reported 14 and 28 days for anagen and refractory telogen for WT mice, respectively, while competent anagen may vary in the range of 0–60 days. It is reported in Murray et al. (2012) that such a time scale cannot be reached: '…in order for the proposed model to yield stochastic excitations at a rate in agreement with observations we find that Tc (competent telogen time) in the model must be of the order of 105 days…' 'A notable feature of our simulations is that competent telogen times must be of the order of 105 days such that the frequency of stochastic excitations across a population of follicles is comparable to the population scale patterns measured by Plikus et al. However, at the single follicle scale Plikus et al. have measured competent telogen times in the range of 0–60 days. In fact, when we used these much shorter competent telogen times the simulations are dominated by stochastic excitations in a manner inconsistent with population scale measurements from wild-type mice (data not shown).'

However, we would like to point out that Chen et al. (2014) use the FHN model and reach biologically observed time scales, as shown in Chen et al. (2014), with a specific choice of parameters, the authors reach TPA (propagating anagen) 15 days, TR (refractory telogen) 30 days and TC (competent telogen) 10 days. However, we fail to reproduce their results without further information about parameter values. Moreover, to reach such a state it needs to sacrifice the sensitivity of telogen length - in particular, that of competent telogen length - with respect to the change of activator level.

2. The dynamics of the FHN model are either periodic or excitable. In the periodic regime there is no competent telogen phase, while in the excitable regime sustained excitability relies on noise or exogenous activating factors, which makes the competent telogen too sensitive to certain core parameters. Our model, however, in its periodic dynamics provides competent telogen with variable length, and noise and/or external factors affect the length of different phases instead of determining the system moves on to the next cycle. In the FHN model, when the nullcline w˙=0 intersects with the other nullcline v˙=0 in its first piece, that is vss < v0, which is the scenario shown in Appendix 2—figure 31A, the system becomes an excitable medium. This is the scenario discussed in Murray et al. (2012), and to end the competent telogen, the noise has to be strong enough to pass over the threshold (Appendix 2—figure 31B). For such an excitable scenario, the mean time spent in competent telogen is estimated by Murray et al. (2012). When the noise is not strong enough, the system will stay in competent telogen for a long time. In contrast, when the inhibitor nullcline w˙=0 intersects the activator nullcline v˙=0 in its second piece (Appendix 2—figure 31C), that is v0<vss<v1, the system will deterministically undergo cycles, since there is no more threshold to pass over that relies on the help of noise. Appendix 2—figure 31D shows the deterministic trajectory (black line) which clearly goes in cycles, unlike the stochastic trajectory in Appendix 2—figure 31B which has to wait in competent telogen for a long time waiting for excitation.

To study a HF system with all four phases (PARC) using the FHN model, it has to be in the excitable regime and depend on noise and/or external factors to get the system excited. In contrast, most simulations of our model are done in the periodic regime, which also provide the competent regime that is missing from a periodic FHN model. In addition, the length of the competent regime in the periodic regime of our model can vary a lot, allowing noise greatly to adjust the length of different phases (Figure 1C).

3. Subjected to changes in the activator/inhibitor production rate (I/J), in FHN model refractory telogen reacts differently to experimental observations, while competent telogen acts too sensitively. Plikus (2012) showed that anagen, refractory and competent telogen lengths change associated to activator and/or inhibitor level changes. Such an activator/inhibitor level change can be simulated via changes of activator/inhibitor production rate I/J, or the positive feedback in the activator pathway α1 — the latter will be discussed later.

Regarding change of I, we did not find published data. For change of J, as is pointed out in Murray et al. (2012) or shown in Chen et al. (2014), 'a notable feature of the (FHN) model is that a decrease in the parameter J does not yield the reduced refractory telogen times observed by Plikus et al.'

We reproduced the phase time to I/J relation using FHN model based on parameter values provided in Murray et al. (2012), shown in Appendix 2—figure 32 for different values of α1.

When the positive feedback in the activator pathway α1 is not too low, refractory telogen length (dashed line, TR) barely changes with activator/inhibitor production rate change I/J (Appendix 2—figure 32B,C,E,F), which is stated in Murray et al. (2012) and against experimental observations. For small α1, refractory telogen length increases as activator production rate increases (Appendix 2—figure 32A) and decreases as inhibitor production rate increases (Appendix 2—figure 32D), also against biological observations.

Competent telogen length (dash-dot line, TC) reacts dramatically to activator/inhibitor production rate changes: it decreases from the order of 104-5 down to 0 quickly as I increases (Appendix 2—figure 32A,B,C) or as J decreases (Appendix 2—figure 32D,E,F). The change rate of TC is affected by the value of α1: with bigger α1 (Appendix 2—figure 32C,F) the change rate is slightly milder comparing to extremely abrupt change when α1 is small (Appendix 2—figure 31A,D). In fact, indicating by Appendix 2—figure 32, only within a small window of (I,J) could TC stay in a reasonable scale. Although the reaction trends coincide with experimental observations – TC decreases as I increases or J deceases; however, it might be too sensitive to changes of I and/or J, comparing to a reasonable change of TPA. If we restrict to the small window of (I,J) where TC behaves reasonable, then TPA reacts too insensitive to changes in I and/or J.

4. The FHN model is unable to reach an extended anagen scenario. Our model gives two equilibrium regimes besides the excitable regime (Figure 1C): one relates to extended telogen achieved with high inhibitor or low activator level (Appendix 2—figure 5A), which can also be achieved in FHN model; however, the other equilibrium regime relates to extended anagen achieved with low inhibitor or high activator level (Appendix 2—figure 5B), this scenario may explain phenomena such as human scalp hair anagen, which can last for years, or the long hair observed in angora rabbits. The FHN model is unable to simulate such an extended anagen scenario.

5. In order to reach both short anagen and telogen in the FHN model, it needs many parameter values changed, while in our model we can easily achieve this state by changing one or two parameters. Thus, our model is a convenient tool to model the interactions between dorsal and ventral domains. In vivo experiments support fast hair cycle dynamics in the ventral domain, more specifically, ventral HFs are short in both anagen and telogen lengths. While in both the FHN and our model we discover parameters that result in short anagen and long telogen (or vice versa), in our model we also have parameters that can be modified to alter anagen and telogen length, RtotA and RtotI (Appendix 2—figure 8, ,9),9), or hmax (Appendix 2—figure 25). In contrast, in the FHN model we do not know if by adjusting one or two parameters we could obtain such a fast cycling HF type. This directly affects the simulation of the interaction between ventral and dorsal domains.

Below are some similarities between the FHN model and our model at a single HF level:

1. Subjected to changes in the activator/inhibitor production rate (I/J), anagen length reacts similar to what is predicted by our model. Despite the noncoincidence of refractory and competent telogen to experimental observations and our model predictions, in FHN model, anagen length increases as activator production rate increases (Appendix 2—figure 32A-C) and decreases as inhibitor production rate increases (Appendix 2—figure 32D-F), which qualitatively coincides with our model predictions. Moreover, the change rate of anagen length to I or J is not greatly affected by the value of α1.

2. Anagen/telogen length change subjected to the positive feedback in the activator pathway α1 coincides with experimental observations. In the FHN model, if the anagen/telogen length change is adjusted via α1 instead of the production rates I/J, the simulation results coincide with experimental observations. However, since in our model there is no direct cross-talk between the activator and inhibitor pathways, there is no way to directly compare the two models regarding such a positive feedback effect.

As summarized in Murray et al. (2012), 'increased positive feedback in the activator dynamics results in the observed phenomena of faster activation wavefronts, shorter refractory and competent telogen times, unchanged anagen time, increased spontaneous initiation rates and the emergence of target patterns at the population scale.’

We reproduce the simulation results of phase length to positive feedback (α1) relation under different values of I and J, shown in Appendix 2—figure 33.

As α1 increases, anagen length (solid line, TPA) increases slightly while refractory telogen length (dashed line, TR) clearly decreases; both coincide with experimental observations.

Depending on the base values of I and J, competent telogen length shows totally different reactions to changes in α1. When J is not high compared to I (Appendix 2—figure 33A-C), the competent telogen length mostly stays at 0, that is, the system dynamics are periodic, absent of competent telogen. On the other hand, when J is relatively high to I, that is, the inhibitor level is high in the system (Appendix 2—figure 33D), competent telogen increases dramatically as α1 decreases. That is, in a system where inhibitor level (J) is high, as the positive feedback in the activator pathway (α1) gets weaker, the mean time spent for the system waiting to get excited increases extremely fast. However, there are certain values of I and J that allow the system - in particular, competent telogen reacts in a similar way to what is observed experimentally or predicted in our model (Appendix 2—figure 33E,F). In these scenarios, competent telogen length decreases at a rate comparable to that of refractory telogen as α1 increases.

In Appendix 2—figure 34, we present the nullclines and trajectories in v-w phase corresponding to values of I, J and α1 in Appendix 2—figure 33.

HF wave

Main differences between our model and the FHN model at the HF wave level:

1. In the context of the heterogeneous two-domain model (i.e. dorsal and ventral domains), the dorsal head-to-tail asynchrony resulted from dorsal-ventral interactions breaks down within approximately 2–3 cycles in our model, which is consistent with the biological observations. Similar breakdowns take many more cycles to achieve in the FHN model. Experimental observations show ventral-dorsal interactions leading to dorsal head-to-tail asynchrony in the first 2–3 cycles, which quickly degrade into spontaneous initiations yet possessing bi-lateral symmetry.

Our model simulated similar patterns (Figure 2F) - the dorsal head-to-tail asynchrony resulted from dorsal-ventral interactions lasts for approximately 2–3 cycles, after which it turns into spreading waves with centers located on the ventral-dorsal borders showing bilateral , together with spontaneous initiations. Using the FHN model, we see the similar results; however, the dorsal-ventral interaction always lasts for more cycles (Appendix 2—figure 35, dorsal-ventral interactions persists for five cycles). We tried different parameter values but the dorsal-ventral interaction always lasts for more than three cycles, as is observed experimentally or predicted by our model.

A possible reason behind this might be that while our model allows us to modulate anagen and telogen lengths by adjusting RtotA and RtotI to approximate ventral HF behavior, that is shorter anagen and telogen, the FHN model has no easy way to modulate them. Therefore, we have to set HFs in ventral domain to have shorter telogen in sacrifice of slightly longer anagen compared to dorsal HFs.

2. In the wave-breaker simulations, we are able to investigate two scenarios using our model: a physical region cut-off from the HF wave domain (like eyes, ears), or a region permeable to signal molecules with disabled HF growth (by wound, for example). The FHN model is unable to distinguish these two scenarios.

Please refer to Appendix 2-Wave breaker, Figure 7A and Appendix 2—figure 18.

3. Experiments also imply ventral HFs are shorter than dorsal HFs, combined with our single HF analysis that within a region, short HF have shorter anagen and telogen (Appendix 2—figure 25), in the future extension of our study, we will investigate the interaction between dorsal and ventral domains characterized by different HF length. The FHN model does not permit for this.

4. In a homogeneous domain with HFs near the head having temporal phase advancement built-in, both models reproduce the head-to-tail hair cycle asynchrony during early cycles, which will degenerate into random initiations with no bilateral symmetry (Appendix 2—figure 14, ,15,15, Appendix 2—video 3, Appendix 2—figure 36). However, the early head-to-tail steady pattern fades quickly in the FHN model while in ours it degenerates gradually. Simulations show that in our model simulations, even after 10 cycles, the wave may still resemble the early head-to-tail pattern. This is because we are sitting in a periodic state with non-zero competent telogen length, and noise is still able to initiate waves allowing the system show excitable properties, and the memory of an early wave pattern can persist through later waves.

On the other hand, in the FHN model, the early head-to-tail pattern may be lost as early as the 2nd cycle. Appendix 2—figure 36 shows the time course of a typical simulation by the FHN model, in this simulation, the head-to-tail pattern persists for the first three cycles, while stochastic initiation takes over starting in the fourth cycle. The reason for such a quick loss of the early memory in the FHN model is because for the early head-to-tail pattern to survive, the system should possess a periodic property; however, the system is in an excitable regime where stochastic effects dominate, hence the stochastic initiations quickly wash away the head-to-tail pattern.

Below are some similarities between the FHN model and our model at HF wave level:

1. In a homogeneous domain, if all HFs are synchronized in telogen, they have equal probability of anagen initiation, in both models we see stochastic initiations dominate and there is no bilateral symmetry (Appendix 2—figure 11, 12, Murray et al., 2012).

2. In a homogeneous domain, when HFs near the head have higher built-in levels of activator signaling, both models show highly successive waves, that is target-like pattern (Appendix 2—figure 13, Figure S4 in Murray et al., 2012).

Supplementary figures

Appendix 2—figure 1.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig1.jpg
Geometry of a single hair follicle.

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

Appendix 2—figure 2.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig2.jpg
Illustration of the hair growth cycle phases.

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

Appendix 2—figure 3.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig3.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig4.jpg
Illustration of 2D and 3D HF model.

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

Appendix 2—figure 5.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig5.jpg
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 Δ+ 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig6.jpg
Effects of inhibitor ligand level and noise on anagen and telogen lengths.

(A, B) Mean and variances of anagen and 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 anagen and 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 anagen and telogen durations and increases their variances. (C, D) Mean and variances of anagen and 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig7.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig8.jpg
Dependence of anagen and telogen length on receptor levels.

In noise-free simulations we adjusted the levels of total available activator and inhibitor receptors (RtotA and RtotI). This allowed differential HF sensitivity to the same ligand level changes. (A, B) Comparative changes in anagen and telogen duration under the conditions of low (solid bars) and high levels (checkered bars) of RtotA=RtotA=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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig9.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig10.jpg
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, anagen (green line) is shortened and telogen (red line) is lengthened; moreover, as the inhibition becomes stronger (larger a), anagen length quickly drops to a minimum while 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 anagen, telogen and cycle lengths in a qualitative similar way as no pathway incorporated (Figure 1C), although the excitable regime is being switched from βI[-9,0.5] down to βI[-13,-4].

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

Appendix 2—figure 11.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig11.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig12.jpg
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 and55 for parameter values, except λA=λI=0.035.

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

Appendix 2—figure 13.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig13.jpg
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 and55 for parameter values, except λA=λI=0.02 and Δ+=1.8 in the 3D simulation 13B.

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

Appendix 2—figure 14.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig14.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig15.jpg
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 and55 for parameter values, except λA=λI=0.035.

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

Appendix 2—figure 16.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig16.jpg
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 and55 for parameter values, except Δ+=3, λA=λI=0.05, αI1=22 and βA(t) given by Appendix 2-Equation (16).

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

Appendix 2—figure 17.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig17.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig18.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig19.jpg
Increasing Δ+ greatly increases telogen length, slightly increases anagen length.

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

Appendix 2—figure 20.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig20.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig21.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig22.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig23.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig24.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig25.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig26.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig27.jpg
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×100 uniform domain, with two values of Δ+=2.6 or 2.7.

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

Appendix 2—figure 28.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig28.jpg
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 RtotA=RtotI=15 for the dorsal while RtotA=RtotI=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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig29.jpg
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 RtotA=RtotI=15 for the dorsal while RtotA=RtotI=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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig30.jpg
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 RtotA=RtotI=15 for the dorsal while RtotA=RtotI=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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig31.jpg
Periodic and excitable regimes in the FHN model.

(A) Nullcline w˙=0 (red dash-dot line) intersects with v˙=0 (green dash line) in its first piece (vss < v0). 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 w˙=0 (red dash-dot line) intersects with v˙=0 (green dash line) in its second piece (v0 < vss <v1). 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig32.jpg
In the FHN model, changes in anagen length (solid line, TPA), 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig33.jpg
In the FHN model, changes in anagen length (solid line, TPA), 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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig34.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig35.jpg
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.

An external file that holds a picture, illustration, etc.
Object name is elife-22772-app2-fig36.jpg
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

Supplementary tables

Appendix 2—table 1.

Parameter values in the reaction-diffusion system: Appendix 2-Equations (1–4). (1) In Figure 1C, βI ranges from −11 to 2, and Figure 1B is obtained under βI=0; both are noise-free, thus λA=λI=0. (2) In Appendix 2—figure 6, low, medium and high inhibitor L correspond to βI=-6,-3,0, respectively, and weak, medium and strong noise correspond to λA=λI=2,4,6 – the same for both multiplicative and additive noise. (3) In Appendix 2—figure 8 and and9,9, low R and high R correspond to RtotA=RtotI=R=6 or 60, respectively, while low, medium and high inhibitor L correspond to βI=-6,-3,0, respectively, which is the same as in Appendix 2—figure 6. Appendix 2—figure 8 and and99 again show noise-free results, hence λA=λI=0. Also see Appendix 2—table 2 for supplementary parameter values.

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

Reaction - Diffusion System
Parameters1D2D3D
DorsalVentral
Preservation of ReceptorsRtotA10101550
RtotI
DiffusionDA10.50.5
DI
ReactionkonA0.3
koffA0
konI0.5
koffI0
DegradationkdegA2
kdegI46.57
Concentration ThresholdΔ+1.52.252.6
Δ-0
NoiseλA2/4/63.2See Appendix 2—table 2.
λI
OthersβA4
βI−11 22.3

Appendix 2—table 2.

Parameter values in the reaction-diffusion system in 3D simulations, supplementary to Appendix 2—table 1.

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

Appendix 2—table 3.

Values of temporal and geometric parameters.

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

Time and geometry
Parameters1D2D3D
dt0.10.01
dxNA1.52
dyNA2
dz0.10.2
hI-1
dI0.50.6
dII0.8
hmax-6
hmin−1.8
XNA9060/100
YNA100
Z−10

Appendix 2—table 4.

Activity strength of activator/inhibitor ligands (L), antagonists (A), receptors (R) and ligand-bond receptors (LR) in different phases within a full HF growth cycle.

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

ActivatorInhibitor
Region IRegion IIRegion IRegion II
(Bulge)(sHG, Matrix, DP)(Bulge)(sHG, Matrix, DP)
Late
Competent
Telogen
(late C)
L3(3,0,4)2(2,0,3)
A2(2,0,1)0(0,0,3)
R2(2,0,4)3(3,0,3)
LR4(4,0,4)1(1,0,3)
Propagating
Anagen
(P)
L2(0,4,4)2(0,3,3)
A2(0,2,3)0(0,2,2)
R2(0,3,4)3(0,3,3)
LR4(0,4,4)1(0,3,3)
Autonomous
Anagen
(A)
L2(0,4,4)2(0,4,4)
A2(0,4,4)0(0,2,2)
R2(0,3,4)3(0,3,3)
LR2(0,4,4)3(0,3,3)
Catagen
(Cat)
L2(2,0,2)2(2,0,4)
A2(2,0,4)0(0,0,1)
R2(2,0,4)3(3,0,3)
LR1(1,0,1)3(3,0,3)
Refractory
Telogen
(R)
L2(2,0,1)2(2,0,4)
A2(2,0,4)0(0,0,1)
R2(2,0,4)3(3,0,3)
LR1(1,0,1)4(4,0,4)
Early
Competent
Telogen
(early C)
L
Similar to refractory telogen


Similar to late competent
telogen
A
R
LR

Appendix 2—table 5.

Parameter values of ligand productions (Appendix 2-equations (9, 10)) and HF growth (Appendix 2-equations (11, 12)).

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

Ligand productionHF growth
Parameters1D2D3DParameters1D2D3D
αA+,162040H+0.2
αA+,0300m+4
αA-1.51.53H-2
αI0−7.8−13−32m-2
αI161020

Supplementary appendix 2 - video captions

Appendix 2—video 1.

3D simulation of dorsal HF waves, where 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 symmetry. Also see Appendix 2—figure 12.

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

Appendix 2—video 2.

3D simulations of dorsal HF waves, where HFs near the head have higher built-in levels of activator signaling.

Here we see highly successive steady waves characterized as the ‘target-like wave’ pattern, also against experimental observations. Also see Appendix 2—figure 13B.

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

Appendix 2—video 3.

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 cycle. However, this asynchrony fails to degrade quickly, and persists during multiple consecutive cycles. This differs from rapid HF growth pattern evolution observable in mice. Also see Appendix 2—figure 15.

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

Appendix 2—video 4.

3D modeling of dorsal-ventral HF wave interaction.

Initial head-to-tail asynchrony in the dorsal domain was introduced by assigning a group of HFs on one end with a one-time increase in activator ligands. We also block ventral HF growth until t = 14 to simulate the temporal-advancement in dorsal domain. Ventral HFs has faster cycle than dorsal HFs. Simulations reproduce head-to-tail asynchrony of first two cycles, which degrades into spontaneous waves possessing bilateral symmetries in later cycles, similar to what is observed from experiments. Also see Figure 2F.

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

Appendix 2—video 5.

Regions like eyelids, physical breaks in the skin can break homogenous wave spreading.

The black region is ‘cut-off’ from the homogeneous outside region, i.e., no HF growth is allowed in this region and no molecules are allowed to diffuse into it. As the HF wave propagates, the black region breaks the wave spreading, creating a distortion-like effect. Also see Figure 7A.

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

Funding Statement

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Funding Information

This paper was supported by the following grants:

  • http://dx.doi.org/10.13039/501100003725National Research Foundation of Korea 2016R1C1B1015211 to Ji Won Oh.
  • http://dx.doi.org/10.13039/501100003725National Research Foundation of Korea 2014R1A5A2010008 to Ji Won Oh.
  • http://dx.doi.org/10.13039/100000054National Cancer Institute T32-CA009054 to Hye-Lim Lee.
  • http://dx.doi.org/10.13039/100000001National Science Foundation DGE-1321846 to Christian Fernando Guerrero-Juarez.
  • http://dx.doi.org/10.13039/100000057National Institute of General Medical Sciences GM055246 to Christian Fernando Guerrero-Juarez.
  • http://dx.doi.org/10.13039/100000069National Institute of Arthritis and Musculoskeletal and Skin Diseases R01-AR061552 to Krzysztof Kobielak.
  • http://dx.doi.org/10.13039/100000069National Institute of Arthritis and Musculoskeletal and Skin Diseases R01-AR056439 to Bogi Andersen.
  • http://dx.doi.org/10.13039/501100000925National Health and Medical Research Council 1023371 to Kiarash Khosrotehrani.
  • http://dx.doi.org/10.13039/100000001National Science Foundation DMS 1161621 to Qing Nie.
  • http://dx.doi.org/10.13039/100000001National Science Foundation DMS 1562176 to Qing Nie.
  • http://dx.doi.org/10.13039/100000057National Institute of General Medical Sciences P50-GM076516 to Qing Nie.
  • http://dx.doi.org/10.13039/100000057National Institute of General Medical Sciences R01-GM107264 to Qing Nie.
  • http://dx.doi.org/10.13039/100000065National Institute of Neurological Disorders and Stroke R01-NS095355 to Qing Nie.
  • http://dx.doi.org/10.13039/100000069National Institute of Arthritis and Musculoskeletal and Skin Diseases R01-AR067273 to Maksim V Plikus.
  • http://dx.doi.org/10.13039/100000875Pew Charitable Trusts 00029641 to Maksim V Plikus.
  • NIH NIAMS R01-AR069653 to Maksim V Plikus.

Additional information

Competing interests

The authors declare that no competing interests exist.

Author contributions

QW, Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Methodology, Writing—original draft, Writing—review and editing.

JWO, Conceptualization, Data curation, Validation, Investigation, Methodology, Writing—original draft, Writing—review and editing.

H-LL, Investigation, Methodology.

AD, Investigation, Methodology.

TP, Data curation, Software, Validation, Methodology.

RR, Investigation, Methodology.

CFG-J, Investigation, Methodology.

XW, Investigation, Methodology.

RZ, Investigation, Methodology.

XC, Investigation, Methodology.

JL, Investigation, Methodology.

MAF, Investigation.

SCJ, Investigation.

ARR, Investigation.

BV, Investigation.

KP, Visualization.

XW, Investigation, Methodology.

NMM, Investigation.

JMP, Investigation.

J-HC, Investigation.

HL, Investigation.

JMDL, Conceptualization, Investigation, Methodology.

EK, Investigation.

JCK, Writing—original draft.

MK, Writing—original draft.

JF, Writing—original draft.

ZY, Writing—original draft, Writing—review and editing.

KKo, Writing—original draft.

BA, Writing—original draft, Writing—review and editing.

KKh, Conceptualization, Investigation, Methodology, Writing—original draft, Project administration, Writing—review and editing.

QN, Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Writing—original draft, Project administration, Writing—review and editing.

MVP, Conceptualization, Supervision, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing.

Ethics

Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#2012-3054 and #2013-3081) of the University of California, Irvine.

Additional files

10.7554/eLife.22772.010

Supplementary file 1.

Dataset 1: Putative activator genes (tabs #1, #2) and putative inhibitor genes (tabs #3, #4) available from a whole skin microarray dataset.

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

10.7554/eLife.22772.011

Supplementary file 2.

Dataset 2: Putative activator and inhibitor genes displaying domain-specific expression patterns at all hair cycle time points on whole skin RNA-seq profiling.

Genes are grouped into individual tabs based on (i) their activator or inhibitor expression profile and (ii) their specificity to one or several skin domains.

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

10.7554/eLife.22772.012

Supplementary file 3.

Dataset 3: Differentially expressed genes specific to refractory telogen dorsal skin, telogen ear skin and cartilage/muscle complex.

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

Major datasets

The following dataset was generated:

Plikus MV,Nie Q,Oh JW,Wang Q,Lee H,Peng T,2016,A multi-scale model for hair follicle reveals heterogeneous skin domains drive rapid spatiotemporal hair growth patterning,https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE85039,Publicly available at the NCBI Gene Expression Omnibus (accession no: GSE85039)

The following previously published dataset was used:

Lin KK,Kumar V,Geyfman M,Chudova D,Ihler AT,Smyth P,Paus R,Takahashi JS,Andersen B,2009,Expression profiling of mouse dorsal skin during hair follicle cycling,https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11186,Publicly available at the NCBI Gene Expression Omnibus (accession no: GSE11186)

References

  • Al-Nuaimi Y, Baier G, Watson RE, Chuong CM, Paus R. The cycling hair follicle as an ideal systems biology research model. Experimental Dermatology. 2010;19:707–713. doi: 10.1111/j.1600-0625.2010.01114.x. [PMC free article] [PubMed] [Cross Ref]
  • Al-Nuaimi Y, Goodfellow M, Paus R, Baier G. A prototypic mathematical model of the human hair cycle. Journal of Theoretical Biology. 2012;310:143–159. doi: 10.1016/j.jtbi.2012.05.027. [PubMed] [Cross Ref]
  • Al-Nuaimi Y, Hardman JA, Bíró T, Haslam IS, Philpott MP, Tóth BI, Farjo N, Farjo B, Baier G, Watson RE, Grimaldi B, Kloepper JE, Paus R. A meeting of two chronobiological systems: circadian proteins Period1 and BMAL1 modulate the human hair cycle clock. The Journal of Investigative Dermatology. 2014;134:610–619. doi: 10.1038/jid.2013.366. [PubMed] [Cross Ref]
  • Bernard BA. The human hair follicle. A Bistable Organ? Experimental Dermatology. 2012;21:401–403. doi: 10.1111/j.1600-0625.2012.01457.x. [PubMed] [Cross Ref]
  • Botchkarev VA, Botchkareva NV, Nakamura M, Huber O, Funa K, Lauster R, Paus R, Gilchrest BA. Noggin is required for induction of the hair follicle growth phase in postnatal skin. FASEB Journal : Official Publication of the Federation of American Societies for Experimental Biology. 2001a;15:2205–2214. doi: 10.1096/fj.01-0207com. [PubMed] [Cross Ref]
  • Botchkarev VA, Komarova EA, Siebenhaar F, Botchkareva NV, Sharov AA, Komarov PG, Maurer M, Gudkov AV, Gilchrest BA. p53 involvement in the control of murine hair follicle regression. The American Journal of Pathology. 2001b;158:1913–1919. doi: 10.1016/S0002-9440(10)64659-7. [PubMed] [Cross Ref]
  • Botchkarev VA, Sharov AA. BMP signaling in the control of skin development and hair follicle growth. Differentiation. 2004;72:512–526. doi: 10.1111/j.1432-0436.2004.07209005.x. [PubMed] [Cross Ref]
  • Brugger SM, Merrill AE, Torres-Vazquez J, Wu N, Ting MC, Cho JY, Dobias SL, Yi SE, Lyons K, Bell JR, Arora K, Warrior R, Maxson R. A phylogenetically conserved cis-regulatory module in the Msx2 promoter is sufficient for BMP-dependent transcription in murine and Drosophila embryos. Development. 2004;131:5153–5165. doi: 10.1242/dev.01390. [PubMed] [Cross Ref]
  • Candille SI, Van Raamsdonk CD, Chen C, Kuijper S, Chen-Tsai Y, Russ A, Meijlink F, Barsh GS. Dorsoventral patterning of the mouse coat by Tbx15. PLoS Biology. 2004;2:E3 doi: 10.1371/journal.pbio.0020003. [PMC free article] [PubMed] [Cross Ref]
  • Castellana D, Paus R, Perez-Moreno M. Macrophages contribute to the cyclic activation of adult hair follicle stem cells. PLoS Biology. 2014;12:e1002002 doi: 10.1371/journal.pbio.1002002. [PMC free article] [PubMed] [Cross Ref]
  • Chen CC, Murray PJ, Jiang TX, Plikus MV, Chang YT, Lee OK, Widelitz RB, Chuong CM. Regenerative hair waves in aging mice and extra-follicular modulators follistatin, dkk1, and sfrp4. The Journal of Investigative Dermatology. 2014;134:2086–2096. doi: 10.1038/jid.2014.139. [PMC free article] [PubMed] [Cross Ref]
  • Chen CC, Wang L, Plikus MV, Jiang TX, Murray PJ, Ramos R, Guerrero-Juarez CF, Hughes MW, Lee OK, Shi S, Widelitz RB, Lander AD, Chuong CM. Organ-level quorum sensing directs regeneration in hair stem cell populations. Cell. 2015;161:277–290. doi: 10.1016/j.cell.2015.02.016. [PMC free article] [PubMed] [Cross Ref]
  • Cheng CW, Niu B, Warren M, Pevny LH, Lovell-Badge R, Hwa T, Cheah KS. Predicting the spatiotemporal dynamics of hair follicle patterns in the developing mouse. PNAS. 2014;111:2596–2601. doi: 10.1073/pnas.1313083111. [PubMed] [Cross Ref]
  • Choi YS, Zhang Y, Xu M, Yang Y, Ito M, Peng T, Cui Z, Nagy A, Hadjantonakis AK, Lang RA, Cotsarelis G, Andl T, Morrisey EE, Millar SE. Distinct functions for wnt/β-catenin in hair follicle stem cell proliferation and survival and interfollicular epidermal homeostasis. Cell Stem Cell. 2013;13:720–733. doi: 10.1016/j.stem.2013.10.003. [PMC free article] [PubMed] [Cross Ref]
  • Chu EY, Hens J, Andl T, Kairo A, Yamaguchi TP, Brisken C, Glick A, Wysolmerski JJ, Millar SE. Canonical WNT signaling promotes mammary placode development and is essential for initiation of mammary gland morphogenesis. Development. 2004;131:4819–4829. doi: 10.1242/dev.01347. [PubMed] [Cross Ref]
  • Clavel C, Grisanti L, Zemla R, Rezza A, Barros R, Sennett R, Mazloom AR, Chung CY, Cai X, Cai CL, Pevny L, Nicolis S, Ma'ayan A, Rendl M. Sox2 in the dermal papilla niche controls hair growth by fine-tuning BMP signaling in differentiating hair shaft progenitors. Developmental Cell. 2012;23:981–994. doi: 10.1016/j.devcel.2012.10.013. [PMC free article] [PubMed] [Cross Ref]
  • Cotsarelis G, Sun T-T, Lavker RM. Label-retaining cells reside in the bulge area of pilosebaceous unit: implications for follicular stem cells, hair cycle, and skin carcinogenesis. Cell. 1990;61:1329–1337. doi: 10.1016/0092-8674(90)90696-C. [PubMed] [Cross Ref]
  • Enshell-Seijffers D, Lindon C, Kashiwagi M, Morgan BA. beta-catenin activity in the dermal papilla regulates morphogenesis and regeneration of hair. Developmental Cell. 2010;18:633–642. doi: 10.1016/j.devcel.2010.01.016. [PMC free article] [PubMed] [Cross Ref]
  • Fessing MY, Sharova TY, Sharov AA, Atoyan R, Botchkarev VA. Involvement of the Edar signaling in the control of hair follicle involution (catagen) The American Journal of Pathology. 2006;169:2075–2084. doi: 10.2353/ajpath.2006.060227. [PubMed] [Cross Ref]
  • Festa E, Fretz J, Berry R, Schmidt B, Rodeheffer M, Horowitz M, Horsley V. Adipocyte lineage cells contribute to the skin stem cell niche to drive hair cycling. Cell. 2011;146:761–771. doi: 10.1016/j.cell.2011.07.019. [PMC free article] [PubMed] [Cross Ref]
  • Foitzik K, Lindner G, Mueller-Roever S, Maurer M, Botchkareva N, Botchkarev V, Handjiski B, Metz M, Hibino T, Soma T, Dotto GP, Paus R. Control of murine hair follicle regression (catagen) by TGF-beta1 in vivo. FASEB Journal : Official Publication of the Federation of American Societies for Experimental Biology. 2000;14:752–760. [PubMed]
  • Greco V, Chen T, Rendl M, Schober M, Pasolli HA, Stokes N, Dela Cruz-Racelis J, Fuchs E. A two-step mechanism for stem cell activation during hair regeneration. Cell Stem Cell. 2009;4:155–169. doi: 10.1016/j.stem.2008.12.009. [PMC free article] [PubMed] [Cross Ref]
  • Guha U, Gomes WA, Samanta J, Gupta M, Rice FL, Kessler JA. Target-derived BMP signaling limits sensory neuron number and the extent of peripheral innervation in vivo. Development. 2004;131:1175–1186. doi: 10.1242/dev.01013. [PubMed] [Cross Ref]
  • Halloy J, Bernard BA, Loussouarn G, Goldbeter A. Modeling the dynamics of human hair cycles by a follicular automaton. PNAS. 2000;97:8328–8333. doi: 10.1073/pnas.97.15.8328. [PubMed] [Cross Ref]
  • Higgins CA, Petukhova L, Harel S, Ho YY, Drill E, Shapiro L, Wajid M, Christiano AM. FGF5 is a crucial regulator of hair length in humans. PNAS. 2014;111:10648–10653. doi: 10.1073/pnas.1402862111. [PubMed] [Cross Ref]
  • Hodgson SS, Neufeld Z, Villani RM, Roy E, Khosrotehrani K. Transgenic flash mice for in vivo quantitative monitoring of canonical wnt signaling to track hair follicle cycle dynamics. The Journal of Investigative Dermatology. 2014;134:1519–1526. doi: 10.1038/jid.2014.92. [PubMed] [Cross Ref]
  • Hsu YC, Pasolli HA, Fuchs E. Dynamics between stem cells, niche, and progeny in the hair follicle. Cell. 2011;144:92–105. doi: 10.1016/j.cell.2010.11.049. [PMC free article] [PubMed] [Cross Ref]
  • Ito C, Saitoh Y, Fujita Y, Yamazaki Y, Imamura T, Oka S, Suzuki S. Decapeptide with fibroblast growth factor (FGF)-5 partial sequence inhibits hair growth suppressing activity of FGF-5. Journal of Cellular Physiology. 2003;197:272–283. doi: 10.1002/jcp.10369. [PubMed] [Cross Ref]
  • Javier AL, Doan LT, Luong M, Reyes de Mochel NS, Sun A, Monuki ES, Cho KW. Bmp Indicator mice reveal dynamic regulation of transcriptional response. PLoS One. 2012;7:e42566 doi: 10.1371/journal.pone.0042566. [PMC free article] [PubMed] [Cross Ref]
  • Kandyba E, Leung Y, Chen YB, Widelitz R, Chuong CM, Kobielak K. Competitive balance of intrabulge BMP/Wnt signaling reveals a robust gene network ruling stem cell homeostasis and cyclic activation. PNAS. 2013;110:1351–1356. doi: 10.1073/pnas.1121312110. [PubMed] [Cross Ref]
  • Kandyba E, Kobielak K. Wnt7b is an important intrinsic regulator of hair follicle stem cell homeostasis and hair follicle cycling. Stem Cells. 2014;32:886–901. doi: 10.1002/stem.1599. [PMC free article] [PubMed] [Cross Ref]
  • Kimura-Ueki M, Oda Y, Oki J, Komi-Kuramochi A, Honda E, Asada M, Suzuki M, Imamura T. Hair cycle resting phase is regulated by cyclic epithelial FGF18 signaling. The Journal of Investigative Dermatology. 2012;132:1338–1345. doi: 10.1038/jid.2011.490. [PubMed] [Cross Ref]
  • Kobielak K, Pasolli HA, Alonso L, Polak L, Fuchs E. Defining BMP functions in the hair follicle by conditional ablation of BMP receptor IA. The Journal of Cell Biology. 2003;163:609–623. doi: 10.1083/jcb.200309042. [PMC free article] [PubMed] [Cross Ref]
  • Kobielak K, Stokes N, de la Cruz J, Polak L, Fuchs E. Loss of a quiescent niche but not follicle stem cells in the absence of bone morphogenetic protein signaling. PNAS. 2007;104:10063–10068. doi: 10.1073/pnas.0703004104. [PubMed] [Cross Ref]
  • Legrand JM, Roy E, Ellis JJ, Francois M, Brooks AJ, Khosrotehrani K. STAT5 activation in the dermal papilla is Important for Hair Follicle Growth Phase Induction. Journal of Investigative Dermatology. 2016;136:1781–1791. doi: 10.1016/j.jid.2016.04.014. [PubMed] [Cross Ref]
  • Leishman E, Howard JM, Garcia GE, Miao Q, Ku AT, Dekker JD, Tucker H, Nguyen H. Foxp1 maintains hair follicle stem cell quiescence through regulation of Fgf18. Development. 2013;140:3809–3818. doi: 10.1242/dev.097477. [PubMed] [Cross Ref]
  • Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323 doi: 10.1186/1471-2105-12-323. [PMC free article] [PubMed] [Cross Ref]
  • Lin KK, Kumar V, Geyfman M, Chudova D, Ihler AT, Smyth P, Paus R, Takahashi JS, Andersen B. Circadian clock genes contribute to the regulation of hair follicle cycling. PLoS Genetics. 2009;5:e1000573 doi: 10.1371/journal.pgen.1000573. [PMC free article] [PubMed] [Cross Ref]
  • Lindner G, Botchkarev VA, Botchkareva NV, Ling G, van der Veen C, Paus R. Analysis of apoptosis during hair follicle regression (catagen) The American Journal of Pathology. 1997;151:1601–1617. [PubMed]
  • Liu S, Leask A. CCN2 modulates hair follicle cycling in mice. Molecular Biology of the Cell. 2013;24:3939–3944. doi: 10.1091/mbc.E13-08-0472. [PMC free article] [PubMed] [Cross Ref]
  • Lustig B, Jerchow B, Sachs M, Weiler S, Pietsch T, Karsten U, van de Wetering M, Clevers H, Schlag PM, Birchmeier W, Behrens J. Negative feedback Loop of wnt signaling through upregulation of Conductin/Axin2 in colorectal and liver tumors. Molecular and Cellular Biology. 2002;22:1184–1193. doi: 10.1128/MCB.22.4.1184-1193.2002. [PMC free article] [PubMed] [Cross Ref]
  • Matsumura H, Mohri Y, Binh NT, Morinaga H, Fukuda M, Ito M, Kurata S, Hoeijmakers J, Nishimura EK. Hair follicle aging is driven by transepidermal elimination of stem cells via COL17A1 proteolysis. Science. 2016;351:aad4395 doi: 10.1126/science.aad4395. [PubMed] [Cross Ref]
  • Maurer M, Handjiski B, Paus R. Hair growth modulation by topical immunophilin ligands: induction of anagen, inhibition of massive catagen development, and relative protection from chemotherapy-induced alopecia. The American Journal of Pathology. 1997;150:1433–1441. [PubMed]
  • Mesa KR, Rompolas P, Zito G, Myung P, Sun TY, Brown S, Gonzalez DG, Blagoev KB, Haberman AM, Greco V. Niche-induced cell death and epithelial phagocytosis regulate hair follicle stem cell pool. Nature. 2015;522:94–97. doi: 10.1038/nature14306. [PMC free article] [PubMed] [Cross Ref]
  • Murray PJ, Maini PK, Plikus MV, Chuong CM, Baker RE. Modelling hair follicle growth dynamics as an excitable medium. PLoS Computational Biology. 2012;8:e1002804 doi: 10.1371/journal.pcbi.1002804. [PMC free article] [PubMed] [Cross Ref]
  • Oshimori N, Fuchs E. Paracrine TGF-β signaling counterbalances BMP-mediated repression in hair follicle stem cell activation. Cell Stem Cell. 2012;10:63–75. doi: 10.1016/j.stem.2011.11.005. [PMC free article] [PubMed] [Cross Ref]
  • Paladini RD, Saleh J, Qian C, Xu GX, Rubin LL. Modulation of hair growth with small molecule agonists of the hedgehog signaling pathway. The Journal of Investigative Dermatology. 2005;125:638–646. doi: 10.1111/j.0022-202X.2005.23867.x. [PubMed] [Cross Ref]
  • Paus R, Stenn KS, Link RE. The induction of anagen hair growth in telogen mouse skin by cyclosporine A administration. Laboratory Investigation; a Journal of Technical Methods and Pathology. 1989;60:365–369. [PubMed]
  • Paus R, Müller-Röver S, Botchkarev VA. Chronobiology of the hair follicle: hunting the “hair cycle clock” Journal of Investigative Dermatology Symposium Proceedings. 1999;4:338–345. doi: 10.1038/sj.jidsp.5640241. [PubMed] [Cross Ref]
  • Paus R, Foitzik K. In search of the “hair cycle clock”: a guided tour. Differentiation. 2004;72:489–511. doi: 10.1111/j.1432-0436.2004.07209004.x. [PubMed] [Cross Ref]
  • Picelli S, Faridani OR, Björklund AK, Winberg G, Sagasser S, Sandberg R. Full-length RNA-seq from single cells using Smart-seq2. Nature Protocols. 2014;9:171–181. doi: 10.1038/nprot.2014.006. [PubMed] [Cross Ref]
  • Plikus M, Wang WP, Liu J, Wang X, Jiang TX, Chuong CM. Morpho-regulation of ectodermal organs: integument pathology and phenotypic variations in K14-Noggin engineered mice through modulation of bone morphogenic protein pathway. The American Journal of Pathology. 2004;164:1099–1114. doi: 10.1016/S0002-9440(10)63197-5. [PubMed] [Cross Ref]
  • Plikus MV, Zeichner-David M, Mayer JA, Reyna J, Bringas P, Thewissen JG, Snead ML, Chai Y, Chuong CM. Morphoregulation of teeth: modulating the number, size, shape and differentiation by tuning bmp activity. Evolution & Development. 2005;7:440–457. doi: 10.1111/j.1525-142X.2005.05048.x. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV, Chuong CM. Complex hair cycle domain patterns and regenerative hair waves in living rodents. The Journal of Investigative Dermatology. 2008a;128:1071–1080. doi: 10.1038/sj.jid.5701180. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV, Mayer JA, de la Cruz D, Baker RE, Maini PK, Maxson R, Chuong CM. Cyclic dermal BMP signalling regulates stem cell activation during hair regeneration. Nature. 2008b;451:340–344. doi: 10.1038/nature06457. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV, Widelitz RB, Maxson R, Chuong CM. Analyses of regenerative wave patterns in adult hair follicle populations reveal macro-environmental regulation of stem cell activity. The International Journal of Developmental Biology. 2009;53:857–868. doi: 10.1387/ijdb.072564mp. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV, Baker RE, Chen CC, Fare C, de la Cruz D, Andl T, Maini PK, Millar SE, Widelitz R, Chuong CM. Self-organizing and stochastic behaviors during the regeneration of hair stem cells. Science. 2011;332:586–589. doi: 10.1126/science.1201647. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV. New activators and inhibitors in the hair cycle clock: targeting stem cells' state of competence. The Journal of Investigative Dermatology. 2012;132:1321–1324. doi: 10.1038/jid.2012.38. [PMC free article] [PubMed] [Cross Ref]
  • Plikus MV, Chuong CM. Macroenvironmental regulation of hair cycling and collective regenerative behavior. Cold Spring Harbor Perspectives in Medicine. 2014;4:a015198 doi: 10.1101/cshperspect.a015198. [PMC free article] [PubMed] [Cross Ref]
  • Rivera-Gonzalez GC, Shook BA, Andrae J, Holtrup B, Bollag K, Betsholtz C, Rodeheffer MS, Horsley V. Skin adipocyte stem cell self-renewal is regulated by a PDGFA/AKT-signaling axis. Cell Stem Cell. 2016;19:738–751. doi: 10.1016/j.stem.2016.09.002. [PMC free article] [PubMed] [Cross Ref]
  • Schneider MR, Schmidt-Ullrich R, Paus R. The hair follicle as a dynamic miniorgan. Current Biology : CB. 2009;19:R132–142. doi: 10.1016/j.cub.2008.12.005. [PubMed] [Cross Ref]
  • Sharov AA, Fessing M, Atoyan R, Sharova TY, Haskell-Luevano C, Weiner L, Funa K, Brissette JL, Gilchrest BA, Botchkarev VA. Bone morphogenetic protein (BMP) signaling controls hair pigmentation by means of cross-talk with the melanocortin receptor-1 pathway. PNAS. 2005;102:93–98. doi: 10.1073/pnas.0408455102. [PubMed] [Cross Ref]
  • Sharov AA, Sharova TY, Mardaryev AN, Tommasi di Vignano A, Atoyan R, Weiner L, Yang S, Brissette JL, Dotto GP, Botchkarev VA. Bone morphogenetic protein signaling regulates the size of hair follicles and modulates the expression of cell cycle-associated genes. PNAS. 2006;103:18166–18171. doi: 10.1073/pnas.0608899103. [PubMed] [Cross Ref]
  • Sick S, Reinker S, Timmer J, Schlake T. WNT and DKK determine hair follicle spacing through a reaction-diffusion mechanism. Science. 2006;314:1447–1450. doi: 10.1126/science.1130088. [PubMed] [Cross Ref]
  • Stenn KS, Paus R. Controls of hair follicle cycling. Physiological Reviews. 2001;81:449–494. [PubMed]
2017; 6: e22772.

Decision letter

Valerie Horsley, Reviewing editor
Valerie Horsley, Yale University, United States;

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "‘A multi-scale model for hair follicles reveals heterogeneous domains driving rapid spatiotemporal hair growth patterning"’ for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors and the evaluation has been overseen by Fiona Watt as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We found that your manuscript describes an interesting phenomenon that mouse skin is composed of different anatomic regions that interact as heterogeneous excitable domains. The authors used thorough analysis of the whole body hair cycle dynamics coupled with novel mathematical modeling to explain the emergence of certain hair cycle patterns in mouse skin. Extensive mRNA profiling analysis offers some hints as to the molecular mechanisms behind some of the observed features. This work opens up a new angle to analyze and understand the nature of hair growth on the whole body scale, an important perspective not appreciated before.

Major comments to be addressed in a revised manuscript:

1) In general, the reviewers thought that these results should be placed in a broader context with the authors' own work and that of others. First, what is the difference between this manuscript, Plikus et al., 2011 and the interaction between follicles that was recently explored by Chen et al., 2015? The authors make some comparisons in the Results section (subsection “‘Model reveals skin is a heterogeneous excitable medium.”’) but this type of logic needs to come into the experimental design for the introduction and results. Second, the focus on Wnt and Bmp signaling seems redundant based on previous work from several labs. Third, the authors should discuss how a similar approach could benefit the study of other tissues?

2) Often the text is written in a confusing manner. For instance, "‘we scaled the model both into two (Appendix 2—figure 4A) and three dimensions (Appendix 2—figure 4B)"’ does not make sense to me. This lack of clarity makes Figure 1 and Figure 2 almost impossible to interpret.

3) Figure 2A and B do not clearly show "‘neck-to-tail asynchrony"’ in hair morphogenesis. Do the authors mean that the center is different from the left to right regions?

In Figure 2F the authors demonstrated a scenario of how bilaterally symmetric hair growth pattern could emerge. Even with the assumptive parameters they described in the text, it is not clear why the spread of anagen from dorsal-ventral only occurs in the posterior part of the ventral skin, such as from t=50 to t=62. Is it assumed that the anterior- posterior parts of the ventral skin are also different? Whether or not this is the case here, the authors should expand to explain this phenomenon because it affects the outcome of the entire model. In this manuscript the authors didn't address the anterior- posterior difference of the skin, but given the nature of this study, it would be helpful to clarify why the authors choose to do so.

4) Figure 3 uses a Wnt reporter to examine hair cycling but it is unclear whether the luciferase signal is skin specific.

The result in Appendix 1 is highly interesting. It suggests the faster cycling nature of the chin skin is intrinsic, instead of relying on neighboring tissues in the chin like the cartilage in the ear. To drive the point home, it would help to do the Flash analysis on the chin graft and compare the dynamics of Wnt signaling in graft compared to native chin skin. Also this result should be moved to main figure because it serves as a strong justification to do subsequent RNA-seq analysis.

5) The implication that the ventral skin initiates hair cycling is interesting and not fully tested experimentally in this manuscript.

6) The hair follicle is certainly a good model to begin to examine these problems and the combination of experiment and modeling is very powerful. Could it apply to other large spread out tissues such as muscle, blood, bones, etc? Maybe this can be discussed. It is very interesting that the model predicts a sort of self-sufficiency of the regenerative process, in which, like a perpetuum mobile, it can keep going potentially without "‘higher-order"’ control or even trigger from a systemic signal. Maybe this aspect can be emphasized.

7) Given the amount of work the authors put in Figure 4, the message is surprisingly little, other than the description of very complicated and a little confusing dynamic expression patterns. The authors do have multiple genetic models used in this study that can help them sort out the functions of some of the signal pathways in defining the different anatomic domains. For instance in their K14-Noggin, K14-BMP4, K14Cre-Wnt7b fl/fl, K14-Wnt7a, K5trTA-teto-Dkk1 models, do the chin skin cycling frequencies change, do the ventral skin anagen duration change? The authors could choose at least 1 representative model of Wnt and BMP signal pathways to do this.

8) In Figure 6H, the percentages of the ear follicles going into anagen in these transgenic models needs to be quantified. It will help address how dominant these signal pathways are in maintaining the quiescence of the ear skin.

Major comments to be addressed in a revised manuscript:

1) In general, the reviewers thought that these results should be placed in a broader context with the authors' own work and that of others. First, what is the difference between this manuscript, Plikus et al., 2011 and the interaction between follicles that was recently explored by Chen et al., 2015? The authors make some comparisons in the Results section (subsection “‘Model reveals skin is a heterogeneous excitable medium.”’) but this type of logic needs to come into the experimental design for the introduction and results. Second, the focus on Wnt and Bmp signaling seems redundant based on previous work from several labs. Third, the authors should discuss how a similar approach could benefit the study of other tissues?

We thank the reviewers for their suggestions. In this submission, we made changes in the Introduction and Discussion to clarify that this work is the first to examine hair regeneration across most of the mouse’s body skin, and its findings represent a significant advance over previous knowledge on group-level hair growth behavior in the dorsal skin. From a modeling perspective, this is the first multi-scale hair follicle model that incorporates both molecular dynamics and physical growth into one model consisting of three dimensions in space, and that faithfully predicts cyclic behaviors at the level of both a single hair follicle and of hair populations collectively. Through both experiments and modeling, our work revealed that:

a) The physiological pace of hair regeneration differs among distinct anatomical sites. It is prominently faster in the chin and ventral skin relative to dorsal skin, and very slow in the ear skin. Thus, the skin as a whole appears to function as a complex regenerative landscape with regions of fast, slow, and very slow hair renewal. Our new data in this revision (Appendix 1—figure 23) shows that this behavior produces fur coat with variable hair density, which likely serves an adaptive role. Prior works on dorsal skin did not consider that hair regeneration dynamics prominently differ between body sites. We now highlight this in the Introduction and in the Discussion.

b) WNT/BMP activator/inhibitor signaling pair modulates hair regeneration in all skin regions. We believe this finding is important and non-redundant with prior works. Our current work shows that WNT/BMP signaling is a universal “‘molecular language”’ for hair growth on the entire body, rather than only on a specific body site, such as the dorsal skin shown in previous works.Indeed, the choice of dorsal skin as a model for studying hair growth in most prior works is primarily due to convenience. By studying hair growth in other, “‘less convenient”’ body sites we observe that the mechanism driven by WNT/BMP signaling pair is general. In this revision we took additional care in further strengthening this important aspect of our work through providing new data on skin-wide WNT and BMP reporter activity (Figure 4G and 4I; Appendix 1—figure21), as well as new data on the changes in hair growth patterns in WNT and BMP mutant mice (Figure 4J and 4K; Appendix 1—figure22).

c) Several additional conclusions stem directly from the findings above. First, since all skin regions share a WNT/BMP “‘language”’, hair-to-hair communications naturally arise across the anatomic boundaries. An alternative to this mechanism would be distinct, site-specific signals and full autonomy of regional hair growth. Because WNT/BMP-based communication is possible across the boundaries, it enables novel hair growth behavior not obvious from prior works – fast cycling skin regions (such as chin skin) function as a kind of pacemaker for skin-wide hair growth. We emphasize this in the Introduction and in the Discussion.

d) Second, we found that non-skin tissues with prominent WNT/BMP expression profiles can act as a signaling macro-environment for hair follicles and strongly influence their regeneration pace. In this study we show that in the mouse ear, WNT-low/BMP-high cartilage/muscle complex slows hair growth in the ear skin to a virtual standstill. This finding expands the repertoire of tissues with a signaling macro-environment function to include closely positioned anatomic structures beyond typical skin cell populations. We highlight this in the Discussion.

e) We posit that some of the newly found hair regeneration features can have analogs in other tissues and organs. For instance, dominant anatomically defined pacemakers are common in the electrically coupled muscle-based tissues. Prominently, cardiac tissue contains pacemaker nodes that generate dominant electric rhythms to drive contractile activity of the entire heart. A similar pacemaker exists in the stomach (gastric pacemaker zone), which controls its coordinated peristalsis. Other actively regenerating organs, such as the gastro-intestinal tract and bone marrow likely contain anatomic regions of faster and slower regeneration and, conceivably, they can be coupled to work in coordination. Knowledge learned from the skin system can guide the search for regenerative landscapes in these and other organs. Because coordination principles observed in the skin may be universal for large-scale patterning systems, the likelihood of them operating in other organs is substantial despite prominent anatomical differences. We highlight this in the Discussion.

2) Often the text is written in a confusing manner. For instance, "‘we scaled the model both into two (Appendix 2—figure 4A) and three dimensions (Appendix 2—figure 4B)"’ does not make sense to me. This lack of clarity makes Figure 1 and Figure 2 almost impossible to interpret.

We apologize for the confusion. In this revision we included an expanded description of the modeling, both in the main text (Results section) and in the legends for the main Figures 1 and and22 and Appendix 2—figure4. We hope that this helps to improve the clarity of the modeling aspects of the study.

3) Figure 2A and B do not clearly show "‘neck-to-tail asynchrony"’ in hair morphogenesis. Do the authors mean that the center is different from the left to right regions?

Indeed, as noted by the reviewers, the neck-to-tail asynchrony in dorsal hair morphogenesis at P0-P1 is fairly subtle. Additionally, subtle differences are observed between the sides and the center of the dorsal skin. We now noted this explicitly in the revised text (Results section). Our data on main Figures 3 and and44 shows that these initial spatio-temporal asynchronies in dorsal hair morphogenesis become effectively “‘wiped out”’ by the dominant chin/ventral to dorsal hair growth spreading at the onset of the second anagen. Nonetheless, we pointed out these initial asynchronies.

In Figure 2F the authors demonstrated a scenario of how bilaterally symmetric hair growth pattern could emerge. Even with the assumptive parameters they described in the text, it is not clear why the spread of anagen from dorsal-ventral only occurs in the posterior part of the ventral skin, such as from t=50 to t=62. Is it assumed that the anterior- posterior parts of the ventral skin are also different? Whether or not this is the case here, the authors should expand to explain this phenomenon because it affects the outcome of the entire model. In this manuscript the authors didn't address the anterior- posterior difference of the skin, but given the nature of this study, it would be helpful to clarify why the authors choose to do so.

We agree with the reviewers’ criticism that the modeling result in the original Figure 2F is not fully consistent with the experimentally observed hair growth patterns. In the resubmission, we critically scrutinized our model through comprehensive exploration of parameters, and found that this feature was the result of a particular set of parameter values chosen for our original simulations. In this revision, we tested a wider range of parameter values and found that dorsal/ventral spreading is not an exclusive feature of the posterior skin, but that it can occur along the entire length of the dorsal/ventral boundary. The updated set of parameter values includes: (a) earlier onset of ventral hair follicle morphogenesis, (b) lower signaling threshold for anagen initiation (Δ+=2.6), (c) lower activator/inhibitor signaling noise (λAI=0.02), and (d) smaller time step length (dt=0.01). Updated results are now shown in the main Figure 2F and the parameters used are shown in the Appendix 2(subsection “‘2D and 3D modeling details”’.1 and Appendix 2—tables 1 to 3. Furthermore, for consistency, we re-run other simulations and updated Appendix 2—figures7, ,12,12, 13B, ,1515 to 18 and 28 to 30. In all cases we obtained qualitatively similar results that do not change our original conclusions. Lastly, we also made three additional changes in the presentation of the modeling results in the Appendix 2 file:

1) The effects of the maximum hair follicle growth length hmax on wave propagation behaviors were presented in a supplemental table in the original submission. They are now shown in Appendix 2—figure27.

2) The effects of total available receptors RAand RI on the ventral-dorsal wave interactions were shown in the movies in the original submission. They are now shown in Appendix 2—figure28 to 30.

3) New Appendix 2—table 2 is added to accompany Appendix 2—table 1. It details the parameter values used in in each simulation.

4) Figure 3 uses a Wnt reporter to examine hair cycling but it is unclear whether the luciferase signal is skin specific.

The original study on Flash WNT reporter mice by Hodgson et al., 2014 carefully examined the origin of the luminescence signal and concluded that it is skin specific. In their work, authors state that '[…]Upon dissection, no signal beyond WT level was detected in organs other than the skin at both short and long exposure times, although luciferase expression could be detected at RNA levels, thus indicating that transgene activity at high and detectable level is exclusive to the skin without internal organ signal contamination […]'. We thus interpret Flash reporter signal in the current study as skin-specific. In the revision, we point out the skin-specific nature of Flash signal and cite the above reference.

The result in Appendix 1 is highly interesting. It suggests the faster cycling nature of the chin skin is intrinsic, instead of relying on neighboring tissues in the chin like the cartilage in the ear. To drive the point home, it would help to do the Flash analysis on the chin graft and compare the dynamics of Wnt signaling in graft compared to native chin skin. Also this result should be moved to main figure because it serves as a strong justification to do subsequent RNA-seq analysis.

We agree that grafting results support intrinsically faster underlying dynamics of WNT and/or BMP signaling in the chin skin. As suggested, we did consider Flash skin grafting experiments; however, our Flash mice are on a mixed genetic background; this precludes chin-to-dorsal grafting between individual Flash animals. In addition, due to the highly traumatic nature of an autologous transplantation procedure from chin region, we were unable to obtain IACUC approval to conduct chin-to-dorsal grafting experiments in the same animal.

Instead, we performed the following additional experiments to strengthen the link between domain-specific hair cycle dynamics and BMP/WNT signaling:

a) We performed additional chin and dorsal WT skin grafting experiments to SCID donors to strengthen the findings that chin skin has faster intrinsic hair cycle dynamics. New grafting results are now shown in main Figure 3.

b) We performed in-depth analysis of the domain specific RNA-seq data, with a focus on early third telogen. This new analysis is now shown on main Figure 3C to 3F. Compared to dorsal skin, chin and ventral skin shows prominent expression differences for the BMP pathway (less BMP ligands and more BMP antagonists) and WNT pathway (less WNT antagonists). Because the RNA-seq was performed on carefully dissected skin, our analysis further confirms that the observed BMP and WNT signaling differences are intrinsic to these skin regions.

c) We also examined lacZ expression patterns in Axin2-lacZ WNT reporter and BRE-gal BMP reporter mice during second telogen (main Figure 4G to 4I and Appendix 1—figure21). Consistent with Flash WNT reporter data, we observe faster increase in Axin2-lacZ reporter activity in the chin and ventral skin as compared to dorsal skin. Dynamic lacZ signal localizes to dermal papillae of telogen hair follicles, where it functions as a marker of telogen competence (Plikus et al., 2011). We see many more lacZ+ dermal papillae in the chin and ventral skin during early second telogen, on day P36. For BRE-gal reporter mice, we see fewer lacZ+ hair follicles (bulge region expression) in the chin and ventral skin as compared to dorsal skin (Figure 4I). These new results strengthen our findings that telogen hair follicles acquire low-BMP and high-WNT signaling status faster in the chin and ventral skin compared to dorsal skin.

d) We also compared second telogen skin of chin, ventral and dorsal sites by histology, immuno-histochemistry and FACS. Consistent with the Gene Ontology analysis on the RNA-seq data (main Figure 3E), we observe site-specific differences in the distribution of dermal adipose tissue and panniculus carnosus skeletal muscle (Appendix 1—figure19), as well as macrophages (Appendix 1—figure20). We speculate that these regional differences in skin composition may contribute to regional differences in WNT and BMP activities, which underlie faster chin and ventral hair cycling. Admittedly, an in-depth follow-up study will be necessary to identify and verify the major site-specific cellular sources for WNT and BMP ligands and antagonists. At the same time, our current study establishes clear site-specific hair cycle dynamics and promising regional differences in skin microanatomy, warranting follow-up experiments.

5) The implication that the ventral skin initiates hair cycling is interesting and not fully tested experimentally in this manuscript.

We speculate that generally faster hair cycle in the chin and ventral skin helps to achieve higher fur density, which can aid with increased insulation and protection against mechanical abrasion of the ventrum in low gait animals, such as mice. In this revision, we examined club hair density distribution between the three regions in adult, 4-month old mice. Indeed, we found that chin and ventral telogen hair follicles consistently had more club hairs (three or more) as compared to dorsal telogen hair follicles (mostly two or one club hair). This new data is now shown on the Appendix 1—figure23, and the relevant speculations are added in the Discussion. We also acknowledged in Discussion that the rate of hair regeneration is only one of the factors determining fur density, and that the latter can also be affected by the rate of hair fallout (aka exogen), a parameter not directly tested in this work.

6) The hair follicle is certainly a good model to begin to examine these problems and the combination of experiment and modeling is very powerful. Could it apply to other large spread out tissues such as muscle, blood, bones, etc? Maybe this can be discussed. It is very interesting that the model predicts a sort of self-sufficiency of the regenerative process, in which, like a perpetuum mobile, it can keep going potentially without "‘higher-order"’ control or even trigger from a systemic signal. Maybe this aspect can be emphasized.

We thank the reviewers for this excellent suggestion. In the revised Discussion, we added that skin with its many cycling hair follicles is just one, albeit prominent example, of many organs that implements large-scale coordination of its physiological activities. As mentioned above, electrically coupled heart and stomach display large-scale coordination of their contractile activities from anatomically defined pacemaker regions. Other organs, prominently the intestine share many attributes with skin – analogous to skin, intestines contain many repetitive stem cell-rich units along its villi. It is conceivable that higher-order coordination of regeneration, analogous to that seen in haired skin occurs between intestinal villi.

Related to hair regeneration, several levels of higher-order coordination exist:

a) Coordination within hair populations occurs via self-organization between individual hair follicles (Plikus et al., 2008 and Plikus et al., 2011). This can be further modulated by the local signaling macro-environment, including adipose cells, immune cells, etc.

b) Coordination between distinct anatomical regions, such as between chin, ventral and dorsal skin sites described in the current study.

c) Systemic hormonal and neuronal coordination. Indeed, large amount of literature exists on hair growth coordination by humoral factors, including but not limited to prolactin and sex hormones, and sympathetic nerves.

In a given species, fur renewal can be regulated by several of the above mechanism working either in parallel or hierarchically, when one mechanism overrides another. In the revised Discussion, we provide succinct speculation of these ideas.

7) Given the amount of work the authors put in Figure 4, the message is surprisingly little, other than the description of very complicated and a little confusing dynamic expression patterns. The authors do have multiple genetic models used in this study that can help them sort out the functions of some of the signal pathways in defining the different anatomic domains. For instance in their K14-Noggin, K14-BMP4, K14Cre-Wnt7b fl/fl, K14-Wnt7a, K5trTA-teto-Dkk1 models, do the chin skin cycling frequencies change, do the ventral skin anagen duration change? The authors could choose at least 1 representative model of Wnt and BMP signal pathways to do this.

We agree with this point. For this revision, we moved most of the correlative RNA-seq analysis into the supplement. Furthermore, we performed more in-depth analysis on the early second telogen RNA-seq (main Figure 4C to 4F). This new analysis shows more clearly that chin and ventral skin express lower levels of WNT antagonists and BMP ligands and higher levels of BMP antagonists, consistent with the possibility that these regions are less refractory. Following on this correlation, we carefully examined lacZ staining patterns in Axin2-lacZ WNT reporter and BRE-gal BMP reporter mice (main Figure 4G to 4I and Appendix 1—figure21). This new data supports the notion that both chin and ventral telogen skin acquires low-BMP and high-WNT signaling status (i.e. less refractory status) faster as compared to dorsal skin. Lastly, following reviewers’ suggestion we examined skin-wide hair growth patterns in Krt14-Wnt7a and Krt14-Bmp4 mutant mice (main Figure 4J and4K and Appendix 1—figure22). Consistent with the speculated role for WNT pathway as the skin-wide hair cycle activator we observe loss of ventral-to-dorsal hair cycle initiation dominance and ectopic hair growth initiation sites appear in the dorsal skin of Krt14-Wnt7a mice at the start of the second (Appendix 1—figure22) and third anagen (Figure 4J). In analogy, consistent with the speculated role for BMP pathway as the skin-wide hair cycle inhibitor we observe stalled ventral-to-dorsal hair cycle spreading, formation of premature hair growth boundaries and generally asymmetric, jagged hair growth patterns in Krt14-Bmp4 mice at the start of third anagen (Figure 4K).

8) In Figure 6H, the percentages of the ear follicles going into anagen in these transgenic models needs to be quantified. It will help address how dominant these signal pathways are in maintaining the quiescence of the ear skin.

Phenotype quantification and additional images of individual ear skin samples from WT, Krt14-Noggin and Krt14-Wnt7a mice have been added in the new Appendix 1—figure28. Both mutants show statistically significant increase in the number of spontaneously cycling ear hair follicles.


Articles from eLife are provided here courtesy of eLife Sciences Publications, Ltd