CIK cell preparation method and system

Abstract

The invention relates to the technical field of biological information processing, in particular to a CIK cell preparation method and system.The CIK cell preparation method comprises the following steps that cell density data of a time sequence is collected, an autocorrelation function value and a partial autocorrelation function value are calculated, a first-order difference sequence and a second-order difference sequence are extracted, ADF statistics and KPSS statistics are calculated according to the difference sequences, and the CIK cell density data is obtained. And comparing the statistical magnitude range to judge the stationarity of the time sequence, and generating stationary cell density sequence data. According to the method, through extraction and optimization of the lag correlation value and the partial correlation value, the accuracy of model parameters is enhanced, reliable support is provided for dynamic correction and conjoint analysis of the future cell proliferation rate, abnormal points are accurately recognized based on analysis of the difference mean value and the weighted accumulation value, the anomaly detection efficiency is remarkably improved, and the method is suitable for large-scale popularization and application. Dynamic adjustment of culture conditions is realized by combining abnormal point indexes and drifting rate calculation of environmental parameters, the influence of environmental fluctuation on cell proliferation quality is reduced, and preparation quality and efficiency are guaranteed.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics processing, and particularly to a method and system for preparing CIK cells. Background Art

[0002] The technical field of bioinformatics processing includes technologies related to the acquisition, storage, analysis, and management of biological data based on computer technology. The core content of this technical field includes the processing of biological information through a computer system to achieve the analysis and interpretation of biological sample data to support multiple application directions such as medicine, pharmaceuticals, genomics, etc. Specifically, bioinformatics processing technology covers multiple aspects such as the interpretation of genomic data, protein sequence analysis, modeling and simulation of cell biology data, and data management and calculation related to medical diagnosis. It is an important field of the cross-integration of modern biology and information technology.

[0003] Among them, the CIK cell preparation method refers to a method for preparing cytokine-induced killer cells through specific technical means. This patent theme focuses on technical matters such as maintaining cell activity and optimizing induction factors during the in vitro amplification and preparation process of CIK cells. It adopts an operation process and processing method based on cell biology experimental techniques, specifically including the separation and culture of peripheral blood mononuclear cells, multi-step induction by optimizing cytokine combinations, regulating key parameters during cell proliferation and differentiation, and separating and purifying the amplified cells. Finally, a CIK cell preparation method system that meets application requirements is formed.

[0004] The prior art lacks effective processing of the non-stationarity of time series data, resulting in the interference of data fluctuations on model construction and unstable analysis results. The selection of prediction models lacks systematic optimization of lag correlation and is difficult to accurately capture complex dynamics. The anomaly detection method relies on static criteria and has insufficient sensitivity, unable to detect abnormal situations in real time. The adjustment of culture parameters is based on empirical rules and lacks real-time data support, resulting in adverse effects of environmental fluctuations on cell quality. The above deficiencies limit the accuracy and consistency of the cell amplification process and may affect the final preparation effect and its reliability in subsequent applications. Summary of the Invention

[0005] The purpose of the present invention is to solve the drawbacks existing in the prior art, and to propose a method and system for preparing CIK cells.

[0006] To achieve the above purpose, the present invention adopts the following technical scheme: A method for preparing CIK cells, comprising the following steps:

[0007] S1: Collect the cell density data of the time series, calculate the values of the autocorrelation function and the partial autocorrelation function, extract the first-order and second-order difference sequences, calculate the ADF statistic and the KPSS statistic based on the difference sequences, compare the ranges of the statistics to judge the stationarity of the time series, and generate the smoothed cell density sequence data;

[0008] S2: Based on the smoothed cell density sequence data, extract the lag autocorrelation value and the lag partial autocorrelation value, calculate the Bayesian information criterion value and the Akaike information criterion value for different lag orders, select the lag order corresponding to the minimum of the two criterion values, calculate the Ljung-Box statistic and the confidence interval of the error residual sequence, and generate the optimal ARIMA model parameter set;

[0009] S3: Based on the optimal ARIMA model parameter set, recursively calculate the predicted values of the cell proliferation rate at future time points, dynamically correct the upper and lower bounds of the confidence interval according to the prediction sequence, and conduct a joint analysis of the predicted cell density and proliferation rate to generate the dynamic predicted values of cell proliferation;

[0010] S4: Based on the dynamic predicted values of cell proliferation, extract the mean difference between the predicted value and the actual value at the current time point, calculate the moving weighted value and the cumulative sum value for the mean difference, combine the upper and lower bounds of the distributions of the two values to determine the trend fluctuations and outliers, and mark the positions of the outliers as cell anomaly marker data;

[0011] S5: Based on the cell anomaly marker data, extract the time index values of the outliers and the corresponding culture environment parameters, perform a sliding window operation on the index values, calculate the parameter drift rates in the previous and subsequent time periods, compare the drift rate values, and generate a set of dynamically adjusted culture parameters.

[0012] The smoothed cell density sequence data includes the values of the autocorrelation function, the partial autocorrelation function, the first-order difference sequence, the second-order difference sequence, the ADF statistic, and the KPSS statistic. The optimal ARIMA model parameter set specifically includes the lag order, the Bayesian information criterion value, the Akaike information criterion value, the Ljung-Box statistic of the error residual sequence, and the confidence interval. The dynamic predicted values of cell proliferation include the predicted values of the cell proliferation rate at future time points, the upper and lower bounds of the confidence interval of the prediction sequence, and the predicted cell density values. The cell anomaly marker data specifically includes the positions of the outliers, the time index values of the outliers, and the culture environment parameters. The set of dynamically adjusted culture parameters includes the time index values, the parameter drift rates in the previous and subsequent time periods, and the culture parameter adjustment values.

[0013] As a further solution of the present invention, the specific steps for obtaining the smoothed cell density sequence data are as follows:

[0014] S111: Based on the collected time - series cell density data, calculate the autocorrelation function values and partial autocorrelation function values of the sequence. Through cumulative analysis of the lag parameters of the autocorrelation function and the lag parameters of the partial autocorrelation function, extract the key parameter sets of the first - order and second - order differences, and generate difference - sequence data;

[0015] S112: Calculate the ADF statistic and KPSS statistic for the difference - sequence data. By comparing the significance probability of the ADF statistic with the critical value of the KPSS statistic, judge the stationarity of the time series, and adopt the time - series stationarization processing rule to generate stationary time - series data;

[0016] S113: Calculate the stationary cell - density characteristics of the sequence for the stationary time - series data, using the formula:

[0017]

[0018] Calculate to obtain the stationary cell - density sequence data;

[0019] where, S t represents the stationary cell - density sequence, D t,i represents the i - th component of the stationary difference data, W t,i represents the adjustment weight corresponding to the lag time t, R t,i represents the factor used to enhance the stationarity of the time series, n represents the total number of data points in the sequence, and t represents the current time point.

[0020] As a further solution of the present invention, the steps for obtaining the optimal ARIMA model parameter set are specifically as follows:

[0021] S211: Based on the stationary cell - density sequence data, extract the autocorrelation values and partial autocorrelation values under different lag orders. For multiple orders within the lag - order range, calculate the corresponding lag - correlation parameter sets. By comparing the fluctuation amplitudes in the lag - correlation parameters, screen the significant lag - order characteristics and generate a set of lag - characteristic values;

[0022] S212: Call the set of lag - characteristic values, calculate the Bayesian information criterion value and Akaike information criterion value for multiple lag orders, mark the lag order corresponding to the minimum of the two criterion values as the candidate lag order, and generate a lag - error residual sequence by calculating the error - residual sequence corresponding to the candidate lag order;

[0023] S213: For the lag - error residual sequence, use the formula:

[0024]

[0025] Calculate the Ljung-Box statistic Q value and the confidence interval boundary, adjust the residual distribution in combination with the candidate lag order, and generate an optimal ARIMA model parameter set;

[0026] Among them, Q represents the Ljung-Box statistic, n represents the number of sample data points, m represents the upper limit of the lag order, and ρ k represents the sample autocorrelation coefficient at lag k, and k represents the current lag order.

[0027] As a further solution of the present invention, the steps for obtaining the cell proliferation dynamic prediction value are specifically as follows:

[0028] S311: Based on the optimal ARIMA model parameter set, recursively calculate the predicted value of the cell proliferation rate at future time points, dynamically adjust the upper and lower bounds of the confidence interval for each prediction point according to the prediction errors at multiple time points, and generate adjusted confidence interval data;

[0029] S312: Call the adjusted confidence interval data, perform error boundary analysis on the predicted value of the cell density at each prediction time point according to the current upper and lower bounds of the confidence interval, and combine the predicted cell proliferation rate to perform error analysis and data correction to generate corrected predicted sequence data;

[0030] S313: Perform joint analysis on the corrected predicted sequence data, refer to the mutual influence between the cell density and the proliferation rate, and use the formula:

[0031]

[0032] Calculate the cell proliferation dynamic prediction value at each time point to generate the cell proliferation dynamic prediction value;

[0033] Among them, P(t) represents the cell proliferation prediction value at time point t, and D i represents the corrected cell density value, R i represents the associated proliferation rate, ΔC i represents the rate change amount at the corresponding time point, β represents the correction coefficient, and σ t represents the predicted standard error, t represents the current time point, and i represents the time point sequence, within the range from t - 5 to t.

[0034] As a further solution of the present invention, the steps for obtaining the cell abnormal marking data are specifically as follows:

[0035] S411: Based on the cell proliferation dynamic prediction value, calculate the difference between the predicted value and the actual value at the current time point, perform mean operation and standard deviation calculation on the difference value sequence, extract the difference mean and the corresponding fluctuation characteristic value, and generate the difference mean parameter;

[0036] S412: Call the difference mean parameter, calculate the sliding weighted value for the difference value sequence, generate dynamic sliding weighted distribution data through iterative calculation adjusted by the weight coefficient and the dynamic time window, and accumulate and calculate the weighted value accumulation of multiple time points based on the distribution data to generate the cumulative sum value parameter;

[0037] S413: Based on the cumulative sum value parameter, combine the upper and lower bounds of the sliding weighted distribution data to calculate the abnormal fluctuation interval, using the formula:

[0038]

[0039] Compare and analyze the cumulative sum value parameter with the abnormal fluctuation interval, mark the time points exceeding the fluctuation interval as abnormal points, and generate cell abnormal marking data;

[0040] Among them, A t represents the abnormal marking value at the current time point, W t represents the sliding weighted value, ΔM t represents the fluctuation amount of the difference mean, C i represents the standardized component of the difference cumulative value, n represents the number of cumulative time points within the sliding window, t represents the current time point, and i is the time point index within the sliding window.

[0041] As a further solution of the present invention, the acquisition steps of the dynamically adjusted culture parameter set are specifically as follows:

[0042] S511: Extract the abnormal point time index value and the corresponding culture environment parameters in the cell abnormal marking data, divide the culture parameters in the time periods before and after the abnormal point into windows according to the time sequence, calculate the parameter mean and change rate within multiple time windows, and generate a parameter drift rate sequence;

[0043] S512: Call the parameter drift rate sequence, for the drift rate within each time window, perform distribution determination by setting the upper and lower bounds of the drift rate, calculate the distribution boundary value of the abnormal time period, extract the time windows and the corresponding drift rates exceeding the distribution range, and generate abnormal drift rate data;

[0044] S513: Based on the abnormal drift rate data, combine the parameter drift rate of the normal distribution to dynamically adjust the abnormal window, using the formula:

[0045]

[0046] Calculate the adjusted culture parameter value of each abnormal window to generate a dynamically adjusted culture parameter set;

[0047] Among them, P c represents the dynamically adjusted culture parameter value, R i represents the drift rate of the i-th time window, Wi Let \(w\) represent the drift weight, \(k\) represent the number of windows, \(\alpha\) represent the adjustment coefficient, \(\Delta R\) represent the difference between the abnormal drift rate and the normal drift rate, and \(i\) be the window index.

[0048] A CIK cell preparation system for performing the above-mentioned CIK cell preparation method, the system comprising:

[0049] The cell density stabilization module collects time-series cell density data, calculates the autocorrelation function value and the partial autocorrelation function value, extracts the first-order difference sequence and the second-order difference sequence, calculates the ADF statistic and the KPSS statistic, compares the statistic range to judge the stationarity of the time series, selects the time series that meets the conditions, and generates a stabilized cell density sequence;

[0050] The model parameter optimization module extracts the lag autocorrelation value and the lag partial autocorrelation value based on the stabilized cell density sequence, calculates the Bayesian information criterion value and the Akaike information criterion value, selects the lag order corresponding to the minimum value, calculates the Ljung-Box statistic and the confidence interval of the error residual sequence, and generates an optimal model parameter set;

[0051] The cell dynamic prediction module recursively calculates the predicted value of the cell proliferation rate at future time points based on the optimal model parameter set, corrects the upper and lower bounds of the confidence interval for the predicted cell proliferation rate, and jointly calculates the cell density and the proliferation rate at future time points to generate a cell proliferation dynamic prediction value;

[0052] The anomaly detection and marking module calculates the mean difference between the predicted value and the actual value at the current time point based on the cell proliferation dynamic prediction value, calculates the sliding weighted value and the cumulative sum value, combines the upper and lower bounds of the two-value distribution to judge the trend fluctuations and anomaly points in the time series, marks the positions of the anomaly points, and generates cell anomaly marking data;

[0053] The culture parameter adjustment module extracts the time index value corresponding to the anomaly point and the culture environment parameters based on the cell anomaly marking data, performs a sliding window operation on the time index value, calculates the drift rate of the culture environment parameters within the sliding window, compares the drift rate values in the previous and subsequent time periods, extracts the difference range, and generates a dynamic adjustment value of the culture parameters.

[0054] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0055] In the present invention, through the smoothing process of the time-series cell density, the interference of data fluctuations on the analysis results is reduced, and the accuracy of subsequent predictions is improved. The extraction and optimization of the lag correlation value and the partial correlation value enhance the accuracy of the model parameters, providing reliable support for the dynamic correction and joint analysis of the future cell proliferation rate. Based on the analysis of the difference mean and the weighted cumulative value, abnormal points are accurately identified, significantly improving the efficiency of anomaly detection. By combining the anomaly point index and the drift rate calculation of environmental parameters, the dynamic adjustment of the culture conditions is realized, reducing the impact of environmental fluctuations on the cell proliferation quality and ensuring the preparation quality and efficiency. Brief Description of the Drawings

[0056] Figure 1 It is a schematic diagram of the work flow of the present invention;

[0057] Figure 2 It is a flow chart of the acquisition steps of the smoothed cell density sequence data of the present invention;

[0058] Figure 3 It is a flow chart of the acquisition steps of the optimal ARIMA model parameter set of the present invention;

[0059] Figure 4 It is a flow chart of the acquisition steps of the dynamic prediction value of cell proliferation of the present invention;

[0060] Figure 5 It is a flow chart of the acquisition steps of the cell anomaly marker data of the present invention;

[0061] Figure 6 It is a flow chart of the acquisition steps of the dynamic adjustment culture parameter set of the present invention. Detailed Description of the Invention

[0062] In order to make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0063] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as limiting the present invention. In addition, in the description of the present invention, the meaning of "a plurality of" is two or more unless otherwise specifically defined.

[0064] Example 1

[0065] Please refer to Figure 1 , the present invention provides a technical solution: a method for preparing CIK cells, comprising the following steps:

[0066] S1: Collect cell density data of time series, calculate the values of autocorrelation function and partial autocorrelation function, extract the first-order and second-order difference sequences, calculate the ADF statistic and KPSS statistic according to the difference sequences, compare the range of the statistics to judge the stationarity of the time series, and generate the stationary cell density sequence data;

[0067] S2: Based on the stationary cell density sequence data, extract the lag autocorrelation value and lag partial autocorrelation value, calculate the Bayesian information criterion value and Akaike information criterion value for different lag orders, select the lag order corresponding to the minimum of the two criterion values, calculate the Ljung-Box statistic and confidence interval of the error residual sequence, and generate the optimal ARIMA model parameter set;

[0068] S3: Based on the optimal ARIMA model parameter set, recursively calculate the predicted values of cell proliferation rate at future time points, dynamically correct the upper and lower bounds of the confidence interval according to the prediction sequence, and conduct a joint analysis of the predicted cell density and proliferation rate to generate the dynamic prediction value of cell proliferation;

[0069] S4: Based on the dynamic prediction value of cell proliferation, extract the mean difference between the predicted value and the actual value at the current time point, calculate the sliding weighted value and cumulative sum value for the mean difference, combine the upper and lower bounds of the distributions of the two values to determine the trend fluctuation and outliers, and mark the positions of the outliers as cell anomaly marker data;

[0070] S5: Based on the cell anomaly marker data, extract the time index value of the outliers and the corresponding culture environment parameters, perform a sliding window operation on the index value, calculate the parameter drift rates in the front and back time periods, compare the drift rate values, and generate the dynamically adjusted culture parameter set.

[0071] The stationary cell density sequence data includes the values of autocorrelation function, partial autocorrelation function, first-order difference sequence, second-order difference sequence, ADF statistic, KPSS statistic. The optimal ARIMA model parameter set specifically includes the lag order, Bayesian information criterion value, Akaike information criterion value, Ljung-Box statistic of the error residual sequence, confidence interval. The dynamic prediction value of cell proliferation includes the predicted values of cell proliferation rate at future time points, the upper and lower bounds of the confidence interval of the prediction sequence, the predicted cell density value. The cell anomaly marker data specifically includes the positions of the outliers, the time index values of the outliers, culture environment parameters. The dynamically adjusted culture parameter set includes the time index value, the parameter drift rates in the front and back time periods, the culture parameter adjustment value.

[0072] Please refer to Figure 2 , the specific steps for obtaining the stationary cell density sequence data are as follows:

[0073] S111: Calculate the autocorrelation function value and partial autocorrelation function value of the sequence based on the collected time-series cell density data. By performing cumulative analysis on the lag parameters of the autocorrelation function and the partial autocorrelation function, extract the key parameter sets of the first-order and second-order differences, and generate difference sequence data;

[0074] First, form an observation sequence set with the collected cell density data distributed according to the time series. Then, perform autocorrelation analysis on the cell density values at each time point in the sequence set. Through the formula Calculate the autocorrelation coefficient of the k-th lag, where is the sequence mean, and R k represents the correlation of the k-th lag. At the same time, the calculation of the partial autocorrelation coefficient needs to use the recursive method. Starting from 1 and increasing the lag order step by step, gradually eliminate the influence of the correlations of other orders. When extracting the first-order and second-order difference sequences, through the difference formula ΔX t = X t - X t-1 and Δ 2 X t = ΔX t - ΔX t-1 Calculate the difference values of each point. After completion, form a difference sequence set and generate difference sequence data.

[0075] S112: Calculate the ADF statistic and KPSS statistic for the difference sequence data. By comparing the significance probability of the ADF statistic with the critical value of the KPSS statistic, judge the stationarity of the time series, and adopt the time series stationarization processing rule to generate the stationarized time series data;

[0076] First, perform an ADF test on each value of the difference sequence data to check whether there is a unit root in the time series. The ADF test is carried out through the model formula to calculate the t statistic of the test value γ. Compare the statistic value with the critical value to judge whether to reject the unit root hypothesis. The KPSS test verifies whether the sequence is stationary by decomposing the residual sum of squares of the time series. Its statistic is calculated through the formula where S t is the partial sum sequence, and σ 2 is the variance of the sequence residuals. Finally, by comparing the significance probability p of the ADF statistic with the critical value q of the KPSS statistic, if the p value is less than the significance level and the q value does not exceed the critical range, then judge that the sequence is stationary, and generate the stationarized time series data after its stationarization processing.

[0077] S113: Calculate the stationarized time series data, extract the cell density stationarization characteristics of the sequence, using the formula:

[0078]

[0079] The smoothed cell density sequence data is calculated;

[0080] Among them, S t represents the smoothed cell density sequence, D t,i represents the i-th component of the smoothed difference data, W t,i represents the adjustment weight corresponding to the lag time t, R t,i represents the factor used to enhance the stationarity of the time series, n represents the total number of data points in the sequence, and t represents the current time point.

[0081] Formula:

[0082]

[0083] The advantage of the formula is that by introducing the lag weight coefficient W t,i and the smoothing factor R t,i to enhance the feature extraction ability of the stationary sequence, and at the same time, through the normalization of the sum of squares of the denominator, the sensitivity to outliers is optimized.

[0084] Detailed explanation of the formula and the derivation process of the formula calculation:

[0085] Perform item-by-item calculation on the smoothed time series data D t,i , introduce the lag weight coefficient W t,i and the smoothing factor R t,i . Assume D t,1 = 5, D t,2 = 4, D t,3 = 3, W t,1 = 0.6, W t,2 = 0.4, W t,3 = 0.2, R t,1 = 1.2, R t,2 = 1.1, R t,3 = 1.0, substitute into the formula:

[0086]

[0087] The result shows that the eigenvalue S t of the smoothed cell density sequence is 0.842, indicating that the feature extraction process of the sequence has been completed, and the value can be further used for the construction of the next prediction model.

[0088] Please refer to Figure 3 , the specific steps for obtaining the optimal ARIMA model parameter set are as follows:

[0089] S211: Based on the smoothed cell density sequence data, extract the autocorrelation values and partial autocorrelation values at different lag orders. For multiple orders within the lag order range, calculate their corresponding lag correlation parameter sets. By comparing the fluctuation amplitudes in the lag correlation parameters, screen out the significant lag order features and generate a lag feature value set.

[0090] First, through comparative analysis of the autocorrelation values and partial autocorrelation values at each order within the lag order range, screen out the significant differences at the corresponding order. After calculating the autocorrelation value for each order, calculate the partial autocorrelation value using the following formula: where, PACF k is the partial autocorrelation coefficient at lag order k, a k,i is the AR model parameter estimate. Through the result sets of the above two formulas, according to the significance screening conditions of the lag correlation parameters, extract the lag order features and mark the lag correlation parameters at the corresponding order, and then generate a lag feature value set.

[0091] S212: Call the lag feature value set, calculate the Bayesian Information Criterion (BIC) value and the Akaike Information Criterion (AIC) value for multiple lag orders. Mark the lag order corresponding to the minimum of the two criterion values as the candidate lag order. By calculating the error residual sequence for the lag feature values corresponding to the candidate lag order, generate a lag error residual sequence.

[0092] First, for each order in the lag feature value set, calculate the Bayesian Information Criterion value using the following formula: where, BIC is the Bayesian Information Criterion value, e t is the error term, k is the number of parameters, n is the number of sample data points. Subsequently, calculate the Akaike Information Criterion value, and the formula is: AIC = 2k - 2ln(L) where, AIC is the Akaike Information Criterion value, L is the likelihood function. By comparing the calculation results of BIC and AIC, mark the lag order corresponding to the minimum value as the candidate lag order. Finally, through error residual sequence analysis of the candidate lag order, calculate the following formula: where, e t is the error residual, X t is the true observed value, is the model predicted value, and generate a lag error residual sequence by integrating the error term.

[0093] S213: For the lag error residual sequence, use the formula:

[0094]

[0095] Calculate the Ljung-Box statistic Q value and the confidence interval boundary, and adjust the residual distribution in combination with the candidate lag order to generate an optimal ARIMA model parameter set.

[0096] Among them, Q represents the Ljung-Box statistic, n represents the number of sample data points, m represents the upper limit of the lag order, ρ k represents the sample autocorrelation coefficient at lag k, and k represents the current lag order.

[0097] Formula:

[0098]

[0099] The advantage of the formula is that by introducing the joint calculation of the lag correlation autocorrelation coefficient and the total number of sample points, the time series characteristics of the error residual distribution can be effectively captured, so as to more accurately evaluate the rationality of the model lag parameters.

[0100] Detailed explanation of the formula and the derivation process of formula calculation:

[0101] In the formula, Q represents the Ljung-Box statistic, n represents the number of sample data points, which is set to 100, m is the upper limit of the lag order, which is set to 10, ρ k is the sample autocorrelation coefficient at lag k, and is calculated through the formula:

[0102]

[0103] The autocorrelation coefficient values of each order are calculated as follows: Assume ρ1 = 0.15, ρ2 = 0.12, and so on until ρ 10 = 0.05, and substitute these parameters into the original formula:

[0104]

[0105] Expand the calculation to get:

[0106]

[0107] Calculate term by term to get the result:

[0108] Q = 100·102·0.015 = 153;

[0109] This result shows that:

[0110] The value of the Ljung-Box statistic, 153, falls within the boundary range of the confidence interval, indicating that the lag error residual distribution conforms to the white noise characteristics. The rationality of the lag order parameter is further verified through this statistic, and the optimal ARIMA model parameter set is generated in combination with the error distribution.

[0111] Please refer to Figure 4 for the specific steps to obtain the predicted value of cell proliferation dynamics:

[0112] S311: Based on the optimal ARIMA model parameter set, recursively calculate the predicted values of cell proliferation rate at future time points. According to the prediction errors at multiple time points, dynamically adjust the upper and lower bounds of the confidence interval for each prediction point to generate adjusted confidence interval data;

[0113] Extract the prediction terms and error terms of each lag parameter in the ARIMA model. By performing step-by-step recursive iteration on each lag parameter, calculate the predicted values of cell proliferation rate at future time points respectively. For each prediction point, first update the autocorrelation coefficient of the lag order, use the previously calculated residual sequence error value as one of the input values for the current prediction, calculate the cumulative amount of error and recursively update the initial value of the proliferation rate. Then call the original data of the upper and lower bounds of the confidence interval. By analyzing the error distribution characteristics of the predicted value at the current time point, adjust the upper and lower bound values of the confidence interval. By analyzing the relationship between the change range of the upper and lower bounds and the dynamic adjustment value, recalibrate the confidence interval range for each future time point to form adjusted confidence interval data. Through this adjustment process, control the uncertainty of the predicted value within a reasonable range to ensure the stability and practicality of the adjusted data.

[0114] S312: Call the adjusted confidence interval data. Based on the current upper and lower bounds of the confidence interval, perform error boundary analysis on the predicted values of cell density at each prediction time point. Combine with the predicted cell proliferation rate to perform error analysis and data correction to generate corrected predicted sequence data;

[0115] First, extract the predicted cell density values and proliferation rate values at each time point. For the cell density values, analyze the deviation ranges from the upper and lower bounds of the confidence interval, divide the deviation ranges into several intervals, calculate the corresponding standard errors for each interval, and screen out the time points with excessive errors by setting the change threshold of the residuals. Then, according to the error ranges of these time points, dynamically correct their proliferation rates, and use the corrected rate as the input parameter for the next time point. At the same time, for the correction terms of the proliferation rate, combine with the adjustment amount of the density data to calculate the corrected prediction error range, and perform step-by-step correction on all time points through multiple rounds of iteration to form corrected predicted sequence data.

[0116] S313: Conduct a joint analysis on the corrected predicted sequence data. Referring to the mutual influence between cell density and proliferation rate, use the formula:

[0117]

[0118] Calculate the dynamic predicted values of cell proliferation at each time point to generate dynamic predicted values of cell proliferation;

[0119] where, P(t) represents the predicted value of cell proliferation at time point t, D i represents the corrected cell density value, Ri Represents the associated proliferation rate, ΔC i Represents the rate change at the corresponding time point, β represents the correction coefficient, σ t Represents the predicted standard error, t represents the current time point, and i represents the time point sequence, within the range from t - 5 to t.

[0120] Formula:

[0121]

[0122] The advantage of the formula is that by introducing the interaction term of the cumulative proliferation rate and density change at historical time points, the predicted value can dynamically reflect the comprehensive influence of historical characteristics on the current time point, and the prediction accuracy is improved through the correction adjustment of the standard error term.

[0123] Detailed explanation of the formula and the derivation process of formula calculation:

[0124] D i = 80, 85, 90, 95, 100 (cell density data, measured by a density sensor);

[0125] R i = 1.2, 1.25, 1.3, 1.35, 1.4 (proliferation rate, calculated by a cell proliferation detection module);

[0126] ΔC i = 0.05, 0.04, 0.03, 0.02, 0.01 (rate change of proliferation, calculated by a first-order difference sequence);

[0127] β = 0.1 (adjustment coefficient, adjusted by error distribution analysis);

[0128] σ t = 0.02 (standard error, calculated by prediction error analysis).

[0129] Substitute into the formula:

[0130]

[0131] Expand the summation:

[0132]

[0133] Average value:

[0134]

[0135] Correction term:

[0136] 1 + β·σ t = 1 + 0.1·0.02 = 1.002;

[0137] Predicted value:

[0138] P(t) = 90.346 · 1.002 = 90.526;

[0139] This result indicates that the predicted dynamic value of cell proliferation at the current time point is 90.526, indicating that the proliferation rate and density at the current point have reached high accuracy after correction and can be used as the core result of the dynamic prediction model.

[0140] Please refer to Figure 5 , and the steps for obtaining cell abnormal marker data are specifically as follows:

[0141] S411: Based on the predicted value of cell proliferation dynamics, calculate the difference between the predicted value and the actual value at the current time point. By performing mean operation and standard deviation calculation on the difference value sequence, extract the mean difference and the corresponding fluctuation characteristic value, and generate the mean difference parameter;

[0142] By analyzing the difference between the predicted value and the actual value at each point in time, extract the difference value sequence, perform mean operation and standard deviation calculation on the difference value sequence, extract the mean difference and the fluctuation characteristic value, where the difference value is the point-to-point difference between the predicted value and the actual value. By separately calling the cell density time series data and the proliferation rate sequence data for the actual value and the predicted value at each time point, and according to the difference calculation formula (predicted value - actual value), obtain the difference value sequence; the mean difference is obtained by summing all the points in the difference value sequence and dividing by the total number of samples n, and the formula is The fluctuation characteristic value is obtained by calculating the standard deviation of the difference value sequence, and the formula is The standard deviation reflects the degree to which the data in the difference value sequence deviates from the mean. In the example, if the predicted values are [10, 15, 20] and the actual values are [12, 14, 22], then the difference value sequence is [-2, 1, -2], the mean difference is -1, and the fluctuation characteristic value is Generate the mean difference parameter.

[0143] S412: Call the mean difference parameter, calculate the sliding weighted value for the difference value sequence, generate the dynamic sliding weighted distribution data through weight coefficient adjustment and iterative calculation of the dynamic time window, and generate the cumulative sum value parameter by cumulatively calculating the weighted value cumulative amount at multiple time points based on the distribution data;

[0144] By setting the sliding window size m to a fixed value, perform weighted operation on the data at each time point and the data of the previous m - 1 time points after assigning weights, where the weight parameter is set according to a decreasing function of the distance of the time point from the current time point, and the decreasing function can be set as where k represents the distance of the time point; the formula for the sliding weighted value is The cumulative sum value is obtained by accumulating the sliding weighted values of each sliding window, and the formula is In the example, if the difference value sequence is [-2, 1, -2], the sliding window m = 2, and the decreasing function is Then the sliding weighted values are [(-2 + 1) / 1, (1 - 2) / 2], and the cumulative sum value is (-1 + -0.5) = -1.5, generating the cumulative sum value parameter.

[0145] S413: Based on the cumulative sum value parameter, combined with the upper and lower bounds of the sliding weighted distribution data, calculate the abnormal fluctuation interval, using the formula:

[0146]

[0147] Compare and analyze the cumulative sum value parameter with the abnormal fluctuation interval, mark the time points exceeding the fluctuation interval as abnormal points, and generate cell abnormal marking data;

[0148] Among them, A t represents the abnormal marking value at the current time point, W t represents the sliding weighted value, ΔM t represents the fluctuation amount of the difference mean, C i represents the standardized component of the difference cumulative value, n represents the number of cumulative time points within the sliding window, t represents the current time point, and i is the time point index within the sliding window.

[0149] Formula:

[0150]

[0151] The benefit of the formula is that by introducing the sliding weighted value W t and the fluctuation amount of the difference mean ΔM t , combined with the standardized component C i of the cumulative value, the accuracy of abnormal point marking is improved;

[0152] Detailed explanation of the formula and the derivation process of formula calculation:

[0153] The sliding weighted value W t is calculated from the formula in the previous paragraph, and the fluctuation amount of the difference mean ΔM t = |difference mean t - difference mean t-1 |, and the standardized component of the cumulative value In the example, if the sliding weighted value W t = -1, the fluctuation amount of the difference mean ΔM t = 2, the cumulative value is [-2, 1, -2], and the difference mean is -1, then

[0154] C i= [-2 / -1, 1 / -1, -2 / -1] = [2, -1, 2];

[0155] Summation

[0156] The formula value is

[0157] This result indicates that the calculated abnormal marker value A t = 0.63. The time points within the abnormal fluctuation range can be marked as abnormal points, generating cell abnormal marker data.

[0158] Please refer to Figure 6 , and the specific steps for obtaining the culture parameter set by dynamic adjustment are as follows:

[0159] S511: Extract the abnormal point time index values and corresponding culture environment parameters in the cell abnormal marker data. Divide the culture parameters in the time periods before and after the abnormal points into windows according to the time series, calculate the parameter means and change rates within multiple time windows, and generate a parameter drift rate sequence;

[0160] Locate the time periods before and after the abnormal points through the time index values, extract the culture environment parameter sets within each time period, divide these sets into multiple windows in chronological order, calculate the parameter means within each window, where the parameters include culture temperature, carbon dioxide concentration, nutrient solution concentration, etc. The culture temperature is directly obtained through a sensor, the carbon dioxide concentration is recorded and output by an online monitoring device, and the nutrient solution concentration is recorded through laboratory quantitative analysis. For window division, a fixed time span is used for sliding, and the end point of the previous window is used as the starting point of the next window. Through the formula: where, P i represents the parameter value at the i-th time point, n represents the number of time points within the window. After calculating the means within each window, calculate the change rate using the means of adjacent windows through the formula: Extract the change rates of all windows as the parameter drift rate sequence, and extract the time periods of abnormal change points by gradually comparing the parameter change rates of each window to generate the parameter drift rate sequence.

[0161] S512: Call the parameter drift rate sequence. For the drift rate within each time window, perform distribution determination by setting the upper and lower bounds of the drift rate, calculate the distribution boundary values of the abnormal time periods, and extract the time windows and corresponding drift rates that exceed the distribution range to generate abnormal drift rate data;

[0162] For the drift rate of each time window in the sequence, upper and lower bounds are set to determine the abnormal range of parameters. The upper and lower bounds are set according to the parameter change records. For example, the upper and lower bounds of the drift rate are set to the average value of the normal change range ± twice the standard deviation. The parameters are obtained through long-term cultivation environment data monitoring, and its standard deviation is calculated by the formula: where μ is the mean drift rate, P i is the drift rate of a single time window, n is the number of windows. By comparing the drift rate with the upper and lower bounds one by one, the time windows that exceed the upper and lower bounds are screened out, and the drift rates corresponding to these windows are recorded as abnormal drift rate values. When extracting these values, the number of consecutive times the drift rate exceeds the upper and lower bounds and the distribution characteristics within each time window are considered. For example, if the cumulative number of exceedances is greater than 2 times, it is judged as an abnormal interval. The time windows of all abnormal drift rates and the corresponding drift rate data are sorted out and output to generate abnormal drift rate data.

[0163] S513: Based on the abnormal drift rate data, combined with the parameter drift rate of the normal distribution, the abnormal windows are dynamically adjusted using the formula:

[0164]

[0165] Calculate the adjusted cultivation parameter values for each abnormal window to generate a dynamically adjusted cultivation parameter set;

[0166] where P c represents the dynamically adjusted cultivation parameter value, R i represents the drift rate of the i-th time window, W i represents the drift weight, k represents the number of windows, α represents the adjustment coefficient, ΔR represents the difference between the abnormal drift rate and the normal drift rate, and i is the window index.

[0167] Formula:

[0168]

[0169] The benefit of the formula is that by introducing the drift weight and the difference between abnormal and normal, the cultivation parameter values can be dynamically adjusted to ensure the stability and continuity of the cultivation environment.

[0170] Detailed explanation of the formula and the derivation process of formula calculation:

[0171] R i : The drift rate value obtained by formula calculation. For example, the drift rate of a certain time window is calculated to be 0.15;

[0172] W i : Set according to the weight parameter of the time window. For example, the weight of the window close to the abnormal point is 0.8, and the weight of the window far from the abnormal point is 0.5;

[0173] α: The adjustment coefficient is determined by the sensitivity analysis of the change in the drift rate in historical data. For example, it is adjusted to 1.2 according to the trend of the drift rate change;

[0174] ΔR: It is calculated by the difference between the abnormal drift rate and the average of the normal drift rates. For example, if an abnormal drift rate is 0.2 and the average of the normal drift rates is 0.1, then ΔR = 0.1;

[0175] k: The total number of time windows, assumed to be 5 windows.

[0176] Substitute into the calculation:

[0177] 1. Calculate the weighted sum of the drift rates:

[0178] Weighted sum = (0.15·0.8) + (0.12·0.6) + (0.18·0.7) + (0.1·0.5) + (0.08·

[0179] 0.4) = 0.12 + 0.072 + 0.126 + 0.05 + 0.032 = 0.4;

[0180] 2. Calculate the total weight sum:

[0181] Weight sum = 0.8 + 0.6 + 0.7 + 0.5 + 0.4 = 3.0;

[0182] 3. Calculate the first partial result:

[0183]

[0184] 4. Calculate the second partial result:

[0185]

[0186] 5. Combine the results:

[0187] P c = 0.1333 + 0.024 = 0.1573;

[0188] This result indicates that the value of the culture parameter after dynamic adjustment is 0.1573. Comparing it with the range of the benchmark normal parameters shows that this value is within the reasonable adjustment range and can be used as the dynamically adjusted value of the culture parameter at the current time point, obtaining a dynamically adjusted set of culture parameters.

[0189] A CIK cell preparation system, which is used to execute the above CIK cell preparation method. The system includes:

[0190] The cell density stabilization module collects time series cell density data, calculates the values of the autocorrelation function and the partial autocorrelation function, extracts the first-order difference sequence and the second-order difference sequence, calculates the ADF statistic and the KPSS statistic, compares the statistic ranges to judge the stationarity of the time series, selects the time series that meet the conditions, and generates a stabilized cell density sequence;

[0191] The model parameter optimization module, based on the stabilized cell density sequence, extracts the lag autocorrelation value and the lag partial autocorrelation value, calculates the Bayesian information criterion value and the Akaike information criterion value, selects the lag order corresponding to the minimum value, calculates the Ljung-Box statistic and the confidence interval of the error residual sequence, and generates an optimal model parameter set;

[0192] The cell dynamics prediction module, based on the optimal model parameter set, recursively calculates the predicted values of the cell proliferation rate at future time points. For the predicted cell proliferation rate, it corrects the upper and lower bounds of the confidence interval, and jointly calculates the cell density and the proliferation rate at future time points to generate the predicted values of cell proliferation dynamics;

[0193] The anomaly detection and marking module, based on the predicted values of cell proliferation dynamics, calculates the mean difference between the predicted value and the actual value at the current time point, calculates the sliding weighted value and the cumulative sum value, combines the upper and lower bounds of the two value distributions to judge the trend fluctuations and anomaly points in the time series, marks the positions of the anomaly points, and generates cell anomaly marking data;

[0194] The culture parameter adjustment module, based on the cell anomaly marking data, extracts the time index value corresponding to the anomaly point and the culture environment parameters, performs a sliding window operation on the time index value, calculates the drift rate of the culture environment parameters within the sliding window, compares the drift rate values in the previous and subsequent time periods, extracts the difference range, and generates the dynamic adjustment value of the culture parameters.

[0195] The above are only the preferred embodiments of the present invention, and do not limit the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention without departing from the technical solution content of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for preparing CIK cells, characterized in that: The following steps are involved: S1: Collect cell density data of time series, calculate the autocorrelation function value and partial autocorrelation function value, extract the first-order and second-order difference sequences, calculate the ADF statistic and KPSS statistic based on the difference sequence, compare the statistical range to judge the stationarity of the time series, and generate stable cell density series data; S2: Based on the stabilized cell density sequence data, extract the lagged autocorrelation value and the lagged partial autocorrelation value, calculate the Bayesian information criterion value and the Akaike information criterion value for the differentiated lag order, select the lag order corresponding to the minimum of the two criterion values, calculate the Ljung-Box statistic and confidence interval of the error residual sequence, and generate the optimal ARIMA model parameter set; S3: Based on the optimal ARIMA model parameter set, recursively calculate the predicted value of cell proliferation rate at future time points, dynamically correct the upper and lower bounds of the confidence interval according to the predicted sequence, jointly analyze the predicted cell density and proliferation rate, and generate a dynamic prediction value of cell proliferation; S4: Based on the cell proliferation dynamic prediction value, extract the difference mean between the prediction value and the actual value at the current time point, calculate the sliding weighted value and the cumulative sum value for the difference mean, determine the trend fluctuation and abnormal point by combining the upper and lower bounds of the distribution of the two values, and mark the abnormal point position as cell abnormality mark data; S5: Based on the cell abnormality mark data, extract the time index value of the abnormal point and the corresponding culture environment parameter, perform a sliding window operation on the index value, calculate the parameter drift rate of the previous and next time periods, compare the drift rate values, and generate a dynamically adjusted culture parameter set.

2. The method for preparing CIK cells according to claim 1, characterized in that: The stabilized cell density sequence data includes autocorrelation function values, partial autocorrelation function values, first-order difference sequences, second-order difference sequences, ADF statistics, and KPSS statistics. The optimal ARIMA model parameter set specifically includes lag order, Bayesian information criterion value, Akaike information criterion value, Ljung-Box statistic of error residual sequence, and confidence interval. The dynamic prediction value of cell proliferation includes the predicted value of cell proliferation rate at future time points, upper and lower bounds of predicted sequence confidence intervals, and predicted cell density values. The cell abnormality marking data specifically includes abnormal point positions, abnormal point time index values, and culture environment parameters. The dynamically adjusted culture parameter set includes time index values, drift rates of parameters in previous and next time periods, and culture parameter adjustment values.

3. The method for preparing CIK cells according to claim 2, characterized in that: The steps for obtaining the stabilized cell density sequence data are specifically as follows: S111: Based on the collected time series cell density data, the autocorrelation function value and the partial autocorrelation function value of the sequence are calculated, and the key parameter sets of the first-order and second-order differences are extracted by cumulatively analyzing the lag parameters of the autocorrelation function and the lag parameters of the partial autocorrelation function to generate differential sequence data; S112: Calculate the ADF statistic and the KPSS statistic for the differential sequence data, judge the stationarity of the time series by comparing the significance probability of the ADF statistic with the critical value of the KPSS statistic, and use the time series stationarity processing rule to generate stationarized time series data; S113: Calculate the stabilized time series data to extract the cell density stabilized feature of the sequence using the formula: The stabilized cell density sequence data is calculated; Among them, S t represents the stable cell density series, D t,i represents the ith component of the stabilized difference data, W t,i represents the adjustment weight corresponding to the lag time t, R t,i Represents the factor used to enhance the stationarity of the time series, n represents the total number of data points in the series, and t represents the current time point.

4. The method for preparing CIK cells according to claim 3, characterized in that: The steps for obtaining the optimal ARIMA model parameter set are specifically as follows: S211: based on the stabilized cell density sequence data, extracting autocorrelation values ​​and partial autocorrelation values ​​under different lag orders, calculating the corresponding lag correlation parameter set for multiple orders within the lag order range, screening significant lag order features by comparing the fluctuation ranges in the lag correlation parameters, and generating a lag feature value set; S212: calling the set of lagged eigenvalues, calculating the Bayesian information criterion value and the Akaike information criterion value for multiple lagged orders, marking the lag order corresponding to the minimum of the two criterion values ​​as the candidate lag order, and generating a lagged error residual sequence by calculating the error residual sequence of the lagged eigenvalues ​​corresponding to the candidate lag order; S213: For the lagged error residual sequence, the formula is used: Calculate the Ljung-Box statistic Q value and confidence interval boundaries, adjust the residual distribution in combination with the candidate lag order, and generate the optimal ARIMA model parameter set; Where Q represents the Ljung-Box statistic, n represents the number of sample data points, m represents the upper limit of the lag order, and ρ k Represents the sample autocorrelation coefficient of lag k, where k represents the current lag order.

5. The method for preparing CIK cells according to claim 4, characterized in that: The step of obtaining the dynamic prediction value of cell proliferation is specifically as follows: S311: recursively calculating the predicted value of the cell proliferation rate at a future time point based on the optimal ARIMA model parameter set, dynamically adjusting the upper and lower bounds of the confidence interval of each prediction point according to the prediction errors at multiple time points, and generating adjusted confidence interval data; S312: calling the adjusted confidence interval data, performing error boundary analysis on the predicted value of cell density at each prediction time point according to the upper and lower bounds of the current confidence interval, performing error analysis and data correction in combination with the predicted cell proliferation rate, and generating corrected prediction sequence data; S313: Jointly analyzing the corrected predicted sequence data, referring to the mutual influence of cell density and proliferation rate, using the formula: Calculate the dynamic prediction value of cell proliferation at each time point to generate the dynamic prediction value of cell proliferation; Among them, P(t) represents the predicted value of cell proliferation at time point t, D i represents the corrected cell density value, R i represents the associated proliferation rate, ΔC i represents the rate change at the corresponding time point, β represents the correction coefficient, σ t Represents the standard error of the prediction, t represents the current time point, and i represents the time point sequence, ranging from t-5 to t.

6. The method for preparing CIK cells according to claim 5, characterized in that: The steps for obtaining the abnormal cell marker data are specifically as follows: S411: Based on the cell proliferation dynamic prediction value, the difference between the prediction value at the current time point and the actual value is calculated, and the difference mean and the corresponding fluctuation characteristic value are extracted by performing mean operation and standard deviation calculation on the difference value sequence to generate a difference mean parameter; S412: calling the difference mean parameter, calculating the sliding weighted value for the difference value sequence, generating dynamic sliding weighted distribution data through weight coefficient adjustment and iterative calculation of the dynamic time window, and accumulating and calculating the weighted value accumulation of multiple time points according to the distribution data to generate a cumulative sum parameter; S413: Based on the cumulative sum parameter, combined with the upper and lower bounds of the sliding weighted distribution data, calculate the abnormal fluctuation range using the formula: Compare and analyze the cumulative sum value parameters with the abnormal fluctuation range, mark the time points beyond the fluctuation range as abnormal points, and generate cell abnormality marking data; Among them, A t Represents the abnormal mark value at the current time point, W t Represents the sliding weight value, ΔM t represents the fluctuation of the mean difference, C i Represents the normalized component of the difference cumulative value, n represents the cumulative number of time points in the sliding window, t represents the current time point, and i is the time point index in the sliding window.

7. The method for preparing CIK cells according to claim 6, characterized in that: The steps for obtaining the dynamically adjusted training parameter set are specifically as follows: S511: extracting the abnormal point time index value and the corresponding culture environment parameter in the abnormal cell mark data, dividing the culture parameters in the time period before and after the abnormal point into windows according to the time series, calculating the parameter mean and change rate in multiple time windows, and generating a parameter drift rate sequence; S512: calling the parameter drift rate sequence, for the drift rate in each time window, performing distribution determination by setting the upper and lower bounds of the drift rate, calculating the distribution boundary value of the abnormal time period, extracting the time window and the corresponding drift rate that exceed the distribution range, and generating abnormal drift rate data; S513: Based on the abnormal drift rate data and in combination with the parameter drift rate of the normal distribution, the abnormal window is dynamically adjusted using the formula: Calculating the adjustment training parameter value of each abnormal window and generating a dynamic adjustment training parameter set; Among them, P c represents the dynamically adjusted cultivation parameter value, R i represents the drift rate of the i-th time window, W i represents the drift weight, k represents the number of windows, α represents the adjustment coefficient, ΔR represents the difference between the abnormal drift rate and the normal drift rate, and i represents the window index.

8. A CIK cell preparation system, characterized in that: According to any one of claims 1 to 7, the CIK cell preparation method, the system comprises: The cell density stabilization module collects time series cell density data, calculates the autocorrelation function value and the partial autocorrelation function value, extracts the first-order difference sequence and the second-order difference sequence, calculates the ADF statistic and the KPSS statistic, compares the statistical range to determine the stability of the time series, selects the time series that meets the conditions, and generates a stabilized cell density sequence; The model parameter optimization module extracts the lagged autocorrelation value and the lagged partial autocorrelation value based on the stabilized cell density sequence, calculates the Bayesian information criterion value and the Akaike information criterion value, selects the lag order corresponding to the minimum value, calculates the Ljung-Box statistic and confidence interval of the error residual sequence, and generates an optimal model parameter set; The cell dynamic prediction module recursively calculates the predicted value of the cell proliferation rate at a future time point based on the optimal model parameter set, corrects the upper and lower bounds of the confidence interval for the predicted cell proliferation rate, and jointly calculates the cell density and proliferation rate at the future time point to generate a cell proliferation dynamic prediction value; The abnormality detection and marking module calculates the mean difference between the predicted value and the actual value at the current time point based on the dynamic prediction value of cell proliferation, calculates the sliding weighted value and the cumulative sum value, combines the upper and lower bounds of the two value distributions to determine the trend fluctuations and abnormal points in the time series, marks the position of the abnormal points, and generates cell abnormality marking data; The culture parameter adjustment module extracts the time index value and culture environment parameters corresponding to the abnormal point based on the cell abnormality mark data, performs a sliding window operation on the time index value, calculates the drift rate of the culture environment parameter within the sliding window, compares the drift rate values ​​of the previous and next time periods, extracts the difference range, and generates a dynamic adjustment value of the culture parameter.

Images