by Hyojung Kang, Rachel Orlowsky, and Gregory J. Gerling (University of Virginia)
As presented at the 2017 Winter Simulation Conference
Abstract
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
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.
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.
Results
An example simulation of one end organ through the four stages of the hair cycle is shown in Figure 2. Based on the conceptual deletion/addition policies between hair cycles developed in the previous study (Marshall 2016), we have generated bottom-up probabilities as shown in Table 2 (next page). In the transition from First Telogen to Anagen, the chance a terminal branch is removed is 0.9 if it does not contain a Merkel cell and the chance is 0.5 if there is a Merkel cell. During the transition, a Merkel cell is deleted with a probability of 0.9. In the transition from Anagen to Catagen, Heminodes are more likely to be pruned if their constituents are empty. For example, a heminode is pruned with a probability of 0.2 if it does not contain any terminal branches Merkel cells, with a probability of 0.1 if it has terminal branches but not Merkel cells, and with a probability of 0.05 if it contains both terminal branches and Merkel cells. In the transition from Catagen to Second Telogen, heminodes follow a similar proportional deletion policy. During the transition, heminodes lose at least one terminal branch with a probability of 0.8. The number of terminal branches to be deleted is between 1 and 3 with the equal probability. A Merkel cell is deleted with a probability of 0.9.
Figure 2. Example simulation of an end organ through the four stages of the hair cycle. For instance, in the transition from First Telogen to Anagen the Merkel cells (MC) decrease from 12 to 1, the terminal neurite branches (TB) decrease from 12 to 5, while the heminodes (HN) remain at a constant number of five.
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).
Table 3. Validation of the baseline model.
Observational data (Mean) |
Model outputs (mean [95 % CI]) |
|||
Merkel Cells |
Heminodes |
Merkel Cells |
Heminodes |
|
First Telogen |
13 |
5 |
13 [12.4, 17.3] |
5 [5.3, 5.5] |
Anagen |
0 |
5 |
0 [0, 0.7] |
5 [5.2, 5.5] |
Catagen |
25 |
4 |
24 [24.1, 24.8] |
4 [3.9, 4.2] |
Second Telogen |
15 |
3 |
14 [13.8, 14.2] |
3 [3.2, 3.9] |
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.
Discussion
As noted in the introduction, this work sought to employ a DES approach 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. The choice of a 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 model appears to well replicate the observed data and speculate on a set of probabilistic values that might be driving the biological processes, at this level of abstraction. Simulating different combinations of rate evolutions and incorporation rules allowed us to represent dynamics in end organ arbors for comparison to spontaneous hair cycle observations and helped fill the gaps left from limitations in observational research studies. For example, the entity-based bottom-up approach allowed for detailed traceability of end organs and each or their components through every stage of the hair cycle, which is restricted in observational research. Also, using the DES model, we were able to determine refined probabilities of addition and deletion for each constituent part of the end organ over each of the hair cycle stages. That said, several assumptions were made in building the DES model. First, it was assumed that observed mean values in the hind limb spontaneous hair cycle dataset are accurate enough to make logical policies. Second, heminodes with the smallest number of Merkel cells were assumed most likely to be removed. On the other hand, Merkel cells in large Merkel cell clusters were assumed more likely to be deleted. Lastly, it was assumed that Merkel cells are randomly added to heminodes.
The applications of a DES approach to problems in biological research have been limited compared to its applications in other areas (Hunt 2009), such as healthcare more broadly. This study demonstrates that DES can be used as a promising add-on to commonly used research practices in biological research. By testing potential biological rules before moving to mouse models, both time and money might be well resourced. Various what-if scenarios tests and sensitivity analysis using DES models can also help design additional experiments by informing critical factors closely related to emergent dynamics. For instance, this model seeks to afford parameter experimentation features. The incentive behind such a model is to study possible confounding experimental factors, such as timing and end organ conditions, that could potentially influence the behavior of an end organ and its constituents’ behavior through the hair cycle. Specifically, experimentation can transpire through simulation of those irregularities in combination with the base model probability rules to determine the actual robustness of the determined probabilistic rules.
Author Biographies
HYOJUNG KANG is a Research Assistant Professor in the Department of Systems and Information Engineering at the University of Virginia. She holds a PhD in Industrial Engineering and Operations Research from Penn State University. Her research interests include developing and applying hybrid simulation models to provide implementable solutions for complex systems, particularly healthcare systems. Her email address is hkang@virginia.edu.
RACHEL L. ORLOWSKY is an Undergraduate Student in the Departments of Systems and Information Engineering and Biomedical Engineering at the University of Virginia. Her interests lie in computational modeling applied to biological domains. Her e-mail address is rlo4kc@virginia.edu.
GREGORY J. GERLING is an Associate Professor in the Departments of Systems and Information Engineering and Biomedical Engineering at the University of Virginia. He earned his PhD from the University of Iowa. His research interests include computational neuroscience, haptics, human factors, human-machine interaction. His chief interests lie in understanding the complexity that underlies the tactile sensory system, which informs our percept of touch and the design of haptic devices. His email address is gg7h@virginia.edu.