Motile Living Biobots Self-Construct from Adult Human Somatic Progenitor Seed Cells
Anthrobots were collected in Pluristrainer Mini's with a 40‐micron pore size (Fisher Scientific #431 004 050) and fixed with 4% paraformaldehyde at room temperature for 30 min. Following phosphate buffered saline (PBS) washes, blocking and permeabilization were performed for 1 h at room temperature on a rocker in a blocking buffer consisting of phosphate‐buffered saline with 10% normal goat serum, 1% bovine serum albumin (BSA), and .15% triton x‐100. Anthrobots were then incubated with mouse anti‐acetylated tubulin (Sigma‐Aldrich #T7451) primary antibodies at 1:250 dilution factor in blocking buffer for 24 h at 4 °C on a rocker. The primary antibodies were labeled with Alexa Fluor 647 donkey anti‐mouse (Thermo Fisher Scientific #A31571) secondary antibodies, at 1:500 dilutions in blocking buffer, for 1 h at room temperature on a rocker. Lastly, Anthrobots were incubated with Alexa Fluor 594‐conjugated mouse anti‐ZO‐1 (Thermo Fisher Scientific #339 194) at a 1:100 dilution in blocking buffer for 24 h at 4 °C on a rocker. Anthrobots were mounted on glass‐bottom 96‐well plates in ProLong Glass Antifade Mountant with NucBlue (Thermo Fisher Scientific #P36981). Neuronal tissues were fixed, blocked and stained using the same protocol, except by using Beta III Tubulin (Tuj1) (Abcam #ab18207) as the primary antibody for staining the hiNSCs. Images were collected using a Leica SP8 Fluorescence Lifetime Imaging Microscopy (FLIM) with a 25x water immersion objective. Z‐stack step size = 3 micron unless otherwise specified.
To find whether there were any unique morphological types like was the case with movement types, each spheroid was processed through a custom‐made analysis pipeline (code attached). First, a 3D model of the bot was created in R where the points corresponding to the body and the cilia were clearly indicated. For this, first the cilia channel was isolated from the Laser Induced Fluorescence (LIF) images of the bots; these cilia‐only images were then run through CiliaQ[ 39 ] using the RenyiEntropy algorithm for detection. These binarized cilia were then imported into the code, along with the points that comprised the body. These “body points” were extracted using the “body channel” of the LIF images by first running the pixels through a logistic transform then thresholding the pixels based on the signal to noise ratio, calculated by using a median filter and comparing the points before (“signal + noise”) to after (“only signal”).
Due to the large volume of body points, to reduce the points to a manageable amount, we first found the outlines of each slice of the spheroid by using a concave hull using the Concaveman package (version 1.1.0).[ 40 ] The cilia points were then projected on the nearest body points by simply choosing the nearest one by distance to get the shadow of the cilia on the body.
The structural index variable Cilia Points was calculated by counting the number of unique projected points on the body. Then, the dbscan package (version 1.1‐10)[ 41 ] function in R was used to find the clusters of cilia. The number of points that fell outside of clusters with this definition were defined as Noise Points.
Following this step, we computed the spanning ellipsoid of the body points by using the “ellipsoidhull” function from the cluster package (version 2.1.3).[ 42 ] The Max Radius variable was calculated directly by the function, and Aspect was defined as the ratio of the largest radius to the shortest radius, all quantities computed by the function. Finally, we used the “ashape3D” function from the alphashape3d package (version 1.3.1)[ 43 ] to generate a 3D alpha hull of the body points, and used the mesh to get the surface area of each spheroid.
Cilia Points/Area was defined as the Cilia Points variable divided by the calculated surface area. Similarly, the Shape Smoothness was defined as the ratio of the volume of the 3D alpha hull to the volume of the spanning ellipsoid. Finally, we found the center of the bot by finding the sum of the centroids of each triangle that makes up the alpha hull weighted by area of the triangle. Polarity was defined as the norm of the vectors from the center to each cilia point divided by the sum of the norm of each vector. The Cilia Distribution Homogeneity was defined as 1 – D statistic of the two sample Kolmogorov‐Smirnov test, where sample A is the 1st nearest neighbor (1NN) distances for the cilia, and sample B the 1NN distances if the same number of cilia points were distributed close to uniformly but randomly across the surface of the bot. After these analyses were carried out, we ran the dataset through a Principal Components Analysis with centering and scaling. Afterwards, a hierarchical clustering was carried out on the resultant dataset with the Ward.D2 method and the resulting classification was plotted as above. In total, 350 bots were put through the pipeline and into the following PCA and included Movers, Nonmovers, Linear and Circulars. Further details could be seen in the code.
To get a confidence interval for the absolute value of the loadings, we bootstrapped the loading value by sampling 250 bots from the 350 that we have, 10 times. We then took the loadings value for the 1st and the 2nd primary component for all 8 variables and calculated the mean and 95% Confidence Interval for the loadings for the PC in question. If there was overlap between the CI of the loadings, they were assigned the same rank, otherwise they were assigned different ranks. Ranks were relative to the “highest” contributor of the rank; i.e, for PC1, Shape Smoothness had overlap with Max Radius, but it also had overlap with Cilia Distribution Homogeneity. However, Cilia Distribution Homogeneity did not overlap with Max Radius. Max Radius had the highest upper limit of the CI among the 3 variables in question, and thus Shape Smoothness was co‐ranked #1 along with Max Diameter, but Cilia Distribution Homogeneity was not.
After finding trends among behavioral and morphological data, we decided to see if there was any potential overlap between the two. To eke out any possible correlation, we first chose to use categories of spheroids behaviorally orthogonal to movers, the non‐movers. The goal was to observe whether there is any overlap between the morphology of movers and non‐movers. Similarly, we had four potential behavioral types that could overlap with our morphological clusters. Eclectics could not be included in the analysis since they are an aggregate of multiple inconsistent patterns and highly uncommitted to their behavior (Figure 2G), and thus cannot be used in a bot‐level analysis (instead of period). Circulars and Linears on the other hand were highly committed behavioral types that were orthogonal to each other (had little to no interconversion on the Markov plot) and prototypes of two extreme movement types with high variability between them. Consequently, the Curvilinear behavioral subtype was also not included since it lacked orthogonality with both Circulars and Linears due to shared traits between them. The morphological indices for each spheroid were calculated as outlined in the methods for the previous sections, and then clustered with Circulars, Linears, Nonmovers and Movers together. To measure the significance of the overlap, if any, between clusters, we decided to use a Fisher test to compute whether the proportion of a certain behavior per cluster type was different from the others. We ran the test twice, once to see if there were any significant differences in number of nonmovers per cluster, and once to compute the difference in the ratio between circulars and linear per cluster.
It showed that the proportion of nonmovers in Cluster 1 versus Clusters 2 and 3 were significantly different with an average p = 2.6*10e‐6 and 3.5*10e‐8 respectively Cluster 2 and 3 also had a statistically significant difference in number of nonmovers (p = 0.01) which is understandable, since Cluster 2 had no non‐movers. For Circular/Straight, Cluster 2 versus 3 were significantly different with p = 0.00011, and Cluster 1 had no Circulars nor Straights.
The generated tracks were analyzed alongside Z‐stacks of designated Anthrobot from a confocal microscope to see if their morphology was connected to their movement. ImageJ was used to compile the slices of the Anthrobot so that a 3D model could be generated and rotated to render a transformation that visually matches a random frame of the Anthrobot from the timelapse. This random selection could be done as the Anthrobot, despite moving around, did not tend to roll and therefore generally maintained the same orientation throughout a timelapse. Additionally, they often moved with a specific side that always faced forward that was designated as a heading. To see if biases in cilia patterns to one side or lack thereof on an Anthrobot affected its movement this heading would serve as the axis along which a plane of symmetry would be extended to bisect the bot. This plane would be defined by 3 points on the bot along this axis, one placed at the centermost point of the axis within the bot, one on the part of the bot that most visually served as the heading in the video and one on the opposite point of the bot from the heading.
After realizing that polarity could play a key role in determining movement type, we decided to see whether the symmetry across the movement axis any trends had compared to the other axes. To calculate this, we followed the procedure used in both Figure 3 and Figure 4 to get representations of Cilia on the body of the bot, then project these representations onto the plane of symmetry defined by the points obtained in the Motility orientation alignment section. The side of the plane (movement axis) each cilia point belonged to was noted using the sign of the dot product of the normal of the plane and the vector to the cilia point. Finally, to better distinguish whether the cilia distribution played a role in movement type (linear vs. circular) we created the Bilateral Symmetry index. This index was modified from the Chamfer distance, and was calculated as the sum of the median/ mean of the distances between all points in set A and the closest point in set B and the median/ mean of the distances between all points in set B and the closest point in set A. To calculate the index, the cilia points were projected onto the plane defined by the three points in the section below. Set A and Set B then became the points projected from one or the other side, respectively, after which the modified Chamfer index was calculated for the two sets using the createTree() function of the SearchTrees package. The statistics used to calculate the asymmetricity between both sides were the difference in points between the two hemispheres, the difference in points/ total cilia points, the median and the mean modified Chamfer distance. In the end, they were visualized and clustered using a PCA to see trends (see code).
In the case of Figure 4E, instead of calculating the asymmetry statistics after getting the equation of the plane, we then used the Rodrigues’ rotation formula to rotate the normal (and thus the plane) with the fixed intersection being the center of the bot. Rotations of 45, 90 and 135 degrees were used yielding 4 axes (in the form of plane equations) including the movement axis. To eliminate the z‐axis, we used a PCA to get the rotation matrix to convert the projected cilia points from 3D to 2D and calculated the Chamfer distance using the formula described at https://github.com/UM‐ARM‐Lab/Chamfer‐Distance‐API , except that we did not square point distances. This statistic was calculated along for the cilia of all linear and circulars and put into a paired Wilcoxon rank‐sum test with an alternative hypothesis of “greater” and “less” for circulars and linear respectively to see if the movement axis was “more asymmetrical” or “less asymmetrical” respectively. For the body we did the same procedure, with the exception that our statistic now involved finding the distance of the body points from the center of the bot (once again segregated into two hemispheres with the dot product). Then, we used the KS test to calculate a D‐statistic which had greater values the more dissimilar the two distance distributions for both hemispheres were. Just like the cilia we then used a paired Wilcoxon rank‐sum test with an alternative hypothesis of “greater” and “less” for circulars and linear respectively to see if the movement axis was “more asymmetrical” or “less asymmetrical” respectively.
We followed a previously established protocol for creating the neuronal cultures[ 28 ] which is summarized as follows. A 150 cm dish was first coated with 0.1% gelatin for 20 min and then aspirated off before seeding mouse embryonic fibroblasts (ATCC #SCRC‐1008) in mouse embryonic fibroblast (MEF) growth media (89% DMEM GlutaMAX, 10% Fetal Bovine Serum (FBS), and 1% Anti‐anti). Once the MEFs were confluent, they were inactivated by adding 20 mL of MEF growth media containing 500 µL of 10 µg mL−1 mitomycin C (Sigma #M4287) and incubating for 2 and 3 h at 37 °C. After incubation, the MEF growth media + mitomycin C media was replaced with hiNSCs at a density of 1/10 of a confluent target vessel in 25 mL of hiNSC growth media (77.6% Knockout DMEM, 20.20% KnockOut Serum Replacement (KOSR), 1% GlutaMAX, 1% Anti‐anti, 0.18% 2‐mercaptoethanol with 0.1% of 20 ng mL−1 bFGF). The day after seeding the hiNSCs required a media change where all the old media was aspirated off, and 25 mL of fresh hiNSCs growth media was added. Media changes were performed every other day until the hiNSCs were 80%–85% confluent. 3 h before performing the differentiation, the destination vessels were first coated with .1 mg mL−1 poly‐d‐lysine (PDL) (enough to coat the bottom of the wells) for 1 h at room temp, and then the PDL was aspirated before adding in 10 ug mL−1 laminin in DPBS (enough to coat the bottom) for 2 h at 37 °C. In the differentiation, the hiNSCs first went through one DPBS wash before adding TrypLE Select for 3–5 min to detach the cells from the plate. The cells were then collected and spun down for 3 min at 500 g then resuspended in neurobasal differentiation media (96% Neurobasal Media, 2% B‐27 supplement, 1% GlutaMAX, 1% Anti‐anti). The hiNSCs were seeded at a concentration of 100 000 cells cm−2.
Once the hiNSCs were in differentiation, there was a media change the day preceding their differentiation and then every other day from there.
To better explain the relationship between bot trajectory and movement in certain environments, we tried to relate the scratch edge to the actual movement of the bot. The steps taken before analysis involved i) creating a background image prototype from the.czi recording, ii) verifying the quality of the background image and saving as .png file, iii) tracking the bot, iv) extracting the coordinates of the scratch, and v) checking whether tracking was correctly carried out by generating a video with a beacon on the bot.
This procedure was carried out on 30+ files and yielded 20 usable datasets, which were whittled down to 17 after a manual check of tracking quality and excluding videos where bots never touched the scratch wall. Using the coordinates of the scratch walls and the tracking of the bot, we used the Rbioformats (version 0.0.74)[ 44 ] and Revision (version 0.6.2)[ 45 ] package tools to test 1) whether bots were more in contact with the scratch when they have a higher rotational tendency and 2) when moving on tissue, whether faster bots tend to cover more area, i.e., explore better. The lm() function was used to model the data after calculation of proportion of bot on tissue, instantaneous angular velocity and linear speed as variables. Before each model was approved, diagnostics were run on the model using the DHARMa package which included analysis of the residuals.
Afterwards, in order to take a better look at the nuances of interactions between bots and scratches, we decided to limit the data and remove any videos that had bots with very low rotational tendency (<0.33) since they were not stable enough in their rotational behavior. Bots with very high rotational tendency (>0.7) were removed since they were prone to skidding instead of interacting with the scratch walls. Finally, bots whose tracking videos were not optimal i.e., they frequently went backward or circularly in the scratch were removed since they did not have consistent forward movement that could be correlated with the scratch wall. After all these removals, our dataset ended up with 13 examples of scratch‐bot interactions which could be effectively analyzed. The Gyration was simply the Rotational Tendency values renamed. The Scratch‐Trajectory Similarity metric was calculated as the larger absolute value of the correlation between the heading angle of the trajectory and the heading angles of the scratch from the surface perpendicular to the bot. These correlation values were then modelled using the lm() function with an expectation of a quadratic relationship for Gyration. The specifics can be seen in the attached code.
Traversal videos of bots moving along a scratch within a neuron plate were processed via Adobe Illustrator to see if the bot faithfully followed the edge of the scratch. The first method of processing aligned the center of the scratch at a horizontal line parallel to the bottom of the screen and placed a point on the center of the bot at each frame of the video as well as straight above and below this point on the edges of the scratch. Lines were made to connect each respective type of point for both of the edges of the scratch and the position of the bot. The output of this for further analysis was a set of coordinates of the end of each line derived from rendering these series of lines as an vector file and exported as text.
To investigate whether these “bridges” were actually akin to neurons, we decided to analyze the pixel densities of various areas on and surrounding the bridge. In order to prevent confusions regarding this process with regard to intensity of color, we binarized the image on ImageJ. If the automatic thresholding did not visually appear similar to the raw image, we adjusted the threshold manually. We ended up with six areas of interest: the neurons above the bridge, below the bridge, to the left but adjacent to the left but far, and to the right, both adjacent and far. These areas were defined relative to a FIJI ROI box on the neuronal bridge which tried to encompass the width of the bridge and the height close to the narrowest point of the scratch channel that we would interact with. A line of one bridge length or lower if the image size required smaller lines to fit the boxes was used in the vertical and horizontal directions (called hereafter as “bridge length”). The above and below bridge measurements were taken by placing the bounding box one vertical bridge length from the box on the neuronal ridge. The adjacent areas on both sides were defined as 1 horizontal bridge length away from the bridge in the scratch. The far areas were 1 bridge length beyond the adjacent areas. The far and adjacent boxes were (vertically) adjusted so they overlaid the scratch as much as possible (Figure S9, Supporting Information). Finally, we used Analyze>Histogram in FIJI to get the size of the box (which was constant) and the number of pixels of scratch tissue (in black) and calculated the proportion. We then used an unpaired two sample T‐test with unequal standard deviation to calculate the significance of the difference, if any.
For all analyses in the paper, the p‐value to symbol correspondence was that a range of 0 to 0.0001 corresponded to ****, 0.0001 to 0.001 corresponded to ***, 0.001 to 0.01 corresponded to **, 0.01 to 0.05 corresponded to * and 0.05 to 1 corresponded to ns. Additionally, all significance tests were evaluated at an alpha value of 0.05. Unless otherwise specified, the alternative hypothesis was always two‐sided for t‐tests. For all statistical analyses listed below we used the Rstudio computational/ statistical software.
For Figure 2, we analyzed tracks from 197 bots for 5 h, collected across 47 timelapse videos (each video featuring 4–5 bots). In the pre‐processing step, we omitted data that is within one bot length (≈100 um) from the edge of the vessel to prevent edge effect as a confounding factor, which yielded a final of 42 235 individual 30 s periods. After cross‐entropy clustering these periods, we used a t‐test to analyze cluster‐specific differences in active periods (with cluster 1 having 6004 periods, cluster 2 with 6700, cluster 3 with 3436 and cluster 4 with 2384), which were further analyzed as shown on the figure. There were 23 711 inactive periods that were excluded from this downstream analysis.
For Figure 3, the pre‐processing step involved binarizing cilia versus body masses of 350 bots (each represented in 3D via confocal Z‐stack images) through CiliaQ[ 38 ] as described in the methods above. The data obtained from these 350 bots were further clustered into 3 groups with sizes of 125, 24 and 201 for clusters 1,2 and 3 respectively. To check which of the 8 variables that were used to compute the PCA were significant for each cluster, we ran a two‐sided, two‐sample t‐test on all pairs of clusters, for all 8 variables.
For Figure 4, we used the same full set of data for 350 bots as Figure 3, which were still clustered into 3 groups with sizes of 125, 24 and 201 for clusters 1,2 and 3 respectively. Of these 350, it had specific information on type of movement for 28 bots (“displacers” – circulars and linears). Thus, of these 350, we then focused on 28 bots, 15 circular and 13 straight. These were analyzed for various metrics of asymmetry (difference between the two hemispheres in number of cilia, or chamfer distance of the hemispheres of cilia etc.) along the movement axis and the axis that was 90 degrees offset from the movement axis as described in methods above. We then measured the change in the chamfer distance asymmetry statistic from the 90‐degree offset to the movement axis (asymmetry difference = asymmetry ≈90 degree offset – asymmetry around the movement axis) for circulars (n = 15) and linears (n = 13) with a two‐sided one sample t‐test for each with a significant result for circulars (p = 0.0482) but not linears (p = 0.1116).
For Figure 5 during preprocessing, we excluded videos where the bot never touched the scratch wall as described in the methods above. A t‐test for the slope was run on the relationship between bots’ rotational tendency and proportion of bot on tissue which yielded a significant (p = 0.017, slope 1.15, n = 17) result. Similarly, when we compared the relationship between bot linear speed and the proportion of bot on tissue using a t‐test again, we received a significant result (p = 0.031, slope 0.0082, n = 17). For a subset of these 17 bots (dataset constrained to non‐stalling bots with rotational tendencies between 0.33 and 0.7 and viable tracking videos as described in methods above), we initially tried to fit a linear model of the relationship between bot rotational tendency and scratch‐trajectory similarity metric. The residuals of this analysis were not centered around a mean of 0 but rather followed a visibly quadratic trend (Figure S9, Supporting Information). This suggested a quadratic model would be a better fit for the relationship between the two. We ran the t‐test for the significance of this quadratic relationship that was significant (p = 0.006, n = 13).
For Figure 6, we ran a two‐sample t‐test pairwise with each category to characterize the pixel density of these various areas (gap closure site, native tissue, sites adjacent and distal to the gap closure sites) for all bridges where the connectivity to both sides of the scratch was maintained throughout the 3‐day experiment, which was 50% of the total N = 10. Difference between gap closure site and native tissue is insignificant (p = 0.37), while the difference between the gap closure site and both adjacent and distal scratch sites are significant (w/ p = 0.006 and p = 0.005, respectively).
This work is partially funded by a sponsored research agreement between Tufts University and a company called Astonishing Labs; co‐author Levin is a scientific co‐founder of Astonishing Labs.
P.S., and B.G.C. contributed equally to this work. G.G., and M.L. performed experimental design and data interpretation. G.G., P.S., B.C., H.L., B.S., S.G. performed data acquisition and analysis. All co‐authors Wrote the manuscript. M.L. performed funding acquisition.
The authors thank Santosh Manicka for helpful technical discussion of Multivariate Classification, Roger Kamm and David Kaplan for in vitro bot inoculation discussion, Benjamin Sweigart for help with statistical analysis of change in motility as a function of varying starting seeding conditions, and Nik Davey for cell culture maintenance and imaging help with the neuronal scratch negative control. The authors also thank Doug Blackiston, Joshua Bongard, Caitlin Grasso and Sam Kriegman for high‐level scientific discussions on the biobots field, Cindy Zhu, Zoe Weiner, and Serena Meng for image acquisition automation and processing help, and Julia Poirier for helpful comments on the manuscript. M.L. gratefully acknowledges support via grant 62212 from the John Templeton Foundation, and a Sponsored Research Agreement from Astonishing Labs. The confocal imaging portion of this work was supported by the NIH Research Infrastructure grant NIH S10 OD021624.
Gumuskaya G., Srivastava P., Cooper B. G., Lesser H., Semegran B., Garnier S., Levin M., Motile Living Biobots Self‐Construct from Adult Human Somatic Progenitor Seed Cells. Adv. Sci. 2024, 11, 2303575. 10.1002/advs.202303575
The data that support the findings of this study are available from the corresponding author upon reasonable request.