Genome-scale transcriptional dynamics and environmental biosensing

The Dynomics microfluidic device was designed for straightforward experimental setup, reliable trap filling and cell retention, and optimal fluorescent signal from each spotted microcolony (Fig. 1 B–D). The single media inlet–outlet device requires only two fluidic connections after cell spotting and chip bonding. The media inlet channel feeds a total of 2,176 4-$\mathrm{\mu }$m-tall cell traps. Trap shape and spacing allow a 6,144 Society for Biomolecular Sciences (SBS) density pin pad to deposit cells into the back of the trap, where they grow toward the tapered opening interfacing with 50-$\mathrm{\mu }$m-tall minor media channels. These minor channels branch off of a larger 230-$\mathrm{\mu }$m-tall major media channel manifold, which eliminates the possibility of cell trap cross-contamination. Once spotted cells have reached confluence, inducer compounds can be pulsed in at user-specified frequencies with the dynamic response of each strain measured down to a 4-min temporal resolution (Fig. 1E).

Using the Dynomics platform with a previously developed GFP E. coli promoter library (22), 1,807 unique E. coli promoters were screened against nine heavy metals (Cu(II), Zn(II), Fe(III), Pb(II), Cd(II), Cr(VI), Hg(II), As(III), Sb(III)) at environmentally relevant concentrations (SI Appendix, Table S3). Screening experiments lasted 7 to 14 d, with cells exposed to a different heavy metal every 24 h. Promoters responsive to each metal can be identified through a combination of clustering and fold-change analysis. A high-level view of the 1,807 promoter time traces (Fig. 2A) and subsequent clustering (Fig. 2B) reveals distinct classes of transcriptional responses to a single 4-h zinc exposure. In Fig. 2B, clusters 1 and 2 include promoters that are up- and down-regulated, respectively, in the presence of zinc, but return to baseline expression levels within 15 h after zinc removal. Clusters 3 and 4 include promoters that are up- and down-regulated, respectively, but with slower dynamics. Gene ontology (GO) enrichment analysis (SI Appendix, Fig. S10) suggests that from these four clusters, genes associated with cellular stress are up-regulated (cellular detoxification, cellular response to toxic substance, and antibiotic catabolic processes) while genes involved in differing metabolic and biosynthesis can be either up- or down-regulated.

Individual responsive strains for each metal were identified, based on their fold-change response (Fig. 2D) to daily 4-h metal exposures (Fig. 2C). Fold-change measurements highlight the promoters displaying the strongest response to each metal. Subsequent investigation of the most responsive strains (Fig. 2E) quantitatively elucidates dynamical properties, such as amplitude, relaxation time, and response speed, all of which are important factors for their use in the study of gene expression regulation and continuous biosensing applications. While many of the identified sensing strains, such as zntA (23) or cueO (24), have well-documented metal interactions (SI Appendix, Tables S5 and S6), others are less studied or poorly annotated, particularly members of E. coli “y-ome” (25). Overall, these analyses demonstrate the utility of this platform as a screening tool for dynamic environmental-response phenotypes in a strain library. However, specific metal discrimination based on fold change alone is difficult to interpret due to promoter nonspecificity, cross-talk, noise, and low-amplitude responses.

To better discriminate between E. coli’s responses to the heavy metals used in our screening, we trained and tested two types of machine-learning models on the Dynomics data (26). The first model, known as extreme gradient boosted trees (XGBoost), is a popular decision tree ensemble-based classifier known for its ability to learn nonlinear models (27). The second one, known as a long short-term memory recurrent neural network (LSTM-RNN), is a DNN (28) selected because of its ability to effectively utilize sample sequence history to classify time series data, a property not shared by XGBoost.

Both classification algorithms outperformed random guessing of the majority class (no toxin) on the standardized experiments’ feature set, with the LSTM-RNN performing the best overall (SI Appendix, Figs. S8 and S9). As seen by examining the diagonal elements in the confusion matrix in Fig. 3A, the LSTM-RNN was able to distinguish both biotic and xenobiotic metal-spiked water from pure water with a high level of reliability.

The LSTM-RNN found iron and copper to be easily detectable biotic metals, which is not surprising given their importance to E. coli cellular function (24, 29). Cadmium was the most readily detected xenobiotic metal with the LSTM-RNN classifier, although it was sometimes confused with zinc. E. coli are known to use the same sensing and transport systems to capture and export excess amounts of these two metals, which possess the same number of valence electrons (23, 30). Most classification errors occurred during the 10 to 40 min at the start or end of the induction periods, when the LSTM-RNN occasionally had difficulty determining the exact time that each metal was added or removed from the media (Fig. 3B). This is most pronounced with the prediction of lead, for which the classifier incorrectly predicted no toxin for 48% of time points where lead was present. This is largely due to the weak promoter responses induced by 0.03 ppm lead, which is only double the Environmental Protection Agency (EPA) maximum contaminant level. In lead exposures with poor prediction, time points at the start of the 4-h induction window are misclassified as no toxin, while lead is accurately predicted near the end of this window (SI Appendix, Fig. S10). While past studies have used machine-learning frameworks to assign cells to chronologically distinct phenotypes based on their transcriptomes (31), we believe this is a different instance of a multiclass classifier successfully leveraging genome-wide transcriptional dynamics in live cells to predict exposure of a biological organism to an environmental stressor.

At present, a major obstacle to making scientific conclusions from machine-learning results is the black box problem: As an algorithm’s ability to model complex phenomena grows, its decision-making processes become more and more obscured from its operators (32). Recently, explainable artificial intelligence techniques have been employed to explain the decision making of machine-learning algorithms in the life sciences (33–35), while contributions from coalitional game theory have led to the development of a mathematically consistent method for understanding the decision-making process of any AI classifier (36, 37).

Taking advantage of these recent advances, we trained a Shapley additive explanations (SHAP) learner on both our XGBoost and LSTM classifiers (36, 38). The SHAP algorithm scores a strain’s impact on the classifier’s predictions by calculating Shapley values from cooperative game theory. Shapley values are the mathematically unique way to divide game payout between players who have collaborated with each other to achieve a common goal, assuming basic rules of fairness (39). A major advantage of SHAP is that Lundberg and Lee (36) demonstrated that it is an umbrella method that mathematically unifies several commonly used feature attribution frameworks, including LIME, Layerwise Relevance Propagation, and DeepLIFT. Viewing both SHAP values (impact on classifier output) and feature values (data fed to the classifier) with respect to time offers insight into how the classifier operates in real time (Fig. 3C). The causes of misclassification are made clearer, as SHAP dynamics reveal that the predictive impact of a strain often varies within an induction window, particularly at its start and end. Furthermore, we see how some promoters, such as zntA, positively contribute to the detection of multiple metals, which causes the classifier to rely on promoters with less-pronounced responses and lower SHAP value magnitudes to distinguish the exposed metal, explaining some misclassification instances. Finally, promoters that may not have been identified as responsive using fold-change analysis because of subtle, low-amplitude, and noisy responses can be identified via XAI. While these responders may not serve as stand-alone biosensor strains, they provide promising targets for future sensor engineering efforts. These insights highlight the ability of the LSTM-RNN classifier to compile the influence of many strains, prominent and subtle, to make an accurate prediction of the metal exposure.

The SHAP algorithm also highlights similarities and differences between how the LSTM-RNN and XGBoost make decisions. Fig. 4A shows the 15 promoters with the highest mean impact on the model, plus the promoterless strain U139, which is included as a negative control. Both methods rely heavily on the metal-sensing promoter zntA for the detection and discrimination of multiple metals, especially cadmium and zinc. Beyond zntA, XGBoost relies heavily on single strains to detect single metals, in a manner comparable to human attention patterns. The LSTM-RNN, on the other hand, utilizes many strains of moderate influence in a combinatorial fashion; this tendency to find a different representation from that of the human visual system has been noted in other works (40). These trends are also seen when looking at the top 15 promoters for each individual metal class (SI Appendix, Fig. S11).

The ability of the explained classifiers to identify promoters involved in metal response serves as a valuable scientific tool, suggesting potential pathways and genes for further investigation. This value is highlighted by looking at a subset of the 10 most-impactful promoters individually for cadmium and iron inductions (Fig. 4B). These summary plots illustrate how the two classifiers make similar decisions through different methods. In the case of cadmium, zntA plays a significant role for both classifiers, while different sets of genes involved in ion transport or amino acid synthesis are identified for each. Most notably, the metE and metB promoters which are involved in methionine synthesis, an amino acid known to chelate cadmium (41), are identified by XGBoost, while the LSTM-RNN uses only the metE regulator, metR, for detection. Similarly with iron, we see XGBoost rely on members of the arginine synthesis, argA and argC, while the LSTM-RNN relies on different promoters that are involved in other metabolic or biosynthetic processes.

Given the severe impact of heavy metals on human health (42) and the persistence of water quality issues in the United States (43), we sought to deploy the Dynomics platform as a real-time water quality biosensor. To verify that this device was functional on waters of varying ion compositions, we conducted experiments with media made from municipal water samples from San Diego, Seattle, Chicago, Miami-Dade, and New York City with added cadmium. Fig. 5A shows the LSTM-RNN classifier predictions for cadmium exposures on each city’s water supply. While there is some misclassification of cadmium as zinc, there are few instances of incorrectly predicting the presence of a toxin versus water, even with largely different water compositions between cities. Additionally, the mean absolute SHAP values for each city correlated strongly with those for laboratory Milli-Q water (${\mathrm{R}}^2$ = 0.853), indicating that water composition did not affect gene response. zntA was the best predictor of cadmium presence across all water compositions (SI Appendix, Fig. S13).

The Dynomics device was also exposed to samples collected from the Gold King Mine spill in August 2015. Fig. 5B shows the predictions of the LSTM-RNN classifier on samples from the spill, collected from the San Juan River. The classifier predictions are output as multiclass, multilabel probability vectors. As the sample was introduced onto the device, the probability of uncontaminated water decreased significantly while the probabilities of the other metals increased. The metal with the highest probability, iron, was also the most abundant metal in the samples, as measured by inductively coupled plasma mass spectrometry (ICP-MS) (SI Appendix, Table S2). Despite the classifier not being trained on combinations of metals or at the concentrations present in these samples, the ability to reliably report the presence of the most prominent metal and, to a lesser degree, the less abundant metals suggests the broad applicability of this platform for heavy metal detection.

Our group has previously described the microfabrication techniques used to pattern SU-8 photoresist onto a silicon wafer to create the mold for our device (46). A poly-dimethylsiloxane (PDMS) device was made from the wafer by mixing 77 g of Sylgard 184 and pouring it on the wafer centered on a level $5\prime \prime$ $\times$ $5\prime \prime$ glass plate surrounded with an aluminum foil seal. The degassed wafer and PDMS were cured on a flat surface for 1 h at 95 °C.

The E. coli promoter library (22) was arrayed using the Singer ROTOR Stinger (Singer Instrument Co. Ltd.) attachment from 96-well density formatted agar plates onto four 1,536-density formatted agar plates to match the layout of the cell traps on the microfluidic device. At the time of experimental setup, the four 1,536-density agar plates were combined onto one 6,144-density agar plate using the Singer ROTOR and grown for 2 h at 37 °C before being transferred to the device.

A PDMS device cleaned with 70% ethanol and adhesive tape was aligned to a custom fixture compatible with the Singer ROTOR. Both the fixture and a clean $4\prime \prime$ $\times$ $3\prime \prime$ glass slide sonicated with 2% Hellmanex III were exposed to oxygen plasma. Cells were spotted from the previously arrayed 6,144-density agar plate to the aligned PDMS device using the Singer ROTOR spotting robot. The device and glass slide were bonded together and cured at 37 °C for 2 h.

Microfluidic experiments were performed on a custom optical assembly described in SI Appendix. Continuous imaging occurred every 10 min, imaging both the transmitted light and GFP fluorescence channels. Cells were grown in the device on LB media with kanamycin, 0.075% Tween-20, and 50 mM methyl $\alpha$-d-mannopyranoside until traps were filled to confluence. The media were then switched to a heavy-metal–trace-free minimal media (HM9) minimal media described in SI Appendix, Table S1, which was based on a previous study (47) and optimized for microfluidic E. coli growth with minimal traces of metals. Cells were grown on HM9 for 48 h before inducing with heavy metals. Heavy metal inductions occurred once a day for 4 h with HM9 media flowing on chip for the remaining 20 h. Quintuplicate inductions of each metal were performed in a random order across multiple experiments, with each experiment lasting 7 to 14 d. A total of 2,176 time traces were collected from each experiment (Fig. 2 A–C). Extracted time traces were normalized to remove device background fluorescence and strain background fluorescence (SI Appendix). Detailed methods on experimental setup and data collection can be found in SI Appendix, Table S3.

Water samples were obtained from the Department of Water Management at the City of Chicago in Chicago, IL; the Alex Orr Water Treatment Plant in Miami, FL; the New York City Department of Environmental Protection and Bureau of Water Supply in Corona, NY; the Seattle Public Utilities Water Quality Laboratory in Seattle, WA; and the Alvarado Water Treatment Plant in San Diego, CA. HM9 media for each city water experiment were prepared by diluting 5$\times$ HM9 concentrate made from Milli-Q water with the water obtained from each city. The microfluidic device was initially grown on LB media with kanamycin, 0.075% Tween-20, and 50 mM methyl $\alpha$-d-mannopyranoside until traps were filled to confluence and then switched to HM9 made with city water for the remainder of the experiment. Cadmium diluted in the HM9 city water media was used to perform inductions as described in SI Appendix.

Water was collected from Mexican Hat, UT in August 2015 when the Gold King Mine spill plume reached the collection point in the San Juan River. Samples were stored in 0.5% HCl acid until tested. HM9 media were prepared by diluting 5× HM9 concentrate made from Milli-Q water with filtered San Juan River samples. The pH was adjusted to 7.05. The metal concentrations of the HM9 San Juan River samples were tested by ICP-MS at the Environmental and Complex Analysis Laboratory (ECAL) at University of California, San Diego. Four-hour inductions were performed as described in SI Appendix.

We transformed our 18 standardized experiments’ time points into a first derivative-based feature for the training and testing feature sets. To optimize the classifiers, extensive Bayesian optimization searches were used to find optimal hyperparameter combinations (48). Throughout our hyperparameter searches, we used leave-one-out cross-validation on a per-experiment basis and appropriate overfitting-prevention strategies to ensure that any resultant classifier would generalize to future datasets. All classifiers were evaluated using the $F_1$-macro scoring metric. The $F_1$-macro score, which is the per-class average of the harmonic mean of precision and recall, was especially well suited because of our dataset’s large multiclass imbalances, with water making up ∼86% of the final feature set (49). Finally, all generalization evaluations were performed by recording the results of using leave-one-out cross-validation with early stopping and then taking the mean prediction across the cross-validation’s output.

Preprocessed, labeled machine-learning features, the corresponding library strain position records, and the relevant metadata and code for our experiments are available on the University of California San Diego Biodynamics Laboratory website (http://biodynamics.ucsd.edu/downloads).