Organ-like automatic culture monitoring method and system

Abstract

The invention provides an organoid automatic culture monitoring method and system, and the method comprises the steps: building an initial model of fluid shear stress distribution according to adhesion strength distribution of matrigel and a culture dish surface measured by an atomic force microscope or a single molecule force spectrometer, and obtaining initial shear stress data; by controlling the Reynolds number and the capillary number within a specific range, the relative influence of the inertia effect and the viscosity effect on flowing is adjusted, the liquid flowing speed is adjusted, and the flow speed adjustment amount is output; under the adjusted pressure control parameters, microscopic morphology characteristics of the culture medium and the organoid interface are obtained, and whether the risk of local matrigel stripping exists or not is judged by comparing and analyzing the image sequence; and iteratively optimizing the pressure control process according to the effect of accurate pressure control and by referring to real-time data monitoring feedback until the local stripping risk is eliminated to an acceptable level.

Description

Technical Field

[0001] The present invention relates to the field of monitoring technologies, and particularly to an automated culture monitoring method and system for organoids. Background Art

[0002] In the research of the automated culture monitoring system for organoids, the design of the culture medium replacement strategy faces a key technical problem: the difference in fluid shear stress caused by the non-uniformity of the adhesion strength of the basement membrane matrix. When replacing the culture medium, the liquid handling module needs to precisely control the liquid in the culture dish to ensure the integrity of the organoids. However, due to the difference in the adhesion strength between the basement membrane matrix and the surface of the culture dish in different regions, uneven shear stress will be generated when the liquid flows. This non-uniformity may cause the peeling of the basement membrane matrix in local areas, thereby affecting the structural stability of the organoids. During the culture medium replacement process, the pressure control scheme of the liquid handling module needs to be dynamically adjusted according to the real-time monitored fluid shear stress data. However, this adjustment faces two main challenges. First, the local difference in the adhesion strength of the basement membrane matrix is difficult to accurately measure in real time by conventional sensors, resulting in insufficient feedback information for pressure control. Second, the complexity of liquid flow and the instability of fluid mechanics make it difficult to meet the requirements of the structural integrity of organoids in terms of the accuracy and response speed of pressure regulation. When the risk of local peeling appears, the liquid handling module needs to precisely adjust the pressure in a specific area without affecting other areas. This adjustment needs to consider the inertial effect, viscous effect of the fluid, and the mechanical properties of the basement membrane matrix. However, due to the complexity of the fluid mechanics model and the non-linear mechanical behavior of the basement membrane matrix, it is difficult to achieve an ideal state in terms of the accuracy of pressure control. This inaccurate control may cause the organoids to be subjected to excessive shear stress in local areas, thereby triggering local damage or deformation of the structure and affecting their long-term culture effect. Summary of the Invention

[0003] To solve the above technical problems, the present invention provides an automated culture monitoring method for organoids, mainly including:

[0004] According to the adhesion strength distribution of the basement membrane matrix and the surface of the culture dish measured by an atomic force microscope or a single molecule force spectrometer, establish an initial model of the fluid shear stress distribution and obtain initial shear stress data;

[0005] For the initial shear stress data, adopt the fluid mechanics simulation method under non-Newtonian fluid to conduct fluid mechanics simulation, and calculate the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation;

[0006] By controlling the Reynolds number and the capillary number within a specific range, adjust the relative influence of the inertial effect and the viscous effect on the flow, and then adjust the liquid flow velocity, and output the flow velocity adjustment amount;

[0007] Using the Bernoulli equation, the adjusted flow velocity is converted into a change in the pressure control parameter. According to the established relationship between flow velocity and pressure, the pressure control parameter for the processed liquid is dynamically adjusted to obtain the adjusted pressure control parameter.

[0008] Under the adjusted pressure control parameter, the microscopic morphological characteristics of the interface between the culture medium and the organoid are obtained, and by comparing and analyzing the image sequence, it is judged whether there is a risk of local Matrigel peeling.

[0009] If there is a local peeling risk, then according to the rheological properties of the Matrigel, the constitutive equation or boundary condition of the fluid shear stress distribution is corrected, and the change trend of the shear stress in each region during the culture medium replacement process is recalculated.

[0010] Through an adaptive algorithm, with the change trend of the shear stress in each region during the recalculated culture medium replacement process as the input of the algorithm, the pressure adjustment amount of the local region is determined, and the pressure adjustment amount of the local region is applied to the hydrodynamic model as an additional boundary condition.

[0011] According to the effect of precise pressure control and referring to the real-time data monitoring feedback, the pressure control process is iteratively optimized until the local peeling risk is eliminated to an acceptable level.

[0012] Furthermore, based on the adhesion strength distribution between the Matrigel and the culture dish surface measured by an atomic force microscope or a single-molecule force spectrometer, an initial model of the fluid shear stress distribution is established to obtain the initial shear stress data, including: collecting the force data between the Matrigel and the culture dish surface by an atomic force microscope, filtering out the sampling noise interference through Fourier transform, and smoothing the sampled adhesion strength curve to obtain a continuous adhesion strength distribution curve. Measuring the surface morphology of the Matrigel by a force spectrometer, using a Gaussian filter to denoise the surface morphology data, and calculating the surface roughness distribution map from the denoised data. For the continuous adhesion strength distribution curve, a radial basis function neural network is used to construct the mapping relationship between the adhesion site density and the surface roughness distribution map, and the surface morphological characteristic data of the Matrigel surface is obtained from the mapping relationship. According to the surface morphological characteristic data of the Matrigel surface, a hexagonal grid unit is established, the grid unit size is set to one-tenth of the Matrigel surface characteristic scale, and the fluid boundary condition parameters are obtained from the grid unit nodes. The Navier-Stokes equation containing the viscous term and the pressure term is solved by the finite difference method, and the flow field velocity distribution data is obtained from the solution result. For the flow field velocity distribution data, the kernel function calculation method in smoothed particle hydrodynamics is used to obtain the stress tensor, and the fluid shear stress data is obtained from the projection of the main direction of the stress tensor. According to the fluid shear stress data, a probability density distribution function is established, and the maximum entropy principle is used to determine the shear stress threshold interval, and the shear stress data is normalized within the threshold interval to obtain the initial shear stress distribution data.

[0013] Furthermore, for the initial shear stress data, a hydrodynamic simulation method under non-Newtonian fluid is used for hydrodynamic simulation. By solving the Navier-Stokes equation and the continuity equation, the fluid viscosity coefficient is calculated, including: constructing a power-law fluid constitutive equation τ = K·γ n , where τ is the shear stress, K is the consistency coefficient, γ is the shear rate, n is the flow index, calculating the medium viscosity change curve from the shear rate data distribution, and obtaining the fluid motion state value by calculating the Reynolds number through the characteristic length and fluid density. A tetrahedral mesh is established for the petri dish geometry, and the side length of the mesh element is taken as one-tenth of the petri dish characteristic size. The flow rate and pressure are set from the inlet boundary condition, and the pressure gradient is set from the outlet boundary condition. The octree mesh refinement algorithm is used to encrypt the hydrodynamic boundary layer region to obtain the mesh node distribution. A momentum conservation equation is constructed according to the medium viscosity change curve, and the equation includes an inertial term, a pressure term, a viscous term, and a body force term. A relationship function between fluid density and pressure is established, and the turbulent viscosity coefficient is determined from the fluid motion state value. The Galerkin finite element method is used to discretize the momentum conservation equation, a trial function space is constructed by selecting linear basis functions, and the node displacement is calculated through a sparse matrix solver to obtain the fluid pressure field distribution. The velocity gradient is calculated for the fluid pressure field distribution, the spatial derivative is calculated using the central difference scheme to obtain the velocity field distribution, and the stress tensor is constructed from the velocity field distribution. The shear stress component is calculated according to the projection of the stress tensor principal axis direction, and the shear stress component is sampled in the time domain, and the shear stress change trend is obtained through the sampled data sequence. It can also be implemented in this way, obtaining the data required for the principal axis direction angle and projection calculation method. Through the principal axis direction angle and projection calculation method, the stress tensor value is processed to obtain the shear stress component. The shear stress component is sampled in the time domain to obtain the time domain sampling points and form a data sequence value. Through the data sequence value, a linear regression algorithm is used to obtain the shear stress change trend line. If the slope of the change trend line exceeds the preset threshold, the pressure field distribution is adjusted through the fluid characteristic parameters, and the spatial derivative value is recalculated. According to the adjusted spatial derivative value, the velocity field distribution is updated to obtain a new shear stress change trend line.

[0014] Furthermore, by controlling the Reynolds number and the capillary number within a specific range, the relative influence of the inertial effect and the viscous effect on the flow is adjusted, and then the liquid flow velocity is adjusted, and the flow velocity adjustment amount is output, including: calculating the fluid viscosity coefficient value according to the shear stress change trend, constructing the Reynolds number Re = ρvL / μ from the fluid density and the characteristic length, and calculating the capillary number Ca = μv / σ through the surface tension coefficient σ and the fluid velocity v, and setting the Reynolds number range between 100 and 200. A stress tensor is constructed for the fluid viscosity coefficient value The principal stress component is calculated from the stress tensor eigenvalue, and the viscous effect contribution value is obtained through the principal stress component and the pressure gradient △p. The inertial force is calculated according to the change of the fluid pressure gradient Calculate the change in kinetic energy ΔEk from the fluid density and velocity curve, and obtain the contribution value of the inertial effect from the change in kinetic energy and the inertial force. Construct a characteristic parameter λ = Fi / Fv for the ratio of the contribution values of the viscous effect and the inertial effect, and use a proportional-integral controller to adjust the inlet flow rate. The controller gain coefficient is optimized through the Jacobian matrix. Adjust the rotational speed of the peristaltic pump according to the controller output signal, obtain the flow rate change value from the rotational speed change curve, and calculate the flow velocity adjustment amount through the flow rate change value and the channel cross-sectional area. Verify whether the Reynolds number and the capillary number meet the range constraints for the flow velocity adjustment amount. If it exceeds the range, recalculate the controller output signal until the flow velocity adjustment amount that meets the constraint conditions is obtained.

[0015] Furthermore, using Bernoulli's equation, convert the output flow velocity adjustment amount into the change amount of the pressure control parameter. According to the established relationship between the flow velocity and the pressure, dynamically adjust the pressure control parameter of the processed liquid to obtain the adjusted pressure control parameter, including: substituting the flow velocity adjustment amount value and the fluid density parameter into Bernoulli's equation p1 + ρgh1 + 1 / 2ρv1 2 = p2 + ρgh2 + 1 / 2ρv2 2 , construct the potential energy term from the liquid level height difference h1 - h2, calculate the velocity ratio v1 / v2 through the cross-sectional area ratio A1 / A2, and calculate the total pressure change value from the potential energy term and the dynamic pressure term. Calculate the friction loss amount h = λ(l / d)(v 2 / 2g) for the pipe resistance coefficient λ, construct the pressure loss distribution curve p = f(x) from the fluid pressure drop value, and obtain the actual pressure adjustment amount through the total pressure change value and the pressure loss correction. Construct the pressure sensor response curve y = kx + b according to the actual pressure adjustment amount, calculate the slope k and intercept b of the response curve using the least squares method, and divide the pressure control interval from the slope of the response curve. Establish the Kalman filter state equation for the pressure sensor data, set the process noise covariance Q and the observation noise covariance R, and filter the pressure data through the Kalman gain K. Use an adaptive neural network ANN to establish the mapping between pressure and flow velocity. The input layer of the network includes pressure and flow velocity parameters, the hidden layer uses the hyperbolic tangent activation function, and the output layer corresponds to the pressure control parameter. Conduct online verification according to the pressure control parameter. If the pressure control parameter exceeds the set range, return to the adaptive neural network for retraining until the adjusted pressure control parameter that meets the control accuracy requirements is obtained.

[0016] Furthermore, obtain the pressure sensor data, use a Kalman filter to establish a state equation to obtain the filtered pressure data, establish the mapping relationship between the filtered pressure data and the flow rate through an adaptive neural network, and determine whether the pressure control parameters exceed the set range. If they exceed, return to the neural network for retraining until the pressure control parameters that meet the control accuracy requirements are obtained. Generate an initial pressure distribution curve from the pipeline flow rate data, and obtain real-time sensor data from the initial pressure distribution curve. Process the real-time sensor data using a Kalman filter, and construct a dynamic state equation to obtain the smoothed pressure data. Establish a flow rate distribution mapping relationship through the smoothed pressure data, and determine whether the flow rate distribution conforms to the preset range. If the flow rate distribution exceeds the preset range, adjust the mapping relationship through an adaptive neural network to obtain the updated flow rate data. Calculate the pressure adjustment coefficient based on the updated flow rate data, and determine the change amount of the control parameter from the pressure adjustment coefficient. Compare the change amount of the control parameter with the preset threshold. If it exceeds, return to adjust the mapping relationship until stable control parameters are obtained. Exemplarily, the pressure data collected by the pressure sensor shows along the pipeline distribution curve that the pressure loss decays non-linearly with the increase of distance, and its distribution function can be described as p(x) = p0e^(-kx), where p0 is the initial pressure value of 10 Pa and k is the attenuation coefficient of 2. Use a Kalman filter to process the pressure sensor data, and set the state equation as x(k) = Ax(k - 1) + Bu(k) + w(k), where the state transition matrix A is 95, the control matrix B is 05, and the covariance Q of the process noise w(k) is set to 02. The observation equation is z(k) = Hx(k) + v(k), the observation matrix H is 1, and the covariance R of the observation noise v(k) is set to 01. Iteratively filter the pressure data through the Kalman gain K = 8. The mean square error of the filtered pressure data is reduced from 1 Pa to 03 Pa, and the dynamic response time is 20 milliseconds. Based on the filtered pressure data, use an adaptive neural network to establish the mapping relationship between pressure and flow rate. The network structure is 2 neurons in the input layer (pressure, flow rate), 10 neurons in the hidden layer, and 1 neuron in the output layer, which is the pressure control parameter. The training data includes 200 sets of corresponding relationships between pressure and flow rate. Use the Levenberg-Marquardt algorithm for optimization, set the learning rate to 005, and the number of training rounds to 2000 times. The mean square error converges to 0005. Compare the pressure control parameter output by the neural network with the set target value. When the deviation exceeds 05 Pa, automatically trigger the neural network for retraining, adjust the weight matrix W and the bias vector b, and the number of retraining rounds is 500 times until the deviation is controlled within 03 Pa. During the liquid processing process, the stability of the pressure control parameter directly affects the operation efficiency of the system. Through the above method, the fluctuation range of the pressure control parameter is controlled within ±05 Pa, meeting the system accuracy requirements.

[0017] Furthermore, under the adjusted pressure control parameters, obtain the microscopic morphological characteristics of the interface between the culture medium and the organoid, and judge whether there is a risk of local Matrigel peeling by comparing and analyzing the image sequence, including: according to the adjusted pressure control parameters, obtain the original image sequence of the interface between the culture medium and the organoid from a confocal microscope, set the sampling time interval to 5 seconds, the acquired image size to 1024×1024 pixels, and remove image noise through a bilateral filter with a spatial standard deviation of 5 and a value domain standard deviation of 0.1 to obtain the filtered image. Extract the fluorescence intensity distribution for the filtered image, calculate the interface thickness distribution from the gray gradient map, perform multi-scale decomposition on the image through three-layer discrete wavelet transform, and select the high-frequency coefficients to construct the interface morphological feature map. Construct a deep neural network composed of a convolutional layer, a pooling layer, and a fully connected layer according to the interface morphological feature map, with the network input being an image block of 256×256 pixels, extract the interface features through max pooling and batch normalization, and quantify the interface undulation degree from the phase distribution map. Perform opening and closing operations on the interface features using a circular structuring element with a radius of 3 pixels, extract the connected regions from the morphological processing results, and identify the interface peeling regions through a region labeling algorithm to obtain the peeling risk distribution map. Calculate the statistical value of the peeling area for the peeling risk distribution map, set the peeling area threshold to 5% of the total area, and if the peeling area exceeds the threshold in three consecutive frames of images, output a peeling risk warning signal.

[0018] Furthermore, if there is a local peeling risk, then according to the rheological properties of the Matrigel, correct the constitutive equation or boundary conditions of the fluid shear stress distribution, and recalculate the shear stress change trend in each region during the culture medium replacement process, including: calculate the local strain rate γ according to the distribution of the Matrigel peeling region, extract the rheological characteristic parameters of the Matrigel from the mechanical response curve, and construct a viscoelastic parameter group composed of the elastic modulus G0 and the viscosity coefficient η0 through the generalized Maxwell equation σ(t) = G0γ(t) + η0dγ(t) / dt to obtain the corrected constitutive equation. Calculate the local stress transfer coefficient for the corrected constitutive equation, construct the stress-strain relationship σ = f(ε) from the material strength index, and solve the stress equilibrium equation through tetrahedral element discretization to obtain the stress distribution boundary conditions. Construct the three-dimensional stress tensor σ according to the corrected constitutive equation and boundary conditions ij , from the principal value equation |σ ij - λδ ijSolve the principal stress components when | = 0, and predict the stress transfer path through a deep neural network composed of a convolutional layer and a fully connected layer. Divide the stress distribution law into three intervals: low stress area, medium stress area, and high stress area. Set the strain rate threshold to twice the initial strain rate, and obtain the stress change curve from the continuously sampled stress data sequence. Calculate the stress gradient according to the stress change curve, calculate the stress time derivative using the central difference formula, and solve the stress spatial distribution through the forward difference formula to obtain the corrected shear stress change trend. Conduct a finite element simulation verification for the corrected shear stress change trend, extract the node displacement and stress distribution from the simulation results, and compare with the experimental measurement data to verify the correction results.

[0019] Furthermore, through an adaptive algorithm, using the recalculated shear stress change trend in each region during the medium replacement process as the algorithm input, determine the pressure adjustment amount in the local area, and apply the pressure adjustment amount in the local area as an additional boundary condition to the hydrodynamic model, including: construct a local peeling index function I = Σ(τi - τ0)2 / N according to the shear stress change trend, where τi is the real-time shear stress and τ0 is the critical shear stress, extract the pressure adjustment target value from the stress gradient distribution, establish a state-action value function through a deep Q-learning network to obtain the pressure adjustment parameter. Calculate the local intervention intensity coefficient k = △p / p0 for the pressure adjustment parameter, extract the oscillation frequency and attenuation coefficient from the transient response curve, and solve the Bellman equation through dynamic programming to obtain the local pressure adjustment amount. Construct a pressure step response function h(t) according to the local pressure adjustment amount, calculate the disturbance propagation characteristics from the pressure fluctuation amplitude, and discretize the Bernoulli equation through the central difference formula to obtain the boundary correction value. Construct a proportional-integral controller for the boundary correction value, set the proportional coefficient kp and the integral time Ti, and optimize the control parameters through the least squares method to obtain the initial control parameter group. Establish a rolling horizon prediction model according to the initial control parameter group, set the prediction horizon length to 5 control cycles, and update the controller parameters through online parameter identification to obtain the optimized control parameter group. Use the optimized control parameter group for pressure closed-loop control, calculate the control deviation from the flow field stability index, and eliminate the pressure fluctuation disturbance through feedback correction to obtain a stable pressure adjustment result.

[0020] Furthermore, according to the effect of precise pressure control and referring to the real-time data monitoring feedback, iteratively optimize the pressure control process until the local peeling risk is eliminated to an acceptable level, including: construct a monitoring and evaluation index R = Σwixi according to the pressure control effect and real-time feedback data, where the weight coefficient wi corresponds to three indexes: pressure control accuracy, stress distribution uniformity, and interface stability, and predict the shear stress distribution through a long short-term memory network to obtain the model update parameter. Correct the elastic modulus G0 and viscosity coefficient η0 in the constitutive equation for the model update parameter, and construct a root mean square error function from the data monitoring accuracy The optimization parameter group is obtained by optimizing the prediction accuracy through the gradient backpropagation algorithm. According to the optimization parameter group, the fluid shear stress distribution function τ(x,t) is reconstructed, the control error sequence is calculated from the control process accuracy, and the step size is iteratively calculated by the steepest descent method. The control parameters are adjusted. The pressure response characteristics are calculated for the control parameters, the overshoot and the adjustment time are extracted by using the step response curve, and the dynamic compensation parameters are obtained by updating the controller gain coefficient through online parameter identification. An adaptive controller is constructed according to the dynamic compensation parameters, the control action intensity is adjusted from the real-time monitoring data, and the pressure is accurately regulated through the proportional-integral-derivative control algorithm to obtain the optimized control result. The convergence condition is judged for the optimized control result. The prediction error threshold is set to 5%, and the critical value of the stripping risk index is 0.1. If the convergence condition is satisfied for 5 consecutive cycles, the iteration is completed; otherwise, return to continue the optimization.

[0021] The present invention provides an organoid automated culture monitoring system, which mainly includes:

[0022] An initial model establishment module, which is used to establish an initial model of the fluid shear stress distribution according to the adhesion strength distribution between the matrix gel and the surface of the culture dish measured by an atomic force microscope or a single-molecule force spectrometer, and obtain the initial shear stress data;

[0023] A fluid mechanics simulation module, which is used to perform fluid mechanics simulation on the initial shear stress data by using the fluid mechanics simulation method under non-Newtonian fluid, and calculate the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation;

[0024] A flow velocity adjustment module, which is used to adjust the relative influence of the inertial effect and the viscous effect on the flow by controlling the Reynolds number and the capillary number within a specific range, thereby adjusting the liquid flow velocity and outputting the flow velocity adjustment amount;

[0025] A pressure control parameter adjustment module, which is used to convert the output flow velocity adjustment amount into the change amount of the pressure control parameter by using the Bernoulli equation, and dynamically adjust the pressure control parameter of the processed liquid according to the established relationship between the flow velocity and the pressure to obtain the adjusted pressure control parameter;

[0026] A microscopic morphology feature acquisition module, which is used to acquire the microscopic morphology features of the medium-organoid interface under the adjusted pressure control parameter, and judge whether there is a risk of local matrix gel peeling by comparing and analyzing the image sequence;

[0027] A shear stress distribution correction module, which is used to correct the constitutive equation or boundary condition of the fluid shear stress distribution according to the rheological properties of the matrix gel if there is a local peeling risk, and recalculate the shear stress change trend in each region during the medium replacement process;

[0028] A pressure adjustment amount determination module, which is used to determine the pressure adjustment amount of a local area through an adaptive algorithm, taking the change trend of the shear stress in each area during the medium replacement process recalculated as the algorithm input, and applying the pressure adjustment amount of the local area as an additional boundary condition to the hydrodynamic model;

[0029] A pressure control optimization module, which is used to iteratively optimize the pressure control process according to the effect of precise pressure control and with reference to the real-time data monitoring feedback until the local peeling risk is eliminated to an acceptable level.

[0030] The technical solution provided by the embodiment of the present invention may include the following beneficial effects:

[0031] The present invention discloses an organoid automated culture monitoring method. The method first measures the adhesion strength between the matrix gel and the culture dish through an atomic force microscope to establish an initial shear stress model. Then, non-Newtonian hydrodynamic simulation is used to calculate the shear stress distribution during the medium replacement process. By adjusting the flow rate and pressure parameters, the shear stress is controlled within a safe range. High-resolution imaging technology is used to monitor the interface morphology in real time to determine whether there is a peeling risk. If there is a risk, the model is corrected and recalculated, and the local pressure adjustment amount is determined through an adaptive algorithm. Finally, iterative optimization is performed until the peeling risk is eliminated and the model converges. The present invention can effectively prevent the peeling of the matrix gel during the organoid culture process, ensure the culture quality and experimental repeatability, and is of great significance for biomedical research. Description of the Drawings

[0032] Figure 1 It is a flowchart of an organoid automated culture monitoring method of the present invention. Detailed Embodiments

[0033] Next, the technical solutions in the embodiments of the present invention will be clearly and detailedly described in conjunction with the drawings in the embodiments of the present invention. The described embodiments are only a part of the embodiments of the present invention.

[0034] Such as Figure 1 , a specific organoid automated culture monitoring method in this embodiment may specifically include:

[0035] S101. Measure the adhesion strength distribution between the matrix gel and the surface of the culture dish by using an atomic force microscope or a single-molecule force spectrometer, and establish an initial fluid shear stress model to obtain initial shear stress data.

[0036] S1011. Collect the interaction force data between the Matrigel and the surface of the culture dish using an atomic force microscope, perform noise reduction processing, and generate a continuous adhesion strength distribution curve. Scan the surface topography of the Matrigel using a single-molecule force spectrometer and generate a surface roughness distribution map after Gaussian filtering. Then, combine a neural network to construct an adhesion characteristic mapping relationship to extract the surface morphological characteristics of the Matrigel. The atomic force microscope uses a four-sided pyramidal probe tip with a resolution of 20 nanometers. The sampling frequency is set to 200 Hz, and each sampling point is measured 32 times. The noise signal above 100 Hz is filtered out through Fourier transform to generate a smooth adhesion strength curve. The scanning area of the single-molecule force spectrometer is set to 10 microns by 10 microns, with a step size of 50 nanometers. The topography data is processed using a Gaussian filter with a filter window of 5 pixels by 5 pixels and a standard deviation of 1.5 pixels. The root mean square value of the surface roughness is calculated to be approximately 85 nanometers. Subsequently, a radial basis function neural network is used, based on a multi-scale Gaussian kernel function, with 128 hidden layer nodes, a learning rate of 0.01, and trained for 5000 iterations. The mapping error is controlled within 3%, and thus the surface morphological characteristic data of the Matrigel is output.

[0037] S1012. After constructing hexagonal grid cells based on the surface morphological characteristic data of the Matrigel and setting the fluid boundary conditions, solve the Navier-Stokes equation using the finite difference method to obtain the flow field velocity distribution. Then, use the smoothed particle hydrodynamics method to calculate the stress tensor and extract the fluid shear stress data. Finally, generate the initial shear stress distribution data through normalization processing based on the maximum entropy principle. The characteristic scale of the Matrigel surface is approximately 500 nanometers. The side length of the grid cell is set to 50 nanometers, and the computational domain is divided into approximately 40,000 cells. The boundary layer region is encrypted to 8 layers and a no-slip boundary condition is applied. When solving, the viscosity of the culture medium is taken as 1.5 mPa·s, the density is 1.03 g / cm³, the inlet velocity is 100 μm / s, and the Reynolds number is approximately 0.5 to ensure a laminar flow state. The momentum equation and the continuity equation are discretized using a second-order central difference scheme. In the smoothed particle hydrodynamics calculation, the cubic spline function is selected as the kernel function, the smoothing length is 100 nanometers, and each particle considers the contributions of approximately 60 neighboring particles. The average angle between the main direction of the stress tensor and the flow direction is 15 degrees. The maximum entropy principle aims to maximize the Shannon entropy, constrains the probability density integral to be 1 and matches the experimental mean, determines the shear stress threshold interval to be 0.5 Pa to 2.0 Pa, and the normalized data approximately follows a lognormal distribution.

[0038] S101 ensures the reliability of the initial shear stress data through the combination of high-precision measurement and multi-scale modeling. It can be understood that the optimization of grid division and boundary conditions effectively improves the model's ability to characterize the adhesion characteristics of the Matrigel, avoiding shear stress prediction deviations caused by data noise or calculation errors.

[0039] S1013. In some alternative embodiments, if the atomic force microscope measurement is affected by probe wear, the adhesion strength data can be corrected by calibrating the probe shape parameters, or the step size can be adjusted to 30 nm during the force spectroscopy scan to improve the topography resolution, which can be flexibly set according to the experimental conditions.

[0040] S1014. In the flow field calculation, if the viscosity of the culture medium varies significantly with the shear rate, the power-law model can be introduced to correct the viscous term of the Navier-Stokes equation, so as to more accurately reflect the characteristics of non-Newtonian fluids. The specific parameters are determined by technicians according to the properties of the culture medium and are not overly limited. The initial shear stress data obtained through the above steps can fully reflect the mechanical properties of the interface between the Matrigel and the culture dish, helping to achieve stable monitoring and optimization adjustment during the organoid culture process.

[0041] S102. Perform a hydrodynamic simulation of non-Newtonian fluids for the initial shear stress data. Calculate the flow field characteristics of each region during the culture medium replacement process by solving the kinetic equation and obtain the shear stress change trend. At the same time, use grid optimization and numerical methods to improve the calculation accuracy to reflect the rheological behavior of the culture medium.

[0042] S1021. Establish a power-law fluid constitutive model based on the initial shear stress data and calculate the viscosity change characteristics of the culture medium, and analyze the fluid motion state through characteristic parameters to determine the influence of turbulence. The power-law model is defined as the shear stress equal to the consistency coefficient multiplied by the flow index power of the shear rate. The consistency coefficient is initially set to 2.1 mPa·s, and the flow index is experimentally fitted to 0.75. As a non-Newtonian fluid, the culture medium has the shear-thinning property. When the shear rate increases from 10 per second to 100 per second, the viscosity decreases from 1.8 mPa·s to 0.8 mPa·s, and this change is accurately characterized by the power-law relationship. To further analyze the flow field characteristics, the Reynolds number is calculated using the density of the culture medium of 1.03 g / cm³ and the characteristic length of the culture dish of 10 mm, and the value is 125, indicating that the fluid is in the laminar flow state and the turbulence effect is weak. The turbulent viscosity coefficient is calculated by the Spalart-Allmaras model, effectively capturing the non-linear change of the culture medium viscosity with the shear rate.

[0043] After generating a tetrahedral mesh for the geometric structure of the petri dish and optimizing the node distribution in the boundary layer region, the fluid pressure field and velocity field are calculated through the momentum conservation equation. Then, the stress tensor is constructed based on the characteristics of the velocity field, and the shear stress change data is extracted. The petri dish has a diameter of 35 mm and a height of 10 mm. The mesh is divided using tetrahedral elements, with a base side length of 1 mm. The total number of meshes is approximately 520,000. To improve the boundary layer accuracy, an octree algorithm is used for four-layer encryption, and the minimum mesh size reaches 0.125 mm. The no-slip condition is applied to the solid wall surface, the inlet flow rate is set to 2 ml per minute, and the outlet relative pressure is 0 Pa. The momentum conservation equation includes an inertial term, a pressure term, a viscous term, and a body force term caused by gravity. Among them, the inertial term is discretized using the second-order upwind scheme to enhance stability, the pressure term is processed by second-order central differences, the viscous term is dynamically adjusted according to the viscosity change, and the body force term is included with a gravitational acceleration of 9.8 m / s². The linear relationship between fluid density and pressure is defined by a compressibility coefficient of 4.8×10 -10 per Pa. The Galerkin finite element method is used to discretize the equation. Linear tetrahedral elements are selected. Each node contains three velocity components and a pressure value. The trial function space is constructed by piecewise linear basis functions. The sparse matrix is solved by the conjugate gradient method, and the convergence residual is less than 10 -6 . The calculation results show that the maximum velocity at the inlet reaches 15 mm / s, the pressure field distribution gradually decreases from the inlet to the outlet, and the velocity gradient is calculated by the central difference scheme with a step size equal to the mesh size. The stress tensor is generated and projected onto the principal axis direction to obtain the shear stress components.

[0044] S1023. In some alternative embodiments, if the viscosity of the culture medium changes significantly, a time-dependent term of the shear rate can be introduced through the power-law model to further correct the calculation formula of the viscous term to adapt to the dynamic flow conditions. The specific parameters can be adjusted by technicians according to experimental data and are not overly limited.

[0045] By performing time-domain sampling on the shear stress components, recording once every 0.1 s, and continuously collecting 100 cycles, a shear stress change curve is obtained. The results show that the peak shear stress at the initial stage of the culture medium replacement is 0.15 Pa, and then it shows an exponential decay trend, stabilizing at 0.08 Pa after about 8 s. This change trend reflects the dynamic evolution process of the fluid shear stress. Compared with traditional static simulations, by introducing the power-law fluid characteristics and mesh optimization, the accurate simulation of non-Newtonian fluid behavior is ensured, and the calculation deviation caused by neglecting viscosity changes is avoided.

[0046] S1024. During the hydrodynamic simulation, if it is necessary to further improve the velocity field resolution, the local grid density can be increased in the high-shear region, or the adaptive time step can be used to adjust the solution rhythm to balance the computational efficiency and accuracy. The specific implementation method can be flexibly set according to actual needs. The above steps comprehensively characterize the flow field behavior during the culture medium replacement by combining the non-Newtonian fluid characteristics and high-precision numerical methods. In particular, the capture of the shear stress change trend not only reveals the spatio-temporal distribution of the fluid action but also ensures the stability of the Matrigel structure during the organoid culture process.

[0047] S103. Based on the shear stress change trend in each region during the culture medium replacement, optimize the fluid flow characteristics by regulating the Reynolds number and capillary number within a specific range and output the flow velocity adjustment amount. At the same time, combine the quantitative analysis of viscous and inertial effects to achieve precise control to maintain the stability of the organoid structure.

[0048] S1031. Calculate the fluid viscosity coefficient according to the shear stress change trend and construct the Reynolds number and capillary number, and analyze the influence of viscous effects on the flow field through stress tensor to extract key control parameters. The fluid viscosity coefficient is determined according to the non-linear relationship between shear stress and shear rate. When the shear stress increases from 0.05 Pa to 0.15 Pa, the viscosity coefficient decreases from 1.8 mPa·s to 0.8 mPa·s, reflecting the shear-thinning characteristics of the culture medium.

[0049] The Reynolds number is calculated by the formula of fluid density multiplied by velocity multiplied by characteristic length divided by viscosity coefficient. Selecting the culture dish diameter of 35 mm as the characteristic length, with a density of 1.03 g / cm³ and a flow velocity of 10 mm / s, the Reynolds number is approximately 125. The capillary number is calculated by multiplying the viscosity coefficient by the flow velocity and dividing by the surface tension coefficient. Taking the surface tension of 72 mN / m and a flow velocity of 10 mm / s, the obtained value is 0.011. The target range is set as the Reynolds number from 100 to 200 and the capillary number less than 0.02. To quantify the viscous effect, a stress tensor is constructed using the viscosity coefficient, where the stress components are jointly determined by the velocity gradient and the viscosity coefficient. At the center region of the culture dish where the velocity gradient is the largest, the shear stress component reaches 0.12 Pa and the normal stress component is 0.08 Pa. The principal stress components are extracted through eigenvalue decomposition, and combined with the pressure gradient along the flow direction, from 200 Pa / m at the inlet to 50 Pa / m at the outlet, the contribution value of the viscous effect is calculated to be approximately 0.85, which can clearly reflect the dominant role of viscous force in the flow field and provide a theoretical basis for subsequent adjustment.

[0050] S1032. After calculating the contribution value of the inertial effect based on the fluid pressure gradient and velocity change and optimizing the flow characteristics by combining with the viscous effect, the inlet flow rate is dynamically adjusted through a proportional-integral controller to generate a flow velocity adjustment amount, ensuring that the flow field parameters meet the preset constraints. The inertial force is calculated by the product of the fluid density and the gradient of the velocity vector. When the flow velocity increases from 5 millimeters per second to 15 millimeters per second, the acceleration effect is significant, and the inertial force reaches 0.15 millinewtons. The change in fluid kinetic energy is determined by the relationship between the density and the square of the velocity, with an increase of approximately 0.12 microjoules. Combining with the inertial force, the contribution value of the inertial effect is obtained. The ratio of the viscous effect to the inertial effect is defined as a characteristic parameter, and the calculated value is 0.15, indicating that the viscous force is dominant. To achieve precise control, a proportional-integral controller is used to adjust the inlet flow rate. The proportional coefficient is set to 0.8, and the integral time constant is 2 seconds. The gain coefficient is optimized to 1.2 through Jacobi matrix iteration. The output signal of the controller drives the rotational speed of the peristaltic pump to increase from 60 revolutions per minute to 120 revolutions per minute, and the flow rate increases from 2 milliliters per minute to 4 milliliters per minute. Based on the channel cross-sectional area of 314 square millimeters, the flow velocity is adjusted from 6.4 millimeters per second to 12.8 millimeters per second. At this time, the Reynolds number rises to 160, and the capillary number is 0.018, both within the target range. The response time of the controller is approximately 0.5 seconds, the adjustment accuracy is better than 5%, and the steady-state error is controlled within 2%. This dynamic adjustment method effectively balances the relative influence of inertia and viscous effects, avoiding the abnormal increase in shear stress caused by excessive flow velocity.

[0051] If the Reynolds number or capillary number exceeds the set range, the flow rate adjustment value can be corrected by recalculating the controller output signal until the constraint conditions are met, ensuring that the flow velocity adjustment amount is consistent with the expected flow field characteristics.

[0052] S1033. In some alternative embodiments, if the surface tension of the culture medium fluctuates due to temperature changes, a temperature compensation model can be introduced to correct the capillary number calculation formula, or the real-time nature of the shear stress change trend can be improved by increasing the sampling frequency. The specific implementation can be flexibly set by technicians according to the experimental environment and will not be overly limited.

[0053] The flow velocity adjustment amount generated through the above steps can effectively adapt to the flow field changes during the replacement of the culture medium. Under the conditions of a temperature of 298 Kelvin and a humidity of 60%, 10 consecutive tests show that the adjustment process has good repeatability, with a standard deviation of less than 3%, indicating that this method has high stability and reliability in practical applications.

[0054] S104. Use Bernoulli's equation to convert the flow velocity adjustment amount into a change in pressure control parameters and dynamically adjust the pressure parameters of the liquid handling module. At the same time, ensure the control accuracy through friction loss analysis and filtering optimization to meet the stability requirements of organoid culture.

[0055] Calculate the total pressure change value according to the flow rate adjustment amount and the fluid density parameter, and construct the pressure loss distribution in combination with the pipeline resistance characteristics. Then, generate the adjusted pressure control parameter through sensor data filtering and neural network mapping. The Bernoulli equation is used to describe the energy conservation relationship between the flow rate and the pressure, where the pressure term, the potential energy term, and the dynamic pressure term together determine the total pressure change. Specifically, the height difference between the inlet and the outlet is set to 50 mm, and the cross-sectional area ratio is 2:1. When the flow rate increases from 10 mm / s to 15 mm / s, the potential energy change is about 0.49 Pa, the dynamic pressure change is about 0.12 Pa, and the total pressure change value is calculated to be 0.61 Pa. To further consider the influence of friction in the pipeline, the pipeline resistance coefficient formula is used to calculate the loss amount. The inner wall roughness of the pipeline is taken as 0.02 mm, and the resistance coefficient is 0.038 when the Reynolds number is 125. The pipeline length-diameter ratio is 20, and the friction loss amount is about 0.15 Pa. The pressure loss along the way shows a linear distribution characteristic, and the slope is about 15 Pa / m. After superimposing the local loss, the actual pressure adjustment amount is corrected to 0.82 Pa. This method can effectively quantify the influence of the flow rate change on the pressure distribution.

[0056] S1041. Construct a response curve for the data collected by the pressure sensor and optimize the signal quality through a Kalman filter, and establish a dynamic mapping relationship between the pressure and the flow rate using an adaptive neural network to achieve precise control. The range of the pressure sensor is 0 to 10 Pa, and the linearity is better than 0.1%. The response curve is fitted by the least squares method, the slope is about 1.02, the intercept is 0.05 Pa, and the correlation coefficient is as high as 0.999. Accordingly, the pressure control interval is divided into a low pressure of 0 to 2 Pa, a medium pressure of 2 to 5 Pa, and a high pressure of 5 to 10 Pa. The Kalman filter processes the pressure data through the state equation, sets the process noise covariance to 0.01, the observation noise covariance to 0.05, and the sampling period to 10 ms. After filtering, the standard deviation of the data drops from 0.08 Pa to 0.02 Pa, the signal-to-noise ratio is increased by about 8 dB, and the response time is less than 50 ms. The adaptive neural network adopts a three-layer structure. The input layer includes two variables, namely the pressure and the flow rate. The hidden layer is set with 8 neurons using the hyperbolic tangent activation function, and the output layer is a single pressure control parameter. 100 groups of data are used in the training process, the learning rate is 0.01, and the mean square error is less than 0.001 after 1000 iterations. If the pressure parameter exceeds the set accuracy of 0.1 Pa, the network automatically adjusts the weights and retrains for 5 to 10 rounds until the requirements are met.

[0057] S1042. In some alternative embodiments, if the inner wall roughness of the pipeline changes due to long-term use, the resistance coefficient can be updated by real-time measurement, or a dynamic noise covariance adjustment strategy can be introduced into the filter to adapt to environmental changes. The specific parameters are set by technicians according to the actual situation.

[0058] The pressure control parameters obtained through the above steps performed excellently in 10 consecutive medium replacements, with an average deviation of only 0.05 Pa and a maximum deviation of 0.08 Pa. The repeatability standard deviation was controlled within 3%. This high-precision pressure regulation effectively reduced the impact of fluid shear stress on the Matrigel, ensuring the structural integrity during organoid culture.

[0059] S105. Obtain the microscopic morphological characteristics of the medium-organoid interface under the adjusted pressure control parameters using confocal microscopy or optical coherence tomography technology, and analyze whether there is a risk of Matrigel peeling. At the same time, quantify the interface changes through image processing and deep learning methods to improve the monitoring accuracy.

[0060] S1051. Collect the original image sequence of the medium-organoid interface using confocal microscopy and perform filtering processing, and extract multi-scale features through wavelet transform to construct the interface topography distribution. The confocal microscopy uses a 488-nm laser to excite the fluorescently labeled medium, is equipped with a 60x oil immersion objective lens with a numerical aperture of 1.4, a lateral resolution of 200 nm, an axial resolution of 500 nm, the acquired image size is 1024×1024 pixels, the sampling interval is 5 seconds, the pixel dwell time is 2 microseconds, and the frame averaging is 4 times, with a signal-to-noise ratio exceeding 35 dB. The original image is processed by a bilateral filter with a spatial standard deviation of 5 and a range standard deviation of 0.1. After removing the noise, the peak signal-to-noise ratio is increased to 42 dB, and the edge sharpness retention rate exceeds 95%. Subsequently, a three-layer discrete wavelet transform is performed on the filtered image, and the db4 wavelet basis function is selected to decompose the high-frequency coefficients in the horizontal, vertical, and diagonal directions. The first-layer high-frequency coefficients highlight the microscopic interface details, and the regions with an amplitude greater than 0.5 correspond to the significantly undulating parts. The second-layer and third-layer coefficients reflect the mesoscopic and macroscopic roughness, and the standard deviation quantifies the interface uniformity, generating the interface topography feature map. This multi-scale analysis can effectively capture the characteristic changes of the interface from small fluctuations to large-scale peeling.

[0061] After constructing a deep neural network based on the interface morphology feature map and extracting interface features, a stripping risk distribution map is generated through morphological operations and region labeling to judge the stability of the basement membrane glue. The deep neural network consists of 4 convolutional layers with a 3×3 convolutional kernel, and the number of channels is 32, 64, 128, and 256 in sequence. Each layer is followed by a max pooling layer and a batch normalization layer, with a dropout rate of 0.5. The input is a 256×256 pixel image patch. The network extracts significant features through max pooling, and batch normalization accelerates convergence. 1000 images are used for training, among which 200 contain stripping labels. The cross-entropy loss function guides the optimization, with a learning rate of 0.001. The accuracy of the validation set reaches 92% after 50 rounds of training. The interface features are subjected to opening and closing operations using a circular structuring element with a radius of 3 pixels. The opening operation eliminates noise smaller than 6 pixels, and the closing operation fills holes smaller than 6 pixels. The morphological results are labeled by 8-neighborhood connected regions, and regions with an area smaller than 50 pixels are filtered to generate a stripping risk distribution map. The stripping area often appears as strips or patches with obvious serrated edges. If the stripping area exceeds 5% of the total area for 3 consecutive frames, a warning signal is triggered. Experiments show that the stripping starts from the edge and can expand to 20% of the field of view within 15 seconds. The 5% threshold warning can be used for early intervention to ensure the stability of the culture.

[0062] S1052. Use optical coherence tomography to verify the three-dimensional structure of the interface and assist in the assessment of the stripping risk. Reconstruct tomographic images through a phase-sensitive algorithm to improve the monitoring reliability. The optical coherence tomograph uses a 1310-nanometer swept-source laser, with a scanning range of 100 nanometers, an axial resolution of 5 micrometers, and a scanning depth of 500 micrometers. The phase-sensitive algorithm reconstructs the three-dimensional morphology of the interface through interference signals. The stripping area appears as a dark band, with a width increasing from 20 micrometers to 100 micrometers and a depth of up to 150 micrometers. This technology complements the planar imaging of the confocal microscope, provides structural information in the vertical direction, and ensures the multi-dimensional accuracy of the stripping determination.

[0063] S1053. In some alternative embodiments, if the image noise increases due to environmental light interference, the range standard deviation of the bilateral filter can be adjusted to 0.15, or residual connections can be introduced into the neural network to improve the robustness of feature extraction. The specific parameters can be optimized according to the experimental conditions.

[0064] Through the above method combining high-resolution imaging and intelligent analysis, the impact of culture medium replacement on the organoid interface can be monitored in real time. The generation of the stripping risk distribution map not only intuitively reflects the local stability of the basement membrane glue but also provides a feedback basis for the further optimization of the pressure parameters, effectively reducing the possibility of irreversible damage during the culture process.

[0065] S106. If a local peeling risk is detected, the fluid shear stress distribution model is corrected according to the rheological properties of the matrix gel, and the shear stress change trend during the medium replacement process is recalculated. The model accuracy is optimized through stress transfer analysis and regional division to improve the stability of organoid culture.

[0066] Calculate the local strain rate for the matrix gel peeling area and extract rheological parameters to correct the constitutive equation, and generate the corrected shear stress change trend through stress distribution analysis to reduce the peeling effect. The strain rate in the peeling area increases from the initial 0.1 per second to 0.8 per second. The elastic modulus of the matrix gel is measured to be about 2.5 kPa and the viscosity coefficient is 1.8 Pa·s through the mechanical response curve. The generalized Maxwell equation is used to describe the viscoelastic behavior, where the stress is equal to the elastic modulus multiplied by the instantaneous strain plus the viscosity coefficient multiplied by the time derivative of the strain rate. At a shear stress of 0.5 Pa, the strain increases by 45% within 60 seconds, showing significant creep characteristics. After correcting the constitutive equation, calculate the local stress transfer coefficient, which drops to 0.65 at the peeling edge and remains 0.92 in the intact area, reflecting the spatial difference in stress transfer efficiency. Construct the stress-strain relationship based on the material strength index, solve the stress equilibrium equation through tetrahedral element discretization, with an element size of 50 μm and a total number of meshes of about 86,000. The boundary stress is redistributed due to peeling and decreases by 25%, generating the stress distribution boundary conditions. This method provides a mechanical basis for model correction by quantifying rheological properties and stress transfer, ensuring that the shear stress calculation is close to the actual state.

[0067] S1061. Construct a three-dimensional stress tensor according to the corrected constitutive equation and boundary conditions and solve the principal stress components, and use a deep neural network to predict the stress transfer path to optimize the stress distribution analysis. The stress tensor solves the principal stress components through the principal value equation. The angle between the maximum principal stress direction and the flow direction is about 35 degrees, and the principal stress ratio is 2.8, showing anisotropic characteristics. The deep neural network consists of 5 convolutional layers and 3 fully connected layers, inputs 32×32×32 stress field data, is trained with 1000 sets of samples, with a learning rate of 0.001, and the prediction accuracy reaches 91% after optimization. The stress transfer path prediction results show that the stress gradually decays from the high-stress area to the low-stress area, and the path is significantly affected by the peeling area. Through statistical analysis, the low-stress area is divided into 0 to 0.3 Pa, the medium-stress area is 0.3 to 0.8 Pa, and the high-stress area is greater than 0.8 Pa. The strain rate threshold is set to 0.2 per second, and the area exceeding the threshold accounts for 12%. Continuously monitor the stress data for 100 seconds, with a sampling frequency of 10 Hz. The stress time series shows periodic fluctuations. Use the 5-point central difference to calculate the stress gradient, with a time step of 0.1 s and a spatial step of 20 μm, to generate the corrected shear stress change trend. This analysis method not only reveals the dynamic evolution of the stress field.

[0068] S1062. In some optional embodiments, if the rheological properties of the Matrigel change significantly due to temperature variations, a temperature-dependent term can be introduced to correct the viscosity coefficient, or the real-time nature of the stress change trend can be enhanced by increasing the sampling frequency to 20 Hz. The specific adjustments are determined by the experimental conditions.

[0069] The corrected shear stress change trend is verified through finite element simulation. The node displacements and stress distributions are extracted, and the relative error compared with the experimental data is reduced to 8%, with the displacement deviation less than 5 microns. The shear stress in the high-stress area is reduced by 35%, the maximum stress gradient decreases from 0.15 Pa per micron to 0.08 Pa per micron, the stress distribution becomes more uniform, the fluctuation amplitude is reduced by 52%, the amplitude of the periodic perturbation is reduced by 65%, and the maximum stress in the peeling area is reduced by 45%. The stress concentration coefficient decreases from 2.8 to 1.6. This optimization effectively alleviates the impact of local peeling on the organoid structure.

[0070] S107. Using an adaptive algorithm, the locally adjusted pressure is determined with the recalculated shear stress change trend as the input and incorporated into the hydrodynamic model. Meanwhile, the global flow field stability is optimized to eliminate the pressure perturbations caused by local interventions, ensuring the structural integrity of the Matrigel during the organoid culture process.

[0071] Based on the shear stress change trend, a local peeling index function is constructed and the pressure adjustment target value is extracted. The local intervention parameters are calculated through deep learning and dynamic response analysis to achieve precise control. The local peeling index function is defined as the sum of the squares of the difference between the real-time shear stress and the critical shear stress divided by the number of sampling points. When the real-time shear stress reaches 0.15 Pa, exceeding the critical value of 0.12 Pa, the index value is approximately 0.023, indicating a peeling risk. A deep Q-learning network is used to establish the state-action value function. The network consists of 4 fully connected layers with the number of neurons being 128, 256, 128, and 64 in sequence. The hyperbolic tangent activation function is used, with a learning rate of 0.01 and a discount factor of 0.95. After 1000 rounds of iterative training, the pressure adjustment parameters are generated. This method optimizes the adjustment strategy through reinforcement learning and can dynamically adjust the pressure target value according to the shear stress. The local intervention intensity coefficient is calculated by the ratio of the pressure adjustment amount to the initial pressure. For example, when the adjustment amount is 0.2 Pa, the coefficient is 0.25, which is within a reasonable range. The oscillation frequency of 2.5 Hz and the attenuation coefficient of 0.8 are extracted from the transient response curve, and the Bellman equation is solved using the dynamic programming method to obtain a local pressure adjustment amount of approximately 0.18 Pa, ensuring that the intervention is both effective and not excessive.

[0072] S1071. Construct a pressure step response function for the local pressure adjustment amount and calculate the boundary correction amount, and optimize the control parameters through a rolling horizon prediction model to achieve closed-loop pressure control. The pressure step response function describes the dynamic response of the flow field to pressure changes, showing the characteristics of a second-order system, with an overshoot of about 15%, a regulation time of 0.8 seconds, and the fluctuation amplitude decaying from 0.1 Pa to 0.02 Pa. The Bernoulli equation is discretized using the central difference scheme with a spatial step size of 0.1 mm, and the calculated boundary correction amount is about 0.15 Pa, which is used to compensate for the disturbance of the local intervention on the flow field. The rolling horizon prediction model uses 5 control cycles, that is, 2.5 seconds as the prediction length, and updates the control parameters through online parameter identification. The recursive least squares method is used in the identification process with a forgetting factor of 0.95, the parameter convergence time is less than 0.3 seconds, and the optimized proportional coefficient is about 1.2, and the integral time constant is 0.5 seconds. Under closed-loop control, the overshoot is reduced to 8%, the regulation time is shortened to 0.5 seconds, the pressure fluctuation amplitude is reduced by 75%, the flow field stability index is increased from 0.85 to 0.95, and the root mean square value of the control deviation is reduced to 0.015 Pa. This optimization method effectively balances local regulation and global stability through the combination of real-time feedback and prediction.

[0073] S1072. In some alternative embodiments, if the flow field parameters fluctuate greatly, the number of training rounds of the Q-learning network can be increased to 1500 rounds, or the rolling horizon prediction length can be adjusted to 7 cycles to improve the control robustness, and the specific settings can be flexibly adjusted according to experimental requirements. The optimized pressure adjustment amount is applied to the hydrodynamic model as an additional boundary condition, significantly suppressing the risk of local peeling. Long-term operation tests show that the controller has strong adaptability to flow field disturbances, the pressure distribution uniformity is improved, and the stress concentration phenomenon in the local area is reduced, ensuring the continuous stability of the organoid environment during the culture medium replacement process.

[0074] S108. Update the fluid shear stress distribution model according to the pressure control effect and real-time monitoring data and iteratively optimize the pressure control process, and gradually eliminate the local peeling risk by constructing an evaluation index and a dynamic compensation mechanism until the model prediction error and the peeling risk index meet the preset convergence conditions to ensure the long-term stability of organoid culture.

[0075] Build monitoring and evaluation indicators based on pressure control accuracy, stress distribution uniformity, and interface stability to predict the shear stress distribution, and reconstruct the fluid shear stress model through parameter correction and error optimization to improve the control accuracy. The monitoring and evaluation indicators are defined as the weighted sum of three sub-indicators. The weight of pressure control accuracy is set to 0.4, reflecting the tracking ability of the controller to the target value. The weight of stress distribution uniformity is 0.35, measuring the spatial consistency of the flow field. The weight of interface stability is 0.25, evaluating the integrity of the matrix gel structure. Use a long short-term memory network to predict the shear stress distribution. The network contains 3 hidden layers, with 128 neurons in each layer. It is trained with 1000 sets of historical data, and the prediction accuracy reaches 93%. Modify the constitutive equation according to the prediction results. The elastic modulus is adjusted from 2.5 kPa to 2.2 kPa, and the viscosity coefficient is reduced from 1.8 Pa·s to 1.5 Pa·s to more accurately reflect the rheological properties of the matrix gel. The reconstructed shear stress distribution function shows that the local maximum stress is 0.18 Pa and the minimum value is 0.05 Pa, reflecting the spatial non-uniformity. The root mean square error function calculates the deviation between the predicted and measured values. The initial error of 12% is reduced to 4.8% after optimization by gradient backpropagation, and the correlation coefficient is increased to 0.96. This method combines data-driven and mechanical modeling to ensure that the model gradually approaches the actual flow field state.

[0076] S1081. Calculate the control error for the reconstructed shear stress distribution and optimize the control parameters, and achieve dynamic compensation through pressure response analysis and online identification to improve the pressure control process. The control error sequence is generated based on the shear stress distribution function, with a standard deviation of approximately 0.025 Pa. The steepest descent method is used to calculate the iteration step size, which is determined by the ratio of the error value to the square of the gradient norm and is dynamically adjusted from 0.5 to 0.2 to accelerate the convergence speed. The overshoot and adjustment time are extracted from the pressure step response curve. After optimization, the overshoot is reduced from 15% to 6%, and the adjustment time is shortened from 0.8 s to 0.4 s, indicating a smoother control response. Online parameter identification uses the recursive least squares method with a forgetting factor of 0.95. The controller gain coefficient is adjusted from 1.2 to 0.85, and the integral time is increased from 0.5 s to 0.8 s, improving the adaptability to dynamic disturbances. The adaptive controller adjusts the control intensity according to real-time data. When the peeling risk index exceeds 0.08, the proportional coefficient increases by 20%, the integral action is enhanced by 15%, and the derivative action is weakened by 10% to achieve hierarchical control. The amplitude of pressure fluctuation is reduced by 65%, the control accuracy reaches 0.01 Pa, the prediction error is stable below 4.2% within 5 consecutive cycles, and the maximum value of the peeling risk index is 0.085, meeting the convergence conditions of an error threshold of 5% and a risk critical value of 0.1. This dynamic optimization process effectively suppresses local peeling and maintains the global flow field balance.

[0077] S1082. In some alternative embodiments, if the monitoring data has significant noise, the training data of the long short-term memory network can be increased to 1500 groups, or the forgetting factor can be adjusted to 0.90 to enhance the robustness of parameter identification. The specific parameters can be flexibly set according to actual requirements.

[0078] The optimized pressure control process demonstrates excellent performance during long-term operation, with a 35% reduction in energy consumption, a 45% reduction in response time, and significantly enhanced adaptability of the controller to process fluctuations and external disturbances. For different peeling risk scenarios, the controller automatically switches strategies, using gentle adjustment in the case of weak risks and quickly increasing the action intensity in the case of significant risks to ensure the best balance between control accuracy and system stability.

[0079] The present invention provides an organoid automated culture monitoring system, mainly including:

[0080] An initial model establishment module, used to establish an initial model of the fluid shear stress distribution based on the adhesion strength distribution between the basement membrane matrix and the culture dish surface measured by an atomic force microscope or a single-molecule force spectrometer, and obtain initial shear stress data;

[0081] A fluid mechanics simulation module, used to perform fluid mechanics simulation on the initial shear stress data by using a fluid mechanics simulation method under non-Newtonian fluids, and calculate the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation;

[0082] A flow rate adjustment module, used to adjust the relative influence of inertial effects and viscous effects on the flow by controlling the Reynolds number and the capillary number within a specific range, and then adjust the liquid flow rate to output a flow rate adjustment amount;

[0083] A pressure control parameter adjustment module, used to use the Bernoulli equation to convert the output flow rate adjustment amount into a change amount of the pressure control parameter, and dynamically adjust the pressure control parameter for processing the liquid according to the established relationship between the flow rate and the pressure to obtain an adjusted pressure control parameter;

[0084] A microscopic morphology feature acquisition module, used to acquire the microscopic morphology features of the medium-organoid interface under the adjusted pressure control parameters, and judge whether there is a risk of local basement membrane matrix peeling by comparing and analyzing the image sequence;

[0085] A shear stress distribution correction module, used to, if there is a local peeling risk, correct the constitutive equation or boundary condition of the fluid shear stress distribution according to the rheological properties of the basement membrane matrix, and recalculate the shear stress change trend in each region during the medium replacement process;

[0086] A pressure adjustment amount determination module, which is used to determine the pressure adjustment amount of a local area by means of an adaptive algorithm with the change trend of the shear stress in each area during the medium replacement recalculated as the algorithm input, and apply the pressure adjustment amount of the local area as an additional boundary condition to the hydrodynamic model;

[0087] A pressure control optimization module, which is used to iteratively optimize the pressure control process according to the effect of precise pressure control and with reference to the real-time data monitoring feedback until the local peeling risk is eliminated to an acceptable level.

[0088] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the principle of the present invention, several improvements and supplements can also be made, and these improvements and supplements should also be regarded as the protection scope of the present invention.

Claims

1. An organoid automated culture monitoring method, characterized in that, The method includes: Based on the adhesion strength distribution between the Matrigel and the surface of the culture dish measured by an atomic force microscope or a single molecule force spectrometer, establish an initial model of the fluid shear stress distribution and obtain initial shear stress data; For the initial shear stress data, use the hydrodynamic simulation method under non-Newtonian fluid to conduct hydrodynamic simulation, and calculate the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation; By controlling the Reynolds number and the capillary number within a specific range, adjust the relative influence of the inertial effect and the viscous effect on the flow, and then adjust the liquid flow velocity to output the flow velocity adjustment amount; Utilize the Bernoulli equation to convert the output flow velocity adjustment amount into the change amount of the pressure control parameter, and dynamically adjust the pressure control parameter of the processed liquid according to the established relationship between the flow velocity and the pressure to obtain the adjusted pressure control parameter; Obtain the microscopic morphological characteristics of the interface between the culture medium and the organoids under the adjusted pressure control parameter, and judge whether there is a risk of local Matrigel peeling by comparing and analyzing the image sequence; If there is a local peeling risk, then according to the rheological properties of the Matrigel, correct the constitutive equation or boundary conditions of the fluid shear stress distribution, and recalculate the shear stress change trend in each region during the culture medium replacement process; Through the adaptive algorithm, take the shear stress change trend in each region during the recalculated culture medium replacement process as the input of the algorithm, determine the pressure adjustment amount in the local region, and apply the pressure adjustment amount in the local region as an additional boundary condition to the hydrodynamic model.

2. The method according to claim 1, characterized in that, The step of based on the adhesion strength distribution between the Matrigel and the surface of the culture dish measured by an atomic force microscope or a single molecule force spectrometer, establishing an initial model of the fluid shear stress distribution and obtaining initial shear stress data includes: Receive the data of the force between the Matrigel and the surface of the culture dish collected by the atomic force microscope, and obtain the continuous adhesion strength distribution curve after filtering out the sampling noise interference through Fourier transform; According to the surface topography data of the Matrigel measured by the force spectrometer, perform noise reduction processing on the surface topography data using a Gaussian filter to obtain the surface roughness distribution map; For the continuous adhesion strength distribution curve, use a radial basis function neural network to construct the mapping relationship between the adhesion site density and the surface roughness distribution map to obtain the surface morphological characteristic data of the Matrigel; Establish a hexagonal grid unit according to the surface morphological characteristic data of the Matrigel, set the grid unit size, and solve the Navier-Stokes equation using the finite difference method to obtain the initial shear stress data.

3. The method according to claim 1, wherein The step of for the initial shear stress data, using the hydrodynamic simulation method under non-Newtonian fluid to conduct hydrodynamic simulation, and calculating the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation includes: Construct a power-law fluid constitutive equation using the shear stress data and the shear rate data, and obtain the culture medium viscosity change curve according to the power-law fluid constitutive equation; Establish a tetrahedral grid unit for the geometric structure of the culture dish, determine the side length of the grid unit according to the characteristic size of the culture dish, and perform encryption processing on the hydrodynamic boundary layer region through an octree grid refinement algorithm to obtain the grid node distribution; Construct a momentum conservation equation according to the viscosity change curve of the culture medium. The momentum conservation equation includes an inertial term, a pressure term, a viscous term, and a body force term, and determine the fluid viscosity coefficient from the fluid motion state value.

4. The method according to claim 1, characterized in that By controlling the Reynolds number and the capillary number within a specific range, adjusting the relative influence of the inertial effect and the viscous effect on the flow, and then adjusting the liquid flow velocity, the output flow velocity adjustment amount includes: Obtain the fluid viscosity coefficient, and the fluid viscosity coefficient value is used to construct the Reynolds number and the capillary number; Construct a stress tensor for the fluid viscosity coefficient value, obtain the principal stress components from the stress tensor, and obtain the viscous effect contribution value through the principal stress components and the pressure gradient; Calculate the inertial force according to the pressure gradient change, obtain the kinetic energy change from the fluid density and velocity curve, and obtain the inertial effect contribution value through the kinetic energy change and the inertial force; Use a proportional-integral controller to adjust the inlet flow rate. If the Reynolds number exceeds the range constraint, recalculate the controller output signal until the constraint condition is met to obtain the flow velocity adjustment amount.

5. The method according to claim 1, wherein Using the Bernoulli equation, convert the output flow velocity adjustment amount into a change amount of the pressure control parameter, and dynamically adjust the pressure control parameter of the processed liquid according to the established relationship between the flow velocity and the pressure to obtain the adjusted pressure control parameter, including: Substitute the flow velocity and the fluid density value into the Bernoulli equation, and obtain the total pressure change value through the liquid level height difference and the cross-sectional area ratio; For the total pressure change value, calculate the friction loss amount using the pipe resistance coefficient, and construct a pressure loss distribution curve from the friction loss amount; Obtain the pressure sensor data according to the pressure loss distribution curve, establish a Kalman filter state equation through the pressure sensor data, and obtain the filtered pressure data from the state equation; Use an adaptive neural network to establish the mapping relationship between the filtered pressure data and the flow velocity. If the pressure control parameter exceeds the set range, retrain the adaptive neural network until the pressure control parameter that meets the control accuracy requirements is obtained; it also includes: obtaining the pressure sensor data, establishing a state equation using a Kalman filter to obtain the filtered pressure data, establishing the mapping relationship between the filtered pressure data and the flow velocity through an adaptive neural network, judging whether the pressure control parameter exceeds the set range, and if it exceeds, return to the neural network for retraining until the pressure control parameter that meets the control accuracy requirements is obtained.

6. The method according to claim 1, characterized in that, Obtain the microscopic morphological characteristics of the interface between the culture medium and the organoids under the adjusted pressure control parameter, and judge whether there is a risk of local matrix gel peeling by comparing and analyzing the image sequence, including: Obtain the original image sequence of the interface between the culture medium and the organoids, which is collected by a confocal microscope, and process the original image sequence using a bilateral filter to obtain the filtered image; Perform three-layer discrete wavelet transform decomposition on the filtered image, and construct an interface morphological feature map according to the high-frequency coefficients of the discrete wavelet transform; Input the interface morphological feature map into a deep neural network, and obtain the interface features through the max pooling and batch normalization layers of the deep neural network; Morphological operations are performed on the interface features using a circular structural element, and the peeling risk distribution map is obtained by processing the results of the morphological operations through a region labeling algorithm.

7. The method according to claim 1, characterized in that If there is local peeling risk, then according to the rheological properties of the matrix gel, the constitutive equation or boundary conditions of the fluid shear stress distribution are corrected, and the shear stress change trend in each region during the medium replacement process is recalculated, including: Calculating the local strain rate according to the distribution of the matrix gel peeling region, and obtaining the rheological characteristic parameters of the matrix gel from the mechanical response curve using the strain rate to obtain the corrected constitutive equation; Calculating the local stress transfer coefficient through the corrected constitutive equation, and obtaining the stress-strain relationship from the material strength index using the stress transfer coefficient to obtain the stress distribution boundary conditions; Constructing a three-dimensional stress tensor according to the stress distribution boundary conditions, and solving the principal stress components from the principal value equation using the stress tensor to obtain the stress transfer path prediction result; Dividing the low stress area, medium stress area, and high stress area according to the stress transfer path prediction result, and obtaining the corrected shear stress change trend by obtaining the stress change curve from the continuously sampled stress data sequence using the stress interval.

8. The method according to claim 1, characterized in that, Through an adaptive algorithm, taking the shear stress change trend in each region during the recalculated medium replacement process as the algorithm input, determining the pressure adjustment amount of the local region, and applying the pressure adjustment amount of the local region as an additional boundary condition to the hydrodynamic model, including: Constructing a local peeling index function according to the shear stress change trend, and obtaining the pressure adjustment target value through the sum of squares of the difference between the real-time shear stress and the critical shear stress in the local peeling index function; Calculating the local intervention intensity coefficient for the pressure adjustment target value, where the local intervention intensity coefficient is determined by the ratio of the pressure adjustment amount to the initial pressure, and obtaining the oscillation frequency and attenuation coefficient from the transient response curve using the local intervention intensity coefficient; Constructing a pressure step response function according to the oscillation frequency and attenuation coefficient, and obtaining the boundary correction amount value through the pressure fluctuation amplitude in the pressure step response function; Establishing a rolling horizon prediction model for the boundary correction amount value, where the rolling horizon prediction model obtains an optimized control parameter group through online parameter identification, and uses the optimized control parameter group for pressure control.

9. The method according to claim 1, characterized in that, It also includes: According to the effect of precise pressure control and referring to the real-time data monitoring feedback, iteratively optimizing the pressure control process until the local peeling risk is eliminated to an acceptable level, specifically including: Constructing a monitoring and evaluation index according to the pressure control accuracy, stress distribution uniformity, and interface stability to obtain the evaluation result; Predicting the shear stress distribution of the evaluation result using a long short-term memory network, and obtaining the model update parameters by correcting the elastic modulus and viscosity coefficient in the constitutive equation through the shear stress distribution; Reconstructing the fluid shear stress distribution function according to the model update parameters, calculating the control error sequence from the shear stress distribution function, and obtaining the control parameters by the steepest descent method. The overshoot and settling time are extracted from the step response curve for the control parameters, and the dynamic compensation parameters are obtained by updating the controller gain coefficient through online parameter identification. If the prediction error is less than the preset threshold and the peeling risk index is less than the critical value, the iterative optimization is completed.

10. An organoid automated culture monitoring system, characterized in that, The system includes: An initial model establishment module, configured to establish an initial model of the fluid shear stress distribution according to the adhesion strength distribution of the matrix glue and the culture dish surface measured by an atomic force microscope or a single-molecule force spectrometer, and obtain initial shear stress data; A fluid mechanics simulation module, configured to perform fluid mechanics simulation on the initial shear stress data by using a fluid mechanics simulation method under non-Newtonian fluid, and calculate the fluid viscosity coefficient by solving the Navier-Stokes equation and the continuity equation; A flow rate adjustment module, configured to adjust the relative influence of the inertial effect and the viscous effect on the flow by controlling the Reynolds number and the capillary number within a specific range, thereby adjusting the liquid flow rate and outputting a flow rate adjustment amount; A pressure control parameter adjustment module, configured to convert the output flow rate adjustment amount into a change amount of the pressure control parameter by using the Bernoulli equation, and dynamically adjust the pressure control parameter of the processed liquid according to the established relationship between the flow rate and the pressure to obtain an adjusted pressure control parameter; A microscopic morphology feature acquisition module, configured to acquire the microscopic morphology features of the medium and organoid interface under the adjusted pressure control parameter, and determine whether there is a risk of local matrix glue peeling by comparing and analyzing the image sequence; A shear stress distribution correction module, configured to, if there is a local peeling risk, correct the constitutive equation or boundary condition of the fluid shear stress distribution according to the rheological properties of the matrix glue, and recalculate the shear stress change trend of each region during the medium replacement process; A pressure adjustment amount determination module, configured to determine the pressure adjustment amount of the local region by using an adaptive algorithm with the shear stress change trend of each region during the recalculated medium replacement process as the algorithm input, and apply the pressure adjustment amount of the local region as an additional boundary condition to the fluid mechanics model; A pressure control optimization module, configured to iteratively optimize the pressure control process according to the effect of precise pressure control and with reference to the real-time data monitoring feedback until the local peeling risk is eliminated to an acceptable level.

Images