Single phage high-throughput screening method based on microfluidic-culturomics coupling

Abstract

The application discloses a single phage high-throughput screening method based on microfluidic-cultivation-omics coupling and belongs to the technical field of phage screening. A full-liquid-phase, non-plaque-dependent fluorescent cracking micro-reactor array is constructed and is allocated with orthogonal coding of nutritional stress. Time sequence fluorescent imaging and nonlinear dynamics analysis are performed on each micro-reactor. A multi-dimensional cracking dynamics characteristic vector is extracted and a multi-dimensional response landscape atlas is generated. According to the atlas, hidden cracking events deviating from a standard cracking decay model are identified and pulse dielectrophoresis driving recovery is performed. Finally, multi-round atlas data is collected. An active learning strategy is adopted to infer and generate a stress-induced culture medium combination in the next round with the goal of maximizing information gain. The application breaks through the plaque dependence limitation, realizes closed-loop evolution of cracking behavior panorama characterization and screening path, and improves the discovery efficiency of unconventional phages.

Description

Technical Field

[0001] This invention relates to the field of phage screening technology, and in particular to a high-throughput screening method for single phages based on microfluidics-culture omics coupling. Background Technology

[0002] Bacteriophages, as natural enemies of bacteria, have broad application prospects in the treatment of drug-resistant bacterial infections, food safety control, and microecological regulation. Efficient and comprehensive screening of phages with lytic activity from complex environmental samples is a fundamental prerequisite for the development and application of phage resources. Traditional phage screening methods mainly rely on the double-layer agar plate method, determining lytic activity by observing whether plaques form on the host bacterial colony. This method has the following significant limitations: First, the screening throughput is limited by the plate area and the efficiency of manual operation, making it difficult to meet the screening needs of large-scale samples; second, this method essentially relies on the ability of phages to form visible plaques on the surface of solid culture media. Phages with slow lysis rates, insufficient lysis amounts, or limited diffusion in solid media are easily missed, leading to the omission of a large number of potentially functional phage resources; third, conventional culture conditions are singular and fixed, failing to simulate the diverse metabolic states faced by host bacteria in the natural environment, thus missing phages that only exhibit lytic activity under specific nutrient stress or stress conditions.

[0003] In recent years, microfluidic technology has provided a new platform for high-throughput single-cell / single-virus analysis. By co-encapsulating bacteriophages and hosts in water-in-oil droplets, they can be independently co-cultured in microscale reactors. However, existing microfluidic-based phage screening methods mostly employ homogenized culture conditions, lacking systematic, differentiated design and quantitative characterization of nutrient environmental factors. This makes it difficult to reveal the dependence of phage lysis behavior on host metabolic state, let alone guide the screening pathway accordingly. Furthermore, existing positive screening criteria typically rely on preset fixed fluorescence thresholds, lacking effective mechanisms for identifying and capturing occult lysis events whose lysis kinetics deviate from the norm, resulting in insufficient completeness and diversity of screening results. Therefore, there is an urgent need for a high-throughput phage screening method that can overcome the limitations of traditional plaque dependence, systematically explore nutrient environmental variables, accurately identify atypical lysis signals, and achieve adaptive evolution of the screening pathway. Summary of the Invention

[0004] This invention overcomes the shortcomings of the prior art and provides a high-throughput screening method for single phages based on microfluidics-culture omics coupling.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: This invention discloses a high-throughput screening method for single phages based on microfluidics-culture omics coupling, comprising the following steps: S1. Construct an oil-in-water microdroplet array containing bacteriophages, hosts, and fluorescent probes to form a fully liquid-phase, non-plaque-dependent fluorescent lysis microreactor array. Couple the microreactor array with the microarray chip and assign a unique orthogonal nutrient stress code to each microreactor. S2. Perform time-series fluorescence imaging on the co-culture system in each microreactor, perform nonlinear dynamic analysis on the fluorescence decay signal of each microreactor, and extract a multidimensional fragmentation dynamic feature vector containing the fragmentation acceleration point, the maximum fragmentation rate point, and the plateau phase. S3. Associate the multidimensional lysis kinetics feature vector with the orthogonal encoding of the nutrient stress of the corresponding microreactor to generate a multidimensional response landscape map describing the phage lysis behavior as a function of nutrient conditions. S4. Identify hidden lysis events that deviate from the preset standard lysis decay model as potential positives, and perform pulsed dielectrophoresis-driven recovery on the microreactor based on the nutrient stress orthogonal coding corresponding to the hidden lysis events, so as to obtain unconventional phages that exhibit liquid phase lysis activity in a homogeneous droplet environment. The unconventional phages include non-plaque-forming phages. S5. By gathering the multidimensional response landscape maps accumulated from multiple rounds of screening as prior knowledge, and adopting an active learning strategy with the goal of maximizing information gain, the stress-inducing culture medium combination required for the next round of screening is inferred and generated.

[0006] This invention addresses the technical deficiencies in the prior art and has the following beneficial effects: This invention constructs a fully liquid-phase, plaque-free fluorescent lysis microreactor array, deeply integrating the concept of differentiated nutrient stress in cultinomy with microfluidic droplet technology. It replaces single-point fluorescence determination with full-process analysis of lysis kinetics, overcoming the dependence of traditional methods on plaque formation. This allows for the capture of unconventional phages that cannot form plaques in solid media but possess liquid-phase lysis activity, thus broadening the scope of phage resource discovery. Simultaneously, through the accumulation of multiple rounds of screening data and the closed-loop evolution of the screening path driven by active learning, the screening strategy possesses self-optimization capabilities. It can autonomously approach the enrichment region of occult lysis events within a vast nutrient stress condition space, effectively overcoming the limitations of traditional methods with single, fixed screening conditions and isolated multi-round iterations. This significantly improves the discovery efficiency of rare phages and the intelligence level of the screening system. Attached Figure Description

[0007] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other embodiments can be obtained from these drawings without creative effort.

[0008] Figure 1 This is a flowchart of a high-throughput screening method for single phages based on microfluidics-culture omics coupling. Detailed Implementation

[0009] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0010] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0011] like Figure 1 As shown, this invention discloses a high-throughput screening method for single phages based on microfluidics-culture omics coupling, comprising the following steps: S1. Construct an oil-in-water microdroplet array containing bacteriophages, hosts, and fluorescent probes to form a fully liquid-phase, non-plaque-dependent fluorescent lysis microreactor array. Couple the microreactor array with the microarray chip and assign a unique orthogonal nutrient stress code to each microreactor. S2. Perform time-series fluorescence imaging on the co-culture system in each microreactor, perform nonlinear dynamic analysis on the fluorescence decay signal of each microreactor, and extract a multidimensional fragmentation dynamic feature vector containing the fragmentation acceleration point, the maximum fragmentation rate point, and the plateau phase. S3. Associate the multidimensional lysis kinetics feature vector with the orthogonal encoding of the nutrient stress of the corresponding microreactor to generate a multidimensional response landscape map describing the phage lysis behavior as a function of nutrient conditions. S4. Identify hidden lysis events that deviate from the preset standard lysis decay model as potential positives, and perform pulsed dielectrophoresis-driven recovery on the microreactor based on the nutrient stress orthogonal coding corresponding to the hidden lysis events, so as to obtain unconventional phages that exhibit liquid phase lysis activity in a homogeneous droplet environment. The unconventional phages include non-plaque-forming phages. S5. By gathering the multidimensional response landscape maps accumulated from multiple rounds of screening as prior knowledge, and adopting an active learning strategy with the goal of maximizing information gain, the stress-inducing culture medium combination required for the next round of screening is inferred and generated.

[0012] Specifically, S1 is: A baseline metabolic model of the host bacteria is obtained. Based on the sensitivity coefficients of each metabolic pathway in the baseline metabolic model, multiple stressor factors targeting different metabolic nodes are selected from a pre-set stressor factor library. Based on the orthogonal experimental design principle, a nutritional stress orthogonal coding table composed of different concentration levels of the stressor factors is generated, where each group of codes corresponds to a unique stress-inducing culture medium formulation. The baseline metabolic model of the host bacterium refers to a set of information that reflects the relative activity or flux distribution of various metabolic pathways under standard culture conditions. Specifically, in this embodiment, the relative flux coefficients of pathways such as central carbon metabolism, amino acid synthesis, nucleotide metabolism, and cofactor synthesis in the host bacterium are obtained by searching publicly available genome-scale metabolic network reconstruction results or by using conventional phenotypic microarray assay data in the art. This information is used to construct the baseline metabolic model of the host bacterium. The core output of the baseline metabolic model is a set of sensitivity coefficients for metabolic pathways. The higher the sensitivity coefficient, the more significant the impact of perturbation of that metabolic pathway on the overall physiological state of the host bacterium, thus providing a basis for the targeted selection of stressor factors.

[0013] A pre-defined stressor library refers to a collection and categorized set of compounds that can inhibit or stress specific metabolic nodes in bacteria. Each stressor in the library corresponds to one or more metabolic targets, such as 2-deoxyglucose targeting carbon metabolism, glyphosate targeting amino acid synthesis, physcin targeting fatty acid synthesis, fosfomycin targeting cell wall synthesis, rotenone targeting the respiratory chain, and sodium chloride targeting osmolarity regulation. The target information and commonly used working concentration ranges for each stressor in the library are pre-recorded.

[0014] Using the preset matrix positions of the microarray chip as an index, each group of codes in the nutrient stress orthogonal coding table is sequentially assigned to each matrix position according to the spatial mapping relationship. The stress-inducing culture medium corresponding to each group of codes is pre-fixed at the corresponding matrix position to form a nutrient stress orthogonal coding array with position coding information. In a microfluidic chip, a suspension of phage to be screened, a suspension of host bacteria, and a fluorescent nucleic acid dye are mixed in a preset ratio to form an aqueous phase, which is then injected into a flow focusing structure together with an oil phase containing a surfactant. By adjusting the flow rate ratio of the two phases, water-in-oil droplets are generated, so that the size of the droplets matches the lattice size of the microarray chip. The generated water-in-oil microdroplets are introduced into the surface of the prepared nutrient stress orthogonal coding array, and each microdroplet falls into and is fixed at the corresponding array point. Through the diffusion reconstruction of the aqueous phase inside the microdroplet and the pre-cured stress-induced culture medium on the array point, the nutrient stress orthogonal coding carried by the array point is assigned to the corresponding microdroplet, forming a co-culture system in which each microreactor has a unique nutrient stress condition identifier.

[0015] In practice, several metabolic pathways with high sensitivity coefficients are screened from the baseline metabolic model, and corresponding stressor factors targeting these pathways are selected from a stressor factor library. For example, three stressor factors A, B, and C are selected, each corresponding to a different metabolic target. Subsequently, an orthogonal coding table for nutritional stress is generated based on the principles of orthogonal experimental design. For instance, if two concentration levels, high and low, are set for each factor, a three-factor, two-level orthogonal design can be used to generate an orthogonal coding table containing several sets of experimental conditions. Each code is uniquely determined by the combination of levels of each factor, such as code 1 - A low B low C low, code 2 - A low B high C high, code 3 - A high B low C high, code 4 - A high B high C low. Each code corresponds to a unique stress-inducing culture medium formulation, which is composed of the specific values ​​of each stressor factor at the corresponding concentration level. The advantage of using orthogonal design is that it achieves systematic coverage of the multi-factor, multi-level nutritional stress space with a smaller number of experimental conditions, and the conditions in each group have statistical independence and balanced dispersion.

[0016] After the coding table is generated, if the microarray chip is arranged in a standard row and column format, the i-th group of codes in the coding table is assigned to the matrix point position in the i-th row and i-th column, forming a one-to-one correspondence between the matrix point coordinates and the nutrient stress orthogonal coding. Using microarray spotting technology, the stress-inducing culture medium solution corresponding to each group of codes is added to the corresponding matrix point at preset volume points. After low-temperature vacuum drying or freeze-drying, the culture medium components are pre-solidified on the matrix point surface, forming a nutrient stress orthogonal coding array with positional coding information. Here, positional coding information means that the physical coordinates of each matrix point implicitly contain the nutrient stress orthogonal coding it carries, without the need for additional markings or labels.

[0017] On the microfluidic chip side, the phage suspension to be screened, the host bacterial suspension, and the fluorescent nucleic acid dye are mixed in a preset ratio to form an aqueous phase. The fluorescent nucleic acid dye used is one that can bind to the nucleic acid within bacterial cells and maintain fluorescence when the cell membrane is intact, such as SYTO 9, whose fluorescence intensity decreases with cell lysis and nucleic acid release. The aqueous phase and the oil phase containing surfactant are co-injected into the flow cytometry focusing structure. By controlling the flow rate ratio of the two phases, shear force is generated at the interface between the aqueous and oil phases, generating uniformly sized monodisperse water-in-oil microdroplets. Crucially, the size of the generated microdroplets must match the size of the microarray chip lattice, so that each lattice can only accommodate one microdroplet.

[0018] Finally, the generated water-in-oil microdroplets were introduced into the surface of the prepared orthogonal nutrient stress coding array. Due to the combined effect of microdroplet density and the physical confinement effect of the array surface, each microdroplet fell into and fixed at its corresponding lattice point. The aqueous phase within the microdroplet diffused and reconstructed with the pre-cured stress-induced culture medium on the lattice point, and the pre-cured stress-generating factors redissolved and diffused into the microdroplet, thus completely assigning the orthogonal nutrient stress coding carried by the lattice point to the corresponding microdroplet. This formed a co-culture system where each microreactor had a unique nutrient stress condition identifier. Subsequent detection results of any microreactor could be traced back to the corresponding nutrient stress coding through its lattice point location, achieving end-to-end traceability from culture conditions to lytic phenotype.

[0019] Specifically, S2 is: For each co-cultured microreactor coupled to the nutrient stress orthogonal coding array, multiple frames of fluorescence images are acquired at preset time intervals. The fluorescence intensity value of each microreactor is extracted from each frame of the image, and an original fluorescence decay sequence arranged along the time axis is constructed for each microreactor. At the same time, the nutrient stress orthogonal coding corresponding to each microreactor is recorded. Detrended fluctuation analysis was performed on the original fluorescence decay sequence of each microreactor. By calculating the local fluctuation function of the sequence, non-biological signal trends caused by instrument noise and photobleaching were filtered out to obtain the corrected fluorescence decay sequence. Based on the corrected fluorescence decay sequence, the first-order difference sequence of each sampling point was calculated. Based on the sign change and amplitude jump of the first-order difference sequence, the corrected fluorescence decay sequence is adaptively divided into three dynamic intervals: hysteresis period, burst fragmentation period, and decay plateau period, and a fragmentation stage division scheme containing the boundary time of each interval and the fluorescence decay rate within the interval is generated. Based on the fragmentation stage division scheme, the moment when the first spontaneous fragmentation threshold in the first-order differential sequence exceeds the preset fragmentation threshold is identified as the fragmentation acceleration point, the peak moment of the first-order differential sequence during the burst fragmentation period is extracted as the maximum fragmentation rate point, and the moment when the fluorescence intensity value tends to a steady state during the decay plateau period is identified as the plateau entry point. By splicing together the pyrolysis acceleration point, maximum pyrolysis rate point, plateau inlet point, and average decay rate during the explosive pyrolysis period of the same microreactor, a multidimensional pyrolysis kinetic feature vector of the corresponding microreactor is formed.

[0020] It should be noted that the imaging platform continuously scans and images all microreactors in the array at preset time intervals, such as one frame every 10 minutes, until the fluorescence intensity of the negative control microreactors, where the host bacteria naturally decay, stabilizes. For each frame of fluorescence image acquired, the contour region of each microreactor is identified using an image segmentation algorithm, and the average fluorescence intensity value of the pixels within that region is extracted as the fluorescence intensity value of that microreactor at that time point. Arranging the fluorescence intensity values ​​of the same microreactor across frames along the time axis constructs the original fluorescence decay sequence of that microreactor. Simultaneously, the orthogonal encoding of nutrient stress corresponding to each microreactor is recorded, establishing a data-level correlation between the fluorescence sequence and culture conditions.

[0021] In the original fluorescence decay sequence, in addition to the true fluorescence decay caused by host bacterial lysis, two main non-biological signal interferences are superimposed: one is the electronic noise of the imaging system itself and the short-term fluctuations of the excitation light source, i.e., instrument noise; the other is the photobleaching effect of the fluorescent nucleic acid dye under continuous excitation light, which causes the fluorescence intensity to decrease monotonically over time, i.e., photobleaching trend. If these two interferences are not filtered out, they will directly affect the accuracy of subsequent lysis stage segmentation. This embodiment uses detrended fluctuation analysis to solve this problem. The basic principle is: the original sequence is divided into several local windows, and the change of fluorescence intensity over time is fitted with a local polynomial within each window to obtain a local trend function reflecting the slow non-biological trend within that window. By splicing the local trend functions of each window, a global non-biological trend curve covering the entire time series is obtained. Then, this global trend curve is subtracted from the original fluorescence decay sequence to obtain the corrected fluorescence decay sequence. After correction, the sequence only retains fluorescence intensity fluctuations caused by changes in the physiological state of the host bacteria (mainly the lysis process). Based on this, the fluorescence intensity difference between adjacent sampling points is calculated along the time axis of the corrected fluorescence decay sequence to obtain the first-order difference sequence. The sign and amplitude of the first-order difference sequence directly characterize the instantaneous direction and rate of change of fluorescence intensity.

[0022] Next, the kinetic interval of the corrected fluorescence decay sequence is adaptively segmented based on the first-order difference sequence. Observations reveal that during a typical fragmentation process, the first-order difference sequence exhibits distinct phased changes: initially, the first-order difference value fluctuates slightly around zero, indicating relatively stable fluorescence intensity; this phase corresponds to the fragmentation lag phase. Subsequently, the first-order difference value shows a significant negative jump and remains at a large negative value, indicating that the fluorescence intensity enters a rapid decay phase; this phase corresponds to the burst fragmentation phase. Finally, the first-order difference value returns to near zero, indicating that the fluorescence intensity tends to stabilize; this phase corresponds to the decay plateau phase. Based on this pattern, this embodiment determines the boundary moments of each stage by detecting the sign change and amplitude jump of the first-order difference sequence. Specifically, searching backward from the start of the sequence, when the first-order difference values ​​of multiple consecutive sampling points are all below a preset upper limit for micro-amplitude fluctuations (e.g., the upper limit is set to 0.5 times the sequence standard deviation), these sampling points are classified as hysteresis periods. When the first-order difference value is detected to first exceed this upper limit for micro-amplitude fluctuations and subsequently remain negative, this abrupt change is marked as the boundary moment between the hysteresis period and the burst lysis period. Continuing the search backward, when the absolute value of the first-order difference value falls back below the upper limit for micro-amplitude fluctuations and remains stable, this fallback moment is marked as the boundary moment between the burst lysis period and the decay plateau period. Thus, a lysis stage division scheme is generated, comprising the boundary moments of three intervals and the average fluorescence decay rate within each interval.

[0023] Based on the lysis phase segmentation scheme, four key kinetic characteristic points are further extracted. First, the lysis acceleration point: During the burst lysis phase, the first-order difference sequence is retrieved sequentially from the boundary time, and the sampling time when the first first-order difference value exceeds a preset spontaneous lysis threshold is determined as the lysis acceleration point. The spontaneous lysis threshold can be set based on the distribution of first-order difference values ​​in the negative control microreactor within the same time period; for example, the 99th percentile of the negative control's first-order difference value distribution can be used as the threshold to ensure that the lysis acceleration point reflects the phage's active lysis behavior, which is significantly different from spontaneous lysis. Second, the maximum lysis rate point: The moment of the negative peak of the first-order difference sequence during the burst lysis phase is extracted as the maximum lysis rate point. Third, the plateau entry point: During the decay plateau phase, the starting moment when the fluorescence intensity value tends to a steady state is identified; for example, when the fluorescence intensity fluctuation amplitude of multiple consecutive sampling points is less than 5% of the peak amplitude of the corrected fluorescence decay sequence, the first point of the continuous interval is determined as the plateau entry point. Fourth, the average decay rate during the explosive pyrolysis period: the absolute value of the first-order difference of each sampling point during the explosive pyrolysis period is averaged to obtain the average decay rate of this stage. The pyrolysis acceleration point, maximum pyrolysis rate point, plateau inlet point, and average decay rate during the explosive pyrolysis period extracted from the same microreactor are spliced ​​together in a predetermined order to form the multidimensional pyrolysis kinetic characteristic vector of the microreactor.

[0024] Specifically, S3 is: Using the orthogonal coding of nutrient stress as the primary key, the multidimensional pyrolysis kinetic feature vectors of each microreactor are grouped according to the coding, a mapping table between each orthogonal coding of nutrient stress and the corresponding feature vector is established, and the statistical distribution parameters of each feature vector in each coding group are calculated. The statistical distribution parameters include the central tendency measure and the dispersion measure of the values ​​of each dimension. All orthogonal codes of nutritional stress in the mapping table are hierarchically clustered according to the types and concentration levels of the stressor factors they contain, generating a hierarchical similarity tree structure of nutritional stress conditions, so that codes with similar metabolic pathway perturbation intensities are grouped into the same branch, and each branch is assigned a corresponding stress perturbation pattern label. Based on the hierarchical relationship of the stress disturbance mode labels, cross-condition comparative analysis is performed on the multidimensional lysis kinetic feature vectors between different coding groups. By calculating the significance of the differences between the feature vectors of each dimension between groups, the deviation direction and magnitude of the lysis kinetic behavior under each stress disturbance mode label compared with the baseline culture medium conditions are identified, forming a lysis behavior deviation matrix. Each row of the rift behavior deviation matrix is ​​mapped to a behavior response surface grid point in the nutrient stress condition space. The nutrient stress condition distance defined by the hierarchical similarity tree structure is used as the surface interpolation skeleton. Conformal interpolation is performed between adjacent grid points to generate a multi-dimensional response landscape map with full coverage. The multidimensional response landscape map is analyzed for topological features to extract peak, valley and saddle point regions. Each region is then labeled based on the stress perturbation pattern label to obtain a set of cleavage behavior response patterns associated with specific metabolic pathway perturbations.

[0025] It should be noted that upon entering this step, the system already possesses two datasets: one is the multidimensional fission kinetic feature vector of each microreactor output in step S2, and the other is the nutrient stress orthogonal code assigned to each microreactor in step S1. First, using the nutrient stress orthogonal code as the primary key, the multidimensional fission kinetic feature vectors of all microreactors are grouped according to the code. For example, if this round of screening uses sixteen groups of orthogonal codes, all microreactors are divided into sixteen code groups, and a mapping table is established between each code group and several feature vectors within that group. Since the same code group contains multiple microreactors with repetitive technologies, it is necessary to calculate the statistical distribution parameters of each dimension within the group to characterize the representative features of the fission behavior under that coding condition. This embodiment uses the median as a measure of central tendency and the interquartile range as a measure of dispersion, thereby obtaining the median and interquartile range values ​​of each dimension of the feature vector corresponding to each code group.

[0026] Next, all orthogonal codes for nutritional stress in the mapping table are hierarchically clustered according to the types and concentration levels of the stressor factors they contain, to construct a hierarchical similarity tree structure of nutritional stress conditions. The logic of the clustering operation is that the similarity between nutritional stress conditions depends on the similarity of the metabolic pathways they perturb and the similarity of the perturbation intensity. Specifically, the concentration level values ​​of each stressor factor in each group of codes are first extracted to construct the nutritional stress feature vector of that code. Then, based on the nutritional stress feature vectors of each code, the Euclidean distance between each pair of codes is calculated, and a bottom-up hierarchical clustering algorithm is used to gradually merge the closest codes to generate a tree diagram reflecting the affinity between codes. In this tree structure, codes with similar metabolic pathway perturbation intensities will be preferentially merged into the same branch. By segmenting the dendrogram at an appropriate height, all codes are divided into several clusters. The stress-induced culture media corresponding to the codes within each cluster share similarities in metabolic targets and stress intensity. Each cluster is assigned a unique stress perturbation pattern label; for example, label A represents a coding cluster primarily targeting carbon metabolism pathways with low stress intensity, while label B represents a coding cluster simultaneously targeting carbon metabolism and cell wall synthesis with high stress intensity. Thus, the discrete orthogonal codes for nutrient stress are organized into a hierarchical coding classification system with biological meaning.

[0027] After obtaining the stress perturbation pattern labels, a cross-condition comparative analysis was performed on the multidimensional lysis kinetic eigenvectors of different coding groups to quantify the degree of deviation of each stress perturbation pattern from the baseline culture conditions. The baseline culture conditions refer to standard culture conditions without any stressor factors, and the corresponding multidimensional lysis kinetic eigenvectors serve as the reference baseline for this comparative analysis. Taking a coding group under a specific stress perturbation pattern label as an example, the median of each dimension of the eigenvector of the coding group was compared one by one with the corresponding dimension values ​​of the baseline culture conditions, and the difference in each dimension was calculated. The sign of the difference indicates the direction of deviation; a positive value indicates a delayed time point, and a negative value indicates an earlier time point; the absolute value of the difference indicates the magnitude of the deviation. To determine the statistical significance of the deviation, the Mann-Whitney U test, a standard procedure in this field, was used. One set of eigenvectors was taken from each microreactor under the baseline conditions, and another set was taken from each microreactor within the coding group. The significance level of the two-sided test was calculated. When the significance level is below a preset significance threshold, such as 0.5, the deviation in that dimension is considered statistically significant, and the direction and magnitude of the deviation are recorded in the corresponding position of the lysis behavior deviation matrix; when the difference is not significant, the position is recorded as zero deviation. The above operation is performed on all coding groups under all stress perturbation mode labels, thus forming a lysis behavior deviation matrix with stress perturbation mode labels as rows and each dimension of the multidimensional lysis dynamics feature vector as columns.

[0028] The task of constructing a multidimensional response landscape atlas is to transform the aforementioned lysis behavior deviation matrix from a discrete tabular form into a continuous, visualized surface form. First, each row of the deviation matrix is ​​considered as a grid point on a behavior response surface in the nutrient stress condition space. The coordinate system of the nutrient stress condition space is defined by a pre-order hierarchical similarity tree structure, where the distance between each encoding in the tree structure is mapped to the nutrient stress condition distance between grid points in that space. Using this tree structure as the surface interpolation skeleton, a conformal interpolation algorithm is used to fill the gaps between adjacent grid points, ensuring that the interpolated surface strictly passes through the measured values ​​at existing sampling points and maintains a monotonous and smooth transition between sampling points, ultimately generating a multidimensional response landscape atlas with full coverage. In the atlas, each coordinate point corresponds to a nutrient stress condition or its interpolated inference condition. The height dimension of the atlas represents the comprehensive deviation of the lysis kinetics. A larger height dimension indicates a more significant deviation of the phage lysis behavior from the baseline condition under that nutrient stress condition.

[0029] Finally, topological feature analysis was performed on the generated multidimensional response landscape map to extract the structured information it contained. The analysis identified three types of topological feature regions in the map: peak regions with significantly higher height values ​​than the surrounding areas, indicating nutrient stress conditions where lysis behavior deviates most significantly and may exhibit specific lysis phenotypes; valley regions with significantly lower height values ​​than the surrounding areas, indicating nutrient stress conditions where lysis behavior is inhibited; and saddle point regions where the height gradient direction changes, indicating critical condition zones where the lysis behavior response mode changes. Based on the stress perturbation pattern labels associated with each point in the map, the peak, valley, and saddle point regions were traced and labeled, i.e., the types and concentration levels of stressor factors corresponding to each region were reviewed to identify which metabolic pathway perturbation led to the corresponding lysis behavior response mode. The above labeling results were summarized to obtain a set of lysis behavior response modes associated with specific metabolic pathway perturbations.

[0030] Specifically, S4 is: The multidimensional pyrolysis dynamics feature vector corresponding to the peak region in the multidimensional response landscape map is called as the standard pyrolysis dynamics template. The values ​​of each dimension of each set of standard pyrolysis dynamics templates are set as the center value of the reference interval of the corresponding dimension. The upper and lower bounds of the reference interval are determined based on the dispersion measure of each dimension, and a standard pyrolysis decay model enclosed by the multidimensional reference interval is constructed. The multidimensional pyrolysis kinetic feature vectors of all microreactors generated in this round of screening are traversed. The values ​​of each dimension of each feature vector to be judged are compared with the reference interval of the corresponding dimension in the standard pyrolysis decay model. At least one deviation feature vector whose value falls outside the corresponding reference interval is identified, and the microreactor to which the deviation feature vector belongs is marked as a candidate hidden pyrolysis microreactor. For each candidate hidden pyrolysis microreactor, a known pyrolysis behavior response pattern with the same stress perturbation mode label as the candidate hidden pyrolysis microreactor is extracted from the multidimensional response landscape map. By comparing the angle between the feature vector of the candidate hidden pyrolysis microreactor and the vector of the known pyrolysis behavior response pattern, it is determined whether the candidate hidden pyrolysis microreactor is a variant of the known pyrolysis behavior pattern. If it is not a variant, it is officially marked as a hidden pyrolysis event. Obtain the nutrient stress orthogonal code corresponding to each hidden pyrolysis event. Based on the spatial mapping relationship between the nutrient stress orthogonal code and the microarray chip lattice point position, look up the physical lattice point coordinates of the microreactor where the hidden pyrolysis event is located, and write the physical lattice point coordinates into the target address register of the pulse dielectrophoresis drive sequence. A non-uniform alternating electric field is applied to the microarray chip region pointed to by the target address register. By adjusting the electric field frequency, a difference in dielectric force is generated between the microreactor and the surrounding continuous phase, driving the microreactor to migrate along the electric field gradient direction to a preset recovery channel. A negative pressure is applied at the recovery port to draw the microreactor that has migrated to the recovery port into the collection chamber, thereby obtaining an enriched solution containing unconventional phages with the liquid phase lysis activity.

[0031] Specifically, determining whether a candidate hidden pyrolysis microreactor is a variant of a known pyrolysis behavior pattern by comparing the angle between the feature vector of the candidate hidden pyrolysis microreactor and the vector of a known pyrolysis behavior response pattern involves: Calculate the cosine similarity between the multidimensional pyrolysis kinetic feature vector of the candidate hidden pyrolysis microreactor and the representative vector of the known pyrolysis behavior response mode under the same stress perturbation mode label. When the cosine similarity is lower than the preset mutation determination threshold, the candidate hidden pyrolysis microreactor is determined to be a variant of the known pyrolysis behavior mode.

[0032] Specifically, this step can be divided into two stages: identification of occult lysis events and pulsed dielectrophoresis-driven recovery.

[0033] In the first stage, a standard fragmentation decay model is constructed. In the multidimensional response landscape map, peak regions represent the most typical and significantly deviated response patterns under specific nutrient stress conditions. The set of multidimensional fragmentation kinetic feature vectors corresponding to the peak regions is extracted from the map and used as the standard fragmentation kinetic template. It should be understood that the standard fragmentation kinetic template refers to a known type of fragmentation behavior that has formed a clear clustering pattern in the multidimensional response landscape and can be captured by the map topology structure of step S3. For each set of standard fragmentation kinetic templates, the median of each dimension is taken as the center value of the reference interval for that dimension, and the interquartile range of each dimension is multiplied by a preset expansion coefficient as the half-width of the reference interval, thereby determining the upper and lower bounds of the reference interval for each dimension. The reference intervals of the above dimensions together enclose a hypercube region in multidimensional space, which constitutes the standard fragmentation decay model.

[0034] Subsequently, the multidimensional pyrolysis kinetic feature vectors of all microreactors generated in this round of screening are iterated. For each feature vector to be judged, its values ​​in each dimension are compared with the reference intervals of the corresponding dimensions in the standard pyrolysis decay model, checking whether each of the four dimensions—the pyrolysis acceleration point, the maximum pyrolysis rate point, the plateau initiation point, and the average decay rate during the explosive pyrolysis period—falls within its corresponding reference interval. If at least one dimension falls outside its corresponding reference interval, the microreactor to which the deviating feature vector belongs is marked as a candidate hidden pyrolysis microreactor. The advantage of this dimension-by-dimensional judgment strategy is that even if the overall pyrolysis profile of a microreactor is similar to the standard template, as long as it exhibits an anomaly at a certain key time point or pyrolysis intensity, it will not be masked by the normality of other dimensions, thus maximizing the capture of subtle deviations in pyrolysis behavior.

[0035] Candidate cryptic lysis microreactors still need to undergo variant exclusion screening. A microreactor may be marked as a candidate for cryptic lysis for two possible reasons: first, the bacteriophage within the microreactor genuinely exhibits a novel lysis behavior pattern, i.e., a true cryptic lysis event; second, the bacteriophage within the microreactor is a naturally occurring variant of a known lysis behavior response pattern, i.e., a variant of a known lysis behavior pattern. The significance of distinguishing between these two is that the goal of this invention is to discover bacteriophages with novel lysis behavior types, rather than collecting variants of known types. The specific distinction method is as follows: first, extract known lysis behavior response patterns with the same stress perturbation pattern label as the candidate microreactor from the multidimensional response landscape map. These known lysis behavior response patterns have already been labeled and stored as part of the map in step S3; then calculate the cosine similarity between the feature vector of the candidate microreactor and the representative vector of the known lysis behavior response pattern. Cosine similarity measures the directional closeness of two vectors, with values ​​ranging from -1 to 1. The closer the value is to 1, the more consistent the directions. When the cosine similarity is below a preset mutation threshold, it indicates that the feature vector of the candidate microreactor differs substantially in direction from the known pyrolysis response pattern. The candidate microreactor is then classified as non-mutant, officially labeled as a latent pyrolysis event. If the cosine similarity is above the threshold, the candidate microreactor is classified as a variant of the known pyrolysis behavior pattern and is excluded.

[0036] In the second stage, pulsed dielectrophoresis-driven recovery is performed on the officially tagged cryptic lysis events. First, the nutrient stress orthogonal code corresponding to each cryptic lysis event is obtained. Based on the spatial mapping relationship between the nutrient stress orthogonal code and the microarray chip lattice point position, the physical lattice point coordinates (row and column positions) of the microreactor containing the cryptic lysis event on the microarray chip are retrieved. These physical lattice point coordinates are used as the target address and written into the target address register of the pulsed dielectrophoresis-driven sequence, completing the coordinate transformation from data space to physical space.

[0037] Subsequently, the dielectrophoresis drive is initiated. A non-uniform alternating electric field is applied to the microarray chip region pointed to by the target address register. The non-uniform alternating electric field can be generated by using interdigitated electrodes or parallel electrode arrays integrated on the chip substrate, achieving non-uniformity in the spatial distribution of the electric field by applying phase and amplitude differences to different electrodes. Under the action of the non-uniform alternating electric field, the water-in-oil droplets generate dielectrophoresis force due to the difference in dielectric constant and conductivity between them and the surrounding continuous oil phase. By adjusting the frequency of the applied electric field, the dielectrophoresis force experienced by the droplets can be selectively made to be positive or negative, thereby controlling the direction of droplet migration along the electric field gradient. In this embodiment, the electric field frequency is adjusted so that the target droplets experience a positive dielectrophoresis force pointing towards the preset recovery channel, driving the droplets to migrate from their original array point position along the electric field gradient direction to the inlet of the recovery channel. A negative pressure interface is set at the outlet of the recovery channel, and a weak negative pressure is applied to draw the droplets that have migrated to the recovery port into the collection chamber. After the microdroplets aggregate and demulsify within the collection chamber, they release the bacteriophages and host bacterial remnants contained therein, resulting in an enriched solution containing unconventional bacteriophages with liquid-phase lysis activity.

[0038] Thus, the entire process from digital identification of deviations in lysis behavior to physical recovery of the target microreactor has been completed. Based on the fluorescence lysis kinetics in a homogeneous droplet environment, the precise capture of unconventional bacteriophages, including non-plaque-forming bacteriophages, has been achieved.

[0039] Specifically, S5 is: The multidimensional response landscape map generated in this round of screening is merged with the multidimensional response landscape map accumulated in previous rounds. Using the stress perturbation mode label as the index, the lysis behavior response modes obtained under the same or similar nutrient stress conditions in different rounds are aligned. The drift trend of each dimension of the lysis dynamics feature vector under each stress perturbation mode label is calculated between multiple rounds to form a cumulative response landscape map library containing time-dimensional evolution information. Density clustering is performed on the distribution areas of hidden rift events that have been marked in the cumulative response landscape map library to identify the nutrient stress condition subspace corresponding to the high incidence area of ​​hidden rift events. The boundary contour of the nutrient stress condition subspace is used as the known positive enrichment region. At the same time, the nutrient stress condition regions in the cumulative response landscape map library that have not been covered by any round of sampling are extracted as the unexplored nutrient stress condition candidate set. Using the known positive enrichment region as the query center, calculate the path distance between each candidate nutrient stress condition in the unexplored nutrient stress condition candidate set and the known positive enrichment region on the hierarchical similarity tree structure, and use the reciprocal of the path distance as the proximity weight of the candidate nutrient stress condition to the known positive enrichment region. For each candidate nutrient stress condition, a surrogate prediction model is constructed based on the distribution of pyrolysis behavior response patterns of sampled nutrient stress conditions in the cumulative response landscape map library. The prediction uncertainty of the pyrolysis behavior response pattern output by the surrogate prediction model under the candidate nutrient stress condition is calculated. The prediction uncertainty is multiplied by the proximity weight to obtain the information gain evaluation value of the candidate nutrient stress condition. The information gain evaluation values ​​of all candidate nutritional stress conditions are sorted in descending order. A preset number of candidate nutritional stress conditions with the highest information gain evaluation values ​​are selected. Based on the types and concentration levels of stressor factors contained in each selected candidate nutritional stress condition, the stress-inducing culture medium combination formula required for the next round of screening is generated by reverse decoding.

[0040] It should be noted that the implementation of this step can be divided into four stages: construction of a cumulative response landscape library, delineation of candidate spaces, evaluation of information gain, and reverse decoding of culture medium formulation.

[0041] In the first stage, the multidimensional response landscape maps generated in this round of screening are merged with those accumulated from previous rounds. Each round of screening outputs a multidimensional response landscape map carrying a round identifier. During merging, the stress perturbation mode label is used as an index to align the lysis behavior response patterns falling under the same or similar stress perturbation mode labels in different rounds. That is, the values ​​of each dimension of the lysis dynamics feature vector under the same label in different rounds are arranged in round order to form a cross-round evolution sequence of lysis behavior under that label. Based on this sequence, the drift trend of each dimension value between multiple rounds is calculated. The drift trend can be characterized by the slope of linear regression or the difference between consecutive rounds, reflecting the directional change of phage population lysis behavior under this nutrient stress condition as the screening generations progress. The merged, aligned, and appended map set constitutes a cumulative response landscape map library containing time-dimensional evolutionary information. It not only records the spatial distribution of lysis behavior under various nutrient stress conditions but also records the dynamic evolution trajectory of these behavior patterns as the screening generations progress.

[0042] In the second stage, known positive enrichment regions and unexplored candidate sets are delineated based on the cumulative response landscape map library. First, density clustering is performed on the distribution areas of hidden rift events labeled in all rounds of the map library. Density clustering uses the coordinates of the nutrient stress conditions corresponding to each hidden rift event as input points, employing a density-based spatial clustering algorithm to identify several high-density regions formed by the aggregation of hidden rift events in the nutrient stress condition space. The boundary contours of each high-density region are extracted using methods such as convex hull or kernel density estimation, thus defining the corresponding nutrient stress condition subspace. The union of these subspaces is taken as the known positive enrichment region. The known positive enrichment region represents the range of nutrient stress conditions where hidden rift events are most likely to be discovered under the current knowledge state. Simultaneously, nutrient stress condition regions in the cumulative response landscape map library that have not yet been covered by any round of sampling are extracted. The set of nutrient stress conditions corresponding to these regions constitutes the unexplored nutrient stress condition candidate set. Each element in the candidate set represents a nutrient stress condition that has not yet been experimentally verified.

[0043] In the third stage, information gain is evaluated for each candidate nutrient stress condition in the unexplored candidate set. The information gain is determined by two factors: proximity weight and prediction uncertainty. The proximity weight is calculated using a known positive enrichment region as the query center, measuring the path distance between the candidate nutrient stress condition and the known positive enrichment region on a hierarchical similarity tree structure. This hierarchical similarity tree structure, constructed in step S3, reflects the similarity between nutrient stress conditions. A shorter path distance indicates a closer similarity between the candidate condition and the known positive enrichment region's nutrient stress characteristics, resulting in a higher prior probability of discovering a hidden lysis event under that candidate condition. The reciprocal of this path distance is taken as the proximity weight, giving higher weights to closer candidate conditions. The prediction uncertainty is calculated based on the construction of a surrogate prediction model. The surrogate prediction model employs a conventional Gaussian process regression model. It uses the coordinates of sampled nutrient stress conditions and the corresponding comprehensive deviation of pyrolysis behavior from the cumulative response landscape database as training input and output, respectively. A kernel function defines the covariance structure among the sampled points, establishing a probabilistic mapping from nutrient stress conditions to pyrolysis behavior responses. For any candidate nutrient stress condition, the surrogate prediction model outputs the predicted distribution of the pyrolysis behavior response pattern under that condition, using the variance of this distribution as the prediction uncertainty. A higher prediction uncertainty indicates a greater lack of current knowledge regarding pyrolysis behavior under that candidate condition, and thus a greater potential information gain from sampling it. Multiplying the prediction uncertainty by the proximity weight yields the information gain evaluation value for the candidate nutrient stress condition.

[0044] All candidate nutrient stress conditions are sorted in descending order of their information gain evaluation values, and a preset number of candidate conditions with the highest information gain evaluation values ​​are selected to proceed to the next round of screening. The preset number can be determined based on the dot size of the microarray chip, for example, it can be consistent with the dot size of the nutrient stress orthogonal coding array.

[0045] In the fourth stage, inverse decoding is performed on the selected candidate nutrient stress conditions to restore their information into stress-inducing culture medium formulations that can be directly used for experimental preparation. Each candidate nutrient stress condition corresponds to a coordinate point in the nutrient stress condition space, which has a clear position in a hierarchical similarity tree structure. The coordinate point is used to trace back along the tree structure to the concentration dimension of each stressor factor, extracting the type and concentration level of each stressor factor. For example, if a candidate condition is located as simultaneously perturbing both carbon metabolism and cell wall synthesis pathways at a moderate level, inverse decoding can restore it to a combination of stressor factor A (e.g., 2-deoxyglucose) at a moderate concentration and stressor factor B (e.g., fosfomycin) at a moderate concentration. This combination is then used to generate the specific stress-inducing culture medium combination formulation required for the next round of screening, according to the culture medium preparation rules described in S1. This formulation directly serves as the input for preparing the nutrient stress orthogonal encoding array in step S1, driving the closed-loop execution of the next round of screening.

[0046] Through the iterative process of the above four stages, the screening path can autonomously approach the nutrient stress condition region enriched by cryptic lysis events in each round of screening, so that the discovery efficiency of rare phages can be continuously improved with each round.

[0047] In a preferred embodiment, it further includes: From the multidimensional response landscape map, multidimensional pyrolysis dynamics feature vectors corresponding to different nutrient stress conditions under the same branch are extracted according to the stress disturbance mode label. Each dimension of each multidimensional pyrolysis dynamics feature vector is used as an independent coordinate axis to construct a high-dimensional point cloud data set representing the pyrolysis behavior state. Using the stress condition distance defined by the hierarchical similarity tree structure of nutrient stress conditions as a metric parameter, distance-based nested filter construction is performed on the high-dimensional point cloud data set. At each stage of the filter evolution, the highest-dimensional void structure composed of all simplexes in each stage is identified, and the set of the highest-dimensional void structures in each stage is taken as the high-dimensional simplex skeleton sequence of the rupture behavior state. In the high-dimensional simplex skeleton sequence, the minimum number of simplex variations required for morphological transformation between adjacent simplex skeleton stages through simplex addition or deletion operations is calculated. The ratio of the minimum number of simplex variations corresponding to each morphological transformation to the increment of the metric parameter is defined as the structural sensitivity coefficient of skeleton morphological transformation. Among all homology classes of the high-dimensional topological filter structure, a specific homology class is selected whose generator is composed of simplexes contained in the high-dimensional simplex skeleton sequence. The change trajectory of each dimensional Betti number of the specific homology class is tracked from generation to annihilation, and spliced ​​to form the dimension transfer spectrum of the specific homology class. The dimensional migration spectrum is subjected to mutation detection to locate the critical event of dimensional migration in the specific homology class where the high-dimensional Betty number suddenly drops and the low-dimensional Betty number suddenly increases. The metric parameter value corresponding to the critical event is used as the topological phase transition threshold for qualitative mutation of the cleavage behavior, and is inversely mapped back to the nutrient stress condition space to define the nutrient stress phase transition boundary.

[0048] It should be noted that a high-dimensional point cloud dataset is constructed by using the nutrient stress conditions under the same stress perturbation mode label branch as the horizontal axis and the dimensions of their corresponding multidimensional fragmentation dynamics feature vectors as the vertical axis. Distance-based nested filtering is performed on the point cloud. The core of this method is to use the stress condition distance defined by the hierarchical similarity tree structure in S3 as the metric parameter, gradually increasing the distance threshold from zero to the maximum value covering all point pairs. At each threshold stage, point pairs with distances less than or equal to the threshold are connected by edges, three pairs of connected points are filled with triangles, and k+1 pairs of connected points are filled with a k-dimensional simplex. This forms a nested sequence of simplex complexes that gradually expands with the threshold, i.e., the nested filtering structure. In each stage of the simplex complex, the highest-dimensional void structure is identified. A void refers to a region enclosed by a simplex but not filled inside. For example, three edges connected end-to-end without filling with triangles constitute a one-dimensional annular void, and four triangles enclosing a tetrahedral face without filling the tetrahedron constitute a three-dimensional void. The set of highest-dimensional voids at each stage constitutes a high-dimensional simplex skeleton sequence of the fracture behavior state, which characterizes the key outline of the evolution of point cloud topology with distance parameters.

[0049] Subsequently, specific homology classes are selected from all homology classes generated by the filter flow. A homology class corresponds to a void that is generated, persists, and eventually annihilates at a certain stage of the filter flow. The selection criterion is that the generators of the homology class (i.e., the most basic simplex combination enclosing the void) must consist of simplexes contained in the aforementioned high-dimensional simplex skeleton sequence. This criterion limits the tracking scope to voids that play a structural supporting role in the skeleton evolution, excluding short-lived topological features generated by noise points. For the selected specific homology classes, the changes in the Betti numbers of each dimension are traced throughout their entire process from generation to annihilation, and these are spliced ​​together to form a dimension migration spectrum. The Betti number reflects the topological dimension type of the void: zero-dimensional corresponds to connected components, one-dimensional corresponds to annular voids, and two-dimensional corresponds to spherical voids. The dimension migration spectrum records the transition path of the topological expression dimension of the same void structure during the filter flow evolution.

[0050] Mutation detection is performed on the dimensionality transition spectrum to locate the critical event where a sudden drop in high-dimensional Betty number occurs while a sudden increase in low-dimensional Betty number occurs for the first time. For example, a dimensionality reduction reconstruction where the two-dimensional Betty number drops to zero while the one-dimensional Betty number increases simultaneously. The metric parameter value corresponding to this critical event is the topological phase transition threshold. This threshold is inversely mapped back to the nutrient stress condition space along a hierarchical similarity tree structure. Corresponding nodes are labeled on the tree branches, and connecting the threshold nodes defines the nutrient stress phase transition boundary. On both sides of this boundary, the topological structure of phage lysis behavior undergoes fundamental reconstruction, providing clear conditional navigation for targeted screening. This allows the screening strategy to bypass nutrient stress regions with highly homogeneous lysis behavior and prioritize the deployment of screening resources in the behavioral transition zones near the topological phase transition boundary, further improving the screening accuracy of rare phages.

[0051] In a preferred embodiment, it further includes: Extract the multidimensional lysis dynamics feature vectors corresponding to all sampled nutrient stress conditions from the multidimensional response landscape map, calculate the feature space distance between any two feature vectors, and construct a diffusion affinity matrix that reflects the local nearest neighbor similarity of lysis behavior state based on the feature space distance; The diffusion affinity matrix is ​​normalized row by row to generate a diffusion transition matrix that describes the state transition probability of cleavage behavior between nutrient stress conditions. The eigenvalues ​​and eigenvectors of the diffusion transition matrix are calculated, and the dominant eigenvectors are selected as diffusion coordinates. Each nutrient stress condition is mapped to a low-dimensional diffusion manifold spanned by the diffusion coordinates. The diffusion coordinate sequence of each nutrient stress condition constitutes the potential evolution trajectory of cleavage behavior. On the low-dimensional diffusion manifold, based on the stress condition adjacency relationship defined by the hierarchical similarity tree structure, the diffusion coordinates of adjacent nutrient stress conditions are locally parameterized along the potential evolution trajectory to generate a discrete arc length parameter sequence describing the trajectory geometry. Based on the discrete arc length parameter sequence, the discrete manifold curvature at each point on the potential evolution trajectory is calculated. The discrete manifold curvature characterizes the degree of drastic change in the direction of the rifting behavior at that point as the nutrient stress condition changes, forming a profile of manifold curvature change along the trajectory. In the manifold curvature variation profile, the curvature peaks that form local maxima of discrete manifold curvature are located, and the diffusion coordinates corresponding to the curvature peaks are mapped back to the nutrient stress condition space. The resulting nutrient stress conditions are marked as the brittle nutrient stress nodes that are most sensitive to nutrient disturbance response in terms of pyrolysis behavior.

[0052] It should be noted that a high affinity is assigned only when two eigenvectors are neighbors in the feature space; the affinity between more distant points approaches zero, thus reflecting the local neighborhood similarity of the fragmentation behavior states. After row-normalizing the diffusion affinity matrix, a diffusion transition matrix describing the probability of fragmentation behavior state transitions between nutrient stress conditions is obtained. This matrix essentially defines a Markov random walk process in the nutrient stress condition space. The eigenvalues ​​and eigenvectors of the diffusion transition matrix are solved, and the top few dominant eigenvectors (excluding trivial eigenvectors) are selected as diffusion coordinates. Since the diffusion coordinates capture the main variation directions of the data's intrinsic geometric structure, the mapping positions of each nutrient stress condition on the low-dimensional diffusion manifold spanned by these diffusion coordinates constitute the potential evolutionary trajectory of the fragmentation behavior.

[0053] On a low-dimensional diffusion manifold, based on the adjacency relationship of stress conditions defined by a hierarchical similarity tree structure, the diffusion coordinates of adjacent nutrient stress conditions are locally parameterized along the potential evolution trajectory to generate a discrete arc length parameter sequence, which describes the geometric extension morphology of the trajectory. Based on the discrete arc length parameter sequence, the discrete manifold curvature at each point on the trajectory is calculated. The discrete manifold curvature reflects the degree of abrupt change in the evolution direction of cleavage behavior with changes in nutrient stress conditions; the larger the curvature value, the more drastic the change in the evolution direction of cleavage behavior near that nutrient stress condition. The calculated curvature values ​​at each point are arranged along the trajectory to form a manifold curvature change profile. In the manifold curvature change profile, the curvature peaks that constitute local maxima of the discrete manifold curvature are located. After the diffusion coordinates corresponding to these peaks are mapped back to the nutrient stress condition space, the resulting nutrient stress conditions are the brittle nutrient stress nodes that are most sensitive to nutrient disturbances in the response of cleavage behavior. Near this node, even minor adjustments to nutrient stress conditions can cause a significant shift in the direction of lysis behavior evolution, making it a critical condition region worthy of close attention in screening strategies. By prioritizing the deployment of the next round of screening in the culture medium formulation design around the brittle nutrient stress node, a broader space of lysis behavior diversity can be explored with less screening rounds, improving the exploration efficiency and discovery coverage of the screening strategy.

[0054] The above are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A high-throughput screening method for single phages based on microfluidics-culture omics coupling, characterized in that, Includes the following steps: S1. Construct a fluorescent pyrolysis microreactor array and couple the microreactor array with the microarray chip, while assigning a unique orthogonal nutrient stress code to each microreactor. S2. Perform time-series fluorescence imaging on the co-culture system in each microreactor, perform nonlinear dynamic analysis on the fluorescence decay signal of each microreactor, and extract multidimensional fragmentation kinetic feature vectors. S3. Associate the multidimensional lysis kinetics feature vector with the orthogonal encoding of the nutrient stress of the corresponding microreactor to generate a multidimensional response landscape map describing the phage lysis behavior as a function of nutrient conditions. S4. Identify occult lysis events that deviate from the preset standard lysis decay model as potential positives, and perform pulsed dielectrophoresis-driven recovery on the microreactor based on the nutrient stress orthogonal encoding corresponding to the occult lysis events to obtain unconventional phages. S5. By gathering the multidimensional response landscape maps accumulated from multiple rounds of screening as prior knowledge, and adopting an active learning strategy with the goal of maximizing information gain, the stress-inducing culture medium combination required for the next round of screening is inferred and generated.

2. The high-throughput screening method for single phages according to claim 1, characterized in that, Specifically, S1 is: A baseline metabolic model of the host bacteria is obtained. Based on the sensitivity coefficients of each metabolic pathway in the baseline metabolic model, multiple stressor factors targeting different metabolic nodes are selected from a pre-set stressor factor library. Based on the orthogonal experimental design principle, a nutritional stress orthogonal coding table composed of different concentration levels of the stressor factors is generated, where each group of codes corresponds to a unique stress-inducing culture medium formulation. Using the preset matrix positions of the microarray chip as an index, each group of codes in the nutrient stress orthogonal coding table is sequentially assigned to each matrix position according to the spatial mapping relationship. The stress-inducing culture medium corresponding to each group of codes is pre-fixed at the corresponding matrix position to form a nutrient stress orthogonal coding array with position coding information. In a microfluidic chip, a suspension of phage to be screened, a suspension of host bacteria, and a fluorescent nucleic acid dye are mixed in a preset ratio to form an aqueous phase, which is then injected into a flow focusing structure together with an oil phase containing a surfactant. By adjusting the flow rate ratio of the two phases, water-in-oil droplets are generated, so that the size of the droplets matches the lattice size of the microarray chip. The generated water-in-oil microdroplets are introduced into the surface of the prepared nutrient stress orthogonal coding array, and each microdroplet falls into and is fixed at the corresponding array point. Through the diffusion reconstruction of the aqueous phase inside the microdroplet and the pre-cured stress-induced culture medium on the array point, the nutrient stress orthogonal coding carried by the array point is assigned to the corresponding microdroplet, forming a co-culture system in which each microreactor has a unique nutrient stress condition identifier.

3. The high-throughput screening method for single phages according to claim 1, characterized in that, Specifically, S2 is: For each co-cultured microreactor coupled to the nutrient stress orthogonal coding array, multiple frames of fluorescence images are acquired at preset time intervals. The fluorescence intensity value of each microreactor is extracted from each frame of the image, and an original fluorescence decay sequence arranged along the time axis is constructed for each microreactor. At the same time, the nutrient stress orthogonal coding corresponding to each microreactor is recorded. Detrended fluctuation analysis was performed on the original fluorescence decay sequence of each microreactor. By calculating the local fluctuation function of the sequence, non-biological signal trends caused by instrument noise and photobleaching were filtered out to obtain the corrected fluorescence decay sequence. Based on the corrected fluorescence decay sequence, the first-order difference sequence of each sampling point was calculated. Based on the sign change and amplitude jump of the first-order difference sequence, the corrected fluorescence decay sequence is adaptively divided into three dynamic intervals: hysteresis period, burst fragmentation period, and decay plateau period, and a fragmentation stage division scheme containing the boundary time of each interval and the fluorescence decay rate within the interval is generated. Based on the fragmentation stage division scheme, the moment when the first spontaneous fragmentation threshold in the first-order differential sequence exceeds the preset fragmentation threshold is identified as the fragmentation acceleration point, the peak moment of the first-order differential sequence during the burst fragmentation period is extracted as the maximum fragmentation rate point, and the moment when the fluorescence intensity value tends to a steady state during the decay plateau period is identified as the plateau entry point. By splicing together the pyrolysis acceleration point, maximum pyrolysis rate point, plateau inlet point, and average decay rate during the explosive pyrolysis period of the same microreactor, a multidimensional pyrolysis kinetic feature vector of the corresponding microreactor is formed.

4. The high-throughput screening method for single phages according to claim 1, characterized in that, Specifically, S3 is: Using the orthogonal coding of nutrient stress as the primary key, the multidimensional pyrolysis kinetic feature vectors of each microreactor are grouped according to the coding, a mapping table between each orthogonal coding of nutrient stress and the corresponding feature vector is established, and the statistical distribution parameters of each feature vector in each coding group are calculated. The statistical distribution parameters include the central tendency measure and the dispersion measure of the values ​​of each dimension. All orthogonal codes of nutritional stress in the mapping table are hierarchically clustered according to the types and concentration levels of the stressor factors they contain, generating a hierarchical similarity tree structure of nutritional stress conditions, so that codes with similar metabolic pathway perturbation intensities are grouped into the same branch, and each branch is assigned a corresponding stress perturbation pattern label. Based on the hierarchical relationship of the stress disturbance mode labels, cross-condition comparative analysis is performed on the multidimensional lysis kinetic feature vectors between different coding groups. By calculating the significance of the differences between the feature vectors of each dimension between groups, the deviation direction and magnitude of the lysis kinetic behavior under each stress disturbance mode label compared with the baseline culture medium conditions are identified, forming a lysis behavior deviation matrix. Each row of the rift behavior deviation matrix is ​​mapped to a behavior response surface grid point in the nutrient stress condition space. The nutrient stress condition distance defined by the hierarchical similarity tree structure is used as the surface interpolation skeleton. Conformal interpolation is performed between adjacent grid points to generate a multi-dimensional response landscape map with full coverage. The multidimensional response landscape map is analyzed for topological features to extract peak, valley and saddle point regions. Each region is then labeled based on the stress perturbation pattern label to obtain a set of cleavage behavior response patterns associated with specific metabolic pathway perturbations.

5. The high-throughput screening method for single phages according to claim 4, characterized in that: Each coordinate point of the multidimensional response landscape map corresponds to a nutrient stress condition, and the height dimension of the map represents the overall deviation of the pyrolysis kinetics.

6. The high-throughput screening method for single phages according to claim 1, characterized in that, Specifically, S4 is: The multidimensional pyrolysis dynamics feature vector corresponding to the peak region in the multidimensional response landscape map is called as the standard pyrolysis dynamics template. The values ​​of each dimension of each set of standard pyrolysis dynamics templates are set as the center value of the reference interval of the corresponding dimension. The upper and lower bounds of the reference interval are determined based on the dispersion measure of each dimension, and a standard pyrolysis decay model enclosed by the multidimensional reference interval is constructed. The multidimensional pyrolysis kinetic feature vectors of all microreactors generated in this round of screening are traversed. The values ​​of each dimension of each feature vector to be judged are compared with the reference interval of the corresponding dimension in the standard pyrolysis decay model. At least one deviation feature vector whose value falls outside the corresponding reference interval is identified, and the microreactor to which the deviation feature vector belongs is marked as a candidate hidden pyrolysis microreactor. For each candidate hidden pyrolysis microreactor, a known pyrolysis behavior response pattern with the same stress perturbation mode label as the candidate hidden pyrolysis microreactor is extracted from the multidimensional response landscape map. By comparing the angle between the feature vector of the candidate hidden pyrolysis microreactor and the vector of the known pyrolysis behavior response pattern, it is determined whether the candidate hidden pyrolysis microreactor is a variant of the known pyrolysis behavior pattern. If it is not a variant, it is officially marked as a hidden pyrolysis event. Obtain the nutrient stress orthogonal code corresponding to each hidden pyrolysis event. Based on the spatial mapping relationship between the nutrient stress orthogonal code and the microarray chip lattice point position, look up the physical lattice point coordinates of the microreactor where the hidden pyrolysis event is located, and write the physical lattice point coordinates into the target address register of the pulse dielectrophoresis drive sequence. A non-uniform alternating electric field is applied to the microarray chip region pointed to by the target address register. By adjusting the electric field frequency, a difference in dielectric force is generated between the microreactor and the surrounding continuous phase, driving the microreactor to migrate along the electric field gradient direction to a preset recovery channel. A negative pressure is applied at the recovery port to draw the microreactor that has migrated to the recovery port into the collection chamber, thereby obtaining an enriched solution containing unconventional phages with the liquid phase lysis activity.

7. The high-throughput screening method for single phages according to claim 6, characterized in that, The step of determining whether a candidate hidden pyrolysis microreactor is a variant of a known pyrolysis behavior pattern by comparing the angle between the feature vector of the candidate hidden pyrolysis microreactor and the vector of the known pyrolysis behavior response pattern is as follows: Calculate the cosine similarity between the multidimensional pyrolysis kinetic feature vector of the candidate hidden pyrolysis microreactor and the representative vector of the known pyrolysis behavior response mode under the same stress perturbation mode label. When the cosine similarity is lower than the preset mutation determination threshold, the candidate hidden pyrolysis microreactor is determined to be a variant of the known pyrolysis behavior mode.

8. The high-throughput screening method for single phages according to claim 1, characterized in that, Specifically, S5 is: The multidimensional response landscape map generated in this round of screening is merged with the multidimensional response landscape map accumulated in previous rounds. Using the stress perturbation mode label as the index, the lysis behavior response modes obtained under the same or similar nutrient stress conditions in different rounds are aligned. The drift trend of each dimension of the lysis dynamics feature vector under each stress perturbation mode label is calculated between multiple rounds to form a cumulative response landscape map library containing time-dimensional evolution information. Density clustering is performed on the distribution areas of hidden rift events that have been marked in the cumulative response landscape map library to identify the nutrient stress condition subspace corresponding to the high incidence area of ​​hidden rift events. The boundary contour of the nutrient stress condition subspace is used as the known positive enrichment region. At the same time, the nutrient stress condition regions in the cumulative response landscape map library that have not been covered by any round of sampling are extracted as the unexplored nutrient stress condition candidate set. Using the known positive enrichment region as the query center, calculate the path distance between each candidate nutrient stress condition in the unexplored nutrient stress condition candidate set and the known positive enrichment region on the hierarchical similarity tree structure, and use the reciprocal of the path distance as the proximity weight of the candidate nutrient stress condition to the known positive enrichment region. For each candidate nutrient stress condition, a surrogate prediction model is constructed based on the distribution of pyrolysis behavior response patterns of sampled nutrient stress conditions in the cumulative response landscape map library. The prediction uncertainty of the pyrolysis behavior response pattern output by the surrogate prediction model under the candidate nutrient stress condition is calculated. The prediction uncertainty is multiplied by the proximity weight to obtain the information gain evaluation value of the candidate nutrient stress condition. The information gain evaluation values ​​of all candidate nutritional stress conditions are sorted in descending order. A preset number of candidate nutritional stress conditions with the highest information gain evaluation values ​​are selected. Based on the types and concentration levels of stressor factors contained in each selected candidate nutritional stress condition, the stress-inducing culture medium combination formula required for the next round of screening is generated by reverse decoding.

Images