A Discrete-Event Simulation Approach to Identify Rules that Govern Arbor Remodeling for Branching Cutaneous Afferents in Hairy Skin

In mammals, touch is encoded by sensory receptors embedded in the skin.  For one class of receptors in the mouse, the architecture of its Merkel cells, unmyelinated neurites, and heminodes follow particular renewal and remodeling trends over hair cycle stages from ages 4 to 10 weeks.  As it is currently impossible to observe such trends across a single animal’s hair cycle, this work employs discrete event simulation to identify and evaluate policies of Merkel cell and heminode dynamics. Well matching the observed data, the results show that the baseline model replicates dynamic remodeling behaviors between stages of the hair cycle – based on particular addition and removal polices and estimated probabilities tied to constituent parts of Merkel cells, terminal branch neurites and heminodes.  The analysis shows further that certain policies hold greater influence than others.  This use of computation is a novel approach to understanding neuronal development.

Introduction

The sense of touch is key to behaviors of everyday living such as feeding, social bonding and avoiding bodily harm. In mammals, touch is encoded by sensory receptors embedded in the skin (Kandel 2012). Sensory receptors include cutaneous light touch afferents as well as those signaling information regarding proprioception, chemoreception and pain. Both the sensory receptors and the skin are continually renewing and remodeling to maintain barrier in normal states and after injury (Chung 2010; Marshall 2016; Müller-Röver 2001; Rajan 2003). In hairy skin, Merkel cell nerve endings are clustered into specialized epithelial structures called “touch domes” (Plikus 2008). Mice have hundreds of touch domes in their hairy skin and humans have similar, yet subtly different nerve endings as well.  More broadly, Merkel cell receptors are found throughout both hairy and glabrous skin in mammals, though local receptor and skin structures vary in each instance.  In receptive populations, such afferents help to signal information regarding the edges and curvature of stimuli, among other attributes (Johnson 2001).

The dynamics of the architecture of the Merkel cell-neurite complex is just beginning to be understood (note abstraction given in Figure 1, and Lesniak 2014). A Merkel cell’s connection to or removal from a terminal neurite, and neurites from heminodes, have been observed to follow trends specific to the stages of the hair cycle in the mouse (Marshall 2016). In particular, as mice age, multiple synchronized hair cycles are observed, where the hair of the animal changes over its entire body in a wave-like fashion (Müller-Röver 2001). There are four stages of the spontaneous mice hair cycle: First Telogen: 4 weeks, Anagen: 5-6 weeks, Catagen: 6 weeks, and Second Telogen: 9-10 weeks. After this point, the hair cycle begins to enter a mosaicking phase whereby the hair over the body of the animal does not change in a wave-like fashion but instead hair is lost and regrown at different rates from seemingly random positions over the skin surface. Cutaneous neurons and Merkel cells may engage plasticity mechanisms during hair-follicle regeneration; however, the dynamics and physiological consequences of neuronal plasticity in touch receptors are not entirely understood (Moll 1996; Nakafusa 2006).

Observational research efforts regarding arbor remodeling are restricted as we are currently unable to trace specific end organs through the hair cycle. A modeling approach, in contrast, can allow for detailed traceability of end organs and each of their components through every stage of the hair cycle. Moreover, despite instances of discrete event simulation (DES) models in biological research, there are still very few published examples of comprehensively validated models. As such, our group recently built a computational model to test an initial set of rules governing arbor remodeling mechanisms (Marshall 2016). With its top down approach, population statistics from observed data were used to construct the computational model. Using the population statistics as reference points, end organ constituents were iteratively created and deleted to reflect observed data. Each of the transitions between hair-cycle stages was modeled as a separate simulation. Within each simulation, there were a number of iterations where Merkel cells and/or heminodes could be added or removed according to just four probabilistic policies, which were informed by morphometric data. That effort identified four policies of Merkel cell and heminode dynamics from observational data and evaluated the policies.

While the study (Marshall 2016) contributed to explaining rules that govern Merkel cell-neurite (or arbor) remodeling processes, those computational models have several limitations. For instance, the individual end organ constituents of Merkel cells, terminal neurites and heminodes were not traceable between phases of the hair cycle. In addition, being constrained by a top down approach hindered the formation of arbor remodeling policies central to individual components of the end organ, rather than the end organ itself.  As well, because of the model structure it was not easy to perform what-if scenario tests and determine optimal parameters associated with the governing rules. Finally, there are no easy optimization or experimentation features within this model, restricting the testing of new parameters that were not previously defined. Therefore, we have sought to overcome these limitations by employing a DES approach. In specific, the objective of this study is to demonstrate how a DES modeling approach can help understand underlying mechanisms of changes in arbors over the course of the hair cycle. In contrast to the prior computational modeling efforts, that presented herein takes on a bottom up approach, with encoded probability rules for each specific constituent as the parent end organ goes through the hair cycle. Policies are a combination of conditional and probabilistic rules. Unlike the computational model, each Merkel cell, terminal neurite branch and heminode relating to a particular end organ was examined individually, and rules governing their addition and deletion were reflective of hypotheses regarding biological characteristics. This bottom-up approach, which models end organs, heminodes, and terminal neurite branches as different entities, allowed us to refine the governing policies and estimate optimal parameters that play a role in the remodeling processes. The Merkel cells are modeled as an attribute of terminal branches that affects the policies and parameters.

METHOD

The objective of this work is to use DES to identify principles that specify how arbors change over the course of the hair cycle and to evaluate policies of Merkel-cell and heminode dynamics. Simulating different combinations of rate evolutions and incorporation rules are done to yield predictions of end organ arbors for comparison to spontaneous hair cycle observations.

Abstraction of Biology for Use in Describing the Modeled Inputs, Outputs, and Transitions

We consider the Merkel cell-neurite complex in the hairy skin of the mouse.  While the architecture of the Merkel cell-neurite complex can be even more complicated in terms of its constituent parts, the level of abstraction given Figure 1 is somewhat typical and will be used to describe the inputs and outputs of the model to be built herein.

Based on the prior modeling effort (Marshall 2016), a preliminary set of rules was developed, which have been modified in the present work to come from a bottom-up perspective as shown in Table 1.  In general, we check first the state of a heminode by checking the associated terminal branches and Merkel cells.

Discrete-Event Simulation

Model Description

A DES model was built to reproduce the dynamics occurring during transitions in the hair cycle and determine key parameters that affect the dynamic behaviors. The model consisted of three entity types, each of which represents end organs, heminodes, and terminal branches. Merkel cells were modeled as a binary attribute of terminal branch entities: every terminal branch either has a Merkel cell (1) or no Merkel cell (0). The structure of the arbors was modeled as a hierarchy of two batching layers. Each end organ entity consists of a batch of heminode entities, and each heminode entity consists of a batch of terminal branch entities. As each end organ goes through a complete hair cycle, its composition continuously changes through batching/un-batching  processes with different rules.  The model was built using the Simio Simulation and Scheduling Software Package.

Parameter Estimation and Validation

To develop the baseline model, we used observational data obtained from the previous study (Marshall 2016). The major parameters in First Telogen included the number of heminodes per end organ, the number of terminal branches per end organ, and the number of Merkel cells per end organ, the number of terminal branches per heminode, and the number of Merkel cells per heminode. These parameters were estimated using empirical distributions.  The duration of each process was 4, 2, and 4 weeks between the stages of the hair cycles, respectively.

Table 1. Graphics- and text-based depiction of rules associated with each of three transitions of the hair cycle, where MC = Merkel cell and TB = Terminal branch and where  represents a probability of a component of an end organ x being removed/pruned in stage s x includes a MC, TB, and HN. Term s includes 1T (First Telogen), A (Anagen), C (Catagen), and 2T (Second Telogen).

One of the goals of this study was to estimate the probabilities of the addition and deletion of each Merkel cell, terminal branch, and heminode relating to an end organ between stages. To calibrate those unknown model parameters, three main outcomes were considered: the average number of heminodes per end organ, the average number of Merkel cells per end organ, and the average number of Merkel cells per heminode. The model calibration was conducted in each stage chronologically, one at a time, because probabilities in different stages are independent of each other, but each hair cycle stages population statistics are dependent on the previous stage. The calibrated baseline model was run with 25 replications and validated by comparing its results with the three main outcomes of the observational data. Sensitivity analysis was performed to understand the impact of changes in the deletion/addition probabilities on the main outcomes between Catagen to Second Telogen. We chose the two stages because more remodeling activities involving heminodes, terminal branches, and Merkel cells occur during the last transition of hair cycle compared to other transitions.

The baseline simulation model was validated using the observational data Marshall et al. (2016) collected.   Table 3 compares the mean number of Merkel cells and Heminodes at each stage derived from 20 observational data with those estimated from the DES model with 20 replications. The model outcomes were considerably similar to the actual data, which indicates that the baseline simulation model is a reasonable representation of the arbor system (Banks 2000).

Table 2. Parameter values used in the simulation model in each of three transitions of the hair cycle, where MC = Merkel cell and TB = Terminal branch and where  represents a probability of a component of an end organ x being removed/pruned in stage s x includes a MC, TB, and HN. Term s includes 1T (First Telogen), A (Anagen), C (Catagen), and 2T (Second Telogen).

Figure 3 compares the results of the observed and simulated data. Relying solely on cues related to the hair cycle, model dynamics produced arbors in close agreement with experimental observations. The distributions of heminodes, terminal branches, and Merkel cells per end organ at each hair cycle and patterns of changes were similar between the observational data and simulation outputs. The simulation model also allowed us to trace each end organ over the course of hair cycle. For example, blue and orange dots in panels D, E, and F in Figure 3 represent two different end organ, respectively. However, the dynamic remodeling behaviors between the experimental observation and model outputs were not compared at the end organ level because of the lack of traceability in the observational data.

Figure 3. Comparison between the (A-C) actual observational data and (D-F) DES model data.  Shown are data regarding each of the heminodes, terminal neurite branches and Merkel cells across each of the four stages of the hair cycle.  In panels D-F, blue and orange dots are used, respectively, to trace two separate end organs through their stages of the hair cycle.  This style is tied as well to the data in Figure 4.

The end organs denoted by orange and blue dots marked in panels D-F of Figure 3 are detailed further in Figure 4.  For example, in Figure 4A, end organ 1 (EO 1) starts in First Telogen with twelve terminal branches and five heminodes where the terminal branches are distributed across the heminodes with 4, 4, 2, 2, 0 per each heminode.  Then when end organ 1 moves to the Anagen stage of the hair cycle, its number of terminal branches decreases to five.  These five terminal branches are distributed across the five heminodes with a correspondence of 2, 1, 1, 1, 0.  Furthermore, note that in the transition from Catagen to Second Telogen the number of heminodes decreases from five to four for this end organ.

A sensitivity analysis was performed to understand which rules have a greater influence on end organ remodeling behavior. First, rules associated with terminal branch and Merkel cell removal from First Telogen to Anagen were examined. During the transition, a deletion of a terminal branch depends on  whether the terminal branch is populated with a Merkel cell or not. For non-populated terminal branches, a 0.1 increase in the deletion probability led to a 20.5% difference in the total number of terminal branches, while the same change for the populated terminal branches led to a 1% difference. Also, two rules governing heminode pruning from Catagen to Second Telogen were evaluated. The analysis showed that as the probability of heminode deletion increases by 0.05 for heminodes without Merkel cells, the number of heminodes changed by 6.7% on average. On the other hand, the same change for populated heminodes only resulted in an average of 1.6% difference in the number of heminodes. From these results, we can conclude that the amount of end organ constituents is more sensitive to deletion policies of end organ constituents that are less populated.

(A)

(B)

Figure 4. Tracing the flow of two example end organs through their stages of the hair cycle from First Telogen, Anagen, and Catagen, to Second Telogen. Panel A details changes in terminal branches and Panel B relates to Merkel cells.  The blue color represents end organ 1 (EO 1) and orange (EO 2).  For example, end organ 1 begins with twelve terminal branches and twelve Merkel cells.  They are distributed across five heminodes at 4, 4, 2, 2, 0 respectively.  Then moving to the Anagen stage of the hair cycle, the number of terminal branches decreases to 5 and Merkel cells to 1, with distributions across five heminodes of 2, 1, 1, 1, 0 and 1, 0, 0, 0, 0, respectively.