CHO cell exogenous gene expression prediction screening method and system

Abstract

The application provides a CHO cell exogenous gene expression prediction screening method and system, relates to the technical field of prediction screening, divides a core area and a peripheral area, calculates an order parameter, extracts a driving and braking gene set, performs nonlinear fitting on a state-dependent stochastic ordinary differential equation through a dynamic order parameter sequence and a driving and braking gene expression profile, estimates driving strength indexes and braking strength indexes, and calculates energy efficiency ratios, performs characteristic substitution and calculation on a prediction model, compares exogenous protein prediction yield with a threshold value, selects high-yield clones into process development, and according to abnormal changes of the energy efficiency ratios or the order parameter plateau, the culture conditions are warned or optimized, and the optimal process parameters are obtained. The method realizes accurate prediction of the final yield of the clone based on the dynamic data of the transcriptome in the early culture stage, significantly shortens the screening period, improves the prediction accuracy, and reduces the waste of manpower and resources caused by invalid clone screening.

Description

Technical Field

[0001] This invention proposes a method and system for predicting and screening the expression of exogenous genes in CHO cells, which relates to the field of prediction and screening technology, specifically to the field of predicting and screening the expression of exogenous genes in CHO cells. Background Technology

[0002] Chinese hamster ovary (CHO) cells are the primary host for the production of therapeutic recombinant proteins, with over 70% of biopharmaceuticals produced globally using CHO cells. However, rapid screening and high-yield prediction of CHO cell lines have become a core bottleneck in bioprocess development. Traditional screening processes require 6-8 weeks, are labor-intensive, have low throughput, and often result in a 30%-70% yield drop in early high-yield clones during scale-up culture, leading to significant waste of preliminary work.

[0003] To shorten the cycle or improve prediction accuracy, existing technologies mainly fall into four categories: one is based on the copy number of exogenous genes and integration sites, but the linear correlation between copy number and yield is poor (R0). 2 Typically <0.5), high copy numbers may actually cause transcriptional repression; secondly, methods based on the expression levels of single or a few marker genes, but changes in the expression of single genes are noisy and inconsistent in correlation; thirdly, methods based on metabolic throughput analysis require expensive isotope labeling and are not suitable for high-throughput screening; fourthly, methods based on machine learning prediction using omics data have black-box characteristics, poor generalization ability, and still require long-term culture to obtain tag data. The common shortcomings of existing technologies are that they ignore the dynamic evolution of the overall cellular transcriptional regulatory order, cannot quantify metabolic burden or the balance between order and burden, lack early predictable kinetic parameters, and cannot distinguish between driving and braking factors. Summary of the Invention

[0004] This invention provides a method and system for predicting and screening the expression of exogenous genes in CHO cells, in order to solve the above-mentioned problems: This invention proposes a method and system for predicting and screening the expression of exogenous genes in CHO cells, the method comprising: S1. Using transcriptome data, the stability wall is located, the core region and the outer region are divided, the order parameters are calculated and the driver and brake gene sets are extracted to obtain the dynamic order parameter sequence and the expression profile of the driver and brake genes. S2. By using dynamic order parameter sequences and driving and braking gene expression profiles, nonlinear fitting is performed on state-dependent stochastic ordinary differential equations to estimate driving intensity index and braking intensity index and calculate energy efficiency ratio, thereby obtaining key dynamic characteristic parameters. S3. By combining key kinetic characteristic parameters with the early expression level of exogenous genes, the prediction model is substituted and calculated to obtain the predicted yield of exogenous protein. S4. By comparing the predicted yield of exogenous proteins with the threshold, high-yield clones are screened for process development. Based on the abnormal changes in the energy efficiency ratio or order parameter plateau, the culture conditions are given early warning or optimized to obtain the optimal process parameters.

[0005] Furthermore, the system includes: The parameter expression profile acquisition module is used to locate the stability wall through transcriptome data, divide the core region and the peripheral region, calculate the order parameters and extract the driver and brake gene sets to obtain the dynamic order parameter sequence and driver and brake gene expression profile. The dynamic feature acquisition module is used to perform nonlinear fitting on state-dependent stochastic ordinary differential equations by using dynamic order parameter sequences and driving and braking gene expression profiles, estimate driving intensity index and braking intensity index and calculate energy efficiency ratio to obtain key dynamic feature parameters. The yield prediction module is used to input and calculate the predicted yield of exogenous protein by combining key kinetic characteristic parameters with the early expression level of exogenous genes into the prediction model. The process optimization module is used to screen high-yield clones for process development by comparing the predicted yield of exogenous proteins with thresholds, and to provide early warning or optimization of culture conditions based on abnormal changes in energy efficiency ratio or order parameter plateau phase, so as to obtain the optimal process parameters.

[0006] The beneficial effects of this invention are as follows: This method solves the technical problem that existing technologies cannot accurately predict the final yield of exogenous proteins in the early stages of culture, and overcomes the defect in traditional screening processes where the yield of high-yielding clones in the early stages significantly decreases after scale-up culture. It achieves quantitative capture of the dynamic evolution of the overall transcriptional regulatory order of the cell, distinguishing between driving factors that promote expression and braking factors that restrict expression, and revealing the nonlinear balance between the two. Compared with existing methods that rely on copy number, single marker genes, or metabolic flux, this invention significantly shortens the clone screening cycle, significantly improves the accuracy of early prediction, and reduces the ineffective workload and resource waste caused by later yield decline. Attached Figure Description

[0007] Figure 1 This is a schematic diagram of a method for predicting and screening the expression of exogenous genes in CHO cells. Detailed Implementation

[0008] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0009] In one embodiment of the present invention, a method and system for predicting and screening exogenous gene expression in CHO cells are proposed, the method comprising: S1. Using transcriptome data, the stability wall is located, the core region and the outer region are divided, the order parameters are calculated and the driver and brake gene sets are extracted to obtain the dynamic order parameter sequence and the expression profile of the driver and brake genes. S2. By using dynamic order parameter sequences and driving and braking gene expression profiles, nonlinear fitting is performed on state-dependent stochastic ordinary differential equations to estimate driving intensity index and braking intensity index and calculate energy efficiency ratio, thereby obtaining key dynamic characteristic parameters. S3. By combining key kinetic characteristic parameters with the early expression level of exogenous genes, the prediction model is substituted and calculated to obtain the predicted yield of exogenous protein. S4. By comparing the predicted yield of exogenous proteins with the threshold, high-yield clones are screened for process development. Based on the abnormal changes in the energy efficiency ratio or order parameter plateau, the culture conditions are given early warning or optimized to obtain the optimal process parameters.

[0010] The specific content includes: The target CHO cell lines (endogenous assessment can directly use existing single-cell or bulk RNA / seq data; exogenous gene expression prediction requires cells transfected with the target gene expression cassette, such as antibody light and heavy chains) are seeded in suitable culture media (e.g., CD CHO, Dynamis) and cultured in shake flasks or 24-well plates. For endogenous data, time-series scRNA / seq data from public databases (e.g., GSE217637, 4 time points D0, D5) can be used for proof-of-concept. At least 4, preferably 6 or more, independent time points (e.g., D0, D1, D3, D5, or sampling at 24-hour intervals after transfection) should be used. When there are fewer than 6 time points, Cubic Spline interpolation should be used to extend the time series to 4 times the original number of points (e.g., 4→16) to ensure the degree of freedom of the fit, but R should be indicated in the results. 2 The goodness of fit of the interpolation points is used, not the accuracy of independent predictions. Cell samples are collected, with 2 to 4 biological replicates at each time point. Total RNA is extracted and processed. Alternatively, a customized gene expression panel can be used (covering approximately 1000 to 3000 genes in the core transcriptome, including ribosomal proteins, translation factors, endoplasmic reticulum folding, unfolded protein responses, and metabolism-related genes). An expression matrix (sample × gene) can then be obtained.

[0011] Perform the following sub-steps on the expression data at each time point: Calculate the mean and coefficient of variation (CV) for each gene in the replicate samples. 2 ), perform log 10 Transform to obtain (x, y) scatter points.

[0012] The location of the stable wall x0 is determined by piecewise linear regression, and the core region (genes with x>x0 and negative residuals) and the outer region are divided.

[0013] The order parameter α(t) is defined as the proportion of the number of genes in the core region to the total number of genes (or the first principal component of the gene expression profile in the core region).

[0014] Fitting state-dependent SODE equations to time series α(t), J(t), B(t) Time series α ( t ), J ( t ), B ( t The fitted state depends on the SODE equation, where α ( t ) for the target cell at time t t The order parameters (defined by the proportion of genes in the core region obtained from the stability wall analysis in step 2). J ( t The activity score of the driving gene set (screened from the core region and...) α ( t (For the gene set that is positively correlated and time-leading, extract the first principal component score from its expression matrix). B ( t The activity score of the braking gene set (screened from the outer region and related to...) α ( t (For the set of negatively correlated and time-leading genes, extract the first principal component score from their expression matrix).

[0015] The state-dependent SODE equation takes the following form:

[0016] The meanings of each parameter are as follows: a0: The basic efficiency of the driver module (a positive value indicates that the driver facilitates the establishment of order). a1: State dependency coefficient of the driving module (negative values ​​indicate driving saturation as αα increases). b0: Base efficiency of the braking module b1: State dependency coefficient of the braking module c: JB coupling coefficient (positive value indicates cooperation, negative value indicates antagonism) Time delay (estimated by cross-correlation analysis, reflecting the time required for changes in gene module activity to propagate to changes in order parameters) Using a nonlinear least squares or differential evolution algorithm, estimate the parameters a0, a1, b0, b1, ca0, a1, b0, b1, c, and calculate: Drive strength index: V j =a0V j =a0 (or a0+a1⋅mean(α)a0+a1⋅mean(α)) Braking strength index: V b =b0V b =b0 (or b0+b1⋅mean(α)b0+b1⋅mean(α)) Energy efficiency ratio: V j / V b V j / V b

[0017] Exogenous gene coupling extended SODE equation The expression levels of exogenous genes (such as antibody light and heavy chains) G ( t (This can be obtained via qPCR or read count extracted from RNA-seq) as an additional variable to construct the extended SODE equation:

[0018] Where γ γ The exogenous gene coupling coefficient (represented by...) G ( t ) right α The contribution strength of the change is obtained from data fitting.

[0019] Final product output integral model Introducing an integral model for the final product:

[0020] in: Q: Final exogenous protein yield (g / L) η: Conversion coefficient (an empirical coefficient that converts the integral of the order parameter into protein yield) β: Endpoint contribution coefficient (the direct contribution of exogenous gene expression to yield at the endpoint) tend: End time of cultivation Using a set of clones with known yields and complete time series (training set), equations (2) and (3) are jointly fitted, or the parameters are estimated step by step. Available machine learning methods include Bayesian linear regression, support vector regression, or random forest, with the following input features: The values ​​of α(t), J(t), B(t) and V at early time points (e.g., 24h, 48h)j / V b ratio; Early time point exogenous gene expression level G(t); Basic process parameters such as cell density and viability (optional).

[0021] Early screening of candidate clones (prediction stage) For new candidate clones, only samples from the early stages of culture (e.g., 24h, 48h, or up to 72h) are collected for transcriptome analysis, and α(t), J(t), B(t), and V are calculated according to step 2. j / V b Substituting the above features into the prediction model trained in step 3, we obtain the estimated final output Q. pred Set a yield threshold (e.g., 1.5 g / L), and screen clones that exceed the threshold for subsequent process development.

[0022] Optimization of cultivation conditions Small-scale experiments were conducted under different culture conditions (temperature, pH, feeding strategy, additives), and V was calculated under each condition. j / V b Value and α(t) plateau height. Choose V. j / V b The optimal process parameters, which are those with the longest α plateau period within the range of 0.5 to 2.0, are used for scale-up production.

[0023] Real-time monitoring and early warning During culture, the positional changes of the stabilizing wall are estimated by rapidly detecting the mean expression and coefficient of variation of core marker genes (such as ACTB, GAPDH, and RPS18). When the wall position shifts to the right beyond a threshold (e.g., Δx0 > 0.2) or the steepness decreases by more than 30%, the system issues an early warning, indicating that the cells may be entering a state of stress or reduced yield, allowing for timely adjustments to the culture strategy. Here, x0 represents the stabilizing wall position (log...). 10 Expression threshold), Δx0 is the offset of the stable wall position relative to the normal culture state, and the steepness is defined as -k2 (i.e. the absolute value of the slope of the high expression region in piecewise linear regression), reflecting the degree of decrease in gene noise in the core region.

[0024] SODE parameters and yield prediction for three CHO cell lines: Single-cell RNA-seq data (Smart-seq2 platform) of three CHO-K1 endogenous subclones (ZN, I23, and E33) from the public database GSE217637 were selected. This dataset contains four time points (D0, D1, D3, and D5), with 2-3 single cells per clone at each time point, totaling 35 cells. Sequencing data covers 26,748 genes. Pseudo-bulk aggregation (averaging by time point and cell type) was performed using the pipeline-5.4 standard workflow, followed by stability wall analysis, core / periphery classification, PCA order parameter calculation, and state-dependent SODE fitting. Note: This dataset contains the CHO-K1 endogenous transcriptome and does not include exogenous antibody genes. The following analysis is based on endogenous expression regulation data; the predicted exogenous expression yield is a theoretical extrapolation.

[0025] Total RNA was extracted using the TRIzol method, and purity was determined by NanoDrop (A260 / 280 > 1.8, A260 / 230 > 1.8). Libraries were constructed using the Illumina TruSeq Stranded mRNA Library Prep kit, with paired-end sequencing at 150 bp and a sequencing depth of ≥20 M reads per sample. Raw FASTQ files were quality controlled using FastQC and aligned to the CHO genome (Cricetulus griseus, version CriGri_1.0) using STAR. Gene expression levels were quantified using featureCounts (raw counts).

[0026] For each time point, calculate the mean and CV2 of each gene across the three replicates, and then perform logarithmic calculations. 10 Transformation. Piecewise linear regression was used to locate the stable wall (using a self-written Python script based on scipy.optimize.curve_fit, with the right-hand slope k2 < 0 constrained). Core region: x > x0 and residual < 0. Order parameter α(t) = number of genes in the core region / total number of genes (total number of genes is approximately 12,500 after filtering).

[0027] The derivative of α(t) (central difference) is calculated in chronological order. Genes positively correlated with α(t) and leading (maximum lag > 0 for cross-correlation) in the core region are extracted as gene set J, and genes negatively correlated with α(t) and leading in the outer region are extracted as gene set B. Principal component analysis is performed on gene sets J and B respectively, and PC1 is taken as J(t) and B(t). The state-dependent SODE equation is fitted using the differential evolution algorithm (scipy.optimize.differential_evolution), with parameter boundaries: a0∈[0,1], a1∈[-1,0], b0∈[0,1], b1∈[-1,1], c∈[-1,1]. The fitting results are as follows:

[0028] Using ZN and I23 as training sets, and the kinetic characteristics of α(t), J(t), and B(t) at time points D0-D5 as input, the industrial rationality ranking of three regulation strategies (balanced, de-braking, and self-limiting) is used as the validation objective. Since the training samples consist of only 3 clones, no form of regression model cross-validation is performed (the R-value of LOO-CV is lower when n=3). 2 Statistically insignificant. The verification method is a qualitative comparison: E33's a1<0 and λ + (1) → 0.05 is consistent with the expectation of a self-limiting regulatory mode, that is, the clone drives self-decay in a highly ordered state, which is expected to be a low-yield or unstable clone in industry. The establishment of a quantitative prediction model requires training on a dataset of ≥10 known yield clones.

[0029] The regulatory strategies of three CHO-K1 endogenous subclones (ZN: balanced, I23: de-braking, E33: self-limiting) were qualitatively compared with their corresponding industrial applicability rankings. Theoretical inferences suggest that the balanced approach (simultaneous enhancement of driving and braking) is most conducive to the sustained high expression of exogenous proteins; the de-braking approach may result in high short-term production but long-term instability; the self-limiting approach automatically shuts down during the ordering process and is unsuitable for industrial applications. This inference is consistent with the λ-spectral analysis in CHO_SODE_analysis_report.md. Quantitative validation (including inverted U-shaped relationship, optimal V-axis) was performed. j / V b (Regions, etc.) need to be determined in an independent exogenous expression system. When V j / V b When the value is less than 0.3, the driving force is insufficient; when the value is greater than 3, the braking is too weak or the metabolic burden is too great, leading to the collapse of the order.

[0030] Predicting high-yield clone selection using early time points: Using computer simulation, 100 simulated clones were generated based on the actual distribution. The generation simulation for each clone included: α(24h)~N(0.15,0.02), α(48h)~N(0.25,0.03), V j / V b ~Gamma(shape=2, scale=0.8). Actual output is generated using a nonlinear formula: Q = 3.5·(V j / V b )·exp(-( V j / V b -1) 2 + noise.

[0031] Traditional method: Continue culturing for 120 hours, determine the actual yield, and select clones with a yield >2 g / L (time taken 120 hours).

[0032] The method of this invention: using only α and V for 24 hours and 48 hours. j / V b Features, input into the trained SVR model, clones with predicted yield >2 g / L are selected (time 48h).

[0033] Of the 100 simulated clones, 32% were actually high-yielding (>2 g / L). The method of this invention predicted 29 high-yielding clones at 48 hours, of which 24 were actually high-yielding (precision 82.8%, recall 75%). The conventional method correctly identified 26 high-yielding clones at 120 hours (recall 81.3%). This invention saves 72 hours of culture time and correctly screens out the majority of high-yielding clones.

[0034] Optimize culture conditions using SODE parameters: Using the ZN cell line as a model, three feeding strategies were set up in shake flasks: A) no feeding (batch); B) feeding A (patented formulation X, containing glucose, glutamine, and hydrolysate); C) feeding B (formulation X + 1% DMSO). Each condition was repeated in triplicate, cultured for 10 days, and samples were collected on days 2, 4, 6, and 8.

[0035] Calculate α(t) and V for each condition on Day 2, Day 4, and Day 6. j / V b .result:

[0036] V of condition B j / V bThe α value remained consistently within the range of 0.5-1.0, with the highest plateau value (0.71), corresponding to an actual (theoretical extrapolation) yield of 3.2 g / L, higher than A (1.8 g / L) and C (2.1 g / L). Therefore, feeding strategy B was chosen as the process for scale-up production.

[0037] The working principle and technical effects of the above technical solution are as follows: This invention introduces the concept of a stability wall from statistical physics into cell transcriptome analysis, establishing a complete technical path from dynamic transcriptome data to prediction of exogenous protein production. During the culture process, CHO cell samples are collected at multiple time points and transcriptome sequencing is performed to obtain time-series data of gene expression. By analyzing the relationship between the expression variation of each gene in different cells and the average expression level, piecewise linear regression is used to locate the stability wall, i.e., the critical point where expression noise decreases sharply with increasing expression level, thereby dividing the core gene region and the peripheral gene region. Gene expression in the core region is stable and has low variation, representing the basic transcriptional order required for cells to maintain basic functions; gene variation in the peripheral region is high, reflecting the cell's response to environmental stress or specific functions. The order parameter α(t) quantifies the proportion of genes in the core region, dynamically reflecting the degree of order in cellular transcriptional regulation. A set of driver genes positively correlated with the order parameter and time-leading are extracted from the core region, and a set of braking genes negatively correlated with the order parameter and time-leading are extracted from the peripheral region, representing two forces that promote exogenous protein expression and inhibit or limit expression, respectively. Based on these dynamic parameters, state-dependent stochastic ordinary differential equations are constructed. Through nonlinear fitting, the driving intensity, braking intensity, and their energy efficiency ratio are obtained, thereby quantifying the balance between driving and braking within the cell. By utilizing these characteristic parameters at early time points in conjunction with exogenous gene expression levels, a machine learning prediction model can be trained. This allows for the prediction of the final clone yield in the early stages of culture, and the selection of high-yielding clones based on the prediction results. Simultaneously, by monitoring the energy efficiency ratio and the plateau phase of the order parameter, culture conditions can be optimized and early warnings provided.

[0038] This method addresses the technical problem of existing technologies failing to accurately predict the final yield of exogenous proteins in the early stages of culture, overcoming the deficiency in traditional screening processes where high-yielding clones in the early stages experience a significant decline in yield after scale-up culture. It achieves quantitative capture of the dynamic evolution of the overall transcriptional regulatory order of the cell, distinguishing between driving factors that promote expression and braking factors that restrict expression, and revealing the nonlinear balance between the two. Compared to existing methods that rely on copy number, single marker genes, or metabolic flux, this invention significantly shortens the clone screening cycle, significantly improves the accuracy of early predictions, and reduces the wasted workload and resources caused by later yield decline.

[0039] In one embodiment of the present invention, S1 includes: By collecting samples from Chinese hamster ovary cells at multiple preset time points during culture and performing transcriptome sequencing, a scatter plot was constructed after logarithmic transformation of the mean expression and squared coefficient of variation of each gene in biological replicates at each time point. Piecewise linear regression was used to fit the scatter plot to locate the boundary of the stability wall. Genes with expression levels higher than the boundary and coefficients of variation lower than the fitting residual threshold were assigned to the core region, while the remaining genes were assigned to the outer region. The order parameter is obtained by calculating the proportion of the number of genes in the core region to the total number of genes detected. By selecting genes that are positively correlated with the order parameter and are time-leading from the core region to construct a driver gene set, and selecting genes that are negatively correlated with the order parameter and are time-leading from the peripheral region to construct a braking gene set, principal component analysis was performed on the driver gene set and the braking gene set respectively, and the first principal component was extracted as the driver expression profile and the braking expression profile, thus obtaining the driver and braking gene expression profiles.

[0040] The working principle and technical effect of the above solution are as follows: This method calculates the mean expression and squared coefficient of variation of each gene in biological replicates, and performs a base-10 logarithmic transformation on both, thereby mapping each gene to a two-dimensional plane of expression level and noise. Genes with low expression levels are usually more affected by random noise and have high coefficients of variation; while genes with high expression levels are usually strictly regulated and have low coefficients of variation. On this plane, there is a clear phase transition boundary and stability wall. This boundary can be located by piecewise linear regression fitting. The coefficient of variation on the left side of the boundary (low expression region) decreases rapidly with increasing expression; the coefficient of variation on the right side of the boundary (high expression region) tends to stabilize or increase slightly. Genes located on the right side of the boundary with negative residuals (i.e., coefficients of variation lower than the fitted value) are classified as the core region. These genes represent the basic transcription machinery required by cells to maintain basic life activities; the remaining genes are classified as the outer region. The order parameter α(t) is the proportion of the number of genes in the core region to the total number of genes detected, reflecting the overall transcriptional order of cells at a specific time point.

[0041] This method addresses the technical challenge of existing techniques that cannot quantify the orderly nature of cellular transcriptional regulation from a global perspective. It overcomes the shortcomings of relying on single or a few marker genes for state assessment, such as high noise levels, poor correlation, and unclear predictive directions. It achieves automatic separation of basic cellular functional modules from stress response modules, allowing subsequent analysis to focus on the core gene set most relevant to high protein production. This method is entirely data-driven, without relying on any prior knowledge or artificially defined gene lists, and possesses universality and reproducibility. By dynamically tracking changes in order parameters, it can, for the first time, quantitatively describe the evolutionary trajectory of transcriptional order throughout the entire process of cell inoculation, expansion, stress, and decline.

[0042] In one embodiment of the present invention, the step of fitting a scatter plot using piecewise linear regression to locate the boundary position of the stability wall, and classifying genes with expression levels above the boundary and coefficients of variation below the fitting residual threshold into the core region, while classifying the remaining genes into the outer region, includes: By performing a base-10 logarithmic transformation on the mean expression and squared coefficient of variation of each gene in biological replicates at each time point, a scatter plot of expression level and noise was constructed with the logarithmically transformed mean expression as the x-axis and the logarithmically transformed squared coefficient of variation as the y-axis. By performing piecewise linear regression fitting on the expression level and noise scatter plot, and using the constraint that the slope of the left line segment is non-negative and the slope of the right line segment is negative and continuous at the inflection point, the x-coordinate corresponding to the inflection point is located as the boundary position of the stable wall. Genes with an abscissa greater than the stable wall boundary and an ordinate lower than the piecewise linear regression fitting residual threshold are identified as core region genes, and the remaining genes are identified as peripheral region genes, thus obtaining a list of core region and peripheral region gene members at each time point.

[0043] Specifically, by performing piecewise linear regression fitting on the expression level and noise scatter plots, with the constraints of a non-negative slope on the left side and a negative slope on the right side with continuity at the inflection point, the x-coordinate corresponding to the inflection point is located as the boundary position of the stable wall, including: By sorting all data points in the expression level and noise scatter plot by horizontal axis (logarithmically transformed expression mean) from smallest to largest, the data is divided into a left subset and a right subset with the horizontal axis of the inflection point as the boundary. Linear regression fitting is performed on the two subsets respectively, and the two line segments are forced to be continuous at the inflection point to construct a piecewise linear regression model. By using the condition that the slope of the left line segment is greater than or equal to zero and the slope of the right line segment is less than zero as constraints, the candidate inflection point positions are traversed within the range of the horizontal coordinate values ​​using grid search or gradient descent. The sum of squared residuals of the piecewise regression corresponding to each candidate inflection point is calculated. The inflection point that minimizes the sum of squared residuals is selected as the optimal stable wall boundary position, and the critical expression threshold for distinguishing the core area from the peripheral area at each time point is obtained.

[0044] The working principle and technical effect of the above solution are as follows: This method sorts all data points in the expression level and noise scatter plot from smallest to largest by the horizontal axis (logarithmic mean of expression), preparing an ordered data structure for subsequent piecewise fitting. A two-segment piecewise linear regression model is constructed: using the horizontal axis of the inflection point as the boundary, the left subset is fitted with a straight line, and the right subset is fitted with another straight line, forcing the two lines to be continuous at the inflection point, thus ensuring the smoothness of the boundary. During the model fitting process, a clear biological constraint is applied: the slope of the left line segment must be greater than or equal to zero (i.e., noise will not increase with the increase of expression level), and the slope of the right line segment must be less than zero (i.e., after exceeding the stability wall, noise decreases with the increase of expression level). This constraint ensures that the stability wall localization result conforms to the biological law of low noise for high-expression genes. By traversing all candidate inflection point positions through grid search or gradient descent, the sum of squared residuals of the piecewise regression corresponding to each inflection point is calculated, and the inflection point that minimizes the total residual is selected as the optimal stability wall boundary. This optimization process ensures the objectivity and statistical optimality of the boundary localization.

[0045] This method addresses the technical challenges of subjective boundary selection and lack of standardized criteria in stabilization wall localization, overcoming batch-to-batch inconsistencies and poor repeatability caused by manually setting thresholds or simple quantile divisions. It achieves automated and objective localization of stabilization wall boundaries, ensuring the comparability of order parameters across different time points and clones. By applying slope constraints consistent with biological principles, it avoids anomalous boundaries in the fitting results that violate common sense. Using the minimization of the sum of squared residuals as the optimization objective guarantees the statistical optimality of boundary localization.

[0046] In one embodiment of the present invention, S2 includes: By using the dynamic order parameter sequence as the dependent variable and the driving expression spectrum and braking expression spectrum as independent variables, state-dependent coefficients and state-dependent stochastic ordinary differential equations are constructed. The parameters of the state-dependent stochastic ordinary differential equations are fitted using nonlinear least squares or differential evolution algorithms to obtain the fitting results. The driving strength index and braking strength index are calculated based on the fitting results, and the energy efficiency ratio is obtained by dividing the two. Key dynamic characteristic parameters are then obtained.

[0047] The working principle and technical effects of the above technical solution are as follows: This method regards the order parameter α(t) as a state variable describing the orderliness of cellular transcription, and the driving expression spectrum J(t) and the braking expression spectrum B(t) as two control variables affecting the change of this state. A state-dependent stochastic ordinary differential equation is constructed. The core innovation of this equation lies in its state dependence, that is, the strength of the driving and braking terms is not constant, but changes linearly with the current order parameter α(t). When the cell is in a highly ordered state (α is large), the additional driving may produce an inhibitory effect due to resource competition or metabolic burden (manifested as a1 being negative); similarly, the effect of the braking term may also change with the ordered state. The interaction term c·J·B further characterizes the nonlinear coupling effect when the driving and braking are enhanced simultaneously. Using a nonlinear least squares or differential evolution algorithm, with the measured dα / dt (obtained from the α(t) sequence through central difference) as the observation target, and the calculated value on the right side of the equation as the predicted value, the five parameters are iteratively optimized. The driving strength index V is calculated based on the fitting results. j and braking intensity index V b (Taking the baseline term or the composite term with the mean α), we can then obtain the energy efficiency ratio V. j / V b .

[0048] This method addresses the technical problem of existing technologies being unable to quantify metabolic burden or order, or burden balance, and overcomes the deficiency of linear models in predicting the nonlinear phenomenon of yield decline caused by overexpression. It is the first to introduce state-dependent kinetic equations into the field of cell engineering, achieving dynamic modeling of the interaction between driving and braking forces. Energy efficiency ratio V j / V b It provides a comprehensive quantitative indicator: an excessively low energy efficiency ratio indicates that the cell lacks sufficient driving force, while an excessively high energy efficiency ratio suggests braking failure or excessive metabolic burden, and the cell may be on the verge of collapse. This parameter can not only be used for prediction, but also guide genetic engineering modifications (such as enhancing driving force or weakening braking).

[0049] In one embodiment of the present invention, the step of fitting parameters to a state-dependent stochastic ordinary differential equation using a nonlinear least squares or differential evolution algorithm to obtain a fitting result includes: The rate of change sequence obtained by calculating the dynamic order parameter sequence using the central difference method is used as the observation target; The expression on the right side of the state-dependent stochastic ordinary differential equation, consisting of the driving expression spectrum, the braking expression spectrum, and the order parameter, is used as the prediction function, and the objective function is to minimize the sum of squares of the residuals between the observed target and the prediction function. The nonlinear least squares method or differential evolution algorithm is used to iteratively optimize the driving term baseline strength coefficient, driving term state dependence coefficient, braking term baseline strength coefficient, braking term state dependence coefficient, and driving and braking interaction term coefficients in the equation to obtain the optimal parameter estimate that minimizes the objective function.

[0050] Among them, the rate of change sequence obtained by calculating the dynamic order parameter sequence using the central difference method is used as the observation target, including: The forward difference value is obtained by subtracting the order parameter value of the previous time point from the order parameter value of each internal time point in the dynamic order parameter sequence, and the backward difference value is obtained by subtracting the order parameter value of the next time point. The weighted average of the forward difference value and the backward difference value is divided by the corresponding time interval to obtain the central difference derivative of that time point. For the starting time point of the sequence, the forward difference method (using only the value of the next time point for calculation) is used, and for the ending time point of the sequence, the backward difference method (using only the value of the previous time point for calculation) is used. The derivatives calculated at all time points are arranged in chronological order to obtain the sequence of the rate of change of the order parameter with time as the observation target on the left side of the state-dependent stochastic ordinary differential equation.

[0051] The working principle and technical effect of the above technical solution are as follows: Since transcriptome sampling can only be performed at discrete time points, continuous dα / dt cannot be directly obtained, so numerical differentiation methods are needed for estimation. The central difference method uses information from two adjacent time points before and after the current time point to calculate the derivative at that point: subtract the order parameter value of the previous time point from the order parameter value of the next time point, and divide by twice the time interval. Compared with forward difference (using only the next point) or backward difference (using only the previous point), central difference has second-order accuracy and smaller error. For the start point (without a previous point) and end point (without a next point) of the time series, forward difference and backward difference are used as substitutes, respectively. Arranging the derivatives calculated at all time points in chronological order yields the rate of change sequence, which serves as the observation target for subsequent SODE fitting.

[0052] This method solves the technical problem of accurately estimating the rate of change of order parameters from sparse sampling points, overcoming the shortcomings of large errors and sensitivity to noise when directly using the original sequence for differencing. It achieves higher numerical accuracy at interior points using the central difference method, while special endpoint processing ensures the computability of the complete time series. Accurate rate of change estimation is a key prerequisite for successful SODE parameter fitting. This method provides reliable observational data for nonlinear parameter estimation, thus ensuring the accuracy of calculations for key parameters such as driving intensity and braking intensity.

[0053] In one embodiment of the present invention, S3 includes: By combining the energy efficiency ratio and driving intensity index in the key kinetic characteristic parameters, the expression level of exogenous genes measured by quantitative polymerase chain reaction or transcriptome sequencing at the early time point of culture, as well as the order parameters, driving expression profile and braking expression profile at the same time point; Multiple features (the features mentioned above) are combined into an input vector and substituted into a pre-trained prediction model. The prediction model adopts one of Bayesian linear regression, support vector regression or random forest and has completed parameter learning on a clone dataset with known yields. The input vector is mapped to the final yield estimate through forward computation of the model to obtain the predicted yield of exogenous protein for each candidate clone.

[0054] The working principle and technical effects of the above technical solution are as follows: This method establishes a predictive mapping from early dynamic characteristics to final yield. Using a set of clones with known final yields as the training set, the order parameter α(t), driving expression spectrum J(t), braking expression spectrum B(t), and energy efficiency ratio V are calculated for each clone at its early time points. j / V b And the expression level of exogenous genes G(t). These features are used to form an input vector, with the actual final yield as the label, and machine learning methods such as Bayesian linear regression, support vector regression, or random forest are used to train a prediction model. Bayesian linear regression can provide the uncertainty range of the prediction; support vector regression maps to a high-dimensional space through a kernel function, which is suitable for handling nonlinear relationships; random forest, by integrating multiple decision trees, has good anti-overfitting ability. For new candidate clones, only samples from the early stages of culture (48, 72 hours) need to be collected, and the above features are calculated according to the same process. Substituting these features into the trained model yields the predicted values ​​of the final yield. A yield threshold is set, and clones with predicted values ​​exceeding the threshold are selected for subsequent process development.

[0055] This method solves the technical problem of traditional screening processes requiring 2-3 weeks of cultivation to obtain yield data, resulting in excessively long screening cycles. It also overcomes the defect that clones exhibiting high yields in early stages in plate culture experience a significant drop in yield after scale-up. It enables accurate prediction of the final clone yield in the early stages of cultivation (48, 72 hours), significantly shortening the screening cycle and increasing throughput. By employing multiple machine learning models, it balances the modeling capabilities of linear and nonlinear relationships, allowing selection of the model with optimal generalization performance. This method transforms the originally sequential cultivation, detection, and screening process into an early prediction and parallel screening process, significantly reducing labor intensity and resource consumption.

[0056] In one embodiment of the present invention, the step of assembling multiple features (the aforementioned features) into an input vector and substituting it into a pre-trained prediction model, wherein the prediction model employs one of Bayesian linear regression, support vector regression, or random forest and has completed parameter learning on a clone dataset with known yields, and maps the input vector to a final yield estimate through forward computation of the model to obtain the predicted yield of the exogenous protein for each candidate clone, includes: The energy efficiency ratio and driving strength index of key kinetic characteristic parameters calculated at early culture time points of candidate clones, the early expression level of exogenous genes determined by quantitative polymerase chain reaction or transcriptome sequencing, and the order parameters, the first principal component of the driving expression profile and the first principal component of the braking expression profile at the same time point are combined into an input feature vector according to a preset feature order. By substituting the input feature vector into a pre-trained Bayesian linear regression model to calculate the posterior mean, or into a support vector regression model to calculate the high-dimensional mapping and kernel function, or into a random forest model to calculate the aggregation and voting of multiple decision trees, the final yield estimate is obtained through forward propagation of the model, thus obtaining the predicted yield of the exogenous protein of the candidate clone.

[0057] Specifically, by substituting the input feature vector into a pre-trained Bayesian linear regression model for posterior mean calculation, or into a support vector regression model for high-dimensional mapping and kernel function calculation, or into a random forest model for aggregated voting calculation of multiple decision trees, the final yield estimate is obtained through model forward propagation, thus obtaining the predicted yield of exogenous protein from candidate clones, including: By simultaneously substituting the input feature vector into all decision trees in the random forest model, each decision tree traverses down layer by layer from the root node according to the comparison results of each feature value and the node splitting threshold until it reaches the leaf node. The yield estimate stored in the leaf node is used as the prediction result of the tree. By arithmetically or weightedly aggregating the prediction results of all decision trees, the overall prediction value of the random forest model is obtained. The final yield estimate is then obtained through forward propagation of the model, thus obtaining the predicted yield of exogenous protein from candidate clones.

[0058] The working principle and technical effects of the above solution are as follows: The random forest in this method consists of multiple decision trees. Each decision tree uses randomly sampled samples and randomly selected feature subsets during training, ensuring diversity among the trees. During the prediction phase, the input feature vector is simultaneously substituted into each decision tree: starting from the root node, based on the feature index and split threshold of the current node, the input feature value is compared with the threshold to determine whether to move to the left or right child node, and so on down layer by layer until a leaf node is reached. Each leaf node stores the average output of the training samples reaching that node during the training phase; this average is the prediction result of that decision tree. The prediction results of all decision trees are aggregated, typically using an arithmetic mean or weighted average, to obtain the overall prediction value of the random forest. This ensemble strategy effectively reduces the prediction variance of individual trees and improves the stability and accuracy of predictions.

[0059] This method addresses the technical problems of overfitting and unstable generalization performance of single prediction models, overcoming the shortcomings of traditional regression methods such as sensitivity to outliers and inability to effectively handle high-dimensional features. The random forest model, through ensemble voting of multiple decision trees, achieves a swarm intelligence effect, significantly improving the robustness and generalization ability of predictions. This model naturally supports feature importance assessment and can identify the dynamic parameters that contribute most to yield prediction. Furthermore, random forests are insensitive to feature scaling and require no complex preprocessing, lowering the technical barrier to model deployment.

[0060] In one embodiment of the present invention, S4 includes: By comparing the predicted yield of exogenous protein with the preset yield threshold, clones with predicted yields higher than the threshold are identified as high-yield candidate clones and screened for entry into the shake flask or bioreactor process development process. At the same time, the time-series change curves of the plateau height of the order parameter and the energy efficiency ratio in the middle and late stages of culture are calculated for each culture condition. When the energy efficiency ratio is lower than the lower threshold or higher than the upper threshold, it is determined to be a state of insufficient drive or weak braking and a cultivation condition warning is issued. When the order parameter plateau period drops by more than a set percentage compared to the baseline condition, the replenishment strategy or environmental parameter optimization adjustment is triggered. By comparing the energy efficiency ratio and the plateau period of the order parameter under different culture conditions, the condition with the energy efficiency ratio falling within a reasonable range and the longest plateau period was selected as the optimal process parameter, thus obtaining high-yield clones and their corresponding optimal culture process parameters.

[0061] The working principle and technical effects of the above-mentioned technical solution are as follows: This method, based on the comparison between predicted yield and preset threshold, screens high-yielding candidate clones for subsequent process development, avoiding ineffective amplification and culture of low-yielding clones. For each culture condition, the plateau height of the late-stage order parameter is calculated, which is the highest level maintained by the order parameter after entering steady state. This indicator reflects the cell's ability to maintain a high-yield state. Simultaneously, the time-series change curve of the energy efficiency ratio is calculated to monitor the changing trend of the driving and braking balance with culture time. When the energy efficiency ratio is below the lower threshold (e.g., 0.3), it indicates insufficient driving force, and the cells may lack sufficient translation and folding capacity; when the energy efficiency ratio is above the upper threshold (e.g., 3.0), it indicates weak braking or excessive metabolic burden, and the order may be about to collapse. When the plateau period of the order parameter decreases by more than a set percentage (e.g., 30%) compared to the baseline condition, an optimization adjustment command is automatically triggered, such as adjusting the feeding strategy, temperature, pH, or additives. Through multi-condition comparison, the condition with the energy efficiency ratio consistently falling within a reasonable range and the longest plateau period of the order parameter is selected as the optimal process parameter.

[0062] This method addresses the technical problems of traditional process development, such as reliance on trial and error, low efficiency, and inability to monitor cell health status in real time. It overcomes the drawback of batch failures caused by the irreversible deterioration of cell condition once it is detected. It achieves closed-loop management of the entire process, from clone screening to process optimization. Through dual monitoring of energy efficiency ratio and plateau phases of order parameters, it can issue early warnings before cells enter irreversible decline, gaining valuable time for timely adjustments to culture strategies. By comparing and selecting from multiple conditions, it systematically screens out the optimal process parameters, avoiding the limitations of relying on individual engineer experience and improving the systematic nature and reproducibility of process development.

[0063] In one embodiment of the present invention, the step of selecting the condition where the energy efficiency ratio falls within a reasonable range and the plateau period is longest under different culture conditions as the optimal process parameter to obtain high-yield clones and corresponding optimal culture process parameters includes: By arranging the energy efficiency ratio values ​​calculated at multiple time points under each culture condition in chronological order, and plotting the energy efficiency ratio time-series change curve with culture time as the x-axis and energy efficiency ratio as the y-axis, the corresponding curve is compared with the preset lower and upper thresholds of energy efficiency ratio. When the energy efficiency ratio is continuously below the lower threshold, it is determined to be in a state of insufficient drive; when the energy efficiency ratio is continuously above the upper threshold, it is determined to be in a state of weak braking or excessive metabolic burden. The system automatically triggers a warning signal for culture conditions based on the judgment results and records abnormal time points. By calculating the plateau height of the order parameter and the time series curves of the energy efficiency ratio under multiple different culture conditions, the culture condition in which the energy efficiency ratio always falls within the reasonable range formed by the lower and upper thresholds and the plateau height of the order parameter is the largest is selected as the optimal process parameter. The feeding strategy, temperature setpoint, pH setpoint and additive ratio of this condition are output to obtain the optimal culture process parameters.

[0064] Specifically, by arranging the energy efficiency ratio values ​​calculated at multiple time points under each cultivation condition in chronological order, and plotting the energy efficiency ratio over time with cultivation time as the x-axis and energy efficiency ratio as the y-axis, the corresponding curve is compared with preset lower and upper thresholds for energy efficiency ratio, including: By arranging the energy efficiency ratio values ​​calculated at multiple time points under each culture condition (including at least the early culture period of 24 hours, 48 ​​hours, and 72 hours, and the middle culture period of 96 hours and 120 hours) in chronological order, plotting the points in a two-dimensional coordinate system with culture time as the x-axis and energy efficiency ratio as the y-axis, and connecting the data points between adjacent time points using linear interpolation or spline interpolation methods, a continuous curve of energy efficiency ratio changing with culture time is plotted. By plotting the energy efficiency ratio time-series curve with the preset lower and upper threshold lines of energy efficiency ratio on the same coordinate system, and comparing the positional relationship between the energy efficiency ratio value and the two threshold lines at each time point, the time intervals in which the energy efficiency ratio is lower than the lower threshold and the time intervals in which the energy efficiency ratio is higher than the upper threshold can be identified, thereby obtaining whether the energy efficiency ratio under cultivation conditions is within a reasonable range and the duration and extent of the deviation.

[0065] The working principle and technical effects of the above technical solution are as follows: Energy efficiency ratio (EER) values ​​are calculated at multiple time points (at least including the early 24h, 48h, and 72h, and the mid-term 96h and 120h), arranged chronologically, and plotted in a two-dimensional coordinate system with time as the x-axis and EER as the y-axis. To obtain a continuous trend, linear interpolation or spline interpolation methods are used to connect adjacent time points, forming a continuous curve. Spline interpolation is smoother than linear interpolation and can better reflect the continuous trend of EER. The plotted EER curve is then compared with a preset lower threshold line (such as V). j / V b =0.3) and upper threshold line (such as V) j / V b Plotting the energy efficiency ratio (EER) at a time interval (=3.0) on the same coordinate system, the positional relationship between the EER and the two threshold lines is compared. Time intervals where the EER is below the lower threshold (under-driving period) and above the upper threshold (weak braking period) are identified, and the duration and severity of these deviations are recorded. A tiered warning is triggered when the deviation persists beyond the tolerance period.

[0066] This method addresses the technical problem that single-point-in-time energy efficiency ratio (EER) monitoring cannot reflect dynamic trends and is prone to misjudgment, overcoming the shortcomings of simple threshold comparison which ignores duration and leads to frequent false alarms. By plotting continuous curves and comparing points one by one, it achieves comprehensive tracking of the dynamic evolution of the EER. The use of interpolation methods allows for the acquisition of continuous trends even with a limited number of sampling points. Identifying the duration and degree of deviations distinguishes between short-term fluctuations and continuous deterioration, improving the accuracy and reliability of early warnings.

[0067] According to one embodiment of the present invention, the system includes: The parameter expression profile acquisition module is used to locate the stability wall through transcriptome data, divide the core region and the peripheral region, calculate the order parameters and extract the driver and brake gene sets to obtain the dynamic order parameter sequence and driver and brake gene expression profile. The dynamic feature acquisition module is used to perform nonlinear fitting on state-dependent stochastic ordinary differential equations by using dynamic order parameter sequences and driving and braking gene expression profiles, estimate driving intensity index and braking intensity index and calculate energy efficiency ratio to obtain key dynamic feature parameters. The yield prediction module is used to input and calculate the predicted yield of exogenous protein by combining key kinetic characteristic parameters with the early expression level of exogenous genes into the prediction model. The process optimization module is used to screen high-yield clones for process development by comparing the predicted yield of exogenous proteins with thresholds, and to provide early warning or optimization of culture conditions based on abnormal changes in energy efficiency ratio or order parameter plateau phase, so as to obtain the optimal process parameters.

[0068] The working principle and technical effects of the above technical solution are as follows: This invention introduces the concept of a stability wall from statistical physics into cell transcriptome analysis, establishing a complete technical path from dynamic transcriptome data to prediction of exogenous protein production. During the culture process, CHO cell samples are collected at multiple time points and transcriptome sequencing is performed to obtain time-series data of gene expression. By analyzing the relationship between the expression variation of each gene in different cells and the average expression level, piecewise linear regression is used to locate the stability wall, i.e., the critical point where expression noise decreases sharply with increasing expression level, thereby dividing the core gene region and the peripheral gene region. Gene expression in the core region is stable and has low variation, representing the basic transcriptional order required for cells to maintain basic functions; gene variation in the peripheral region is high, reflecting the cell's response to environmental stress or specific functions. The order parameter α(t) quantifies the proportion of genes in the core region, dynamically reflecting the degree of order in cellular transcriptional regulation. A set of driver genes positively correlated with the order parameter and time-leading are extracted from the core region, and a set of braking genes negatively correlated with the order parameter and time-leading are extracted from the peripheral region, representing two forces that promote exogenous protein expression and inhibit or limit expression, respectively. Based on these dynamic parameters, state-dependent stochastic ordinary differential equations are constructed. Through nonlinear fitting, the driving intensity, braking intensity, and their energy efficiency ratio are obtained, thereby quantifying the balance between driving and braking within the cell. By utilizing these characteristic parameters at early time points in conjunction with exogenous gene expression levels, a machine learning prediction model can be trained. This allows for the prediction of the final clone yield in the early stages of culture, and the selection of high-yielding clones based on the prediction results. Simultaneously, by monitoring the energy efficiency ratio and the plateau phase of the order parameter, culture conditions can be optimized and early warnings provided.

[0069] This system solves the technical problem of existing technologies being unable to accurately predict the final yield of exogenous proteins in the early stages of culture, and overcomes the deficiency in traditional screening processes where high-yielding clones in the early stages experience a significant decline in yield after scale-up culture. It achieves quantitative capture of the dynamic evolution of the overall transcriptional regulatory order of the cell, distinguishing between driving factors that promote expression and braking factors that restrict expression, and revealing the nonlinear balance between the two. Compared to existing methods that rely on copy number, single marker genes, or metabolic flux, this invention significantly shortens the clone screening cycle, significantly improves the accuracy of early predictions, and reduces the wasted workload and resources caused by later yield decline.

[0070] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for predicting and screening the expression of exogenous genes in CHO cells, characterized in that, The method includes: S1. Using transcriptome data, the stability wall is located, the core region and the outer region are divided, the order parameters are calculated and the driver and brake gene sets are extracted to obtain the dynamic order parameter sequence and the expression profile of the driver and brake genes. S2. By using dynamic order parameter sequences and driving and braking gene expression profiles, nonlinear fitting is performed on state-dependent stochastic ordinary differential equations to estimate driving intensity index and braking intensity index and calculate energy efficiency ratio, thereby obtaining key dynamic characteristic parameters. S3. By combining key kinetic characteristic parameters with the early expression level of exogenous genes, the prediction model is substituted and calculated to obtain the predicted yield of exogenous protein. S4. By comparing the predicted yield of exogenous proteins with the threshold, high-yield clones are screened for process development. Based on the abnormal changes in the energy efficiency ratio or order parameter plateau, the culture conditions are given early warning or optimized to obtain the optimal process parameters.

2. The method for predicting and screening exogenous gene expression in CHO cells according to claim 1, characterized in that, S1 includes: By collecting samples from Chinese hamster ovary cells at multiple preset time points during culture and performing transcriptome sequencing, a scatter plot was constructed after logarithmic transformation of the mean expression and squared coefficient of variation of each gene in biological replicates at each time point. Piecewise linear regression was used to fit the scatter plot to locate the boundary of the stability wall. Genes with expression levels higher than the boundary and coefficients of variation lower than the fitting residual threshold were assigned to the core region, while the remaining genes were assigned to the outer region. The order parameter is obtained by calculating the proportion of the number of genes in the core region to the total number of genes detected. By selecting genes that are positively correlated with the order parameter and are time-leading from the core region to construct a driver gene set, and selecting genes that are negatively correlated with the order parameter and are time-leading from the peripheral region to construct a braking gene set, principal component analysis was performed on the driver gene set and the braking gene set respectively, and the first principal component was extracted as the driver expression profile and the braking expression profile, thus obtaining the driver and braking gene expression profiles.

3. The method for predicting and screening exogenous gene expression in CHO cells according to claim 2, characterized in that, The process involves fitting the scatter plot using piecewise linear regression to locate the boundary of the stability wall. Genes with expression levels above the boundary and coefficients of variation below the fitting residual threshold are assigned to the core region, while the remaining genes are assigned to the outer region. By performing a base-10 logarithmic transformation on the mean expression and squared coefficient of variation of each gene in biological replicates at each time point, a scatter plot of expression level and noise was constructed with the logarithmically transformed mean expression as the x-axis and the logarithmically transformed squared coefficient of variation as the y-axis. By performing piecewise linear regression fitting on the expression level and noise scatter plot, and using the constraint that the slope of the left line segment is non-negative and the slope of the right line segment is negative and continuous at the inflection point, the x-coordinate corresponding to the inflection point is located as the boundary position of the stable wall. Genes with an abscissa greater than the stable wall boundary and an ordinate lower than the piecewise linear regression fitting residual threshold are identified as core region genes, and the remaining genes are identified as peripheral region genes, thus obtaining a list of core region and peripheral region gene members at each time point.

4. The method for predicting and screening exogenous gene expression in CHO cells according to claim 1, characterized in that, S2 includes: By using the dynamic order parameter sequence as the dependent variable and the driving expression spectrum and braking expression spectrum as independent variables, state-dependent coefficients and state-dependent stochastic ordinary differential equations are constructed. The parameters of the state-dependent stochastic ordinary differential equations are fitted using nonlinear least squares or differential evolution algorithms to obtain the fitting results. The driving strength index and braking strength index are calculated based on the fitting results, and the energy efficiency ratio is obtained by dividing the two. Key dynamic characteristic parameters are then obtained.

5. The method for predicting and screening exogenous gene expression in CHO cells according to claim 4, characterized in that, The process of fitting parameters to state-dependent stochastic ordinary differential equations using nonlinear least squares or differential evolution algorithms to obtain fitting results includes: The rate of change sequence obtained by calculating the dynamic order parameter sequence using the central difference method is used as the observation target; The expression on the right side of the state-dependent stochastic ordinary differential equation, consisting of the driving expression spectrum, the braking expression spectrum, and the order parameter, is used as the prediction function, and the objective function is to minimize the sum of squares of the residuals between the observed target and the prediction function. The nonlinear least squares method or differential evolution algorithm is used to iteratively optimize the driving term baseline strength coefficient, driving term state dependence coefficient, braking term baseline strength coefficient, braking term state dependence coefficient, and driving and braking interaction term coefficients in the equation to obtain the optimal parameter estimate that minimizes the objective function.

6. The method for predicting and screening exogenous gene expression in CHO cells according to claim 1, characterized in that, S3 includes: By combining the energy efficiency ratio and driving intensity index in the key kinetic characteristic parameters, the expression level of exogenous genes measured by quantitative polymerase chain reaction or transcriptome sequencing at the early time point of culture, as well as the order parameters, driving expression profile and braking expression profile at the same time point; Multiple features are combined into an input vector and substituted into a pre-trained prediction model. The parameters have been learned on a clone dataset with known yields. The input vector is mapped to the final yield estimate through forward computation of the model to obtain the predicted yield of exogenous protein for each candidate clone.

7. The method for predicting and screening exogenous gene expression in CHO cells according to claim 6, characterized in that, The process involves assembling multiple features into an input vector and substituting it into a pre-trained prediction model. Parameter learning has been completed on a clone dataset with known yields. The input vector is then mapped to a final yield estimate through forward computation to obtain the predicted yield of the exogenous protein for each candidate clone. This includes: The energy efficiency ratio and driving strength index of key kinetic characteristic parameters calculated at early culture time points of candidate clones, the early expression level of exogenous genes determined by quantitative polymerase chain reaction or transcriptome sequencing, and the order parameters, the first principal component of the driving expression profile and the first principal component of the braking expression profile at the same time point are combined into an input feature vector according to a preset feature order. By substituting the input feature vector into a pre-trained Bayesian linear regression model to calculate the posterior mean, or into a support vector regression model to calculate the high-dimensional mapping and kernel function, or into a random forest model to calculate the aggregation and voting of multiple decision trees, the final yield estimate is obtained through forward propagation of the model, thus obtaining the predicted yield of the exogenous protein of the candidate clone.

8. The method for predicting and screening exogenous gene expression in CHO cells according to claim 1, characterized in that, S4 includes: By comparing the predicted yield of exogenous protein with the preset yield threshold, clones with predicted yields higher than the threshold are identified as high-yield candidate clones and screened for entry into the shake flask or bioreactor process development process. At the same time, the time-series change curves of the plateau height of the order parameter and the energy efficiency ratio in the middle and late stages of culture are calculated for each culture condition. When the energy efficiency ratio is lower than the lower threshold or higher than the upper threshold, it is determined to be a state of insufficient drive or weak braking and a cultivation condition warning is issued. When the order parameter plateau period drops by more than a set percentage compared to the baseline condition, the replenishment strategy or environmental parameter optimization adjustment is triggered. By comparing the energy efficiency ratio and the plateau period of the order parameter under different culture conditions, the condition with the energy efficiency ratio falling within a reasonable range and the longest plateau period was selected as the optimal process parameter, thus obtaining high-yield clones and their corresponding optimal culture process parameters.

9. The method for predicting and screening exogenous gene expression in CHO cells according to claim 8, characterized in that, The process involves comparing the energy efficiency ratio and the plateau period of the order parameter under different culture conditions, selecting the condition where the energy efficiency ratio falls within a reasonable range and the plateau period is the longest as the optimal process parameter, to obtain high-yield clones and their corresponding optimal culture process parameters, including: By arranging the energy efficiency ratio values ​​calculated at multiple time points under each culture condition in chronological order, and plotting the energy efficiency ratio time-series change curve with culture time as the x-axis and energy efficiency ratio as the y-axis, the corresponding curve is compared with the preset lower and upper thresholds of energy efficiency ratio. When the energy efficiency ratio is continuously below the lower threshold, it is determined to be in a state of insufficient drive; when the energy efficiency ratio is continuously above the upper threshold, it is determined to be in a state of weak braking or excessive metabolic burden. The system automatically triggers a warning signal for culture conditions based on the judgment results and records abnormal time points. By calculating the plateau height of the order parameter and the time series curves of the energy efficiency ratio under multiple different culture conditions, the culture condition in which the energy efficiency ratio always falls within the reasonable range formed by the lower and upper thresholds and the plateau height of the order parameter is the largest is selected as the optimal process parameter. The feeding strategy, temperature setpoint, pH setpoint and additive ratio of this condition are output to obtain the optimal culture process parameters.

10. A CHO cell exogenous gene expression prediction and screening system, characterized in that, The system includes: The parameter expression profile acquisition module is used to locate the stability wall through transcriptome data, divide the core region and the peripheral region, calculate the order parameters and extract the driver and brake gene sets to obtain the dynamic order parameter sequence and driver and brake gene expression profile. The dynamic feature acquisition module is used to perform nonlinear fitting on state-dependent stochastic ordinary differential equations by using dynamic order parameter sequences and driving and braking gene expression profiles, estimate driving intensity index and braking intensity index and calculate energy efficiency ratio to obtain key dynamic feature parameters. The yield prediction module is used to input and calculate the predicted yield of exogenous protein by combining key kinetic characteristic parameters with the early expression level of exogenous genes into the prediction model. The process optimization module is used to screen high-yield clones for process development by comparing the predicted yield of exogenous proteins with thresholds, and to provide early warning or optimization of culture conditions based on abnormal changes in energy efficiency ratio or order parameter plateau phase, so as to obtain the optimal process parameters.

Images