Inferring regulatory networks from experimental morphological phenotypes: a computational method reverse-engineers planarian regeneration
The comprehensive model of planarian regeneration reverse-engineered by our method represents the first quantitative model able to recapitulate regeneration under genetic knock-downs, pharmacological treatments, and surgical manipulations. Unlike conventional arrow diagrams derived from molecular genetic experiments, this system identifies models that not only include necessary components (without which regeneration cannot occur normally), but are fully-specified as a constructive model showing which dynamics are sufficient to give rise to the remarkable pattern homeostasis of planaria. Most models of regeneration are based on generalized mechanisms and do not consider the specific dynamic regulatory mechanisms or network topology necessary to precisely recapitulate the observed patterning phenotypes [92–96]. Meinhardt’s pioneering work on the mechanisms of pattern formation represents the only dynamic models of planarian regeneration proposed to date, based on reaction-diffusion mechanisms and able to recapitulate the head-versus-tail polarity regeneration and midline formation [23, 64, 97, 98]. However, this approach was purely numerical as a proof of the general dynamic mathematical principles, without characterizing any of the regulatory products, and hence accounting only for surgical amputations. Our model was inferred directly from experimental data and includes particular genetic regulatory components able to precisely predict genetic and pharmacological interventions in addition to surgical manipulations. Hence, the models inferred with our method can be used to predict the morphological outcomes in specific genetic knock-downs.
The method can identify those interactions most strongly implied by the dataset, by performing multiple searches and extracting the common pathways found in the resultant set of regulatory networks. Interestingly, the consensus model found in this way includes most of the genetic regulations of head vs. tail planarian regeneration published in the field to date, as well as novel genetic regulations only discovered recently in other model organisms, such as the inhibition of wnt by notum [99]. Furthermore, the method can be used as a generable protocol for automatically finding the less-universal regulatory interactions inferred from the data, and for automatically suggesting additional perturbations for in vivo experimental testing. Importantly, the robustness of the method to infer predictive regulatory networks was validated with a subtraction control test, which successfully produced a regulatory network that not only predicted all the experiments in the dataset used during the search, but also predicted the exact resultant phenotypes from a set of new in vivo experiments that were not part of the search process. In summary, these results validate the capacity of our method to reverse engineer robust regulatory networks with a high predictive power.
Although our method has produced the most comprehensive model of planarian regeneration to date, it contains several limitations. We have restricted our experimental formalization and simulation to 2-dimensional spatial data; thus, the discovered models do not yet address the regulatory mechanisms necessary to specify the dorsoventral axis patterning in planaria [100–103], or the detailed patterning of individual internal organs. In addition, the discovered models are deterministic, and do not account for the stochasticity shown by some partial penetrance phenotypes. Adding a stochastic component to such equations does not represent any technical difficulty; however, the computational requirements of the method to quantify the frequency for each of the possible resultant phenotypes in each experiment would increase by several-fold. Because the basic paradigm is fundamentally very flexible, future work will address these limitations, leading to further improvements in the ability to reverse engineer models that are more complete, including specific modeling of the numerous cellular mechanisms that physically implement such outcomes, such as cell migration, division, differentiation, or apoptosis.
Our approach is broadly applicable to any model system whose experimental procedures and anatomical outcomes can be formalized [104, 105] and can readily be extended to other problems in morphogenesis, including embryonic development or the programmed self-assembly of hybrid systems such as bioinspired robots [106–109]. The models discovered with this method allow the identification of the key mechanisms and the major regulatory products, including those directly perturbed during the experiments as well as-yet unidentified necessary products, explaining the resultant experimental phenotypes. Such models are required for the identification of intervention strategies to produce desired changes in large-scale shape, for birth defects, regenerative medicine, or synthetic bioengineering research. Our method represents a proof of principle towards the use of evolutionary search and quantitative spatial simulation to help constructively understand complex morphological outcomes in embryogenesis, regeneration, and synthetic bioengineering.
We created the input datasets of formalized experiments for the automated search algorithm with the software tool Planform [70]. Planform uses a functional ontology based on mathematical graphs [47], a set of interconnected nodes [110], to unambiguously describe the main characteristics of the morphology, including the overall shape and the location of specific phenotypic regions (head, trunk, and tail regions in the worm). Experimental procedures are described in the ontology as a nested set of basic operations, including amputations, cuts, genetic/pharmacological perturbations, and their parameters. Using the graphical user interface, we created a separated dataset with the phenotypic experiments presented in each of the main publications of head-versus-tail planarian regeneration [72–79], and an additional dataset including all the experiments together.
To apply an automated discovery system capable of finding complex spatial and temporal dynamic networks, we modeled the behavior of gene, protein, and metabolite regulatory network with a system of nonlinear partial differential equations (PDE). Products can act as intercellular signals or be confined intracellularly, decay with time, and be activated or inhibited by other products in the regulatory network. Each product can be regulated by several other products, where interactions can be combined in either a necessary or sufficient fashion.
A regulatory network is made of phenotypic products and signaling products. Phenotypic products represent phenotypic regions in the organism, and as such cannot regulate other products. In this way, morphological features of the phenotype are abstracted as a single product, resulting in inferred regulatory networks centered on the signaling mechanisms and not the molecular details to form specific morphological features. For example, full-body worm phenotypic data are formalized using head, trunk, and tail regions, whereas the inferred networks employ corresponding specific products representing head, trunk, and tail outcomes.
Signaling products can regulate other products, and they can represent the product of specific genes (such as β-catenin or wnt1) inferred from the perturbation experiments in the dataset, or be found de novo as necessary by the search algorithm. In addition, special products are used to model specific aspects of an experiment. In particular, we implemented a wound signal product, which is produced in the area adjacent to a surgical cut during an experiment.
Each equation in the system models the production rate of a product as the linear relation between a production term, a decay term, and a diffusion term. The production term is modeled with a combination of Hill functions, a widely-used nonlinear model of biochemical interactions and genetic regulation [111]. Each Hill function models the activation or repression of a product by another product (including itself). A product can be regulated by several regulatory interactions simultaneously, and these interactions can be grouped in a necessary (both regulators are required to produce the regulated product), sufficient (one regulator is enough to produce the regulated product), or any combination of them. Sufficient interactions are grouped together in a max operator, while necessary interactions are grouped together in a min operator; the set of sufficient interactions is considered as a necessary interaction by itself, and hence it is included inside the min operator. The rate or production is modulated by a production constant, which multiples the result of the combined Hill functions regulation. Products decay in an exponential fashion. Thus, the decay term is modeled with a decay constant multiplying the product current concentration. Intercellular signaling mechanisms are essential in the regulation of developmental and regenerative processes. We modeled the propagation of intercellular signals as a diffusion term in the differential equation, modulated by a diffusion constant. This allows the implementation of products that can propagate intercellularly, carrying signals regulating other products. The diffusion constant of a product can be zero, in which case the product is considered exclusively intracellular.
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): ∂a∂t=ρamin(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+Da∇2a 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 D a is the diffusion constant.
In summary, each product of a regulatory network is defined according to four parameters (production, decay, and diffusion constants, and the initial concentration value), while each regulatory interaction is defined with three parameters (Hill coefficient, dissociation constant, and whether the regulation is necessary or sufficient). The values of all the parameters are automatically inferred by the search algorithm.
We implemented a simulator able to load morphological phenotypes and perform surgical, genetic, and pharmacological experiments formalized with the functional ontology. The simulator takes as input a formalized experiment with the functional ontology and a regulatory network described with a system of PDEs. The simulator outputs the resultant morphology after numerically integrating the PDE system and performing in silico the formalized experiment. Due to the dynamic boundaries of a developmental and regenerative simulated organism, we implemented an Euler finite difference method [112] to integrate the set of PDEs corresponding to a regulatory network.
The simulator performs an experiment in two stages. During the first stage, the original wild-type morphology is loaded into the simulator, the initial product concentrations are set according to the model, and the regulatory network defined in the PDE system is integrated for a fixed amount of time. This first stage allows the dynamical system to converge into a steady state, which will be used for the initial state in the second stage. A formalized morphology and regulatory network is loaded into the simulator by setting the initial concentration value for every product in the regulatory network. The concentration of phenotypic products (head, trunk, and tail in the worm dataset) is initialized according to the corresponding phenotypic regions. For example, the positions corresponding with a head region will be initialized with head product concentration of 1.0 and trunk and tail product concentrations of 0.0, whereas the positions corresponding with a trunk region will be initialized with trunk product concentration of 1.0 and head and tail product concentrations of 0.0. Signaling products are initially set homogenously according to either a numerical parameter for each product between 0.0 and 1.0 stored in the model or indicating a configuration similar to a phenotypic product.
During the second stage of the simulation, the surgical manipulations and genetic and pharmacologic perturbations are performed and the PDE system is integrated for another fix amount of time. The final state of the system is the resultant morphology of the in silico experiment. Surgical manipulations change the boundary of the system and set the concentration of all products outside of the new boundaries to zero. A genetic knock-down (RNAi) will eliminate all the activation regulations of the corresponding product. Pharmacological treatments (octanol) will set the corresponding product diffusion constant to zero, simulating a block of gap junction channels.
To calculate the error (predictive power) of a regulatory network, the resultant phenotypes of the simulated experiments using the network are scored by comparing them with the resultant phenotypes from the physical experiments. For this end, we implemented a distance metric between phenotypes.
Wild type planarians can vary their size by about an order of magnitude due to feeding and starvation [113], a common situation during regenerative experiments where worms may lack the ability to feed. In consequence, we made the distance metric between phenotypes tolerant to small variations between phenotypes. More precisely, the phenotypic metric is invariant in scale, for which the phenotypes are first centered and scaled before comparison. In addition, we included a concentration tolerance parameter (ε) and a radius tolerance parameter (r) within the metric, as defined below, which even out small differences between phenotypes.
The goal of the search algorithm is to find regulatory networks that produce stable phenotypes, and not transient states that are only temporally similar to the resultant phenotypes of the physical experiments. To bias the search towards stable networks, we included a concentration change penalty that is applied when the maximum concentration change in the last time step of a simulation is higher that certain parameter threshold (μ).
We then define a Euclidian distance between two locations a and b that measures the squared averaged distance between a set p of phenotypic products within a tolerance ε: ∥a−b∥ε=∑pϵP(([p]a−[p]b)2−ε)+ where [p]a and [p]b are the concentrations of product p in the locations a and b, respectively.
Then, we define a distance metric between two phenotypes A and B of size w-h as the mean logarithmic minimum distance between every location a of phenotype A and every location b inside a radius r from a of phenotype B: d(A,B)=1w⋅h∑i=1w∑j=1hlog(1+minδ,θ∈(−r,r)∥ai,j−bi+δ,j+θ∥ε)
Finally, we define the error of a regulatory network model M for a set E of n experiments as: error(M,E)=1n∑i=1n(d(ΨeiM,Pei)+(ΔeiM−μ)+) where ΨeiM is the resultant phenotype of simulating experiment ei with the model M, Pei is the resultant phenotype from the physical experiment ei, ΔeiM is the average concentration change in the last time step of simulating experiment ei with the model M, and μ is the penalty concentration change threshold. The error of a network is calculated with the set of experiments formalized in the input dataset, plus an additional experiment with no surgical manipulation or perturbation to assure that a discovered regulatory network maintains the correct wild type morphology in the absence of any perturbation.
Having an automated measurement of the error of a given regulatory network model for a set of formalized experiments, we then implemented an optimization method to search for models that minimize the error. Our method is flexible enough to find the parameters, the topology, and the necessary products of the network.
We employed an evolutionary algorithm [114] approach to search for regulatory networks, where a population of candidate networks evolve in parallel until a network with zero error is found. The initial population comprises random networks with random parameters and regulations between the phenotypic products, the wound product, and the genetic and pharmacological perturbed products from the set of experiments to search.
New regulatory networks are produced from existing ones through crossover and mutation operators. A crossover mixes randomly two networks to produce two new networks. Products that are in common between the two networks are copied to the new networks randomly, each network receiving one of each product, while products not shared are distributed randomly between the two new networks. Products are copied to a new regulatory network together with their regulatory links. If the regulatory product of a copied link does not exist in the new network, it is substituted randomly by another regulatory product.
Mutations alter the regulatory network randomly. Each product or link parameter can be substituted by a random value with 1% probability. Products and links are duplicated with 1% probability. After duplicating a product, a new regulatory link and a new regulator link is created for the new product. After duplicating a link, the regulated and regulator products are chosen randomly. Products and links can be deleted with 1.5% probability, except phenotypic and perturbed products in the experiments, which cannot be deleted. These evolutionary parameters are not optimized; however, a higher probability of deletion with respect duplication is necessary to bias the evolution towards simpler networks and prevent bloating [115].
The evolutionary algorithm stops when a network with zero error is found and the complexity (number of products and links) of the simpler network with zero error have not been decreased for a certain number of generations. This extra evolutionary time is used to simplify the best network found, since the mutation operators are biased towards simpler networks.
Since new regulatory networks in a population can be simulated and evaluated independently, we implemented a parallel version of our evolutionary algorithm in a cluster computer using 256 cores. We used an island distribution approach [116], which improves performance and preserves genetic diversity by using many independently-evolving subpopulations. We used 32 parallel subpopulations with 64 regulatory networks each. For every 250 generations on average, all subpopulations are randomly paired and their regulatory networks are shuffled randomly; this compensates the trend for a single suboptimal regulatory network to saturate a single subpopulation.
We used the deterministic crowding selection method [117] with 75% crossover, 1% parameter change mutation, 1% duplication mutation, and 1.5% deletion mutation. All the parameters in the regulatory network can vary in the range (0,1), except the Hill coefficient, which can vary in the range (1,5). To calculate the regulatory network error, we used an Euclidian distance tolerance ε of 0.1, a distance comparison radius r of 2, and a penalty concentration change threshold of 10–4. We used 250 extra generations in the criteria to stop the algorithm after a network with zero error is found.
The simulation and search method was implemented in C++ using the Standard Library and the Eigen library (http://eigen.tuxfamily.org). Visualizations used the Qt libraries (The Qt Company Ltd.) and the Qwt library (Uwe Rathmann and Josef Wilgen). The software is freely available at http://www.daniel-lobo.com/planarianmodels.