Engineered Living Systems With Self-Organizing Neural Networks: From Anatomy to Behavior and Gene Expression

In summary, neural growth in neurobots significantly impacts their shape and the distribution of MCCs, and this growth could be potentially increased by modulating serotonergic signaling. The degree to which neural expression contributes to trajectory complexity remains elusive. The lack of observed correlation between neurite growth and complexity index points to the potential heterogeneities in cell type expression and connectivity profiles of neurites across neurobots.

Like the structure and function of the central nervous system, transcriptomes are typically viewed as products of a long evolutionary history of selection. In standard organisms, they are further shaped by ongoing neural inputs [52, 53]. What would the transcriptome of a novel construct with a nervous system look like? With this question in mind, we next asked what changes to default biobot transcriptomes, if any, would be induced by the presence of neural tissue. To characterize the transcriptome of neurobots and compare it to its non‐neuronal counterparts (i.e., biobots and sham neurobots), we performed bulk RNA sequencing of their tissue. For each bot type, four biological samples were included (neurobots: NB1‐4, biobots: BB1‐4, sham neurobots: SH1‐4). Due to the small size of the bots, and, therefore, the small amount of RNA, each sample comprised tissue from multiple bots (see Methods).

We found a high degree of correlation between normalized gene expression levels among all samples within each group (Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced, FPKM [54]), indicating reliability and repeatability of the results (Figure 10a). Moreover, we found that gene expression levels in biobots and sham neurobots were much more correlated to one another than to neurobots (Figure 10a, see Methods). Similarly, neurobots could clearly be separated from shams and biobots based on the principal component analysis on the normalized gene expression value (FPKM) of all samples (Figure 10b, see Methods).

We next compared the gene count variability across the samples of biobots, neurobots, and sham neurobots on a gene‐by‐gene basis (Figure 10c–f). We found that neurobots showed a significantly higher variability, quantified by coefficient of variation (CV, see Methods) in their normalized gene counts (FPKM) across neurobot samples, compared to both biobots and shams, and samples of sham neurobots showed a higher variability compared to samples of biobots (Figure 10c). We then compared pairs of groups (e.g., NB and BB) to determine the fraction of genes that had a greater CV in normalized gene count in one group than the same gene in the other group. For the chosen pair of groups, genes were ranked by the mean count value across all pools of both groups, and the CV of each gene's counts across the pools of each group was calculated (Figure S10, see Methods). The CV list was split into 100 bins (percentiles) containing an equal number of genes, and the fraction of genes in the bin for which the CV of the first group was greater than that of the second group was found and plotted.

We found that in all bins, more than half of the neurobot genes showed higher CV compared to the biobot and sham group (Figure 10d,f; cyan horizontal line is at 0.5). For the sham neurobot group, genes in most, but not all bins showed higher CV compared to biobots (Figure 10c). In addition, we found a trend for genes with higher normalized counts showing a higher degree of variability in their expression when comparing neurobots with biobots and shams (Figure 10d,f). This pattern was significantly different from what is expected from chance in most bins (dark blue bins), that is, relative to the average CV calculated across all bins regardless of the order (orange line, see Methods). The overall higher variability seen in neurobots and shams compared to biobots could be due to several factors. First, in both neurobots and shams, the implanted cells are harvested from ∼50 embryos, whereas a biobot is made out of a single embryo. Moreover, higher variability in the implanted bots is expected due to the high variability in the size of the implants in both neurobots and shams. However, these are likely not the only factors involved, as neurobots showed significant differences in their gene count variability compared to shams. Neural differentiation, therefore, likely plays an important role in the increased variability in gene counts seen in neurobots.

Additionally, we found that neurobots included a significantly larger number of genes that were differentially expressed relative to biobots and sham neurobots (Figure 11a,b), whereas biobots exhibited a smaller subset of differentially expressed genes relative to sham neurobots (Figure 11c). Moreover, the number of significantly upregulated genes (p < 0.05, red dots with positive fold change) in neurobots compared to biobots and shams (6774 and 6859 genes, respectively), were much higher than those that were significantly downregulated (red dots with negative log fold change, 3578 in neurobots vs biobots and 4010 in neurobots vs shams), resulting in highly asymmetric volcano plots (Figure 11a,b). This was not the case when comparing shams with biobots (Figure 11c, 1733 upregulated and 1429 downregulated genes). These results are consistent with a gain‐of‐function as a result of neural growth in neurobots.

We next investigated which biological functions or pathways are significantly associated with the differentially expressed genes. We used Gene Ontology (GO) enrichment analysis to annotate genes to biological processes (bp), molecular function (mf), and cellular components (cc). Due to the large number of upregulated genes in neurobots, we focused on the highly overexpressed genes for this analysis (4 log‐fold or more increase in expression, 2445 genes when comparing neurobots to biobots and 2026 genes when comparing neurobots to shams). The most significantly upregulated pathways in neurobots relative to biobots, as well as in neurobots relative to shams, related to nervous system development, and synapse and neuron projection (Figure 12a,b). One of the interesting genes we found up‐regulated specifically in neurobots relative to biobots and shams is Dact‐4, a member of an evolutionarily conserved family of Dishevelled‐binding proteins involved in the regulation of Wnt and TGF‐beta signaling which is expressed in the Spemann organizer [55]. This suggests that the presence of neurons might exert an organizational influence on the surrounding soma; this molecular signature is consistent with known roles of the nervous system to direct cell behavior in cancer suppression [56, 57, 58], regeneration [59], and embryonic morphogenesis [53, 60]. Trans‐synaptic signaling and neurotransmitter receptor activity were also significantly upregulated in neurobots. This included glutamatergic, GABAergic, cholinergic, dopaminergic, serotonergic, and glycinergic receptors. Surprisingly, exclusively in neurobots, we also found significant enrichment in pathways involved in visual perception (Figure 12a).

Although pathways relating to neuron projection and trans‐synaptic signaling were slightly upregulated in sham neurobots compared to biobots, there were far fewer genes in each pathway, and they were less significant in their degree of upregulation compared to those in neurobots (Figure S8a). Further, there were relatively fewer genes downregulated when comparing neurobots with biobots, and shams (63 and 116 genes respectively), and very few genes were significantly downregulated by 4‐log folds when comparing shams with biobots (38 genes). Our enrichment analysis showed that the largest group of downregulated genes in neurobots compared to biobots belonged to the cellular component pathway (extracellular region in Figure S8b). Interestingly, this pathway includes some of the genes expressed in the Xenopus skin, including glycoprotein 2 (gp2), and mucin (muc17) suggesting that the properties of the “skin” of neurobots might be different from those of biobots.

In order to identify functional biological modules of differentially expressed genes, we extracted the largest protein–protein‐interaction (PPI) sub‐network of these genes using the STRING database. We then performed network embedding and clustering using multi‐nonnegative matrix factorization (MNMF) [61] to find specific functional biological modules. The clusters were subsequently enriched using g: Profiler [62]. Based on this analysis, we identified 25 clusters for neurobots vs biobots comparison and 5 clusters in neurobots vs sham neurobots comparison (Spreadsheet S2, Figures S9 and S10). There were not enough upregulated genes between biobots and sham neurobots to allow for this analysis. Similarly, due to the low number of downregulated genes, the network analysis could not be performed at either a 4‐fold or a 2‐fold change threshold level.

Consistent with the findings from the enrichment analysis, we found clusters containing genes critical for the development of the nervous system, cell fate commitment, and wnt signaling pathways (Cluster 3, Figure S9a, see the NB vs BB tab in Spreadsheet S2). A plethora of growth factors (e.g., various FGFs, BDNF, and EGF), and their receptors were revealed by Cluster 17, which also showed enrichment in enzyme‐linked receptor protein signaling pathways (Figure S9b, see the NB versus BB tab in Spreadsheet S2).

Neurobots contained genes encoding various neurotransmitter receptors including glutamate (e.g., gria1‐4), kainate receptors (e.g., grik1,2,3,5), GABAergic (e.g., gabara3,5; gabarb3), and glycinergic receptors (e.g., glra3; glrb), genes encoding voltage gated calcium channels (cacng3,4,5,7,8), as well as those involved in the uptake of neurotransmitters (e.g., slc1a1,2,3). Genes with important roles in synaptic plasticity were also present in neurobots (arc, camk2b, Cluster 15, Figure 13, see the NB vs. BB tab in Spreadsheet S2) [52]. Cholinergic and muscarinic (chrna2,3,4,5,7; chrnb2,3,4, chrm2,4,5), serotonergic (htr1a,b,e; htr2a), and dopaminergic (drd1,2,4) were among other neurotransmitter receptors (Cluster 23, Figure 13b, see the NB vs BB tab in Spreadsheet S2).

Notably, one of the largest clusters (Cluster 1) contained genes encoding various aspects of visual perception, phototransduction, and photoreceptor development (Figure 13c, see the NB vs. BB tab in Spreadsheet S2). Specifically, this cluster included red and violet cone opsins (opn1lw, opn1sw), retinal G‐protein couple receptors (rgr), melanopsin (opn4, opn5), rhodopsin (rho), as well as many other related genes that encode proteins involved in visual processing. In addition to Cluster 1, Cluster 12 (Figure S9c) also included genes related to eye, lens, and retina development, including genes found in major retinal cell types [63], that is, retinal ganglion cells (neurod1,2; pou4f1), and horizontal cells (onecut1, lhx1). Genes found in bipolar cells (unc5d), and amacrine cells (prdm13) were also present in Clusters 21 and 22, respectively (Figure S9d,e, see the NB vs. BB tab in Spreadsheet S2), suggesting that neurobots could potentially sense and process light stimuli.

Many other clusters included significant enrichment in genes encoding various aspects of the nervous system. This included Cluster 5, which contained various synapsins, tubulins, and microtubule‐associated proteins, which are implicated in biological processes such as the synaptic vesicle cycle, neuron development, regulation of neurotransmitter secretion, and synaptic vesicle localization (Figure S9f, see the NB vs. BB tab in Spreadsheet S2). Cluster 9 revealed the presence of voltage‐gated ion channels including various types of sodium and potassium channels as well as voltage‐gated calcium channels (Figure S9g, see the NB vs. BB tab in Spreadsheet S2). Cluster 11 contained genes encoding various G‐protein coupled receptors, as well as those implicated in the modulation of chemical synaptic transmission (Figure S9h, see the NB vs. BB tab in Spreadsheet S2). Cluster 14 revealed the presence of various hormones and neuropeptides and their receptors (Figure S9i, see the NB vs. BB tab in Spreadsheet S2). Cluster 21 included genes important for axonogenesis and neuron projection development (Figure S9d, see the NB vs. BB tab in Spreadsheet S2).

Moreover, we found that there was significant upregulation in genes encoding ECM constituents (Figure 12a, Cluster 20, and Figure S9j, see the NB vs. BB tab in Spreadsheet S2) including collagen (col17a, col4a2), which is the most abundant fibrous protein and constitutes the main structural element of ECM [46], and fibulins [64], which are glycoproteins that are secreted in the ECM and provide mechanical support in connective tissue (fbln1). Our Second Harmonic Generation imaging indeed showed the presence of some collagen fibers within the central cavity of neurobots, although they were only detected at very low levels (Figure S4). Future experiments are required for assessing the presence of fibulins and other potential constituents of the ECM in this space.

As expected, network embedding and clustering analysis of the upregulated genes in neurobots relative to shams similarly revealed overexpression of genes relating to synapse organization, regulation of neurotransmitter and receptor activity, chemical synaptic transmission (within Clusters 4,5), visual perception (within Cluster 3), neuron projection, and perineuronal nets (within Cluster 2, Figure S10). These findings shed light on biological pathways/molecular functions that are innately present in neurobots in the absence of any external manipulations and will inform future work toward building specialized neurobots through selective enhancement of these pathways.

Finally, we tested the hypothesis that neurobots are expressing a more ancient transcriptome as a result of their nascent evolutionary history. We applied a phylostratigraphic analysis for the differentially expressed genes in the different conditions (Figure 14, NB vs SH and NB vs BB). Interestingly, we found that more than 54% of upregulated genes in neurobots fall into the two categories of most ancient genes (“All living organisms” and “Eukaryota”, Figure 14a). By comparison, very few ancient genes are downregulated. In total 279 genes are downregulated in these two strata for the NB vs BB conditions, and 233 for the NB vs SH condition (Figure 14b), while for the upregulated genes, we obtained 941 and 1109, respectively. Therefore, we conclude that the development of neurobots involves a transcriptomic shift toward very ancient genes for neurobots compared to biobots and shams.

In this study, we built and investigated behavioral, anatomical, and transcriptional properties of novel living constructs with incorporated neural tissue. Using Xenopus laevis embryonic cells, we built two types of living constructs: one using ectodermal cells (biobots) as reported in prior studies [25, 48, 65], and another, novel construct, made using ectodermal and neural precursor cells (neurobots, Figure 1a). We showed that neurobots are viable and self‐motile like their non‐neuronal counterparts (Videos S3 and S4), and that the implanted neural precursor cells indeed differentiate into neurons and extend their processes throughout the construct (Figures 1 and 2; Figure S3).

We found that neurobots became significantly larger than biobots and were significantly more elongated (Figure 1b,c). In addition to overall changes in shape and size revealed by our two‐dimensional analysis, neurobots exhibited differences in the 3D distribution of multiciliated cells and in overall body morphology. In the future, a more advanced 3D analysis incorporating the spatial positions of MCCs will help quantify further differences between neurobots and biobots and elucidate how these may correlate with their behaviors.

Neurobots exhibited a large degree of variability in the amount innervation and in neural architecture (Figure 9; Figure S3). The variability in the size and number of the implanted neural precursor cell clumps (due to the manual nature of their creation), could be a major contributing factor to the variability in neural growth and patterning. Indeed, we found a small positive correlation between the relative size of the implanted tissue and the resulting neurite length and number of neural terminals (see Figure 9c). Another source of variability may stem from the need to pool neural precursors from multiple embryos to generate sufficiently large neuronal clumps for manual implantation. Consequently, the individual cells within each clump were unlikely to be genetically identical and may have differed slightly in their developmental stage. Developing automated methods to implant a defined number of cells, ideally derived from the same embryo, will facilitate more accurate quantification of the relationship between implant size and subsequent growth, and will likely reduce variability across neurobots. Efforts are underway to develop robot scientist platforms for the field of synthetic morpho‐engineering.

We found that the majority of neural processes emanating from the implanted neural precursor clumps grew within the central cavity of the bot. However, some processes clearly extended toward the outer epithelium (arrowheads Figure 2; Figure S3). Future experiments will investigate more thoroughly how neurons modulate the activity of outer epithelial cell types, including multiciliated cells, goblet cells, and serotonergic cells. For example, optogenetic activation of neurons combined with imaging of fluorescent microbead flow patterns could be used to test the causal relationship between neural activity and ciliary beating frequency. When paired with genetic knockdown of relevant cell‐surface receptors, such as serotonergic receptors [29], such approaches could help elucidate the mechanisms by which neurons influence movement in neurobots. It is important to note that, unlike in some invertebrates where ciliary beating frequency, and thus cilia‐driven movement, is under neural control [32], the beating frequency of MCCs on tadpole skin has not been shown to be under neural control. Instead, it is thought to be controlled by molecules secreted by other epithelial cells for example, serotonin (through small secretory cells) [29], and may also be controlled by factors known to influence beating frequency of motile cilia in other mucociliary epithelia such as extracellular ATP molecules secreted in the context of injury or inflammation [66, 67, 68], and noxious sensory stimuli [69]. The emergence of neural regulation of ciliary beating frequency in neurobots would therefore represent a novel property.

Despite variation in the shape and innervation patterns of neurobots, we found that all of them exhibited a central cavity largely devoid of cell bodies (as indicated by the absence of nuclear staining) and containing only a few neurites with sparse expression of the presynaptic marker synapsin‐1 (Figures 2 and 3; Figure S3). We confirmed the largely acellular nature of this space through staining with phalloidin, which labels F‐actin present in all eukaryote cells (Figure S4). We speculated that this region may be filled with extracellular matrix‐like structures. The presence of neurites extending in extremely straight courses points to the presence of such a supporting structure (Figure 2b, red arrowhead). Further, our transcriptomic analysis demonstrated a significant upregulation of genes encoding ECM proteins, including both fibrillar collagens and non‐fibrillar components such as fibulins. In order to assess whether this space contains fibrillar ECM proteins, we used Second Harmonic Generation, which is a label‐free method for visualizing fibrillar ECM proteins such as collagen [45]. Although we found evidence for the presence of collagen fibers, these fibers were extremely sparse (Figure S4). Other non‐fibrillar ECM components, such as proteoglycans, laminins, or fibronectins, may occupy much of this space. It will be important for future work to characterize the molecular composition of this region in greater detail.

Neurobots exhibited a diverse range of movement patterns, and these patterns tended to be more complex than those observed in their non‐neuronal counterparts, indicating that neural expression may affect movement either directly, that is, through neural signaling to the motor effectors, or via changes in the expression patterns of motor effectors. Indeed, we observed a negative correlation between the degree of neural expression and the density of multiciliated cells (Figure 9c). Future modeling efforts could help determine the extent to which differences in the spatial distribution of MCCs alone contribute to the observed increase in movement complexity.

Simultaneous recording of neural activity and behavior could be used to assess the potential neural correlates of the behavior. Indeed, our calcium imaging experiments indicated that the implanted neurons were spontaneously active; however, due to technical difficulties in measuring calcium signals in freely moving neurobots, we were not able to make conclusions about neural correlates of observed movements. Specifically, measuring calcium in freely moving neurobots was complicated by their rotational movements, which often resulted in losing track of specific regions of interest (Video S3). Such 3D rotational movements could be suppressed by making flattened, disk‐like neurobots, however, these neurobots tended not to move as much (data not shown). Future experiments in which the 3D movement of the bots can be better controlled without flattening, as well as experimental setups that allow repeated measurements, will be critical for establishing correlations between neural activity and behavior.

Albino, pigment‐less Xenopus laevis embryos were used to generate the body of neurobots intended for calcium imaging, as the lack of pigmentation facilitates visualization of neuronal calcium activity in the interior. Albino embryos exhibit developmental trajectories similar to those of wild‐type embryos [70], but differences in visually evoked behaviors have been reported in albino tadpoles [71, 72]. Thus, possible differences in the pattern of neural growth and behavioral phenotypes of neurobots with and without body pigmentation, especially in the context of visually‐evoked behaviors, may be an interesting subject of future investigation.

Consistent with a role of neural activity in modulating behavior, we found that treatment of biobots and neurobots with the GABAA receptor antagonist PTZ resulted in significantly different outcomes. To our surprise, we found that most biobots decreased their movement complexity with PTZ treatment, indicating the presence of non‐neuronal drug targets. Indeed GABAergic receptors are found on the surface of the mucus‐secreting goblet cells, and treating Xenopus embryos with bicuculline, which is also a GABAA antagonist, was shown to inhibits mucus secretion [73]. It is possible that treatment with PTZ changes the beating frequency of MCCs through changes in mucus secretion. The effect of PTZ on neurobots was significantly different from that on biobots (Figure 8). In fact, the majority of neurobots showed an increase in movement complexity, suggesting that neural activity may contribute to the observed differential effect. Additional experiments are required to identify the targeted cell types (neuronal or non‐neuronal) and to assess whether the effect we observed is caused by GABA signaling or is due to off‐target effects of PTZ. Although PTZ is a noncompetitive antagonist of GABAA receptors, it may also exert indirect effects on other targets at higher concentrations, including glutamate receptors and voltage‐gated ion channels [74]. Future experiments where the activity of all or a specific population of neurons could be modulated pharmacologically or optogenetically could provide important insight on the causal relationship between neural activity and behavior, as well as the underlying mechanisms.

We raised a group of neurobots in zolmitriptan, which is a selective 5‐hydroxytryptamine (5‐HT) 1B/1D receptor agonist, known to increase the degree of ectopic neural growth in Xenopus embryos [51]. Interestingly, three out of four of these neurobots showed an increase in movement complexity when treated with PTZ (filled circles Figure 8b; Figure S7). Moreover, these neurobots showed a tighter correlation between the amount of implanted tissue and the degree of innervation (Figure S7). Interestingly, 5‐HT1B receptors are shown to modulate GABA release [75], and serotonergic signaling is thought to affect the migration pattern of cortical interneurons [76]. Altered serotonergic signaling in neurobots treated with zolmitriptan may therefore impact the expression of GABAergic neurons during their development. Further experiments are needed to discover the mechanisms underlying this effect and whether zolmitriptan treatment results in biased expression of specific neural subtypes.

Our transcriptomics analysis revealed the landscape of differentially expressed genes between neurobots, biobots, and sham neurobots. Overall, neurobots showed a significant upregulation in gene expression compared to biobots and sham neurobots, whereas biobots and shams were more similar to one another (Figure 10a,b and Figure 11). Functional enrichment analysis revealed that neurobots, compared to both biobots and shams, exhibit a high level of enrichment in genes involved in nervous system development, synapse formation, neuron projection, and trans‐synaptic signaling (Figure 12a,b). Genes encoding major neurotransmitter receptors were present in the transcriptome of neurobots. This included glutamatergic, GABAergic, cholinergic, dopaminergic, serotonergic, and glycinergic receptors.

Our gene network analysis resulted in the identification of multiple functional clusters, allowing us to more deeply examine the genes and pathways that are upregulated in neurobots. Notably, we found a large cluster containing genes with important roles in visual perception (Cluster 1, Figure 13c; Figure S9c). This cluster contained genes normally expressed exclusively in Xenopus eyes, including various members of the opsin family such as a retinal G protein‐coupled receptor, various cone opsins, rhodopsin, as well as genes encoding many other proteins implicated in visual processing. This remarkable finding suggests the possible presence of visually evoked behaviors in neurobots. The next and most exciting step will be to test this hypothesis and discover the ways that light could modulate motor output in neurobots. If present, this will be a completely novel emergent behavior. Follow‐up proteomic analyses will be necessary to determine whether the upregulated transcripts are translated into corresponding proteins and to investigate their spatial organization.

We showed that neurobots exhibit a significant increase in the variability of their gene counts compared to their non‐neuronal counterparts (BBs and SHs, Figure 10c–f). Although some of this variability could be attributable to the manual nature of their creation, and/or pooling of implanted cells from multiple embryos, the excess variability seen in neurobots relative to shams suggests that neurons may play an important role in guiding the way cells explore the gene expression landscape. The nervous systems of animals are known to influence the behavior of non‐neural cells and tissues [60, 77], so it is possible that the information processing activity of the neurons, in response to the unique “life experiences” of individuals, or internally‐generated spontaneous signaling, might account for the neurobots exhibiting the largest inter‐individual gene expression variability of the groups. Moreover, studies show that when cells are exposed to novel stressors that they do not have existing homeostatic mechanisms to resolve, they resort to making random changes in the expression levels of many genes. It is possible that the bots we report here are undergoing stressors that evolution did not prepare them for, and may be employing this kind of exploration of gene expression space [78, 79, 80, 81, 82]. Future single‐cell RNA‐seq experiments, combined with methods for implanting defined numbers of cells, ideally derived from a single embryo, will provide critical insight into cell‐to‐cell variability and determine whether increased transcriptional variability in neurobots is confined to specific cell types.

Finally, based on a phylostratigraphic analysis, we show that the majority of upregulated genes in neurobots consist of the most ancient genes (Figure 14), a pattern that differs significantly from that observed in biobots and shams. In all cases, the cells were wild‐type, and no genomic editing, synthetic biology circuits, scaffolds, or drugs were used. These results suggest that novel configurations of cell types can have large‐scale systemic effects on the transcriptome of the resulting multicellular construct and move it toward the gene expression profiles of the evolutionary past.