Revealing non-trivial information structures in aneural biological tissues via functional connectivity

This study was designed and performed under oversight from the Tufts University Animal Care and Use Committee (IACUC). All experimental protocols involving amphibians were reviewed and approved by the IACUC prior to the work beginning, and were certified under protocol number M2020-35 in compliance with institutional, state, and federal ethical standards for animal welfare.

All experiments were conducted using tissue sourced from the amphibian Xenopus laevis. Wild type embryos were collected 30 minutes post-fertilization and raised in 0.1x Marc’s Modified Ringer’s solution (MMR), pH 7.8, until microinjection at the 4-cell stage and animal cap excision at Nieuwkoop and Faber stage 9 [82].

Microinjection of synthetic mRNA was performed at the 4-cell stage using a pulled glass capillary, with each of the 4 cells being injected to ensure ubiquitous expression across the embryo. Synthetic mRNA was synthesized from a linear DNA template using commercially available kits (Life Technologies), which was stored at –80 ∘C until used. Directly prior to injection, cohorts of healthy wild type embryos were transferred to a laser etched petri dish containing 3% Ficoll solution. The 4 individual cells of each embryo were then injected with a pulled glass capillary, delivering approximately 500 ng of mRNA in 50nL of volume to each cell. After healing for 1 hour, the embryos were washed twice in 0.1x MMR, pH 7.8, to remove the Ficoll solution, and any damaged embryos were discarded before moving the dish to a 14 ∘C incubator. Two mRNA’s were co-injected in the reported work; GCaMP6s, a reporter of calcium activity [83,84], and the intracellular domain of Notch (Notch ICD), which is known to inhibit multiciliated cell induction in developing frog epidermis [53,85,86]. Multiciliated cells were molecularly inhibited in the current study as the presence of these motile structures causes the mucociliary organoid to move during observation, complicating image analysis [87–89].

At Nieuwkoop and Faber stage 9, the animal cap of each embryo was removed to generate epidermal organoids. Cohorts of injected embryos were transferred to a Petri dish containing 0.75x MMR, lined with 1% agarose to reduce cell/tissue adherence. Using a pair of sharpened microsurgery forceps, the vitelline membrane of each embryo was removed, and the animal cap (the central portion of the pigmented top of each embryo) was surgically excised and inverted in the dish. These explants are known to develop into irregular epidermis if untreated [41–43,88,90]. Following excision, the remainder of the embryos were discarded, and the tissue was allowed to heal into a spheroid over the course of 3 hours at room temperature. Following healing, the developing tissue moved to new dishes containing 0.75x MMR and 5 ng/μl gentamicin, lined with 1% agarose, and placed back at 14 ∘C. After an additional 24 hours of development, the animal caps were placed under a glass cover slip for 3 hours at room temperature, generating continuous compression, which resulted in a permanent flattened tissue which improved optical measurements. Following compression, the explants were kept at 14 ∘C for 5-6 further days of development until imaging, at which point the tissue had differentiated into a modified epithelial organoid.

All calcium imaging was performed on an Olympus BX-61 microscope equipped with a Photometrics CoolSNAP DYNO CCD camera and CoolLED pE-300 light source. Individual organoids were placed in a depression slide containing 0.75x MMR under a 4x objective. Images were captured using a FITC filter at a rate of 1 frame every 5 seconds, across a total 20 minutes of observation. Capture rate was determined by pilot studies which identified the minimum time scale to record calcium flashes in individual cells, while also minimizing exposure to illumination to avoid photobleaching and/or phototoxicity. For the first 20 minutes, basal rates of calcium activity were recorded. After 20 minutes the image capture was paused, and a pulled glass needle with a tapered tip diameter of 10-15μm was used to place a puncture near the center of the organoid. Tip diameter was chosen to minimize overall damage to the organoid, and the depth of the wound traversed the entire width of the tissue. Immediately following injury, image capture was reinitiated, and proceeded for an additional 20 minutes of observation. Each organoid was imaged, and injured, individually before being transferred to a new dish, separate from the samples awaiting processing. Between each observation period, the glass depression slide was washed with distilled water, cleaned with a Kimwipe (Kimtech Science), and loaded with fresh 0.75x MMR to avoid sample contamination across trials. Organoids were imaged across two successive days of development, corresponding to Nieuwkoop and Faber embryonic stages 37-40. All images were captured in tiff format and combined into AVI video files for computational analysis using the FIJI software package [91].

A series of preprocessing steps were performed to transform the videos into time series of calcium intensities for each identified cell in the organoids over time. The puncture event caused extreme movement of the organoids as well as a high intensity flash of calcium across the entire organoid obscuring cellular boundaries, making it difficult to track cells throughout the course of the entire video. Full videos were therefore separated into two distinct videos: pre- and post- puncture event. The end of the pre- puncture video was aligned with the time of puncture and the start of the post- puncture video was aligned with the end of immediate high intensity flash that the puncture caused. Image registration was performed on each pair of videos to correct for rotational movement of the organoids, improve the quality of video with a flatfield correction, and do any necessary video cropping. This process was carried out using Advanced Normalization Tools (ANTs) software [58]. Motion correction aligns cells in the organoid throughout time such that a segmentation algorithm can be applied to the time-average of all frames in the video to obtain cell boundaries for the entire series. Cellpose [59], a generalized deep learning model, was used for cellular segmentation. Hyperparameters of the Cellpose model were tuned for each video based on visual inspection (cell diameter = 15, cell threshold = − 2 . 0, flow threshold = 0 . 8, resample = False). Pixel intensities within each cell boundary were extracted and averaged at each frame to produce a time series of calcium intensities. These steps result in two time series arrays (pre- and post- puncture) of size # cells × # frames for each organoid (S1_Fig, S2_Table).

Utilized here was an information-theoretic pipeline to infer pairwise statistical dependencies between individual cells. Information theory has been previously discussed as a general framework for inferring systems-level structures in complex, biological systems [33,92]. Two signal preprocessing steps were applied to emphasize underlying structures in the data and allow for more meaningful inferences: global signal regression and transformation into a conditional entropy rate.

Global signal regression attempts to remove global artifacts by regressing out the mean signal across all cells. The transformation into a conditional entropy rate is a little more involved; inspired by Daube et al., [60], the calcium data for every cell was transformed into time series of the instantaneous local entropy rates. This transforms the raw calcium series into a feature-series that highlights those parts of the signal that we think are relevant to cell-cell interactions, in the style of [57]. Intuitively, this transformation highlights those moments where the cell’s behavior is deviating highly from the trend defined by its own historical dynamics. These deviations could come from two places: intrinsic randomness in X’s own dynamics, or from perturbation by another cell Y, whose activity informs on X’s own activity.

Formally, for a given cell, X, at every time t, the information content of the observation xt is given by the local entropy:

Where P(xt) is the probability of observing X = x. The local entropy (also called the Shannon information content, or surprisal) quantifies how much information about the state of X is disclosed by the observation of xt. This information can be decomposed into two components:

where i(xt−1;xt) is the information about xt that could be learned by observing the immediate past xt−1 (sometimes called the local active information storage [93,94]), while h(xt|xt−1) is the remaining information that could not be learned by observing the past (sometimes called the local conditional entropy rate [93]). All local entropies were estimated using Gaussian estimators and computed using the JIDT [95] and IDTxl [96] packages.

It is important to stress, that, in contrast to Daube et al., [60], we do not interpret this transformation into local entropy rates as “whitening" the data, in the sense of removing autocorrelation while preserving the same information. While the local entropy rate signal is less autocorrelated than the raw signal, this is not necessarily guaranteed to be the case for all data. Instead, we interpret it as a feature, highlighting those moments when the signal is deviating from the expected trend.

Finally, after transformation, excessively noisy frames associated with recording artifacts were deleted. Frames where the absolute value of the change in local entropy rate were greater than two times that standard deviation were classified as outliers and removed. The classifications were manually checked by visual inspection as well, to ensure only artifact frames were removed.

Undirected networks for each time series were generated based on instantaneous correlation (functional connectivity) [26–28,75]. Nodes of these networks are identified cells and edges are functional connections between each pair of cells in the network computed as the Gaussian mutual information between the pair’s signal time series. Gaussian mutual information was computed based on the identity:

where r is the Pearson correlation coefficient between X1 and X2 [97]. Mutual information was chosen as the transformation because unlike the Pearson correlation coefficient, is strictly non-negative, a key desiderata when analyzing networks. Edges were retained only if the p-value associated with the mutual information computation was greater than or equal to α≤10−3 (Bonferroni-corrected against the number of possible edges in the network). Each network was Bonferroni-corrected independently, making the corrected significance threshold 10−3∕((N2−N)∕2), where N is the number of nodes in a given network. Thresholding a functional connectivity network of this type remains controversial due to trade-offs for it or lack thereof. Thresholding unstructured networks can bias the resulting network towards a more complex topology [98], however, unthresholded statistical networks can include false-positive edges, creating an illusion of greater integration. Similarly, how to best infer the structure of the network is an open debate. Here we opted for a classic approach based on descriptive statistics, other options are available. In particular, approaches utilizing generative models have recently become a topic of intense research [76,77,80]. Ultimately, the decision here was thresholding on the functional connectivity networks and classic statistics on the networks. Interested researchers should defer to the particular demands of the study under question to determine what approach is suitable to the required network inference and analysis.

The co-fluctuation analysis was done following the method described in Zamani Esfahlani et al., 2020 [62]. Briefly, each pair of nodal time series was z-scored and multiplied together elementwise to construct an edge time series, where the value of the series at a given time reflects the degree to which those two edges were co-fluctuating together or in opposite directions. Then, the root sum squared deviation from the mean was computed framewise to identify how global co-fluctuations are distributed throughout the duration of the recording (see [63], for more details on high-amplitude co-fluctuations). The instantaneous co-fluctuation bears a strong resemblance to the pointwise mutual information [93,99], another time-resolved measure of dependency between dynamic variables, although the interpretations and meaning of the signs differ. Continuing with the theme of analytic flexibility, future researchers should consider whether the instantaneous co-fluctuation/edge time series or the instantaneous, pointwise mutual information makes the most sense for their particular analysis.

Code and videos for the analyses described herein is publicly available at: https://github.com/caitlingrasso/bio-connectivity.git.

This research was sponsored by the Defense Advanced Research Projects Agency (DARPA) under Cooperative Agreement No. HR0011-180200022, Grant No. DOD060 (to JB, ML, and SIW), the Lifelong Learning Machines funcprogram from DARPA/MTO. The content of the information does not necessarily reflect the position or policy of the government, and no official endorsement should be inferred. Approved for public release; distribution is unlimited. This research was also supported by the Allen Discovery Program through the Alfred P. Sloan Foundation Matter-to-Life program Grant No. G-2021-16495 (to JB and ML). Additionally, this material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. 1842491 (to CG). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. We also acknowledge support from Army Research Office contract No. W911NF-23-1-0327 (to JB). Finally, we gratefully acknowledge the support of Grant 62212 from the John Templeton Foundation (to ML). The opinions expressed in this publication are those of the author(s) and do not necessarily reflect the views of the John Templeton Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Code and videos for the analyses described herein is publicly available at: https://github.com/caitlingrasso/bio-connectivity.git.