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

Building neurobots with predictable nervous system architecture and in large numbers is one of the major remaining challenges of this study. Future automation and standardization efforts will enable higher throughput and consistency, allowing repeated measurements to be performed on neurobots with both identical and distinct nervous system architectures. The impact of pharmacological, optical, and other types of stimulation on the neurobot behavior could be assessed. Additionally, the creation and neuroanatomical characterization of large numbers of identical bots will allow for the discovery of frequently emerging patterns and motifs, thereby shedding light on the potential space for nervous system architectures in novel living constructs whose precise layout has not been shaped by selection for this specific behavioral configuration.

This study establishes a model system and experimental roadmap to increase our understanding of the plasticity of evolutionarily determined hardware of living beings to adapt on developmental (not evolutionary) timescales and to provide interfaces to bioengineered living constructs that may provide novel control capabilities for useful synthetic living machines and shed light on the origins of novelty in the evolution of nervous systems.

All experiments were approved by Tufts University Institutional Animal Care and Use Committee (IACUC) under the protocol number M2023‐18.

Biobots were constructed as described previously [25], by excising tissue from the animal hemisphere of a Nieuwkoop and Faber stage 9 Xenopus laevis embryo (animal cap).

To construct neurobots and sham neurobots we excised 40–50 such animal caps and let them sit with external surface facing up in 60 mm petri dishes filled with a calcium and magnesium free solution (50.3 mM NaCl, 0.7 mM KCl, 9.2 mm Na2HPO4, 0.9 mm KH2PO4, 2.4 mm NaHCO3, 1 mm edetic acid (EDTA), pH 7.3), and coated with 1% agarose made in the same solution. After about 30–40 min the cells were fully dissociated. The dissociated cells were transferred to a deep 60 mm petri dish containing 0.75 Marc's Modified Ringer (MMR) solution, using a P200 pipette, taking as little liquid as possible. For constructing sham neurobots, the dissociated cells were immediately reaggregated and formed into clumps (see below). For constructing neurobots, we dispersed the dissociated cells as far as possible by moving the solution in the dish using a P1000 pipette. The cells were left still in the dish for ∼3–4 h. To reaggregate cells, the dish containing dissociated cells was placed on a shaker, and cells were thereby brought together in the middle of the dish. They were then allowed to reaggregate for approximately 1 h, at which point various clumps of cells formed spontaneously.

Using a P1000 pipette, clumps of ∼>10 cells were moved into the wells of an agarose‐coated 6‐well plate. Larger clumps were broken into smaller ones that could fit inside the body of the bot, that is animal caps excised from a new set of embryos (see Figure 1a). Depending on the size of the clumps, one or more were implanted.

Next, a new set of animal caps was dissociated from a second batch of embryos from a later fertilization (at late blastula, early gastrula stage), and were placed with the external surface facing down individually in the wells of the same 6‐well plate. The excised animal caps slowly formed a bowl and eventually closed up within approximately 10–15 min. Clumps of neural precursor cells (or non‐neuronal clumps in the case of shams) were placed inside this bowl (Figure 1a (ii)) before it was closed using fine forceps, and enough time was allowed for the animal cap to fully close before moving the dish to the incubator. Bots were housed in an incubator set to 14°C, and experiments were performed in a laboratory maintained at approximately 18°C.

Bots were fixed overnight at 4°C in 4% paraformaldehyde (Thermo Fisher Scientific) with 0.25% gluteraldehyde (Electron Microscopy Sciences) individually in 96‐well plates. The next day, they were washed three times at room temperature in PBS‐Triton X(PBT, 0.1% Triton X‐100 (Sigma) in PBS‐/‐ (Gibco)) for at least 15 min and then incubated in the 10% Casblock (invitrogen) dissolved in PBT for at least 1 h. They were then transferred into the solution containing primary antibodies and Hoescht (33342, Thermo Scientific). The plate was sealed using parafilm and covered in foil for light protection and placed on a shaker in the cold room for 3 days at 4°C. The bots were next washed 3 times in PBT at room temperature and then transferred and incubated overnight at 4°C in the secondary antibodies, with the dish sealed with parafilm and covered in foil. Finally, the bots were washed again 3 times in PBT and either stored in PBS at 4°C or mounted into 15‐slide 18‐well flat dishes (81821, Ibidi) in an antifade mounting medium (Vectashield, Vectorlabs H‐1000‐10) for confocal imaging. Primary antibodies used were (Anti‐acetylated tubulin antibody, Mouse monoclonal, Sigma T7451; Anti‐synapsin‐1 antibody, Guinea pig monoclonal recombinant antibody, Sysy Antibodies 106 308; MAP2 antibody, Rabbit polyclonal antibody, Cell Signaling technologies 4542), and the corresponding secondary antibodies were (Goat anti‐mouse Alexa 594, Invitrogen A32742, Goat anti‐guinea pig IgG (H+L) highly cross‐adsorbed secondary antibody, Alexa Fluor 488, ThermoFisher Scientific A‐11073, Goat anti‐rabbit IgG (H+L) cross‐adsorbed secondary antibody, Alexa Fluor 488, ThermoFisher Scientific A‐11008). For F‐actin staining fixed and permeabilized neurobots were incubated with the F‐actin‐specific stain, Alexa Fluor 647 Phalloidin (Thermo Fisher Cat # A22287) at 1:50 in PBS (‐/‐) overnight at 4○C.

Neurobots were washed 3× with PBS (‐/‐) at room temperature and subsequently mounted in Vectashield.

Confocal images were acquired using the Stellaris 8 microscope (Leica Microsystems), equipped with a Chameleon Vision II (Coherent; 80 MHz repetition rate). Z‐stacks were collected at 2 µm intervals with a 25 ×/0.95 NA or a 40×/ 1.1 NA water‐immersion objectives. Each signal was acquired in series. To ensure deep tissue penetration, Alexa 647‐Phalloidin was excited by two‐photon excitation at 819 nm, with emission collected between 636 and 746 nm. Second harmonic generation (SHG) imaging was performed at 890 nm excitation, with emission detected between 465 and 616 nm, as previously described [45]. Fiji (Image J), the LAS X (Leica Microsystems), and Imaris (Oxford Instruments) programs were used to process the images.

Embryos at the four‐cell stage were microinjected in all four blastomeres with mRNA encoding the genetically encoded fluorescent calcium indicator GCaMP6s. These embryos were used for obtaining clumps of neural precursor cells for implantation. Albino embryos were used as the outer shell of the neurobots in these experiments so that the fluorescence signal from neurons could be visualized more easily as wild‐type embryos are pigmented. We used a custom‐built microscope to measure calcium activity in freely moving neurobots (Video S1). We used Fiji's [83] Descriptor Based Series Registration plugin to correct for the motion of the neurobot (Video S2), then used Suite2p software [84] to identify active units (Figure 5). To avoid movements in the Z‐direction, which resulted in changes in the plane of focus, we created flattened neurobots as described previously [48]. Briefly, on the day after their formation, neurobots were pressed down using a glass coverslip, which was gradually lowered over them as small amounts of MMR were removed from the dish. Neurobots were left under pressure for 3 h, after which MMR was gradually added to the dish resulting in the release of the coverslip.

We used the brush tool in Fiji [83] to fill in the shape of the bot and calculated the area and roundness index defined as the major_axis / minor_axis. We used the same tool to estimate the relative amount of implanted neural precursor tissue by dividing the area of the implanted clumps to the outer shell (Figure 1a). We used Imaris software (Oxford Instruments) to quantify neural expression using confocal stacks acquired from the bots that were stained with antibodies against acetylated α‐tubulin, which labeled neurons and cilia in multiciliated cells. We manually traced neurites using the filament function and exported values corresponding to the total length of neurites (dendrite length sum parameter in Imaris) and the number of terminal points (number of dendrite terminal points parameter in Imaris). For the analysis of the multiciliated cell distribution, we estimated the total number of MCCs by marking the center of each MCC using the Spots tool in Imaris and calculated the total number of multiciliated cells. We then used this value to calculate the MCC density by dividing this number by the total area of the bot.

Videos of bot movements were taken over 30 min under various conditions and tracked with the DLTdv digitizing tool [49] in MATLAB, and the x and y coordinates of the center of mass were calculated. A custom Python code was used to extract various kinematic variables using the time series of the coordinates. We calculate total Euclidean distance travelled, average speed, and average acceleration. Additionally, we calculated the percentage of the well that was traversed by the bots by dividing the space of each well into 0.1 mm bins. We then calculated the percent covered area by dividing the number of unique visited bins by the total number of bins. We calculated a complexity index by first calculating the power spectral density (PSD) of the trajectory time series along the x and y coordinates and identifying peaks in the power. We calculated Welch's power spectral density estimate with a window size of 400 s and an overlap of 1 s between windows for each of the x and y time series. We picked a threshold of 10 pixels2/Hz (∼0.17 mm2/Hz = 0.4 mm/Hz) to detect peaks in the PSD of the x and y coordinates. This threshold was chosen empirically to remove the baseline noise corresponding to the tracking of the center of mass of bots that had an average radius of 0.4 mm. We then defined the complexity index as the total number of unique peaks in the x and y PSDs.

To further confirm the significance of power at the peak frequencies identified by Welch power analysis, we used the wavelet transform and surrogate data. For each time series, we used MATLAB to calculate the continuous wavelet transform and from there calculated the wavelet power spectrum across all time. To assess the significance of dominant frequencies found from the Welch spectral analysis, we generated surrogate time series that preserve the original signal's amplitude spectrum but randomize its phase information in the wavelet space. These surrogates were used to compute a distribution of wavelet power spectra, from which significance thresholds (95th percentile) were derived for hypothesis testing. Only 4 out of 243 peaks identified through Welch analysis across all bots did not pass the significance test in the wavelet space. This small difference is expected because the two methods differ in their frequency resolution: Welch analysis uses uniform resolution across all frequencies, whereas wavelet analysis employs logarithmic spacing, providing finer detail at lower frequencies and coarser resolution at higher frequencies.

We submitted 12 samples (4 samples per biobot type, 5–15 bots per sample: Number of bots per sample: {NB1 = 14, NB2 = 12, NB3 = 12, NB4 = 12}, {BB1 = 16, BB2 = 6, BB3 = 6, BB4 = 5}, {NN1 = 5, NN2 = 6, NN3 = 6, NN4 = 5}) submerged in Trizol (Invitrogen) in 2 mL Eppendorf tubes, to Novogene (Novogene Corporation Inc., Sacramento, CA) for low‐input, high‐lipid, bulk RNA extraction. The first batch submitted to Novogene (NB1 and BB1) contained a larger number of bots to test the efficacy of the low‐input method. The next series of sample batches (NB2‐4, BB2‐4, NN1‐4) contained a smaller number of bots based on the results from the first batch. Equal quantities of RNA were sequenced from each sample using the NovaSeq6000 sequencer, resulting in consistent library size across samples. Clean reads were extracted from FASTQ files, removing reads with adapter contamination, when uncertain nucleotides constitute more than 10 percent of either read (N > 10%), and when low‐quality nucleotides (Base Quality less than 5) constitute more than 50 percent of the read. The index to the reference genome (Xenopus laevis version 10.1) was built using Hisat2 v2.0.5 [85], and clean reads were aligned to the reference. The mapped reads of each sample were assembled using StringTie (v1.3.3b) [86], and FeatureCounts v1.5.0‐p3 [87] was used to count the read numbers mapped to each gene. FPKM of each gene was calculated based on the length of the gene and read counts mapped to this gene. Differential expression analysis was performed using the DESeq2 R package (1.20.0) [88], and the resulting p‐values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate.

The webapp g: Profiler [62] was used to perform functional enrichment analysis of differentially expressed genes. For each comparison, upregulated genes (p‐adjusted < 0.05; log2foldchange > 4) and downregulated genes (p‐adjusted < 0.05; log2foldchange < ‐4) were separately mapped from Xenopus to human symbols using the HGNC Comparison of Orthology Predictions (HCOP) tool [89]. Genes lacking an established gene symbol were removed from analysis. Each gene list was separately queried using g: Profiler across all data sources. The statistical data scope included only annotated genes, and the g: SCS method was used for computing multiple testing corrections for p‐values at a threshold of p < 0.05. The R package ggplot2 (v3.5.1) [90] was used to generate dot plots of gene ontology driver terms from g: Profiler. Driver terms were determined by grouping significant terms into sub‐ontologies based on their relations, then identifying the leading gene sets that give rise to other significant functions in the ontology neighborhood.

For network analysis and clustering, we applied network analysis techniques to discover biological functional modules [91, 92]. By integrating gene expression and interaction data, we extracted PPI for the different biobots in their respective conditions and applied network embedding and clustering techniques similar to those described by Cantini et al. [93] and Pio‐Lopez et al. [94]. Specifically, we used the MNMF algorithm developed by Wang et al. [61] for network embedding and clustering. To create a network for the bots, we started by isolating genes of interest and identifying corresponding human orthologs under various conditions using the HCOP database [95]. We then used the STRING database [96] to extract relevant PPI networks. The clusters identified through our network embedding and clustering method were further analyzed for enrichment using g: Profiler [62].

The normalized gene count variability was compared between groups (BBs, NBs, and SHs) using a MATLAB script and the method summarized in Figure S11. We excluded NB1 and BB1 since they were sent to Novogene for sequencing as a different batch than other NBs (NB2‐4) and BBs (BB2‐4), and we wanted to reduce the impact of inter‐batch variability (see RNA‐sequencing and bioinformatics). We kept all four NN groups as they were sent as a part of the same batch (NN1‐4). For each gene in each pair of groups being compared, the mean count value across all pools of both groups was found, and genes were ranked from greatest to least mean. Genes for which any of the counts across all pools was 0 were discarded. Because the count value of a given gene in a given pool represents the mean value of all the individuals in that pool, the standard deviation of the pools gives the standard error of the means (SE) of the group. The SE is related to the number of individuals per pool (n) and the standard deviation of the individuals within the pools (σ) with the equation SE = σ/sqrt(n). By multiplying the SE by sqrt(n), σ can be calculated. Dividing σ by the mean count value of the pools gives the coefficient of variation (CV). The ranked gene CV lists were then split into 100 bins (percentiles) containing equal numbers of genes from highest to lowest counts. Within each bin, the number of genes with greater CV for the first group than the second was counted and divided by bin size to find the fraction of genes in the bin with greater CV in the first group. These fractions were then plotted as bar graphs in Figure 10c–e, with a blue line marking 0.5. The bin values appeared to vary with gene count percentile, so to determine the statistical significance of each bin's departure from its expected value, a permutation test was used.

For each plot (each pair of groups), the order of the gene pairs was randomly shuffled (keeping pairs together), and new bins were generated. This was repeated 1000 times for random shuffles to produce a distribution of bin fraction values for each bin. The p‐value of each bin was defined as the proportion of the bin fractions from the distribution that were further in absolute value from the distribution mean than the true bin fraction. Bins with p‐values of p < 0.05 were deemed statistically significant and were colored dark blue. in Figure 10d‐f. All other bins were colored light blue. The mean bin value from the distributions was plotted as a yellow line.

We employed the phylostrat package [97] to conduct a phylostratigraphic analysis on transcripts of all three bot types, with Xenopus laevis (taxon ID 8355) designated as the reference species. This software automated several key steps in evolutionary analysis: (1) it built a clade tree using species from the UniProt database and aligned it with the latest NCBI taxonomy; (2) the clade tree was trimmed to maintain a phylogenetically diverse selection of representative species for each phylostratum; (3) a comprehensive protein sequence database was constructed from hundreds of species based on this clade tree, with additional data such as human and yeast proteomes manually added, resulting in 329 species for our study; (4) a similarity search was performed by conducting pairwise BLAST comparisons between the proteins encoded by Xenopus laevis and those of the target species; (5) the best hits were identified, and gene homology was inferred between Xenopus laevis and the target species; (6) each gene was assigned to a phylostratum that corresponded to the oldest clade for which a homolog was identified. Genes specific to Xenopus laevis were classified as orphan genes and placed within the Xenopus laevis phylostratum. The evolutionary stages we focused on include: All living organisms (bacteria, eubacteria), Eukaryota, Opisthokonta, Metazoa, Eumetazoa, Bilateria, Deuterostomia, Chordata, Vertebrata, Gnathostomata, Euteleostomi, Sarcopterygii, Tetrapoda, Anura, Xenopus, and Xenopus laevis. This methodological approach enabled a detailed examination of gene emergence and their evolutionary trajectories across various taxa. By implementing Phylostratr, we systematically mapped the age of the bot genes in the different conditions with a specific phylostrata to understand the distribution of ages of the bots’ overexpressed genes.

We used the upregulated and downregulated genes in neurobots (log2foldchange > 4 and log2foldchange < ‐2 respectively).

Statistical analyses of behavioral and shape data were performed using the non‐parametric Kruskal–Wallis test. When appropriate, multiple comparisons conducted using Tukey's honestly significant difference test. For PTZ experiments, we used a one‐sample t‐test to assess whether the population mean of the relative complexity index was significantly different from zero. Robust linear regression was used to quantify relationships between pairs of neurobots’ physical, neuroanatomical, and behavioral measures; the significance of regression coefficients and Pearson correlation coefficients was assessed using two‐tailed t‐statistics. Significance of peak frequencies identified by Welch analysis was validated using continuous wavelet transforms and phase‐randomized surrogate time series to generate null distributions of wavelet power, with significance assessed at the 95th percentile. In box plots, whiskers showed the non‐outlier extent, + signs depicted outliers, and the top and bottom of the box showed the upper and lower quartiles of the data. The horizontal bar inside the box showed the median. MATLAB (Mathworks, Natick, MA) functions were used for all statistical analyses.

The statistical analysis for differential gene expression data was performed using the R package DESeq2, which employed the two‐sided Wald test. Multiple hypothesis testing corrections were used to obtain adjusted p‐values. For the analysis of gene expression variability, the statistical significance of bin deviations from expected values was assessed using a permutation test. For each group comparison, gene pairs were randomly shuffled 1000 times (with pairs kept intact) to generate null distributions of bin fractions; bin p‐values were defined as the proportion of permuted values more extreme than the observed value, and bins with p < 0.05 were considered significant.

This research was supported by HR0011‐18‐2‐0022, W911NF1920027, awarded by the Department of Defense, and grants from the John Templeton Foundation and Northpond Ventures.

M.L. is a scientific co‐founder and consults for Fauna Systems, a company seeking to commercialize frog cell‐based biobot technology.

We would like to thank Drs. Douglas Blackiston, Patrick McMillen, Jack Lindsay, Pai Vaibhav, Harini Rajendiran, and Russel Gould for their technical help and scientific insights, Meghan Short for statistical analysis consultation, Susan Marquez, Thomas Ferrante, Ramses Martinez, and Kostyantyn Shcherbina for their help with microscopy and engineering, Jeantine Lunshof for philosophical and ethics discussions, and Drs. Donald Ingber, Florian Engert, and Michael Super for their input and support throughout the project. Thanks to Gordon Allen for copyediting, and Julia Poirier and Tomika Gotch for assistance with the manuscript. This research was supported by HR0011‐18‐2‐0022, W911NF1920027 awarded by the Department of Defense, and grants from the John Templeton Foundation (grant #62212) and Northpond Ventures.

RNA‐sequencing data generated during this study are available in the NCBI Gene Expression Omnibus (GEO) public repository under accession number GSE295614. All other data supporting the findings of this study are available from the corresponding authors upon reasonable request.

RNA‐sequencing data generated during this study are available in the NCBI Gene Expression Omnibus (GEO) public repository under accession number GSE295614. All other data supporting the findings of this study are available from the corresponding authors upon reasonable request.