Discovering novel phenotypes with automatically inferred dynamic models: a partial melanocyte conversion in Xenopus

Progress in regenerative medicine requires reverse-engineering cellular control networks to infer perturbations with desired systems-level outcomes. Such dynamic models allow phenotypic predictions for novel perturbations to be rapidly assessed in silico. Here, we analyzed a Xenopus model of conversion of melanocytes to a metastatic-like phenotype only previously observed in an all-or-none manner. Prior in vivo genetic and pharmacological experiments showed that individual animals either fully convert or remain normal, at some characteristic frequency after a given perturbation. We developed a Machine Learning method which inferred a model explaining this complex, stochastic all-or-none dataset. We then used this model to ask how a new phenotype could be generated: animals in which only some of the melanocytes converted. Systematically performing in silico perturbations, the model predicted that a combination of altanserin (5HTR2 inhibitor), reserpine (VMAT inhibitor), and VP16-XlCreb1 (constitutively active CREB) would break the all-or-none concordance. Remarkably, applying the predicted combination of three reagents in vivo revealed precisely the expected novel outcome, resulting in partial conversion of melanocytes within individuals. This work demonstrates the capability of automated analysis of dynamic models of signaling networks to discover novel phenotypes and predictively identify specific manipulations that can reach them.


Results
In silico discovery of a novel phenotype. We employed a stochastic model of the regulation of melanocyte conversion in Xenopus, uncovered by a Machine Learning approach 21 , to computationally discover a novel phenotype and the pharmacological perturbations necessary to produce it. Figure 2A shows the inferred dynamic model, including the entities used to perturb the system in past work 21 , the signaling elements participating in the pathway 20,25 , the incidence of conversion produced by various phenotypes, unknown required components inferred (predicted) by the system, and all the interactions between all the components. The perturbation experiments consist in the application of a precise combination of reagents during Xenopus development, which stochastically results in a wild-type phenotype (correct numbers of round melanocytes in their normal locations) or a completely converted phenotype (embryo covered extensively with highly-arborized melanocytes that invade many internal tissues and blood vessels) (Fig. 1). Interestingly, no intermediate phenotype was observed in the training data set (20 experiments with none, one, or two reagents applied) or the validation data set (seven experiments with two or three reagents applied)-tadpoles converted or not, stochastically but always in an all-or-none manner for each animal 21 . Could this coordination be broken, and if so, how?
We then performed in silico all possible experiments involving the application of up to three reagents, for a total of 576 different experiments, each repeated 100 times. Figure 2B shows, for each experiment, the average distance from the resultant phenotype (degree of conversion) to either of the extreme phenotypes: normal or complete conversion (see methods section for the exact formulation of the distance metric). In green are those experiments (combination of reagents) used in the training dataset to infer the model, in red are those used in the validation dataset, and in blue are all the experiments that had never been performed in vivo. The results show that all except one experiment (red arrow) are below the 0.1 threshold to an extreme phenotype (normal or converted). The experiment above the 0.1 threshold indicates an intermediate phenotype between the previously observed wild-type normal pigmentation and the converted phenotype (i.e., a novel situation in which only some melanocytes within a single animal convert to the metastatic-like phenotype). This phenotype is predicted by the inferred model to result from a combined application of altanserin (an R2 inhibitor), reserpine (a VMAT inhibitor), and VP16-XlCreb1 (constitutively active CREB).
Dynamics of the partially converted phenotype. We first investigated in silico the dynamic characteristics of the predicted partially converted phenotype with the inferred model using the combined application of the drugs altanserin and reserpine, and constitutively active CREB. Figure 3 shows the phase portraits of 100 simulations of the wild type and the treated embryos with the trajectories of the state of two serotonin receptors (5HT-R1 and 5HT-R2) and the degree of conversion through time. The simulation starts in the instable initial state (yellow dot) and, due to the inherent stochasticity of the model, ends in one of several stable attractors, representing a final partial converted state. Without treatment (wild type, Fig. 3A contrast, applying the discovered treatment with three drugs, the embryos stochastically develop with an intermediate conversion level (green dot, 75% of simulations) or a normal conversion level like the wild type (blue dot, 25% of simulations). In contrast, animals resulting from any of several treatments 20,39 that depolarize instructor (GlyR-expressing) cells or perturb their downstream serotonergic signaling exhibit extensive overabundance and uniform coverage with highly arborized melanocytes (red arrowhead). This occurs in an all-or-none manner in some percentage of the animals (frequency depending on the specific manipulation) 21,25 . Sectioning (level of section shown in schematics to the right) reveals the main features of melanoma-like phenotype: over-proliferation, arborization, and invasiveness. (C) Normally round melanocytes dorsal to the neural tube (green arrow) become highly arborized and drop down over the neural tube itself (D, red arrows). Inset panels show blood vessels, normal in c' and covered by invasive melanocytes in d' , as occurs in melanoma. Sections further along the tail likewise show small numbers of round melanocytes in control larvae (E, green arrows) compared to the excess of long, abnormally extended and ectopically localized melanocytes in converted animals (F, red arrows). In converted animals, cells can be seen invading the lumen (G) or neural tissue (H) of the neural tube, also often forming networks (I) as has been observed in vasculogenic mimicry of mammalian melanoma 40   The "partially converted" attractor (green dot) was not observed in any previous simulation or experiment in vivo 21 , representing a novel bifurcation in the signaling network's dynamics due to the specific perturbation caused by the application of the three reagents. In this way, the inferred model predicted a partially converted phenotype at about 25% between wild-type and totally converted.
The effects of the three reagents in the treatment can be observed in the temporal dynamics of the level of conversion and pathway components targeted by the drugs (Fig. 4). Embryos with no treatment (wild type, Fig. 4A) stochastically develop normal levels of conversion (99% of the simulations) or high levels of conversion (1% of the simulations). These embryos always have high levels of VMAT and low levels of CREB. However, the levels of serotonin receptor R2 are stochastically either low (blue line, corresponding to high levels of conversion) or high (all other colors, corresponding to normal levels of conversion). Embryos under the three drugs treatment (Fig. 4B), stochastically develop either the normal level of conversion (25% of the simulations) or the partial conversion phenotype (75% of simulations). The plots show how the inhibition drugs reserpine and altanserin causes VMAT and R2, respectively, to converge to a low level, whereas the constitutively active CREB causes CREB levels to reach a high transient state and then converge to a lower level, but overall staying higher than the wild type.
In vivo validation of the computationally discovered treatment and phenotype. We then sought to validate the partially converted phenotype predicted by the inferred model by applying in vivo the newly discovered combination of reagents. Control embryos exposed to vehicle alone were 99% normal (N = 258, Fig. 5A-A"). Embryos exposed to ivermectin (which depolarizes instructor cells; the effect has been characterized extensively 20,25 ) exhibited 97% conversion (N = 313, Fig. 5B-B"). In contrast, embryos expressing constitutively active CREB (tom-VP16-XlCreb1) and simultaneously treated with a cocktail of Altanserin (R2 blocker, 10 uM) and Reserpine (VMAT blocker, 100 uM) exhibited the predicted never before-seen partial conversion in 22% of the animals (N = 167, Fig. 5C-D"), while 24% were completely normal, and 54% exhibited the standard full conversion phenotype.

Discussion
Here, we have presented the discovery of a novel in vivo phenotype, and the precise intervention required to obtain it, by means of a computational methodology. From a reverse-engineered signaling network able to recapitulate the stochasticity of two phenotypes representing wild-type and extreme levels of melanocyte conversion in Xenopus laevis larvae, we sought to discover a pharmacological intervention that could generate a never-seen partially converted phenotype. We performed in silico all the possible experiments combining one, two, or three reagents and found that only one specific combination that was able to produce an intermediate phenotype. The model dynamics of this intervention showed that this precise combination produced a bifurcation in the dynamical system, generating an attractor corresponding to a phenotype 25% in between the known wild type and totally converted phenotypes. Almost a decade of highly diverse experiments focused on this pathway had not revealed such an outcome 20,25 . While detailed mechanistic dissection of the pathway explained signaling on a single cell level, explanations for the stochastic concordance throughout each animal, the possibility of partial conversion, or the reagents that could achieve it, had all remained unknown without the described automated approach that strongly augmented human analysis of this complex dynamical system.
We validated these predictions in vivo, confirming the in silico discovery of the novel phenotype. When we applied the discovered combination of three reagents, two drugs and constituently active CREB, we observed for the first time the predicted partial converted phenotype-animals in which some melanocytes and melanocyte-free regions were normal, while others were converted and colonizing ectopic sites. It is crucial to note that each of (B) Applying the treatment combining the discovered three reagents (altanserin, reserpine, and constitutively active CREB), the dynamics reveal a bifurcation in the dynamical system and the appearance of a new attractor with a partial conversion level (0.25, green dot) in addition to another attractor with a normal conversion level (0, blue dot). The attractor dot size is proportional to the number of converging trajectories to that final state. these reagents alone causes some incidence of full conversion-the individual reagents have already been validated to work in Xenopus (i.e., they do penetrate the embryo, reach target cells, and affect this pathway sufficiently to fully convert some percentage of each treated sample) 21 ; however it is of an all-or-none character. Moreover, we have previously shown that injection of just a few cells in the embryo (for example, with depolarizing reagents) is sufficient to cause complete conversion 21,25 -the partial conversion here is not due to mosaic expression of the CREB construct, as injected alone, in one cell of the 32-cell stage embryos as above, causes an incidence of ~45% complete conversion. The model predicted the wild-type phenotype within 1% (24% in vivo vs 25% in silico), yet it aggregated the partially and fully converted phenotypes observed in vivo (22% partial, 54% fully converted) into a single partial phenotype (75% partial). This difference in the frequency of partial and fully converted phenotypes between the in silico and in vivo results may be explained by the lack of time information in the training dataset and the observed dynamics of the system (Figs 3 and 4). In particular, a number of trajectories ending in the partially converted phenotype actually are crossing the levels of highly converted phenotypes. This observation may indicate that the full conversion that we observed in vivo may be transient. Further investigations, and methods that enabled the tadpoles to survive to older stages for tracking any subsequent changes in melanocyte status, would be necessary to validate predictions of transient (time-dependent) changes of the phenotypes predicted by the model. Indeed, the training dataset does not contain any time information, and hence the current system does not yet capture realistic timing, which would allow a finer-course analysis of trajectories and comparison with in vivo results; in future work we will extend this platform to include time-series data, so that these transient phenotypes can be better investigated.
Identification of interventions that alter system-level physiology or pattern formation is a crucial nascent field with many implications for biomedicine 32,33 . The complexity of biological regulation makes it very hard to completely and mechanistically understand signaling networks, the prediction of resultant experimental outcomes, and hence the discovery of novel phenotypes. Computational dynamic modeling of complex regulatory networks is essential for the discovery of new capabilities of biological systems and prediction of specific perturbations that result in desired outcomes. Using these models, we can computationally and systematically interrogate the system to find the precise interventions necessary for the development of specific phenotypes 33 . The computational methodology presented here shows how a reverse-engineered signaling system can be used in this way to precisely find a set of reagents that produce a never before seen outcome at the intersection of embryogenesis and metastatic-like conversion. These computational methods are paving the way for the discovery and further understanding and manipulations of multicellular, stochastic, and multidimensional biological phenomena, as well as for the discovery of novel and effective therapies that target cooperative or group aspects of cell function in vivo.

Methods
Animal husbandry. Xenopus embryos were maintained according to standard protocols 34   In normal animals' tails, melanocytes never exit the stripes of axial muscle (do not venture into the fin). (B) Converted animals produced by Ivermectin-induced depolarization exhibit spread-out (arborized) melanocytes all along the flank (B) and around the eyes (B')-they colonize every region of the head. The arborized melanocytes also leave the mid-body and colonize the fin (B"). (C) In contrast, embryos treated with VP16-XlCreb1and exposed to 5HT-R2 and VMAT blockers revealed a never before-seen phenotype where some regions of the animal were converted (C, periocular region of the right eye, red arrowhead), and some were wild-type (C, periocular region of left eye, blue arrowhead). This can be also seen in the transformation of cell shape over the brain (C', red arrowhead) but normal melanocyte morphology elsewhere (C' , blue arrowhead). In the tail (C",C'") we observed some cells invading the fin, but they had normal morphology. We observed a great diversity of different combinations that were intermediate between normal and converted outcomes (only a representative sample is shown here), which otherwise never occurred in the same animal. (D) Furthermore, we noticed a new, cancer-like behavior that had not occurred in any previous experiment. This included the induction of individual nodules, which were both pigmented and unpigmented (D,D'); this occurred in ~10% of the animals that likewise exhibited partial conversion (D", red arrowhead region vs. blue arrowhead region).
Scientific RepoRts | 7:41339 | DOI: 10.1038/srep41339 M2014-79. All experiments were performed in accordance with relevant guidelines and regulations. Although we have not observed light to make any difference to the converted phenotype, all groups were subjected to the same light conditions. Microinjection, drug exposure, sectioning, and image collection. Embryos were injected with mRNA encoding constitutively active CREB (tom-VP16-XlCreb1 36 ) at 1/32 cell stage and then exposed to a mix of Altanserin (5HT-R2 blocker, 10 uM) and Reserpine (VMAT blocker, 100 uM) using the same methods described previously 20,21,37,38 . Sections were made at 40 μ on a Leica vibratome after JB4 embedding. Images were collected using a Nikon SMZ-1500 microscope and each phenotype was visually characterized as having normal melanocyte status, partial conversion, or full conversion; these three classes are discrete and readily distinguished as categorical data; the quantification of the definition of these states was described in refs 20 and 39.
Dynamic signaling model of melanocyte conversion. We used the dynamic signaling model of Xenopus melanocyte conversion presented in ref. 21. The model is based on Hill-kinetics and includes 14 stochastic ordinary differential equations describing the interactions of signaling molecules, pharmacological compounds (drugs, set to 0 when not applied and 1 when applied), and the level of melanocyte conversion (the phenotype, equivalent to the pigmentation level) of a developing embryo. The interaction between signaling molecules and drugs are modeled with a combination of Hill functions, each of them representing the interaction between two products. When multiple interactions regulate a product, they can be grouped as sufficient, necessary, or a combination of them. Sufficient interactions can regulate the production of the target product independently of other sufficient interactions (modeled with a max 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 min operator, similar to a logic AND interaction; hence, inhibition has preference over activation). This scheme allows the modeling of a rich variety of regulatory interaction combinations. In addition, each product includes a decay term (with a specific decay constant) and a Gaussian random noise of zero mean to model the stochasticity of the system.
The following set of differential equations, including all the kinetic parameters, define the complete model (ξ(t) represents a Gaussian random noise of zero mean): This set of equations (defined in the model) together with a set of initial conditions of the signaling products (defined in the model) and the drugs (defined in the experiment) form an initial value problem which is then numerically solved until reaching a steady state. Due to the stochastic nature of the equations, different runs with the same initial conditions (same experiment) can result in different phenotypic outcomes (level of melanocyte conversion) as observed in vivo. The model was reverse-engineered directly from a dataset of 20 experiments combining different pharmacological drugs and resulting in a specific frequency of two different and extreme phenotypes (within each phenotype, the animal-to-animal variation in pigmentation level was not visually detectable): either normal (pigmentation level 0) or high level (pigmentation level 1) of melanocyte conversion, never an intermediate level as where d(E) is the mean distance of experiment E (a combination of reagents) when simulated 100 times with the model and resulting in phenotypes with level of conversion p i (the level of conversion is a dynamic variable in the model, with values from 0 for the wild type to 1 for a total conversion).