Research ArticleDevelopmental Biology

Serotonergic regulation of melanocyte conversion: A bioelectrically regulated network for stochastic all-or-none hyperpigmentation

See allHide authors and affiliations

Sci. Signal.  06 Oct 2015:
Vol. 8, Issue 397, pp. ra99
DOI: 10.1126/scisignal.aac6609

Driving melanocyte proliferation and invasion

Melanocytes play key physiological functions; one of the easiest to see is pigmentation. In frogs, the number, distribution, and shape of melanocytes are determined by a subpopulation of cells called “instructor cells,” which are regulated by changes in membrane potential. Forced depolarization of instructor cells can result in excessive melanocyte proliferation, altered melanocyte cell shape, and abnormal migration of melanocytes into multiple tissues, which results in darkly colored tadpoles through a stochastic all-or-none process; the embryos are either normally pigmented or hyperpigmented. Lobikin et al. unraveled the molecular signaling pathway and physiological circuit that mediates this melanocyte conversion process, and they used computational approaches to explain how this all-or-none, stochastic process can occur.


Experimentally induced depolarization of resting membrane potential in “instructor cells” in Xenopus laevis embryos causes hyperpigmentation in an all-or-none fashion in some tadpoles due to excess proliferation and migration of melanocytes. We showed that this stochastic process involved serotonin signaling, adenosine 3′,5′-monophosphate (cAMP), and the transcription factors cAMP response element–binding protein (CREB), Sox10, and Slug. Transcriptional microarray analysis of embryos taken at stage 15 (early neurula) and stage 45 (free-swimming tadpole) revealed changes in the abundance of 45 and 517 transcripts, respectively, between control embryos and embryos exposed to the instructor cell–depolarizing agent ivermectin. Bioinformatic analysis revealed that the human homologs of some of the differentially regulated genes were associated with cancer, consistent with the induced arborization and invasive behavior of converted melanocytes. We identified a physiological circuit that uses serotonergic signaling between instructor cells, melanotrope cells of the pituitary, and melanocytes to control the proliferation, cell shape, and migration properties of the pigment cell pool. To understand the stochasticity and properties of this multiscale signaling system, we applied a computational machine-learning method that iteratively explored network models to reverse-engineer a stochastic dynamic model that recapitulated the frequency of the all-or-none hyperpigmentation phenotype produced in response to various pharmacological and molecular genetic manipulations. This computational approach may provide insight into stochastic cellular decision-making that occurs during normal development and pathological conditions, such as cancer.


Bioelectrical mechanisms contribute cell-to-cell communication and large-scale pattern formation processes. Slowly changing ion currents and the resulting spatial gradients of plasma membrane resting potential across cell fields affect the signaling pathways that control cell proliferation (1, 2), differentiation (35), and migration (6, 7). Spatiotemporal gradients of transmembrane potential (Vmem) contribute to the regulation of many patterning processes during embryogenesis and regeneration (812). Bioelectric signaling is also implicated in cancer biology (1316), with the recognition of ion channels as oncogenes (1719) and as drug targets (2022), as well as with the identification of bioelectrical states that are associated with neoplastic transformation and tumor prevention (2326).

One aspect of embryogenesis is the coordination of migration, proliferation, and morphogenesis of specific subpopulations of cells. To study the role of bioelectric cues in coordinating these processes, we use the model organism Xenopus laevis, which is a system in which biophysical cues, molecular genetic signaling (signaling by gene products), and morphogenetic outcomes can be readily characterized. A proportion of tadpoles in which instructor cells are experimentally depolarized acquire a hyperpigmented phenotype as a result of altered melanocyte morphology and behavior (6, 27): The melanocytes overproliferate, acquire a highly abnormal arborized cell shape, and invade the neural tube, brain blood vessels, and gut, properties that are characteristic of dysregulation of the normal epithelial-to-mesenchymal transition that occurs during embryogenesis (28, 29). Instructor cells are characterized by the presence of the glycine-gated chloride channels GlyR (also known as GlyCl), are present throughout the somatic tissues of the body, and control melanocyte proliferation, shape, and migration. The hyperpigmentation phenotype results from increased numbers, altered cell shape, and ectopic localization of melanocytes, not changes in melanin production or pigment granule behavior. This conversion process involves the whole body, not only the skin (6, 27), and is induced by a bioelectrical signal, not requiring DNA mutation or chromatin-modifying drug exposure.

Depolarization of the instructor cells stimulates a serotonin-dependent signal, involving several serotonin receptors, the serotonin transporter SERT, and the vesicular monoamine transporter VMAT (6, 30). Induction of hyperpigmentation requires serotonin receptors 5HT-R1, 5HT-R2, and 5HT-R5 (6, 30). These serotonin receptors are heterotrimeric guanine nucleotide–binding protein (G protein)–coupled receptors (GPCRs) that signal through either the second messenger systems of adenylyl cyclase (AC) and adenosine 3′,5′-monophosphate (cAMP) or phospholipase C (PLC) and inositol trisphosphate (31). Because both 5HT-R1 and 5HT-R5 are GPCRs coupled to Gi/o, which inhibits cAMP signaling, and cAMP signaling plays a role in melanoma (32), we used pharmacological tools to explore the involvement of cAMP signaling in control of melanocyte behavior and found that depolarization of instructor cells signaled to melanocytes through cascades involving cAMP production and the transcriptional activity of the cAMP response element–binding protein (CREB), as well as the transcription factors Sox10 and Slug, all of which are implicated in developmental regulation of melanocytes (33) and pathogenesis of melanoma (32, 3438).

We also found that conversion of melanocytes by instructor cells involved melanocyte-stimulating hormone (MSH)–secreting melanotrope cells of the pituitary. To uncover previously unknown targets of long-range bioelectrical signaling, we performed a genome-wide analysis of transcriptional changes in embryos in which the instructor cells were depolarized.

One puzzling aspect of this phenotype is that hyperpigmentation is stochastic, yet binary, at the level of the organism: A given treatment affects only some percentage of the manipulated embryos, and every animal is either hyperpigmented or not. Because currently available reverse-engineering methods for steady-state or time-series concentration data cannot be used to infer stochastic models directly from phenotypes (3944), we developed an automated method based on evolutionary computation to derive a stochastic network model that predicted how divergent body-wide phenotypes can stochastically arise. Our method reverse-engineered the necessary components, regulatory interactions, and parameters of a stochastic ordinary differential equation model directly from training data sets and accurately predicted the stochastic phenotypes resulting from pharmacological manipulation of the embryos. Our approach revealed the molecular steps and circuits by which change in the membrane voltage of dispersed, specific cell populations can affect the behavior of widely distributed cells of other types.


Ivermectin-induced instructor cell signaling is mediated by cAMP, PKA, and CREB

To depolarize instructor cells of 1-day [stage 10; note that all stages are defined as Nieuwkoop and Faber (NF) (45)] Xenopus embryos, we used the GlyR-opening drug ivermectin, which causes the efflux of chloride ions from the GlyR-positive instructor cells (6, 46). Because melanocytes do not express the gene encoding GlyR (6), ivermectin does not directly affect the melanocytes. We measured the percent of the population that became hyperpigmented at stage 45, 7 days of ivermectin exposure (Fig. 1, A to C). Inspection of tadpole morphology, behavior, and survival indicated no overt toxicity or teratogenesis due to hyperpigmentation, consistent with previous studies (6).

Fig. 1 cAMP and CREB are involved in mediating instructor cell signaling.

(A and B) Frogs treated with pharmacological agents at stage 10 are scored for pigmentation phenotypes at stage 45. (C) Percent of the population exhibiting hyperpigmentation in embryos at stage 45 after treatment at stage 10 with the cAMP antagonist 2′5′-dideoxyadenosine (2′5′-DDA) with or without ivermectin (IVM) or the AC activator forskolin. (D) Percent of population exhibiting hyperpigmentation in embryos at stage 45 after treatment at stage 10 with the protein kinase A (PKA) antagonist H89 with or without ivermectin or injected into one cell at the four-cell stage with the mRNA for the indicated CREB with or without ivermectin. A-CREB encodes dominant-negative (DN) CREB; XlCreb1 encodes wild-type (WT) CREB; and VP16-XlCreb1 encodes constitutively active (CA) CREB. In (C) and (D), sample sizes (number of embryos) are noted for each condition. Error bars represent 1 SEM. *P < 0.0001 based on Pearson’s χ2 test.

We used a suppression screen strategy (6, 30, 47) to identify signals involved in mediating the hyperpigmentation phenotype. We applied inhibitors of various aspects of cAMP signaling together with ivermectin and assessed the effect on the frequency of hyperpigmentation. Treating stage 10 embryos with the cell-permeable cAMP antagonist 2′5′-dideoxyadenosine and ivermectin significantly reduced the incidence of hyperpigmentation that occurred with ivermectin treatment alone (87% hyperpigmented with 1 μM ivermectin alone and 67% hyperpigmented with 500 μM 2′5′-dideoxyadenosine and ivermectin) (Fig. 1C). We also exposed stage 10 embryos to ivermectin and forskolin, which directly stimulates AC to increase cAMP, or to forskolin alone and found that forskolin alone led to a significant incidence of hyperpigmentation compared with the occurrence in control animals (Fig. 1C). These results suggested the involvement of cAMP in Vmem-induced changes in melanocyte behavior.

Changes in intracellular cAMP abundance regulate the activity of the cAMP-dependent PKA. To determine whether PKA was involved in the hyperpigmentation phenotype, we exposed stage 10 Xenopus embryos to ivermectin and the PKA antagonist H89 dihydrochloride. Coapplication of H89 resulted in a significant reduction in the incidence of hyperpigmentation compared with embryos exposed to ivermectin alone (Fig. 1D), suggesting that PKA is involved in this signaling process.

cAMP and PKA signaling affect gene expression through CREB. Translocation of activated PKA to the nucleus phosphorylates Ser133 in CREB, thereby activating it (48). To address whether CREB is involved in mediating the hyperpigmentation phenotype in Xenopus, we injected one cell of a stage 5 (32-cell) embryo with synthetic mRNAs encoding either wild-type CREB (XlCreb1), dominant-negative CREB (A-CREB) (49), or a constitutively active version of CREB (VP16-XlCreb1) (50). Injection of the RNA for wild-type CREB had no effect on the percent of hyperpigmented animals, whereas injection of the RNA for constitutively active CREB resulted in hyperpigmentation of ~42% of tadpoles (Fig. 1D). Injection of the RNA for the dominant-negative CREB into one cell of a stage 5 embryo had no significant effect on the percent of hyperpigmented animals when injected alone and compared to control uninjected embryos or when combined with ivermectin and compared with embryos exposed to ivermectin alone (Fig. 1D). Together, these results implicate the involvement of cAMP-PKA-CREB pathway in transducing signals from the depolarized instructor cells to melanocytes.

Hyperpigmentation phenotype involves the pituitary gland

Instructor cells are sparse yet widely distributed throughout the entire embryo and the responding neural crest–derived melanocytes, which can be many cell lengths away, suggesting that the signal may not be direct from the instructor cells to the melanocytes. Noticing that rare embryos with spontaneously missing midline structures never underwent melanocyte conversion after ivermectin expression, we suspected involvement of one of the unpaired (centrally located) endocrine organs. Because melanotrope cells of the pituitary gland secrete αMSH, which promotes the dispersion of melanin granules in the Xenopus melanocytes (51), we hypothesized that the pituitary may also be involved in controlling other aspects of melanocyte behavior.

We individually removed regions containing the pineal gland, the pituitary, or, as a control, a region below the cement gland (Fig. 2, A and B) from tail bud–stage (stage 32) embryos that had been exposed to ivermectin continuously starting at the neurula stage (stage 10). The first melanocytes develop at stage 33 or 34; therefore, we removed the pineal or pituitary glands before the appearance of the first melanocyte. The embryos were reared to free-swimming tadpole stage and scored for hyperpigmentation. Removing the pineal gland had no significant effect on ivermectin-induced hyperpigmentation, whereas removing the pituitary gland decreased hyperpigmentation incidence by ~50% compared with ivermectin-exposed controls (Fig. 2, C and D). These results indicated that hyperpigmentation involved the pituitary and revealed that depolarization-induced hyperpigmentation can be interrupted by removing the pituitary until at least stage 32.

Fig. 2 The pituitary gland is necessary for ivermectin-mediated hyperpigmentation.

(A) Side view showing the main regions of Xenopus tadpole brain. (B and C) Xenopus embryos were treated with ivermectin from neurula stage (stage 10), cut at tail bud stage (stage 32), raised to tadpole stage (stage 45), and scored for hyperpigmentation (HP). Cuts were performed on tail bud stage (stage 32), removing the pineal gland, pituitary gland, or a control region below the cement gland and away from both the pituitary and pineal. (D) Effect of control cuts (n = 54 embryos), pineal cuts (n = 46 embryos), or pituitary cuts (n = 59 embryos) on the percent of ivermectin-induced hyperpigmented tadpoles. Effect of inhibition of MSH with αMSH release–inhibiting factor (MSH-RIF), an MSH agonist SHU 9119, or a 5HT receptor antagonist altanserin on the percent of hyperpigmented tadpoles. Error bars represent 1 SEM. *P < 0.0001, Pearson’s χ2 test; NS, not significant.

Because the pituitary releases αMSH, a known regulator of melanocyte function, we exposed Xenopus embryos to αMSH agonists and antagonists. Nearly all of the embryos exposed to the αMSH agonist SHU 9119 from the neurula to tadpole stages developed a hyperpigmented phenotype (Fig. 2D). Embryos exposed to MSH-RIF and depolarizing ivermectin displayed 32% less hyperpigmentation incidence than did those exposed to ivermectin alone (Fig. 2D).

To determine whether serotonin signaling and αMSH were part of the same pathway, we exposed the embryos to combinations of inhibitors of αMSH and 5HT signaling. Either MSH-RIF or altanserin, a 5HT-R2–specific antagonist, significantly decreased ivermectin-induced hyperpigmentation when applied to embryos individually (Fig. 2D). Embryos exposed to all three pharmacological agents (ivermectin, MSH-RIF, and altanserin) showed a similar proportion of hyperpigmented tadpoles as those exposed to either ivermectin and MSH-RIF or ivermectin and altanserin, indicating that both serotonin and αMSH signaling are in the same pathway.

Instructor cell depolarization results in transient increase in Sox10 expression

Sox10 is a transcription factor that promotes the specification of neural crest progenitors to the melanocyte lineage (33). Slug is a member of the Snail family of transcription factors that are important in the epithelial-to-mesenchymal transition in neural crest cells (37, 52). Depolarization of specific regions of developing Xenopus embryos by forced expression of genes encoding depolarizing ion channels increases the expression of genes encoding several transcription factors, including Sox10 (27). Whereas injection of mRNA encoding depolarizing ion channels into blastomeres results in depolarization of all of the injected cells and daughters thereof, ivermectin depolarizes only the instructor cells. Therefore, we compared the effects of ivermectin and injection of XminK, which encodes a subunit of a depolarizing potassium channel subunit, on Sox10 and Slug expression by quantitative real-time polymerase chain reaction (RT-qPCR) of tail bud–stage embryos (Fig. 3, A and B). We used RT-qPCR of whole embryos to obtain quantitative results that could not be achieved by in situ hybridization in the ivermectin-treated embryos. Either injection of XminK at the one- or two-cell stages or exposure to ivermectin at stage 12 resulted in a similar induction of Slug in the entire tail bud–stage embryo. Both also induced Sox10, but injection of XminK produced a larger increase in Sox10 abundance (Fig. 3B).

Fig. 3 Both localized depolarization and sparse, widely distributed depolarization increase Sox10 transcripts.

(A) Expression of Sox10 in embryos injected with a mixture of XminK and β-gal mRNAs at the four-cell stage that were fixed at tail bud stages. Far left shows injection, and middle shows the section plane for the data shown at the right. (B) Effect of ivermectin treatment or XminK injection on the expression of Sox10 and Slug in tail bud stage embryos (stage 25) as assessed by RT-qPCR. XminK-injected animals were injected into one cell at the four-cell stage. Ivermectin-treated animals were exposed to the drug from neurula stage (stage 10) onward until processing for RT-qPCR at stage 25. Control animals were uninjected and untreated. Red dashed line denotes no fold change compared to control. All experimental treatments resulted in significant increase in expression of both Sox10 and Slug compared to control (P < 0.05, Student’s t test; n = 10 embryos per sample, samples run in triplicate, three biological replicates). (C) Effect of ivermectin exposure started at neurula stage on Sox10 expression as assessed by RT-qPCR in embryos collected at the indicated stages. ST, stage. NF stages 15 to 35, n = 10 embryos per sample, samples run in triplicate, three biological replicates; NF stage 45, n = 5 embryos per sample, samples run in triplicate, three biological replicates. P < 0.05, Student’s t test. (D to F) In situ hybridization for Sox10 performed on stage 15 control (CTRL) embryos (D) or embryos that had been exposed to ivermectin treatment starting at stage 10 (E and F). (F) is an enlargement of the area boxed in (E), and arrowheads indicate positive staining outside the main Sox10-positive area.

Because ivermectin resulted in a lower than expected increase in Sox10 expression in the tail bud–stage embryos, we analyzed Sox10 by RT-qPCR from embryos exposed to ivermectin at stage 10 and then fixed at different stages. This revealed that Sox10 mRNA increased in abundance by stage 15, which occurs ~12 hours after the addition of ivermectin under our experimental conditions, and then decreased at later stages (Fig. 3C). To determine the spatial profile of Sox10 expression, we analyzed stage 15 embryos that had been treated with ivermectin by in situ hybridization. In control embryos, Sox10 expression was limited to symmetrical regions along the neural fold (Fig. 3D). In contrast, ivermectin-treated embryos displayed ectopic Sox10 expression in a punctate pattern throughout the embryo (Fig. 3, E and F), in a pattern similar to that of the GlyR-expressing instructor cells (6). Moreover, in contrast to the consistent expression pattern of Sox10 in controls, the ivermectin-treated embryos exhibited a patchier signal, suggesting a reduction in expression in the regions where Sox10 is normally expressed.

These data indicated that similar to regional depolarization through mRNA injection, ivermectin, which produces abnormally depolarized Vmem only in the instructor cells, can alter gene expression, in particular, producing foci of ectopic Sox10 expression throughout the embryo in a pattern similar to that of instructor cells, but Sox10 induction is transient compared with that achieved by regional depolarization by mRNA injection.

Instructor cell depolarization produces a distinct pattern of gene expression

Given that instructor cell depolarization increased the expression of genes encoding Sox10 and Slug, two neural crest–specific transcription factors, we performed transcriptional microarray analysis of whole embryos at early neurula (stage 15) and tadpole (stage 45) stages, comparing stage-matched controls and embryos treated with ivermectin from stage 10.

We found that ivermectin exposure resulted in 45 transcripts (Fig. 4A and table S1) that were differentially expressed by stage 15 embryos and 517 transcripts that were differentially expressed in stage 45 embryos (Fig. 4A and table S2). Hierarchical clustering of the data revealed that developmental stage, rather than ivermectin, was the main driver of the differences in gene expression (fig. S1). Of the five mRNAs that were differentially expressed in both stage 15 and stage 45 embryos exposed to ivermectin, three are unknown and represent expressed sequence tags (one for UniProt ID Xl.52612 and two for Xl.52879). The two genes that were differentially expressed in embryos of both stages in response to ivermectin exposure were HIG1 (encoding hypoxia-inducible domain family, member 1A) with a 2.4-fold increase in the early embryos relative to the controls and a 4.5-fold increase in late embryos relative to the controls, and ANKRD37 (encoding ankyrin repeat domain 37) with a 2.9-fold increase in the early embryos relative to the controls and 9-fold increase in the late embryos relative to the controls. Too few differentially expressed genes were identified in the early embryos for functional enrichment analysis; nevertheless, our data showed that five mRNAs were altered in abundance within a few hours of depolarization of instructor cells (in the stage 15 embryos) and remained differentially expressed from controls at the second time point (in the stage 45 embryos).

Fig. 4 Human homologs of the genes with increased expression resulting from ivermectin-induced depolarization are associated with neoplasms and cancer.

(A) Venn diagram of the transcripts that were differentially regulated (increased or decreased) in transcriptional microarrays of Xenopus tadpoles treated with ivermectin starting at stage 10 and collected at early (stage 15) and late stage (stage 45). See fig. S1 for the differentially expressed transcripts clustered by stage and condition. See tables S1 and S2 for a list of the transcripts and genes. (B and C) Pathway analysis of the differentially expressed transcripts in stage 15 embryos (B) and stage 45 embryos (C). Proteins are red shapes, diseases are purple boxes, stimulatory regulatory events are indicated by an arrow and a plus sign on the relationship line, inhibitory regulatory events are indicated by a blunt line, and arrows without any sign indicate direct binding of proteins.

We identified sufficient differentially expressed transcripts in the later stage embryos (the same stage at which the hyperpigmentation phenotype is strongly evident) to perform enrichment analysis using the GOrilla database (53) with Homo sapiens as the reference species. We analyzed the transcripts for enrichment of Gene Ontology (GO) for “biological process” and “molecular function,” which identified the immune response, metabolism, and ion regulation (table S3) as the major enriched processes.

Pathway analysis (54) using human homologs of the frog genes revealed that 45 mRNAs that were differentially expressed in response to ivermectin in the early stage 15 embryos were involved in the cell processes related to chromatin remodeling, apoptosis, cell differentiation, and mitosis (data S1). We also used the program Pathway Studio to determine whether the human homologs of the differentially expressed genes were associated with diseases, which revealed that 6 of the 45 differentially expressed genes were associated with neoplasms (Fig. 4B).

Pathway analysis by subnetwork enrichment analysis (54) indicated that 9% of the ivermectin-affected genes (45 of 517) in stage 45 embryos were associated with cancer (Fig. 4C and table S4). The results of the pathway and disease analyses implied that the transcriptional alterations that occurred in response to depolarization of the instructor cells may be similar to those associated with disease states and developmental dysregulation in humans.

Modeling predicts both the large-scale dynamics and the quantitative stochastic response to instructor cell depolarization

The hyperpigmentation phenotype is both bistable at the level of the individual (no partially hyperpigmented “salt and pepper” tadpoles occur) and stochastic at the level of the population (different conditions result in different percentages of hyperpigmented and normal tadpoles within a treated cohort). To better understand this stochastic process, we undertook an automated computational modeling approach to identify a model of a signaling network that could characterize all of the known necessary signaling pathways controlling this behavior as a result of the embryo-wide depolarization that produces embryo-wide hyperpigmentation. To identify a model that accurately predicted hyperpigmentation frequency in the population, (i) we generated a simple arrow diagram network with a minimal set of regulatory connections and (ii) we automated the remaining steps in the process—the simulation of a given mathematical model with a specific set of parameters, the comparison of the results with data obtained experimentally, and the assignment of specific link parameters, new connections, and yet-to-be-discovered nodes within the network to discover new mathematical models (text S1). The model included both intracellular and intercellular signaling events but did not include parameters representing the spatial location of each of the regulatory events and elements; thus, it is a nonspatial, dynamic model.

Using this process, we discovered a mathematical model of ordinary differential equations representing every known and three unknown elements in the network with a defined starting concentration value that was either specified (for the drugs) or inferred by the automated method. Elements in the network could be a signaling molecule (such as cAMP or serotonin), a protein (such as SERT), or a pharmacological compound (such as ivermectin). The changes in their concentrations over time (dynamics) can involve regulatory interactions with other elements in the network, exponential decay, and fluctuations due to random noise (text S1). Once we had defined the elements in the starting network, we defined a set of required regulatory links connecting the elements in the network controlling depolarization through melanocyte activation, as suggested by the functional experiments in this and previous literature (6, 30, 38, 5558). Note that the starting network did not include GlyR, the target of ivermectin, because the automated method could identify and add this element reliably if necessary, and we preferred to start with the minimal number of humanly defined elements and regulatory links in the model. These required links could not be removed during exploration of the network model space, yet their parameters needed to be defined by the automated method. Using an automated iterative simulation and modeling process (described in text S2), we reverse-engineered a network that exhibited bistable and stochastic behavior, resulting in wild-type or hyperpigmented outputs with frequencies matching experimental data (Fig. 5A). To reverse-engineer diverse regulatory phenomena, we used Hill kinetics (59) as a generalization of Michaelis-Menten kinetics and an algorithm based on evolutionary computation (44, 60).

Fig. 5 Inferred dynamic stochastic model recapitulating the observed phenotype frequencies.

(A) Reverse-engineered network that produces hyperpigmentation phenotype at frequencies matching those observed experimentally. Starting network elements were the drugs (text without a shape) and the named nodes and blue interaction lines, representing known elements and relationships in the hyperpigmentation (HP node) pathway. Elements added algorithmically included three unknown required components (nodes a, b, and c) and green interaction lines. The model also described the kinetic parameters (see text S3 for details). Multiple signaling interactions are combined together in a necessary (dashed lines) or sufficient (solid lines) fashion (see Materials and Methods for details). Activating regulatory interactions are indicated with arrows, and inhibiting interactions are indicated with blunt lines. (B) Simulations of the model showing the dynamic changes in concentration of each colored node after running the model to steady state and the proportion of normally pigmented (pigmentation value of 0) and hyperpigmented (pigmentation value of 1) tadpoles at the end of the simulation. Simulations of control (no treatment), the presence of cyanopindolol, an inhibitor of 5HT-R1, and ivermectin and cyanopindolol are shown.

The evolutionary algorithm maintains a population of signaling networks, including the initial set of fixed regulatory links in addition to new randomly generated links and random parameters based on biochemically plausible rules. The algorithm periodically refines new network models by crossing and mutating existing network models and by using random operations to add and delete nodes and links and to alter parameters, which enables the exploration of a large set of possible models. By applying error analysis, the algorithm identified the best models in each evolution, discarding the poor scoring network models and proceeding iteratively for a defined number of generations to identify the network model with the best fit to the experimental data.

We generated a set of data using various pharmacological manipulations with or without ivermectin, similar to those described in Fig. 1, and scored the frequency of hyperpigmentation in each condition (data S2). We used a subset of these pharmacological experiments as a training data set using the evolutionary algorithm and error analysis to identify a network model (Fig. 5A and text S3 for the ordinary differential equations) and a second set as a test data set to explore the predictive power of the model.

The model identified using the training data set contained three yet-to-be-characterized elements (labeled a, b, and c in Fig. 5A). Simulations using the model network parameters and various test conditions converged to a steady state in which the hyperpigmentation was complete (values close to 1) or nonexistent (values close to 0); no simulation produced an intermediate pigmentation value, thus matching the all-or-none phenotype. Simulations of the no-treatment condition resulted in a hyperpigmented phenotype in 3% of the cases (Fig. 5B), which is consistent with the frequency of hyperpigmentation that occurred in the control embryos (data S2). Simulations of the effect of cyanopindolol (an inhibitor of 5HT-R1) resulted in the hyperpigmentation phenotype in 32% of cases (Fig. 5B), which matches within 10% error the frequency we observed (Fig. 6A). Simulations of the effect of combined ivermectin and cyanopindolol treatment produced hyperpigmented phenotypes 89% of the cases (Fig. 5B), which matches within 5% error the frequency that we observed (Fig. 6A).

Fig. 6 Inferred model performance with the experiments in the training data set and the experiments in the validation data set.

(A) Percentage of correct outcomes of the model for each of the experiments used to reverse-engineer the model (training data set). Dashed line indicates 85% performance accuracy. (B) Percentage of correct predictions of the model for a set of experiments not used in model generation (validation data set). See data S2 for details. MTP, methiothepin.

The performance (as the percentage of correct outcomes) of the automatically reverse-engineered network compared to the experimental results in vivo was greater than 85% for all the 20 experiments that we used in the training data set during the network search by the evolutionary algorithm (Fig. 6A). Thus, the inferred model uncovered by this process correctly predicted both the large-scale dynamics and the quantitative stochastic outcomes of the data set. We also assessed the performance of the model by comparing the outcome to data from experiments that were not in the training set. All the outcomes predicted by the model for this test data set correctly predicted the outcome, with the model exhibiting a performance greater than 75% (Fig. 6B).

We examined the dynamics of the model with simulations of different treatments to gain an understanding of the all-or-none behavior. Projecting the trajectories of several experiments in the phase space (the state of different components in the network through time) revealed different dynamic attractors (stable outcomes, in terms of product concentrations and melanocyte properties, to which the system can converge) for the hyperpigmented and nonhyperpigmented phenotypes, along with phase bifurcations (a global change in the system behavior) after applying different treatments. We assessed the effect of no treatment (Fig. 7A), ivermectin (Fig. 7B), or forskolin (Fig. 7C) on the trajectories of two serotonin receptors (5HT-R1 and 5HT-R2) and the resulting degree of hyperpigmentation through time in 100 simulations of the inferred model, starting from the initial state (early embryo, yellow dot) and ending in the final stable state (hyperpigmented and nonhyperpigmented tadpole, red and blue dots, respectively). The attractors defined by the dynamical system, representing the final stable hyperpigmented and nonhyperpigmented phenotypes (red and blue dots, respectively, in Fig. 7, A to C), were always located in any of the two extreme pigmentation states (1, completely hyperpigmented; 0, completely nonhyperpigmented), explaining the bistability observed for all treatments.

Fig. 7 Phase space of the inferred model showing stochastic developmental trajectories and pharmacological treatment bifurcations.

(A to C) Trajectories of the state of two serotonin receptors (5HT-R1 and 5HT-R2) and degree of hyperpigmentation in the inferred model during 100 simulations for each of three representative treatments. “0” represents inactive on the 5HT-R axes; “2” or “2.5” represents full activity. “0” on the hyperpigmentation axis represents normally pigmented outcome; “1” represents a fully hyperpigmented outcome. The initial state of the embryo (yellow dot) is the same for all the simulations for a given treatment. The attractors represent hyperpigmented phenotypes (red dots) and nonhyperpigmented phenotypes (blue dots) and are the stable states to which the system converges. The attractor dot size is proportional to the number of trajectories that converge to that phenotype. Note that in the ivermectin condition, there are two attractors for the hyperpigmented state, a representation of how different trajectories can produce the same phenotypic outcome. (D) Qualitative phase-space trajectories summarizing the dynamics of the inferred model during three different treatments. In the model, each treatment represents a change to the initial concentration of the indicated pharmacological drug, resulting in a bifurcation (a global change in the dynamics of the system) in the signaling network, a shift in the attractors, and different frequencies for the resultant phenotypes. The trajectory linewidths are proportional to the frequency of simulations producing the trajectory.

The phase-space trajectory diagrams revealed two additional properties of the network model: the presence of a separatrix (a boundary between two major modes of behavior in the system) and bifurcations in the state space. In each of the examples shown, the trajectories passed near a separatrix. For example, Fig. 7A shows a single trajectory resulting in the hyperpigmented state (red dot), and the rest of trajectories end in the nonhyperpigmented state; a separatrix must lie between these two types of trajectories, defining the border between the basins (regions of influence) of the two attractors. The presence of a separatrix derives from the dynamics of the network model and can explain how trajectories that start at the same initial conditions can converge to different attractors and, hence, produce the observed stochastic penetrance of the hyperpigmented phenotype. The trajectory diagrams also revealed the creation of bifurcations in the state space due to the pharmacological treatments, which produced a variation in the location and basin of the attractors (note the different positions of the red and blue dots, which are especially evident when comparing Fig. 7, B and C). Thus, by altering the dynamics of the network, pharmacological treatments produced different trajectories to the attractors, which can explain the change in frequencies of the resultant stochastic hyperpigmented and nonhyperpigmented phenotypes.

A qualitative summary of the trajectories and attractors under different treatments highlights the bifurcations produced in the state space, the phenotypes associated with the different attractors, and their relative frequency for each treatment (Fig. 7D). Starting with an identical embryo (identical initial conditions), the system defined by the inferred model can stochastically converge to either a hyperpigmented or a nonhyperpigmented phenotype within a specific frequency, and the frequency can be precisely altered with the application of certain pharmacological treatments (which change the initial conditions of the corresponding drugs in the network model). This analysis revealed how an apparently similar hyperpigmented phenotype can result from two different molecular states: either a high level of 5HT-R1 activity combined with a low level of SERT activity (bottom right in Fig. 7D) or a low level of 5HT-R1 combined with a high level of SERT (top left in Fig. 7D).


Previous studies have examined the role of Vmem in regulating the pigmentation of X. laevis embryos and have suggested serotonin signaling and the involvement of multiple serotonin receptors in translating the bioelectrical signal to a change in melanocyte behavior (6, 30). Here, we present aspects of the molecular mechanism of the pathway and evidence for a physiological circuit that controls melanocyte behavior and ultimately tadpole pigmentation (Fig. 8). Furthermore, we identified a dynamic model that quantitatively predicted the stochastic effects of pharmacologically altering various control points of this process. The reverse-engineered model accounted for the bistable all-or-none hyperpigmentation phenotype and for the stochastic penetrance observed under multiple conditions, correctly reproducing the results from experiments used to identify the model and predicting results from independent experimental conditions

Fig. 8 Schematic pathways for melanocyte control downstream of instructor cell signaling.

Under normal conditions the polarized instructor cells produce a relatively small serotonergic signal (stars) due to reuptake and retention of 5HT by the SERT into the instructor cells. At this relatively low concentration of serotonin, only a high-affinity serotonin receptor 5HT-R5 on the pituitary melanotrope cells may become activated, thereby producing normally pigmented animals. Right: When the instructor cell population is depolarized, SERT exports serotonin, resulting in increased serotonin concentrations in the microenvironment of the melanocytes and the pituitary. At this higher concentration, serotonin may activate 5HT-R5 and 5HT-R2 on both the pituitary melanotrope cells and on melanocyte cells. In the melanotrope cells, activation of 5HT receptors stimulates AC activity (not shown), thereby enhancing MSH release. Serotonin signaling increases the abundance of pro-opiomelanocortin (POMC), a precursor of MSH. MSH binding to melanocortin 1 receptors (MC1R) on melanocytes stimulates AC activity (not shown), increasing cAMP production, PKA activity (not shown), and CREB phosphorylation and activity (not shown), thereby increasing the expression of Sox10 and Slug. By increasing cAMP in melanocytes, serotonin may also contribute proliferation, invasive migration, and the altered morphology of the melanocytes. Although melanocytes can also respond to serotonin directly, this pathway is not shown here because the most parsimonious network that explains all of the results, including the long-range signaling by instructor cells, does not require this link.

At the level of the molecular pathway, our pharmacological and genetic manipulations indicated a cAMP-dependent transcriptional pathway as involved in the depolarization-induced hyperpigmentation phenotype. Our results also showed that the depolarization of a sparse cell population (instructor cells) stimulated a similar increase in transcription factor gene expression as that associated with hyperpigmentation mediated by the regional injection of mRNA encoding the depolarizing KCNE1 ion channel (27), which would be expected to depolarize only the injected cells and their progeny, not necessarily instructor cells. We found that the expression of the genes encoding the transcription factors Sox10 and Slug was increased in embryos after ivermectin-induced depolarization. Sox10 is expressed in premigratory neural crest cells, and its expression becomes gradually restricted to cells in the glia or melanocyte lineages (61, 62). In the melanocyte, Sox10 plays an indispensable role in melanocyte survival, proliferation, and migration (63) and is also expressed in primary and metastatic melanoma (3436, 64). We detected increased transcripts for Slug at the tail bud stage after exposure to ivermectin at stage 15, suggesting a persistent increase in expression. In contrast, we found that transcripts for Sox10 only transiently increased after depolarization, gradually returning to baseline by stage 35. These results indicate that these two transcriptional regulators may serve different functions and are subject to different regulatory inputs.

Microarray analysis of stage 15 and stage 45 embryos that had been treated with ivermectin compared to control embryos revealed differential expression of a relatively small set of genes (44) at stage 15 and a larger set (517) at stage 45. Pathway analysis of homologs of the 517 genes in humans indicated that several of them are implicated in different disease conditions, including melanoma (table S4). For example, ivermectin exposure resulted in an increase in the transcripts for POMC. POMC, a precursor to αMSH and melanocortin, has been linked to a number of metabolic disorders (65) and several cancers, including lung cancer (66) and melanoma (6770). The receptor for melanocortin, a peptide derived from POMC, is associated with depressive disorder and is implicated in serotonergic signaling in human disease (71), consistent with our model’s unification of serotonergic, melanocyte regulatory, and cancer-related pathways. Microarray analysis also suggested that transcripts encoding transmembrane serine 2 protease, which is encoded by a gene commonly expressed in prostate cancer (72), were increased by ivermectin exposure and that transcripts for albumin, a soluble, monomeric protein that comprises about one-half of the blood serum protein and that has been implicated in breast cancer (7375), were decreased. Other responsive genes remain to be identified, because we analyzed pooled embryonic cells for RT-qPCR, which likely obscured changes in transcripts that were limited to small populations of cells or that increased in some cells and decreased in others. Application of region- or cell-specific next-generation sequencing approaches should reveal such differentially regulated transcripts.

One curious aspect of this phenotype is its all-or-none character: Tadpoles were never seen to be partially hyperpigmented: Each tadpole is either normal or is covered with abnormal melanocytes. This type of bistable phenomena has been reported in bacteria (76), in synaptic plasticity (77), in developmental processes, including the epithelial to mesenchymal transition (78), the specification of neuron cell types (79), and the maturation of Xenopus oocytes (80), and in pathological situations, such as cancer (81). Signaling networks with feedback loops can produce irreversible switch-like responses (82), which could explain these all-or-none phenomena at the level of individual cells. Creating switch-like responses is a goal of synthetic biology, and the approach that we took may be useful for developing models that can be used by synthetic biologists to create systems with this all-or-none property.

However, at a population level, the hyperpigmented phenotype is stochastic in that many of the conditions produced less than 100% penetrance of the phenotype. The decision of hyperpigmentation occurred not at the level of individual cells but at the level of the individual. Similar tissue-level decision-making has been reported in development during left-right patterning with an entire region becoming left or right (8387), but the mechanisms that achieve this coordinated response are unknown.

Here, we presented and provided a proof of principle for a computational method for reverse-engineering dynamic signaling networks that can explain both bistable and stochastic resultant phenotypes. Starting with a data set containing the incidence data on hyperpigmentation under 20 different pharmacological perturbations, the automated method discovered a dynamic network that recapitulated the bistability and the stochasticity shown in the set of experiments performed in vivo. The automated method evaluated more than 20 million different networks, comprising 40 billion virtual experiments, to reverse-engineer the dynamic network that could recapitulate the stochastic phenotypes in 20 pharmacological assays. The network model is multiscale, including candidates (the diffusible factors MSH, 5 HT, and cAMP) for the global coordination of cell state within each embryo and intracellular signaling relationships, and has only three unknown nodes. Because the network model connects these unknown nodes to other nodes into the model, candidate approaches, for example, by searching protein-protein interaction databases, or untargeted screening methods could be used to identify these unknown nodes.

The model correctly predicted the outcomes of new experiments that were not used in generating the model, validating the automated computational method and the reverse-engineered network. The model also predicted that cAMP signaling did not occur directly through the serotonin receptors (5HT-R1, 5HT-R2, and 5HT-R5) involved in the hyperpigmentation response (the discovered network does not contain any direct link from the serotonin receptors toward cAMP; Fig. 5A), which signal through G proteins that either inhibit cAMP production or are coupled to inositol trisphosphate and diacylglycerol signaling or protein kinase C. In summary, the method that we presented here can lead to the discovery of dynamic signaling networks that provide testable predictions and mechanisms of complex stochastic resultant phenotypes.

The discovered network is a model of signaling cascades mediating the control of cell behaviors by an initial bioelectric signaling event; it is consistent with proposals that some disease-promoting stimuli, such as ultraviolet exposure, mediate their effects through changes of resting membrane potential (88, 89). Moreover, the model suggests a number of tractable control points for regulating aberrant cell behaviors, for example, the use of serotonin reuptake inhibitors, which are used in the treatment of depression, a disease believed to involve aberrant serotonergic signaling (6, 90). If the discovered model reflects a network that is similar to that controlling such cell behaviors in mammals, then the model may represent a starting point for in silico testing of proposed biomedical interventions.


Animal husbandry

Xenopus embryos were maintained according to standard protocols (91) in 0.1× Marc’s Modified Ringers (MMR; pH 7.8). Xenopus embryos were staged according to Nieuwkoop and Faber (92). All experimental procedures involving Xenopus embryos were approved by the institutional animal care and use committees and Tufts University Department of Laboratory Animal Medicine under protocol M2014-79.


Capped, synthetic mRNAs were dissolved in nuclease-free water (Ambion) and injected into embryos at cleavage stages in 3% Ficoll using standard methods (91). The mRNA injections were made using borosilicate glass needles calibrated to bubble pressures of 50 to 70 kPa in water, delivering 50- to 70-ms pulses. After 30 min, embryos were washed and cultured in 0.1× MMR until desired stages. Constructs used included XlCreb1 (gene ID, BC041206), VP16-XlCreb1 (50, 93), injected into one cell of NF stage 5 (32-cell) embryos, and XminK [also known as KCNE1 and isK (27, 94)], injected into one cell of NF stage 3 (4-cell) embryos.

Drug exposure

Stocks of ivermectin (Sigma) were stored at 10 mM concentration in dimethyl sulfoxide. Embryos were exposed in 0.1× MMR (91) from stage 10 to stage 45 in the presence of 1 μM ivermectin (Sigma), 500 μM 2′5′-dideoxyadenosine (Sigma), 5 μM forskolin, 25 μM MSH-RIF, 500 nM SHU 9119, 50 μM cyanopindolol, 10 μM altanserin, 2.5 μM SB 699551, 5 mM 5HT added to medium, 10 nM methiothepin, 15 μM fluoxetine, 100 μM reserpine, and 20 μM H89 (sources are listed in data S2).

All compounds were obtained from Tocris, unless otherwise noted, and prepared in Millipore (18 megohms) water, unless indicated. All drug treatments were performed using embryos from mixed batches of fertilizations, using at least three biological replicates. Control experiments were performed using embryos in normal media (0.1× MMR).


Before microsurgery, stage 32 embryos were anesthetized in a 0.02% tricaine solution (pH 7.5) in 0.1× MMR. Using a surgical blade (Feather #11), tissue was excised from one of the three regions depending on the treatment. In the first experimental group, a wedge-shaped cut was made between the eye and the cement gland, corresponding to the region of the developing pituitary. In the second experimental group, a single cut removed the entire anterodorsal region of the head (above the eye), corresponding to the region of the developing pineal gland. In the third experimental group, a wedge-shaped cut removed the tissue between the cement gland and stomach, which served as a control for wounding. After tissue removal, animals were kept anesthetized for 30 min, after which they were washed twice with 0.1× MMR. After washing, operated animals were raised at 16°C under a 12-hour light/12-hour dark cycle and were scored for the presence or absence of hyperpigmentation at stage 46.

Expression analysis

In situ hybridization was performed as previously described (95). Xenopus embryos were collected and fixed in MEMFA for 30 min (91). Before in situ hybridization, embryos were washed in phosphate-buffered saline (pH 7.4) plus 0.1% Tween 20 and then transferred to methanol through a 25%/50%/75% series. Probes for in situ hybridization were generated in vitro from linearized templates using a digoxigenin-labeling mix (Roche). Chromogenic reaction times optimized signal/background ratio. Analyses represent consistent patterns from 50 to 60 embryos for each marker. Probes used for in situ hybridization include Sox10 (GenBank accession no. AY149116) and Slug (GenBank accession no. AF368040).

RNA extraction and cDNA synthesis

Embryos (collected n = 10 per Eppendorf tube, five biological replicates) were washed in ribonuclease (RNase)– and deoxyribonuclease (DNase)–free water and homogenized in TRIzol reagent (Life Technologies). Homogenized samples were stored at −80°C for up to 1 month. Total RNA was extracted using TRIzol according to the manufacturer’s instructions (Life Technologies). RNA yield and quality were assessed by spectrophotometry (ND-1000, NanoDrop) and gel electrophoresis, respectively, to assess integrity of 28S and 18S RNA.

Reverse transcription was performed using ThermoScript RT-PCR System (Life Technologies). Each in vitro reverse transcription reaction was performed using 1 μg of total RNA and 50 μg of oligo(dT)20 primers (Life Technologies). RNA and primers were mixed, denatured for 5 min at 65°C, and placed on ice before adding the reaction mix according to the manufacturer’s instructions. Reverse transcription reaction was carried out at 50°C for 45 min. The reaction was terminated by incubating at 85°C for 5 min, followed by RNA degradation using 1 μg of RNase H for 20 min at 37°C. The complementary DNA (cDNA) was stored at −20°C until use. The quality and quantity of cDNA were validated using Advantage 2 PCR kit (Clontech) on cDNA samples using α-tubulin primers (96).

Quantitative real-time PCR

Primers (table S5) were designed using Primer3Plus enhanced web interface (97) for Sox10 (AY149116) and Slug (AF368040). Ornithine decarboxylase (ODC), a widely used endogenous control for Xenopus, was used to normalize target gene expression, and primer sequences have been previously published (98). The PCR specificity was verified by BLAST ( using the National Center for Biotechnology Information X. laevis reference sequence. Desalted primers were obtained from Invitrogen by Life Technologies.

For each primer pair, standard curve primer analysis was performed using serial dilutions of cDNA from stage 25 control embryos [1 (undiluted), 10−1, 10−2, 10−3]. Formation of primer-dimer and amplification specificity were assessed by efficiency and melt curve analysis.

The cDNA from validated RNA was used to perform RT-qPCR assays. For each biological sample, three technical replicates were run in each RT-qPCR experiment. Each treatment contained five biological replicates. Triplicate negative controls lacking template were also run for each cDNA sample for each reaction.

PCRs were assembled manually. Samples were prepared by adding 1 μl of cDNA (diluted 1:5 in ddH2O), 10 μl of 2× Power SYBR Green PCR Master Mix (Applied Biosystems), and 0.5 μl of each primer (diluted to 10 μM) in a final volume of 20 μl. Reactions were incubated in 96-well MicroAmp Optical Reaction plates at 95°C for 10 min followed by 40 cycles at 95°C for 15 s and at 60°C for 1 min in a StepOnePlus qPCR instrument (Applied Biosystems).

The RT-qPCR data were analyzed using the StepOne software v.2.3, and ΔΔCT values were calculated (Applied Biosystems). Fold change of target genes relative to the amount of the control gene ODC was calculated as −2ΔΔCT.

Microarray analysis

Gene expression analysis was performed using samples treated with ivermectin from NF stage 10 onward and collected at two developmental stages: stage 15 (early neurula) and stage 45 (free-swimming tadpole). Embryos were collected in Eppendorf tubes (n = 50 for stage 15 and n = 5 for stage 45) and frozen at −80°C. RNA extraction and microarray analysis were performed by the Beth Israel Deaconess Medical Center Genomics and Proteomics Center (Harvard) according to the standard Affymetrix protocol, using a high-throughput hybridization and scanning system. Microarray hybridization was performed using the Affy 3′ IVT Express Kit. Fragmented and biotin-labeled and amplified RNA was hybridized to the GeneChip Xenopus laevis Genome 2.0 Array as per the manufacturer’s protocol. The Affymetrix GeneChip Xenopus laevis Genome 2.0 Array has 32,400 probe sets representing more than 29,900 X. laevis transcripts.

All microarray data were analyzed using Bioconductor packages in R. The quality of hybridized arrays was assessed using Affymetrix guidelines on the basis of scaling factor, background value, mean intensity of chip, and 3′ to 5′ ratios for spike-in control transcripts. The outlier analysis was performed using unsupervised clustering and principal components analysis. All high-quality arrays were normalized using the MAS5 algorithm developed by Affymetrix. The absent or present calls for the transcripts were calculated using the MAS5 algorithm. The differentially expressed transcripts were identified on the basis of fold change and Affymetrix transcripts calls. Transcripts with increased abundance were considered to be those that changed more than twofold in the experimental group compared to the control group, with present calls in the experimental group. Transcripts with decreased abundance were considered to be those that changed more than twofold in the experimental group compared to the control group, with present call in the control group. Differentially expressed transcripts were identified using customized script in R. The list of differentially expressed genes was annotated using the Affymetrix Xenopus platform.


Functional enrichment for differentially expressed genes was conducted using the GOrilla database (53) (updated on 2 August 2014). X. laevis gene names were uploaded and mapped using the official gene symbol. H. sapiens was used as the reference species for gene ontology assignments. Pathway Studio 9.0 (Ariadne) was used for pathway analysis to construct gene interaction networks. For early embryos, 21 genes were mapped to an annotated pathway in the database (data S1), whereas in the late embryo experiments, 210 genes were successfully mapped to an annotated pathway (data S1) [see (99) for additional details on performing pathway analysis].

Model definition

We modeled dynamical signaling networks as systems of stochastic ordinary differential equations. A signaling network was composed of a set of interconnected elements, including the products perturbed during the experiments, the pharmacological drugs used, an element representing the degree of pigmentation in the tadpole, and any additional unknown products necessary to achieve the correct stochastic behavior of the system and the resultant all-or-none states. Each element in a network is represented with an equation, modeling its production rate as the linear relation between a production term, a decay term, and a noise term. The production term is modeled with a combination of Hill functions, each of them representing the regulatory interaction between two products. When multiple interactions regulate a product, they can be grouped as sufficient, necessary, or a combination thereof. Sufficient interactions can regulate the production of the target product independently of other sufficient interactions (modeled with a maximum operator, similar to a logic OR interaction; hence, activation has preference over inhibition). In contrast, all the necessary interactions need to be coordinated to regulate the production of the target product (modeled with a minimum operator, similar to a logic AND interaction; hence, inhibition has preference over activation). This scheme allows the modeling of many different regulatory interaction combinations.

The following equation illustrates a model of the production of product a as regulated by two necessary products (activator b and inhibitor c) and two sufficient products (activator d and inhibitor e):dadt=ρa min(bη1α1η1+bη1,α2η2α2η2+cη2, max(dη3α3η3+dη3α4η4α4η4+eη4))λaa+ξ(t)where ρa is the production constant, ηi are the Hill coefficients, αi are the dissociation constants in the Hill functions, λa is the decay constant, and ξ(t) is a Gaussian random noise of zero mean (see text S1 for further details).

Model simulation

Given a signaling model network described as a system of ordinary differential equations, an experiment can be simulated in silico by numerically integrating the system with a set of initial conditions. The initial conditions include any pharmacological element used in the experiment as a constant equal to 1.0 in the case of its presence in the experiment or equal to 0.0 in the case of its absence. Because of the stochastic nature of the system of equations, the pigmentation level resultant from a simulation for a given set of initial conditions can vary among different runs (see text S1 for further details).

Model evaluation

To calculate the predictive power of a given signaling network model, each experiment (defining a set of pharmacological treatments) was simulated 100 times. The resultant phenotype of a simulation was considered hyperpigmented if the level of the corresponding product was above 0.9 and nonhyperpigmented if it was below 0.1. Then, we intuitively define the error of a signaling network model M for a set E of n experiments as the mean square error between the in silico and in vivo outcomes for each of the experiments:Error(M,E)=1ni=1n((HiHi)2+(NiNi)2)where Hi and Ni are the resultant frequency of hyperpigmented and nonhyperpigmented phenotypes in the in silico simulation of experiment i, and Hi and Ni are the resultant frequency of hyperpigmented and nonhyperpigmented phenotypes in the in vivo assay of experiment i. In addition, the model error used during the search penalizes phenotypes that are not similar to any of the observed resultant phenotypes in the in vivo experiment, as well as models that do not converge to a steady state. See text S1 for a detailed description of the error calculation.

Automated reverse-engineering method

We implemented a computational method based on (44) to automate the discovery of the topology (regulatory links) and parameters of a signaling network that minimize the error for a given set of training experiments. In addition to the training data set of experiments, the algorithm takes as input a set of regulatory links from known interaction pathways that will constrain the generation of candidate models—these sets of links will be included in any model considered during the search. The method uses an evolutionary algorithm approach, in which a set of signaling networks evolves in silico by crossing over and mutating them and replacing those with worse error (100, 101). After a fixed number of generations, in our case 10,000, the best network is returned by the algorithm. See text S2 for further details.

Inferred model of bistable and stochastic hyperpigmentation

We applied our automated method to a training data set of 20 experiments to reverse-engineer a model represented with a system of ordinary differential equations that could explain the bistability and exact stochasticity reported in the experiments with the least error. The experiments included the application of 15 different pharmacological drugs affecting 10 known products (an experiment can combine several drugs) and the resultant frequencies of hyperpigmented and nonhyperpigmented phenotypes. The search was constrained with a set of 11 links, representing known regulatory interactions between signaling products, and additional 17 links, representing the interactions between pharmacological drugs and signaling products. No parameters were constrained in these predefined links, except the well-known activator or inhibitor character of the interactions between drugs and signaling products. The algorithm took 36 hours in a computer cluster with 22 CPUs to reverse-engineer the model shown in Fig. 5A. As shown in the reversed-engineered model, the algorithm inferred, as necessary, 3 additional uncharacterized products (labeled a, b, and c in Fig. 5A) and 24 regulatory links between products (green links in Fig. 5A), as well as a total of 145 numerical parameters (see text S3 for the complete system of equations and parameters).


Text S1. Model implementation, simulation, and evaluation.

Text S2. Method for reverse-engineering a stochastic model of hyperpigmentation.

Text S3. System of equations and kinetic parameters for the reverse-engineered model.

Fig. S1. Hierarchical clustering of the differentially regulated transcripts.

Table S1. A list of the 45 transcripts differentially expressed by stage 15 in response to ivermectin.

Table S2. A list of the 517 transcripts differentially expressed by stage 45 in response to ivermectin.

Table S3. Enriched GO affected in stage 45 embryos.

Table S4. Cancer-related genes in humans of the homologs of the genes differentially expressed in stage 45 Xenopus tadpole embryos after depolarizing ivermectin treatment.

Table S5. Reference genes and primers for RT-qPCR.

Data S1. Differentially expressed transcripts in early and late embryos and their association with disease or cellular process.

Data S2. Results of the training set and validation set experiments and the performance of the model for each.


Acknowledgments: We thank A. Allen, R. Lubonja, and E. Switzer for the general laboratory assistance and frog husbandry, J. Lmire for molecular biology assistance, S. Doughty for high-performance computation support, D. Marshall for assistance with statistical analyses, and the members of the Levin laboratory and the bioelectricity community for many useful discussions. We thank V. Pai, C. Fields, M. McVey, S. Fuchs, and J.-F. Pare for the comments on a draft of the manuscript. Funding: Computation used a cluster computer awarded by Silicon Mechanics and the Campus Champion Allocation for Tufts University TG-TRA130003 at the Extreme Science and Engineering Discovery Environment, which is supported by NSF grant ACI-1053575. This work was supported by the G. Harold and Leila Y. Mathers Charitable Foundation. Author contributions: M. Lobikin performed inhibitor and mRNA injection experiments. M. Levin and M. Lobikin planned the experiments and interpreted data. D.J.B. performed extirpation surgeries. D.L. performed computational modeling, software development, and in silico data analysis. E.T. performed RT-qPCR. C.J.M. performed bioinformatics analysis. All co-authors wrote the paper together. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data have been deposited to the National Center for Biotechnology Information Gene Expression Omnibus database (accession no. GSE70834, platform GPL10756).
View Abstract

Navigate This Article