Batch fitting method and system for parallel kinetic parameters of high-density SPR microarray

Abstract

This invention provides a parallel batch fitting method and system for dynamic parameters of high-density SPR microarrays, relating to the field of surface plasmon resonance technology. The method includes: iteratively solving a weighted local dynamic fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; acquiring the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip; constructing a spatial neighborhood regularization constraint term based on the physical spatial coordinate information; and performing spatial smoothing correction on the initial binding rate constant and initial dissociation rate constant using the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site. This invention improves the calculation accuracy of the binding rate constant and dissociation rate constant of each microarray site in a high-density SPR microarray chip.

Description

Technical Field

[0001] This invention relates to the field of surface plasmon resonance technology, and in particular to a method and system for batch fitting of parallel dynamic parameters for high-density SPR microarrays. Background Technology

[0002] In high-density SPR microarray chip applications, microarray sites are densely arranged in a matrix on the chip surface. During detection, a fluid sample continuously flows across the chip surface, binding and dissociating with probe molecules at the microarray sites. This reaction process excites surface plasmon resonance signals, generating real-time sensing data that varies continuously over time. The continuous temporal fluctuations reflect the dynamics of molecular interactions. Simultaneously, the continuous flow of the fluid sample across the chip surface results in a spatially continuous concentration distribution of the fluid environment around adjacent microarray sites. By analyzing the real-time sensing data, binding rate constants and dissociation rate constants, reflecting the rates of molecular binding and dissociation, can be calculated. For a massive number of microarray sites, computational load weights and data complexity feature vectors need to be introduced to assess the data processing pressure at each site.

[0003] However, existing technologies have the following technical shortcomings when processing the real-time sensing data of the aforementioned high-density SPR microarray chips: First, existing technologies typically employ a globally uniform weight allocation method when constructing the dynamic fitting objective function. However, due to significant instantaneous fluctuations in the surface plasmon resonance signal during the binding and dissociation stages, and the influence of local fluid disturbances on each microarray site leading to substantial differences in signal change rates, existing technologies fail to measure the instantaneous change rate of the surface plasmon resonance signal point-by-point. This results in the inability to accurately identify local noise interference regions and the difficulty in dynamically allocating weights to each data point. Consequently, the fitting process is susceptible to abnormal fluctuations, reducing the accuracy of the calculation of the binding rate constant and dissociation rate constant.

[0004] Second, existing technologies typically treat each microarray site as an independent entity for isolated calculations when performing batch calculations of kinetic parameters. Since the continuous flow of fluid samples on the chip surface exhibits a continuous spatial concentration distribution, and adjacent microarray sites are in highly similar physical environments, their actual binding and dissociation rate constants should exhibit a smooth transition in spatial distribution. Existing technologies fail to perform global cumulative calculations and spatial smoothing corrections for kinetic parameters within the spatial neighborhood. This leads to unreasonable spatial jumps in the calculated kinetic parameters of adjacent microarray sites, reducing the consistency and reliability of the overall parameter distribution of high-density SPR microarray chips. Summary of the Invention

[0005] This invention provides a method and system for batch fitting of parallel dynamic parameters for high-density SPR microarrays, which improves the calculation accuracy of binding rate constants and dissociation rate constants of each microarray site in a high-density SPR microarray chip.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: Firstly, a parallel batch fitting method for dynamic parameters of high-density SPR microarrays is provided, the method comprising: Step 1: Obtain real-time sensing image data of each microarray site on the high-density SPR microarray chip, extract multi-dimensional time-series features from the real-time sensing image data to obtain the data complexity feature vector; calculate the computational load weight of each microarray site using the data complexity feature vector. Step 2: Based on the computational load weight, a dynamic task scheduling strategy is used to cluster and divide the real-time sensor map data to obtain multiple parallel fitting sub-task clusters with balanced computational load. Step 3: Extract the local residual distribution features and spatial neighborhood heterogeneity features of the sensing map data of each microarray site in the parallel fitting subtask cluster, and construct a local heterogeneity perception model. Step 4: Calculate the adaptive weight coefficients of each microarray site in the dynamic fitting process using the local heterogeneity sensing model; based on the adaptive weight coefficients, construct the weighted local dynamic fitting objective function corresponding to each microarray site. Step 5: Iteratively solve the weighted local dynamics fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; obtain the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip. Step 6: Construct a spatial neighborhood regularization constraint term based on physical spatial coordinate information; perform spatial smoothing correction on the initial binding rate constant and the initial dissociation rate constant through the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site.

[0007] Secondly, a parallel dynamic parameter batch fitting system for high-density SPR microarrays includes: The acquisition module is used to acquire real-time sensing image data of each microarray site on the high-density SPR microarray chip, extract multi-dimensional time-series features from the real-time sensing image data to obtain a data complexity feature vector, and calculate the computational load weight of each microarray site through the data complexity feature vector. The partitioning module is used to cluster real-time sensor image data based on computational load weights and a dynamic task scheduling strategy to obtain multiple parallel fitting subtask clusters with balanced computational load. The module is used to extract the local residual distribution characteristics and spatial neighborhood heterogeneity characteristics of the sensing map data of each microarray site in the parallel fitting subtask cluster, and to construct a local heterogeneity perception model. The computation module is used to calculate the adaptive weight coefficients of each microarray site in the dynamic fitting process through the local heterogeneity sensing model; based on the adaptive weight coefficients, a weighted local dynamic fitting objective function corresponding to each microarray site is constructed. The iterative module is used to iteratively solve the weighted local dynamic fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; and to obtain the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip. The correction module is used to construct a spatial neighborhood regularization constraint term based on physical spatial coordinate information; the initial binding rate constant and the initial dissociation rate constant are spatially smoothed and corrected by the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site.

[0008] Thirdly, a computing device includes: One or more processors; A storage device for storing one or more programs that, when executed by one or more processors, cause the one or more processors to implement the method.

[0009] Fourthly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0010] The above-described solution of the present invention has at least the following beneficial effects: By extracting multidimensional temporal features from real-time sensor data and measuring the instantaneous change rate of surface plasmon resonance signals point by point, the fluctuation amplitude of signals at each microarray site can be accurately quantified at different reaction stages. The objective function is fitted with weighted local dynamics, which enables accurate identification and dynamic adjustment of local noise interference regions. This effectively reduces the interference of abnormal fluctuations on the fitting process and improves the calculation accuracy of binding rate constants and dissociation rate constants at each microarray site of the high-density SPR microarray chip.

[0011] Based on the physical spatial coordinate information of the microarray sites, the dynamic parameters in the spatial neighborhood are calculated and corrected by global accumulation. This fully considers the continuous spatial concentration distribution characteristics of the fluid sample when it flows continuously on the surface of the high-density SPR microarray chip. This correction mechanism effectively eliminates unreasonable spatial jumps in dynamic parameters between adjacent microarray sites, ensures a smooth transition between the binding rate constant and the dissociation rate constant on the chip surface, and improves the consistency and reliability of the overall parameter distribution of the high-density SPR microarray chip.

[0012] By calculating the data complexity feature vector and computational load weight of each microarray site, a dynamic task scheduling strategy is adopted to cluster and divide the real-time sensor image data, and a parallel fitting subtask cluster with balanced computational load is constructed. This mechanism effectively avoids the problem of uneven distribution of computational resources when high-density SPR microarray chips process massive real-time sensor image data in batches, shortens the overall time consumption of batch fitting of parallel dynamic parameters, and improves the efficiency of large-scale data processing. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating a method for batch fitting of parallel dynamic parameters for high-density SPR microarrays provided in an embodiment of the present invention.

[0014] Figure 2 This is a schematic diagram of a parallel dynamic parameter batch fitting system for high-density SPR microarrays provided by an embodiment of the present invention. Detailed Implementation

[0015] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0016] like Figure 1 As shown, this embodiment uses a high-density surface plasmon resonance microarray chip as the detection platform and an antibody-antigen molecule interaction system as the detection object. The specific implementation process is as follows.

[0017] The detection system consists of five parts: an optical detection unit, a microfluidic unit, a temperature control unit, a high-density SPR microarray chip, and a data processing unit.

[0018] The optical detection unit employs a Kretschmann prism-coupled surface plasmon resonance imager. The light source is a near-infrared light-emitting diode with a wavelength of 760 nm. After passing through a collimating lens group and a polarizer, it forms a parallel polarized beam, which is then incident at a specific angle onto a gold film on the back of the chip via a high-refractive-index prism. The gold film, with a thickness of 47 nm, is deposited on the prism coupling surface via vacuum evaporation. The reflected light beam is projected onto a 2D complementary metal-oxide-semiconductor (SMOS) photodetector array via an imaging lens group. The detector array has a pixel resolution of 2048×2048, with each pixel having a physical size of 5.5 μm × 5.5 μm and a frame rate of 30 frames per second. The microfluidic unit includes a high-precision syringe pump, a 6-way injection valve, a polydimethylsiloxane (PDMS) microfluidic channel, and a waste liquid collector. The microfluidic channel has a cross-sectional dimension of 500 μm wide and 50 μm high, with a total channel length of 25 mm. The fluid sample is driven by the syringe pump at a constant flow rate of 10 μL per minute, flowing through the surface of the high-density SPR microarray chip. In the temperature control unit, a semiconductor thermoelectric cooler and a high-precision thermistor are built into the chip base to control the chip surface temperature within a range of 25 degrees Celsius ± 0.1 degrees Celsius, in order to eliminate the influence of temperature fluctuations on the molecular binding-dissociation kinetics rate.

[0019] The high-density SPR microarray chip uses a glass slide with a refractive index of 1.72 as its substrate, with a 47 nm gold film deposited on the surface as the SPR sensitive layer. 400 microarray sites are arranged in a 20x20 matrix on the gold film surface using a microarray spotter. Each microarray site has a diameter of 200 μm, and the center-to-center spacing between adjacent sites is 400 μm. The first 10 sites serve as reference channels, immobilized with control proteins that have no affinity for the target analyte, such as bovine serum albumin. The remaining 390 sites are immobilized with different types of monoclonal capture antibody probes, each with a different immobilization density at different sites, used to detect the corresponding antigen analyte in the fluid sample. The antibody probes were covalently immobilized on the gold membrane surface via a gold-thiol self-assembled monolayer and a chemical coupling method involving N-hydroxysuccinimide ester / ethyldimethylaminopropylcarbodiimide. After immobilization, the remaining active sites on each site were sealed with ethanolamine to prevent non-specific adsorption. The data processing unit was a high-performance computing workstation equipped with a 32-core CPU and a graphics processing unit accelerator card, running a batch fitting program for parallel dynamic parameters. This program implemented the computational logic of each step of the method of this invention using C++ language and the CUDA parallel computing framework.

[0020] The specific implementation of step 1, real-time sensor image data acquisition and load weight evaluation, includes: Step 11, Acquisition of initial optical response signal: A high-density SPR microarray chip was mounted onto the prism coupling surface, and refractive index matching oil was filled between the chip and the prism to eliminate interface reflection. The syringe pump was started, and phosphate buffer (pH 7.4, sodium chloride ionic strength 150 mmol / L) was flowed through the chip surface at a rate of 10 μL / min for 30 minutes to allow each microarray site on the chip to reach a stable baseline response. A 6-way injection valve was then used to switch to a fluid sample containing the target antigen analyte. The initial concentration of the antigen in the fluid sample was 100 nanomoles / L, and the solvent and buffer components were identical. Incident light at a specific angle was irradiated onto the gold film on the back of the chip. The incident light excited surface plasmon waves at the interface between the gold film and the dielectric. The propagation characteristics of these waves were extremely sensitive to changes in the refractive index of the extremely thin layer on the gold film surface. As the fluid sample flows over the surface of each microarray site, the antibody probes at the sites specifically bind to the antigen analytes in the fluid. The gradual accumulation of antigen molecules on the surface of the antibody layer causes a continuous increase in the local refractive index of the gold film surface, thereby changing the SPR resonance condition and causing a corresponding attenuation of the reflected light intensity. The 2D photodetector array continuously acquires reflected light images at a rate of 30 frames per second. Each frame records the spatial distribution of reflected light intensity in the 400 microarray site regions across the entire chip's field of view, obtaining a time-varying sequence of initial optical response signals. For the complete binding-dissociation cycle of 400 sites, lasting 300 seconds in the binding phase and 300 seconds in the dissociation phase in this embodiment, a total of 18,000 frames of 2D image data were acquired.

[0021] Step 12, pixel-level position analysis and extraction of raw temporal response data: Based on the design drawings of the dotting mask used in the dotting of high-density SPR microarray chips, the spatial mapping relationship of each microarray site on the chip to the detector pixel array is established. Specifically, the standard center coordinates and diameter information of each microarray site are extracted from the CAD design file of the dotting mask. This standard coordinate is then registered with the pixel coordinates of the photodetector array through an affine transformation. Pre-made positioning markers at the four corners of the chip are used as control points for the affine transformation during the registration process. The rotation matrix and translation vector from the image coordinate system to the chip coordinate system are solved using the least squares fitting method, ensuring that the registration error is less than 0.5 pixels. The complete solution process for the rotation matrix and translation vector based on least squares fitting is as follows: A two-dimensional coordinate affine transformation model is constructed. This registration scenario involves a rigid planar coordinate transformation, involving only coordinate rotation and translation, without scaling or distortion. A transformation model is established from the image pixel coordinate system to the chip physical coordinate system. Let the standard coordinates of any control point in the chip coordinate system be... The measured coordinates in the image pixel coordinate system are The transformation formula is: ; In the formula, The parameters for fitting the rotation matrix are determined by the rotation angle. Decide, , ; These are the translation vector components in the X and Y directions; These correspond to the positioning control points at the four corners of the chip.

[0022] Step 2: Construct a system of least-squares overdetermined equations. Four control points were selected for this registration, resulting in eight sets of transformation equations. The unknown parameters to be solved are... There are four equations in total, forming an overdetermined system that cannot be solved directly and precisely. The overall fitting error must be minimized using the least squares method. A single control point fitting error function is defined as follows: ; The global total fitting error is the sum of the errors of all control points. The goal is to find the optimal combination of parameters that minimizes the total error. To obtain the minimum value, according to the multivariable function extremum theorem, a necessary condition for minimizing the total error is that the partial derivatives of the total error with respect to the four unknown parameters are all zero, that is: ; By expanding and simplifying the partial derivative formula, the optimal solution for the translation vector components is derived: ; In the formula, The average of the chip coordinates of the four control points; The mean of the pixel coordinates of the four control points is used. Substituting the translation parameters into the error equation, the rotation parameters are further simplified and solved: ; ; Solve for the rotation matrix and translation vector, obtained from the fitting. Parameters, construct a two-dimensional rotation transformation matrix : ; Obtained by fitting Parameters, construct translation vector Finally, the transformation matrix from the image coordinate system to the chip coordinate system is solved. The obtained rotation matrix and translation vector are substituted into the transformation formula, and reprojection calculation is performed on the four control points to obtain the reprojection coordinates of each control point. The pixel deviation between the reprojection coordinates of each control point and the standard coordinates is calculated, the root mean square error is solved, and the parameters are iteratively fine-tuned until the root mean square value of the overall registration reprojection error is less than 0.5 pixels, thus completing accurate registration. After pixel-level position analysis and high-precision least square registration, for each microarray site, the pixel area it covers in each frame of the image is determined. Each microarray site has a diameter of 200 micrometers, corresponding to a circular area with a diameter of about 36 pixels on the detector. The arithmetic mean of the response signal intensity values ​​of all pixels in the circular area is calculated to obtain the single response signal value of the microarray site in that time frame. After performing an averaging operation on all 18,000 frames of images, the original time-series response data of each microarray site is obtained. The original time-series response data of each microarray site is a 1-dimensional sequence containing 18,000 time sampling points, with an interval of 1 / 30 seconds between adjacent time sampling points.

[0023] Step 13, Background subtraction, baseline calibration, and real-time sensor data generation: For each microarray site belonging to the measurement channel, the raw time-series response data of the reference channel site closest to its physical location is selected as the background reference signal. The response signal value of the measurement site at each time point is subtracted from the corresponding time point response signal value of the reference site. This operation can effectively eliminate common interference factors such as overall baseline drift caused by changes in the overall refractive index of the fluid, gradual changes in system temperature, and non-specific adsorption on the gold film surface. The signal after background subtraction is baseline calibrated. The average value of the signal within the last 60 seconds before fluid sample injection is taken as the baseline value of the microarray site. The baseline value is subtracted from all signal values ​​of the entire time series of the site to make the baseline zero. The calibrated signal is the real-time sensing image data of the microarray site, with its vertical axis in resonance units. The unit of measurement, one resonance unit, corresponds approximately to the refractive index change caused by the adsorption of 0.1 nanograms of protein per square millimeter on the gold film surface, representing 1 / 1,000,000 of the refractive index change of the surface binding amount. The binding phase is defined as a continuous 300-second period starting from the 120th second, during which the antigen analyte binds to the antibody probes on the surface of each microarray site in continuous flow, and the sensor signal shows an upward trend. The dissociation phase is defined as a continuous 300-second period after the binding phase ends and the buffer is switched back, during which the bound analyte gradually dissociates from the probe and is carried away by the flowing buffer, and the sensor signal shows a downward trend. The real-time sensor data of each microarray site contains a sequence of response signal values ​​at 18,000 time points, fully covering the three phases of baseline, binding, and dissociation.

[0024] Step 14, Multidimensional Temporal Feature Extraction and Data Complexity Feature Vector Calculation: For the real-time sensing data of each microarray site, three-dimensional time-series features are extracted to form a three-dimensional data complexity feature vector. The first dimension is the signal fluctuation variance. Within the steady-state response interval of the binding phase, specifically the time period from 100 seconds to 250 seconds after the start of the binding reaction, this time period corresponds to 4500 continuous sampling points. Within this interval, the sum of squares of the deviations of the response signal value from the mean signal value within this interval is calculated, and then divided by the number of sampling points minus 1 to obtain the signal fluctuation variance. The calculation formula is as follows: ; In the formula, For the first The signal fluctuation variance of each microarray site; This is the total number of sampling points within the steady-state response interval of the combined stage; For the first The microarray site within the interval The response signal value of each sampling point; For the first The variance is the arithmetic mean of the response signal values ​​of each microarray site within the steady-state interval. The physical meaning of this variance is the degree of random fluctuation of the surface plasmon resonance signal near the theoretical steady-state value during the period when the binding reaction has reached a steady state. A smaller variance indicates a higher signal-to-noise ratio and a clearer, more stable signal at that site; a larger variance indicates a more severe influence from factors such as local fluid disturbances, detector electronic noise, or uneven probe fixation. The second extracted dimension feature is the baseline drift slope. Within the baseline interval of 60 seconds before fluid sample injection, a linear trend fitting is performed on 1800 response signal values ​​in this interval. Specifically, with the sampling time point as the independent variable and the response signal value at each time point as the dependent variable, the slope of the line that minimizes the sum of the squares of the differences between the response values ​​at each sampling point and the corresponding values ​​on the fitted line is calculated. This slope is the baseline drift slope, expressed in resonance units per second. The calculation formula is as follows: ; In the formula, For the first Baseline drift slope of each microarray site; This represents the total number of sampling points within the baseline interval; For the first interval within the baseline The time value corresponding to each sampling point; For the first The number of microarray sites within the baseline interval The larger the absolute value of the baseline drift slope of the response signal value of each sampling point, the more obvious the baseline change is caused by factors such as slow aging of the gold film, slight changes in the ionic strength of the buffer solution, or residual fluctuations in temperature control during the detection process. This requires extra attention in subsequent batch fitting.

[0025] The third feature extracted was the number of residual convergence steps. A first-order Langmuir single-exponential combined-dissociation kinetic model was used to perform a preliminary rapid fitting of the real-time sensor data for each microarray site. The fitting process employed the Levenberg-Marquardt iterative algorithm, recording the sum of squared residuals after each iteration.

[0026] The number of iterations required for the residual sum of squares to fall below 5% of its initial value is taken as the residual convergence step number. If the residual sum of squares does not fall below 5% of its initial value within the preset maximum number of iterations, the residual convergence step number is set to 300. This feature reflects the degree of conformity of the sensor map data of the microarray site to the standard first-order Langmuir dynamic model. The fewer the convergence steps, the more standard the data and the easier the fitting. The more the convergence steps, the more complex non-ideal behavior of the data and the greater the fitting difficulty. The above three feature values, signal fluctuation variance, baseline drift slope and residual convergence step number, are combined to form a 3-dimensional vector, which is the data complexity feature vector of the microarray site. The three components of the vector characterize the overall processing complexity of the sensor map data of the site from three dimensions: signal noise level, system drift degree and dynamic model matching difficulty.

[0027] Step 15, calculate the weighted sum of the load weights: The three components of the data complexity feature vector of each microarray site are subjected to dimensionless normalization, and the three features with different dimensions and value ranges are uniformly mapped to the dimensionless interval from 0 to 1. The specific normalization method is as follows: For the signal fluctuation variance component, the signal fluctuation variance value of the microarray point is divided by the maximum value of the signal fluctuation variance among all 400 microarray points; for the baseline drift slope component, the absolute value of the baseline drift slope of the microarray point is taken and then divided by the maximum value of the absolute value of the baseline drift slope among all 400 microarray points; for the residual convergence step number component, the residual convergence step number of the microarray point is divided by the maximum value of the residual convergence step number among all 400 microarray points. The three normalized feature components are then weighted and summed according to preset weight coefficients. The three weight coefficients are determined based on empirical evaluation of the actual contribution of different types of data complexity features to the fitting difficulty: the weight of signal fluctuation variance is 0.4, as this feature directly reflects the noise level of the sensor image data and has the greatest impact on fitting accuracy; the weight of baseline drift slope is 0.2, as this feature reflects the gradual change of the system and has a secondary impact on fitting accuracy; and the weight of residual convergence step number is 0.4, as this feature directly reflects the difficulty of matching the data with the theoretical model. The scalar value obtained after weighted summation is the computational load weight of that microarray site, ranging from 0 to 1. The larger the computational load weight value, the more difficult the processing of the sensor map data of that microarray site is, and more computational resources should be allocated in subsequent parallel task allocation to ensure fitting accuracy; the closer the computational load weight value is to 0, the higher the quality of the sensor map data of that site is and the simpler to process, and fewer computational resources can be allocated.

[0028] Step 2, specifically the parallel fitting of subtask cluster partitioning for dynamic task scheduling and computational load balancing, includes the following: Step 21, Calculate the cluster load baseline: The computational load weight values ​​of all 400 microarray sites are summarized. The sum of these 400 computational load weight values ​​is then calculated to obtain the total load weight value for all microarray sites. This total load weight is divided by the preset number of parallel computing processes to obtain the cluster load baseline value. The cluster load baseline value is the ratio of the total load weight value of all microarray sites to the number of parallel computing processes, calculated using the following formula: ; In the formula, This is the baseline value for cluster load. The preset number of parallel computing processes; This represents the total number of microarray sites; For the first The computational load weight of each microarray site is determined by the data processing unit, which is equipped with 32 central processing unit cores. Considering that the operating system kernel process and input / output management need to occupy some core resources, 20 cores are actually allocated for parallel dynamic parameter fitting calculation. Therefore, the preset number of parallel computing processes is 20. The calculation method is to divide the total load weight value by 20. The quotient is the cluster load benchmark value. This benchmark value represents the normalized computational load that each parallel fitting subtask cluster should bear under the ideal uniform distribution.

[0029] Step 22, initial task queue partitioning under the descending first-fit allocation strategy: A descending first-fit allocation strategy was adopted to distribute the real-time sensor map data of all 400 microarray sites to 20 initial fitting subtask clusters. The allocation process was as follows: the 400 microarray sites were arranged in descending order of computational load weight value to form a sequence of sites to be allocated. The site ranked first in the sequence had the largest computational load weight, and the site ranked last had the smallest computational load weight. The site with the largest load weight among the currently unallocated sites was taken out from the sequence of sites to be allocated in sequence, and the total allocated load of each of the 20 initial subtask clusters was checked. The real-time sensor data corresponding to the site is assigned to the initial subtask cluster with the smallest total assigned load. After assignment, the calculated load weight value of the site is added to the total assigned load of the cluster. The above four steps of extracting, finding the minimum load cluster, assigning, and updating are repeated until all 400 microarray sites in the site sequence to be assigned are assigned to an initial fitting subtask cluster, ultimately forming 20 initial fitting subtask clusters. Due to the use of a descending first-fit strategy, sites with large load weights are preferentially assigned to clusters with lighter loads. After the initial assignment, the load difference between the clusters has been initially controlled. Each cluster contains real-time sensor data of approximately 18 to 22 microarray sites.

[0030] Step 23, Cross-cluster data migration and load balancing optimization based on load deviation values: Based on the initial allocation, further cross-cluster data migration operations are used to eliminate residual load imbalance. The actual total load of each initially fitted subtask cluster is calculated, which is the sum of the calculated load weights of all allocated microarray sites within the cluster. This actual total load is compared with the cluster load baseline to calculate the load deviation. The deviation is equal to the actual total load minus the algebraic difference of the baseline value. A positive deviation indicates that the cluster is overloaded, and a negative deviation indicates that the cluster is underloaded. The deviation tolerance threshold is set at 15% of the cluster load baseline. When the absolute value of the load deviation of any cluster exceeds the tolerance threshold, a cross-cluster data migration operation is triggered. The specific process of the migration operation is as follows: From the clusters with positive load deviation values ​​that exceed the tolerance threshold, the computational load weight values ​​of each microarray site contained in each cluster are checked one by one. The site whose computational load weight value is closest to the current deviation value is found, and the real-time sensor map data corresponding to the selected site is removed from the overloaded cluster. At the same time, it is added to the cluster with the current negative load deviation value and the largest absolute value of the deviation. After the migration is completed, the total load of the two affected clusters is recalculated and their respective load deviation values ​​are updated. The above deviation detection and cross-cluster data migration operation is repeated. After each round of migration operation, the absolute value of the maximum load deviation value gradually decreases. When the absolute value of the load deviation values ​​of all 20 clusters drops to within 15% of the baseline value, the migration operation terminates, and the maximum actual load difference between the clusters does not exceed 30% of the baseline value, achieving a state of computational load balance. The 20 parallel fitting subtask clusters obtained in the end are the parallel fitting subtask clusters for computational load balance.

[0031] The specific implementation of step 3, the construction of the local heterogeneity perception model, includes the following: Step 31, Calculation of local residual time series: For each microarray site in the cluster of 20 parallel fitting subtasks, the local residual time series of its real-time sensor map data is calculated. Specifically, the response signal values ​​of two consecutive adjacent sampling times in the real-time sensor map data of the microarray site are subtracted, that is, the response value of the later time time is subtracted from the response value of the previous time time. The difference is the time series differential signal value within the adjacent sampling interval. The calculation formula is as follows: ; In the formula, For the first The first microarray site The time-series differential signal value of each sampling interval; For the first The first microarray site Real-time sensor response signal value at each sampling time; The total number of sampling points for the real-time sensor image data. The physical meaning of this time-series differential signal value is: the actual instantaneous change of the surface plasmon resonance response signal within the adjacent sampling interval. Arranging the differential values ​​of all adjacent sampling point pairs in the entire time series in chronological order, we obtain a local residual time series containing approximately 17,999 differential values. This local residual time series reflects the instantaneous fluctuation characteristics of the sensor image signal deviating from the ideal smooth single exponential change trend during the binding and dissociation process. This residual sequence should fluctuate slightly randomly around the value of 0. At microarray sites with poor signal quality, a large instantaneous jump appears in the residual sequence. These jumps correspond to abnormal peaks or dips on the sensor image signal curve.

[0032] Step 32, Local residual distribution feature extraction and spatial neighborhood site determination: For the local residual time series of each microarray site, two higher-order statistical moments, skewness and kurtosis, are calculated as morphological characteristics describing the residual distribution at that site. The skewness value is calculated as follows: the deviation of each difference value in the residual time series from the mean of the entire series is cubed, and the average value is obtained. Then, this average value is divided by the cube of the standard deviation of the series. The calculation formula is as follows: ; In the formula, For the first The residual distribution skewness value of each microarray site; The total length of the local residual time series; For the first The kurtosis is the arithmetic mean of the local residual time series of each microarray site. A positive skewness value indicates that the residual distribution has a long tail in the positive direction, a negative skewness value indicates that the residual distribution has a long tail in the negative direction, and a skewness value close to 0 indicates that the residual distribution is roughly symmetrical. The kurtosis value is calculated by raising the deviation of each difference value in the residual time series to the fourth power of the mean of the entire series, taking the average value, and then dividing this average value by the fourth power of the standard deviation of the series. The calculation formula is as follows: ; In the formula, For the first The kurtosis value of the residual distribution of each microarray site is used. A kurtosis value greater than 3 indicates that the residual distribution is more angular than a normal distribution, with residuals concentrated in a small fluctuation range, but occasionally exhibiting extreme large fluctuations. A kurtosis value close to 3 indicates that the residual distribution is close to a normal distribution. A kurtosis value significantly less than 3 indicates that the residual distribution is relatively flat. Skewness and kurtosis values ​​together constitute the local residual distribution characteristics, used to quantitatively characterize the statistical distribution of noise in the sensor map of each microarray site. Based on the array index position of each microarray site in the 20x20 matrix array, its spatial neighborhood site set is determined. The specific determination rule is: for each microarray site located in the matrix... For each microarray point within the matrix, its spatial neighborhood is formed by selecting eight adjacent points centered on its array coordinates: the point above, the point below, the point to the left, the point to the right, the point to the upper left, the point to the upper right, the point to the lower left, and the point to the lower right. For microarray points located at the matrix edges, adjacent points within the 400-point array range are selected. The neighborhood of the upper edge point contains five points, and the same applies to the lower edge point. The neighborhood of the left edge point contains five points, and the same applies to the right edge point. For points located at the four corners of the matrix, their neighborhood contains three points.

[0033] Step 33, Calculation of spatial neighborhood heterogeneity characteristics: For each microarray site, the spatial neighborhood heterogeneity characteristics are calculated using real-time sensor data from its neighboring sites. Specifically, the real-time sensor data of the microarray site is compared time-by-time with the real-time sensor data of each neighboring site in its neighboring site set. At each sampling time point, the absolute value of the difference between the response signal values ​​of the microarray site and its neighboring sites is calculated. The arithmetic mean of the absolute differences over the entire time series is then calculated to obtain the difference in response signal between the microarray site and its neighboring sites, expressed in resonance units. The calculation formula is as follows: ; In the formula, For the first The microarray site and its first The difference in response signals between neighboring sites; For the first The _ neighboring sites The response signal value at each sampling time; The total number of sampling points for the real-time sensor image data is represented by this difference, which reflects the degree of deviation of the time-averaged response signal between a pair of adjacent microarray sites due to factors such as differences in probe fixation density, local fluid velocity, or differences in probe-analyte affinity at the two sites. The spatial neighborhood heterogeneity characteristic of the microarray site is obtained by calculating the arithmetic mean of the response signal differences between the microarray site and all its neighboring sites, as shown in the following formula: ; In the formula, For the first Spatial neighborhood heterogeneity characteristic values ​​of microarray sites; For the first The total number of spatial neighborhood sites of each microarray site. The physical meaning of this feature is: the degree of deviation of the average sensor map response of the site relative to its immediate spatial neighborhood. The larger the heterogeneity feature value, the more obvious the influence of local fluid concentration field inhomogeneity or probe fixation mass difference on the site, and the lower the confidence weight should be given when fitting dynamic parameters; the smaller the heterogeneity feature value, the more consistent the sensor map response of the site with the surrounding neighborhood sites, and the higher the data reliability.

[0034] Step 34, Construction and pre-training of the multilayer perceptron network: A multilayer perceptron network containing a residual distribution coding layer and a spatial heterogeneity coding layer is constructed as a local heterogeneity sensing model. This model is used to automatically infer the differential weights that should be assigned to each microarray site in the subsequent dynamic fitting process, based on the local residual distribution characteristics and spatial neighborhood heterogeneity characteristics of each microarray site. The specific structure of the multilayer perceptron network is as follows: The first layer is the input layer, which receives a 4-dimensional input vector. The four components of the input vector are: the first component is the skewness value of the microarray site, the second component is the kurtosis value of the microarray site, the third component is the spatial neighborhood heterogeneity characteristic value of the microarray site, and the fourth component is the load weight value of the microarray site. Before being passed to the next layer, the four input components are z-score standardized, that is, each component is subtracted from the mean of the component over all 400 microarray sites and then divided by the standard deviation of the component to eliminate the dimensional and numerical scale differences among the four components. The second layer is the residual distribution coding layer, i.e., the first hidden layer, which contains 16 neuron nodes. Each neuron node is connected to the four nodes of the input layer through fully connected weights.

[0035] Each neuron operates as follows: it multiplies the signal values ​​from the four nodes in the input layer by their respective connection weights, sums them, adds its own bias term, and outputs the sum after a nonlinear transformation using a modified linear unit activation function. The modified linear unit maps negative values ​​to 0 and keeps positive values ​​unchanged. Its function is to introduce nonlinear expressive power, enabling the network to learn the possible nonlinear complex mapping relationship between residual distribution characteristics and target weights. The physical meaning of this layer is residual distribution encoding, which encodes the original residual statistical characteristics and related signal quality indicators into a set of 16-dimensional intermediate latent variables. This set of latent variables captures the inherent discriminative information of different microarray sites in the signal quality dimension. The third layer is the spatial heterogeneity encoding layer, i.e., the second hidden layer, containing 8 neuron nodes. This layer operates in the same way as the first hidden layer: each node multiplies the output values ​​from the 16 nodes in the first hidden layer by their respective connection weights, sums them, adds a bias term, and then outputs the sum after passing it through the modified linear unit activation function. The physical meaning of this layer is spatial heterogeneity encoding. It further compresses and integrates the 16-dimensional signal quality latent variables extracted from the first hidden layer into an 8-dimensional high-level feature representation. During dimensionality reduction, it retains the information dimensions most critical to weight decisions while filtering out redundant and noisy dimensions. The fourth layer is the output layer, containing only one neuron node. This node receives the outputs from the eight nodes in the second hidden layer, performs a fully connected weighted summation, adds a bias term, and maps the result to a numerical range between 0 and 1 using a sigmoid function. The output of the sigmoid function approaches 1 when the input is a large positive value, approaches 0 when the input is a large negative value, and outputs 0.5 when the input is 0. This output value is the adaptive weight coefficient calculated by the local heterogeneity sensing model for this microarray site.

[0036] The parameters of all connection weights and bias terms in the entire multilayer perceptron network were determined through offline pre-training. The pre-training process is as follows: Several sets of sensor map samples covering different signal quality levels and different degrees of spatial heterogeneity were selected from the SPR sensor map dataset accumulated by our laboratory. Each set of samples contains complete real-time sensor map data of one microarray site, and the optimal local fitting weight allocation scheme determined by experienced experimenters through repeated manual fitting and cross-validation is used as training labels. The skewness, kurtosis, spatial neighborhood heterogeneity feature values, and computational load weight values ​​of the samples are used as network inputs, and the manually determined optimal weights are used as the expected outputs. The mean squared error loss function and backpropagation algorithm are used to iteratively update all connection weights and bias term parameters in the network using the mini-batch stochastic gradient descent method until the loss function value no longer decreases on an independent validation sample set. After pre-training is completed, all network parameters are fixed to form a local heterogeneous sensing model that can be directly deployed online.

[0037] Step 4, specifically the calculation of adaptive weight coefficients and the construction of the weighted local dynamics fitting objective function, includes the following: Step 41, forward propagation calculation of adaptive weight coefficients: For each microarray site to be processed in the current batch, the skewness and kurtosis values, spatial neighborhood heterogeneity features, and load weights of that site are combined into a 4-dimensional input vector. These four input values ​​are then z-score standardized and fed into the input layer of a pre-trained, constructed local heterogeneity perception model of a multilayer perceptron. In the residual distribution coding layer, each of the 16 neurons performs fully connected weighted summation and bias addition on the four standardized input values. After transformation by a modified linear unit activation function, 16 intermediate feature values ​​are output. This layer maps the original 4-dimensional feature data into a 16-dimensional latent variable space, capturing the essential differences in signal quality between different microarray sites. In the differential mode, at the spatial heterogeneity coding layer, eight neurons perform a second fully connected weighted accumulation and bias addition on the 16 intermediate feature values ​​from the first hidden layer. After transformation by the modified linear unit activation function, eight high-level feature values ​​are output. This layer further compresses and integrates the 16-dimensional intermediate representation, extracting the few latent variable dimensions that are most critical to the weight decision. At the output layer, a single neuron performs a third fully connected weighted accumulation and bias addition on the eight high-level feature values ​​from the second hidden layer. The accumulation result is mapped to a scalar value between 0 and 1 through a sigmoid function. This scalar value is the adaptive weight coefficient calculated by the local heterogeneity perception model for this microarray site.

[0038] The physical meaning and value rules of the adaptive weight coefficient are as follows: For microarray sites with low signal fluctuation variance, basically symmetrical and close to normal residual distribution, and low spatial neighborhood heterogeneity, the model outputs a larger adaptive weight coefficient, usually between 0.8 and 0.95. During the fitting of dynamic parameters, the measured data points of these high-quality sites will be given a higher contribution weight, and the fitting process will rely more on these reliable data. For microarray sites with high signal fluctuation variance, severely skewed residual distribution or heavy tails, and high spatial neighborhood heterogeneity, the model outputs a smaller adaptive weight coefficient, usually between 0.1 and 0.3. During the fitting process, the interference of abnormal fluctuations of these low-quality sites on the fitting results is effectively suppressed. The weight coefficient plays the role of softly filtering out abnormal data points, which is different from the simple hard rejection in traditional methods. This avoids the information loss caused by rejecting real but low signal-to-noise ratio data points.

[0039] Step 42, Generation of the weighted local response bias sequence: Using adaptive weighting coefficients, the deviation between the real-time sensor data obtained from actual measurements and the predicted values ​​of the theoretical dynamic model at each microarray site is differentially weighted. Specifically, for each sampling time point of each microarray site, the deviation between the actual measured response signal value and the theoretical response signal value predicted by the current first-order Langmuir dynamic model at that time point is calculated. The deviation value equals the measured value minus the predicted value. The sign of the deviation value indicates whether the measured value is higher or lower than the theoretical predicted value. Multiplying this deviation value by the adaptive weighting coefficient of the microarray site yields the weighted local response deviation at that time point. The core effect of this operation is that for a high-quality site with an adaptive weighting coefficient of 0.9, the response deviation at each time point remains almost unchanged, contributing significantly to the objective function. For a low-quality site with an adaptive weighting coefficient of 0.2, the response deviation at each time point is significantly compressed to 1 / 5 of its original value, reducing the interference of abnormal fluctuations at that site on the overall fitting process. Arranging the weighted local response deviation values ​​of all 18,000 sampling time points of each microarray site in chronological order yields the weighted local response deviation sequence of that microarray site.

[0040] Step 43: Construction of the weighted local dynamics fitting objective function Using the weighted local response deviation sequence as the core input, a weighted local dynamic fitting objective function is constructed for each microarray site. This objective function uses the binding rate constant and dissociation rate constant as two independent variables to be optimized. The construction of the objective function is based on the classical first-order Langmuir binding-dissociation dynamics model. During the binding stage, the evolution of the theoretical response signal value over time is determined by the binding rate constant, the analyte concentration in the fluid, the maximum response capacity of the site, and the time variable. Its typical function shape is a curve that asymptotically rises to the maximum response capacity in an exponential manner. The calculation formula is as follows: ; In the formula, For the first Each microarray site in The theoretical response signal value at time t; For the first Maximum response capacity of each microarray site; This is the binding rate constant; It is the dissociation rate constant; The initial concentration of the analyte in the fluid sample; The time variable represents the current sampling time. To determine the duration of the binding phase, during the dissociation phase, the theoretical response signal value is determined solely by the dissociation rate constant, the actual response value at the end of the binding phase, and the time variable. Its typical function form is a curve that gradually decays exponentially to the baseline, starting from the response value at the end of the binding phase. The calculation formula is as follows: ; In the formula, For the first The measured response signal values ​​of each microarray site at the end of the binding phase are used to construct the weighted local dynamics fitting objective function for each microarray site. The specific construction method is as follows: For each sampling time point in the time series of that site, calculate the difference between the actual measured response value and the theoretical response value predicted by the Langmuir model at that time. Multiply this difference by an adaptive weighting coefficient and square it. Summate the weighted squared deviation values ​​of all 18,000 time points to obtain a scalar function value. This scalar value is the weighted local dynamics fitting objective function value for that microarray site under the current binding rate constant and dissociation rate constant values. The calculation formula is as follows: ; In the formula, For the first Weighted local dynamics fitting objective function for each microarray site; For the first Adaptive weighting coefficients for each microarray site. For the first The first microarray site The real-time sensor response signal value at each sampling time; the smaller the objective function value, the better the fit of the current parameter value to the weighted measured data; the larger the objective function value, the worse the fit. Due to the setting of adaptive weight coefficients, this objective function is different from the traditional uniform weighted least squares objective function. Under uniform weighting, all approximately 7,200,000 data points of all 400 microarray sites have the same unit contribution weight in the objective function. Noise and abnormal fluctuations of low-quality sites are equally weighted with the accurate data of high-quality sites and included in the fitting residual. Data from high-quality sites are automatically assigned a higher contribution weight, while the contribution of data from low-quality sites is artificially compressed, so that the optimization process naturally tends to fit those data with high reliability.

[0041] The specific implementation of step 5, including iterative solution of dynamic parameters and acquisition of physical space coordinate information, is as follows: Step 51: Iterative optimization using the Levenberg-Marquardt algorithm For the weighted local kinetic fitting objective function corresponding to each microarray site, the Levenberg-Marquardt algorithm is used for iterative optimization to solve for the optimal binding rate constant and dissociation rate constant of that site in a weighted sense. The detailed process of iterative optimization is as follows: parameter initialization, setting the initial guess value of the binding rate constant to 10000 per mole per second, and the initial guess value of the dissociation rate constant to 0.001 per second. The selection of these two initial values ​​is based on the common order of magnitude range of kinetic parameters of typical medium affinity antibody-antigen molecule interaction systems. The initial value of the Levenberg-Marquardt damping parameter is set to 0.01. The parameter vector to be optimized in the next iteration is The initial iteration parameters are ,in , Initial values ​​of damping parameters The process then enters the main iterative loop. In each iteration, the following five sub-operations are performed: Calculate the objective function value for the current binding rate constant and dissociation rate constant; substitute the current parameters into the Langmuir dynamics model to calculate the theoretical predicted response sequence for this site across all 18,000 sampling time points; calculate the difference between the measured and predicted values ​​for each time point, multiply by an adaptive weighting coefficient, square the result, and sum all weighted squared differences to obtain the current objective function value. The calculation formula is as follows: ; In the formula, For the first The number of microarray sites in the first The objective function value in the next iteration; , For respectively the first The values ​​of the binding rate constant and dissociation rate constant are taken in the next iteration. The gradient vector of the objective function with respect to the two parameters to be optimized is calculated. The first component of the gradient vector is the partial derivative of the objective function with respect to the binding rate constant. While keeping the dissociation rate constant constant constant, a small perturbation is applied above and below the current binding rate constant value. The objective function values ​​corresponding to the two perturbation values ​​are calculated respectively. The difference between the two values ​​is divided by twice the perturbation amount to obtain an approximate value of the partial derivative. The second component of the gradient vector is the partial derivative of the objective function with respect to the dissociation rate constant. The rules for the value of the perturbation amount are as follows: ; The finite difference approximation formula for the two partial derivatives is as follows: ; ; The gradient vector is composed of two partial derivative components: ; In the formula, For the first The gradient vector of the objective function at the next iteration; , To represent the finite-difference perturbations of the binding rate constant and the dissociation rate constant, respectively, an approximate Hessian matrix is ​​constructed. The Hessian matrix is ​​a 2x2 symmetric matrix whose four elements are the second-order partial derivatives of the objective function with respect to the binding rate constant, the mixed partial derivatives of the objective function with respect to the binding and dissociation rate constants, and the second-order partial derivatives of the objective function with respect to the dissociation rate constant. The values ​​of each element are approximated by the difference in the rate of change of the gradient vector components in different directions of the parameter perturbation. The approximate Hessian matrix form is as follows: ; The approximate formula for the difference calculation of each matrix element is as follows: ; ; ; In the formula, For the first Approximate Hessian matrix at the next iteration; This is an approximation of the second-order partial derivative of the objective function with respect to the binding rate constant; This is an approximation of the second-order partial derivative of the objective function with respect to the dissociation rate constant; , The search direction is approximated by the mixed second-order partial derivatives of the objective function with respect to two parameters. The search direction is controlled by the Levenberg-Marquardt damping parameter. When the damping parameter is small, the search direction is mainly determined by the product of the inverse of the approximate Hessian matrix and the gradient vector, approaching the search direction of the Gauss-Newton method, suitable for achieving fast convergence when the parameters are close to the optimal solution. When the damping parameter is large, the search direction is mainly determined by the product of the scaling of the identity matrix and the gradient vector, approaching the search direction of the steepest gradient descent method, suitable for robust approximation when the parameters are far from the optimal solution. The specific search direction is obtained by solving a system of 2x2 linear equations with the approximate Hessian matrix plus the damping parameter multiplied by the second-order identity matrix as the coefficient matrix and the negative gradient vector as the right-hand side. The equation system is as follows: ; In the formula, For the first Damping parameters at the next iteration; It is a second-order identity matrix; For the first The parameter update increment vector for the current iteration is used to solve the system of linear equations. and The parameters are updated and evaluated along the search direction. The current binding rate constant and dissociation rate constant are shifted along the search direction by a step size controlled by both the gradient magnitude and the damping parameter to obtain candidate new parameter values. The objective function value is recalculated using the candidate new parameters. If the new objective function value is less than the current objective function value, the parameter update is accepted, and the binding rate constant and dissociation rate constant are updated to candidate values, while the damping parameter is halved. If the new objective function value is not less than the current objective function value, the update is rejected, the parameters remain unchanged, and the damping parameter is doubled. The formula for calculating the candidate parameters is: The rules for parameter updates and damping adjustments are as follows: ; ; In the formula, , For candidate updated binding rate constants and dissociation rate constants; For the objective function value corresponding to the candidate parameters, the iteration termination condition is determined by checking whether the following two termination conditions are met after each parameter update: the relative decrease in the objective function value in three consecutive iterations is less than one part per million. The formula for calculating the relative decrease is: ; When all three consecutive iterations satisfy When convergence is reached, the iteration loop is determined. If the current iteration count has reached the preset upper limit of iteration count (500 times in this embodiment), the iteration loop is exited when any termination condition is met. The current parameter value at the time of iteration exit is the initial binding rate constant and initial dissociation rate constant of the microarray site under weighted fitting. The above-mentioned Levonburg-Marquardt iterative solution process is executed in parallel in 20 parallel subtask clusters. Each parallel subtask cluster is assigned to an independent central processing unit core of the data processing unit to run on, processing the dynamic parameter fitting task of about 18 to 22 microarray sites. The 20 cores run in parallel, so that the total time for batch solving the dynamic parameters of 400 microarray sites is about 1 / 20 of the time for solving a single site.

[0042] Step 52: Extract the center pixel coordinates based on the dot pattern mask design drawing The high-density SPR microarray chip's dot matrix mask design drawing is loaded into the coordinate analysis program. The microarray point dot layer is identified from the layer information of the design drawing. All geometric entities in this layer are traversed. Each microarray point is represented as a circular boundary data in the GDSII file, containing the center coordinates and radius parameters. For each circular entity, the coordinates of its center in the GDSII drawing coordinate system are extracted as the center pixel coordinates of that microarray point. The origin of the drawing coordinate system is usually located at the alignment mark at the lower left corner of the chip, with the X-axis pointing to the right and the Y-axis pointing upwards, in nanometers. For 400 microarray points arranged in a 20x20 matrix, the center pixel coordinates of each microarray point can be directly calculated from the array arrangement pattern and spacing parameters. Let the... Line number The x-coordinate of the center pixel of the column microarray site is The vertical axis is row index With column index The values ​​for are all between 1 and 20, and the coordinate calculation formula is as follows: ; ; In the formula, For the first Line number The x-coordinate of the center pixel of the column microarray site; For the first Line number The ordinate of the center pixel of the column microarray site; This is the starting X offset of the array, which is the x-coordinate of the center pixel corresponding to the position in the first row and first column; This is the starting Y offset of the array, which is the ordinate of the center pixel corresponding to the position in the first row and first column; The horizontal pixel spacing is the pixel interval between the centers of two adjacent columns of points in the X direction. The vertical pixel spacing is the pixel interval between the centers of two adjacent rows of points in the Y direction. The center pixel coordinates of all 400 points are summed to obtain the set of center pixel coordinates.

[0043] Step 53: Mapping from pixel coordinates to physical space coordinates A mapping ratio between the set of center pixel coordinates and the actual physical size of the chip is established. The center pixel coordinates of the chip boundary reference marks are extracted from the GDSII design drawing of the dot-matrix mask. The actual physical distance between the alignment marks at the four corners of the chip is measured using a high-precision optical coordinate measuring machine or scanning electron microscope. The ratio of the pixel spacing in the drawing to the actual physical spacing of the two diagonal marks in the X direction is calculated to obtain the pixel-to-micrometer mapping ratio coefficient in the X direction. The mapping ratio coefficient in the Y direction is calculated by multiplying the X component of the center pixel coordinates of each microarray point by the X-direction mapping ratio coefficient and the Y component by the Y-direction mapping ratio coefficient to obtain the physical spatial coordinates of that microarray point on the surface of the high-density SPR microarray chip, with the coordinate unit being micrometers. The physical spatial coordinate system has the upper left corner of the chip as the origin, with the positive X-axis pointing horizontally to the right and the positive Y-axis pointing vertically downwards. After mapping, the physical center spacing between adjacent microarray points in this embodiment should be 400 micrometers. The above coordinate mapping operation is performed on all 400 microarray points one by one to obtain the physical spatial coordinates of each point, where the point number ranges from 1 to 400.

[0044] Step 6, specifically the spatial neighborhood regularization and spatial smoothing correction of dynamic parameters, includes the following: Step 61: Constructing Spatial Neighborhood Regularization Constraints Based on Physical Spatial Coordinates Using the physical spatial coordinates of all 400 microarray sites, the Euclidean distance between any two microarray sites is calculated. For any two microarray sites, the Euclidean distance is the square root of the sum of the squares of the differences in their X and Y coordinates, expressed in micrometers. The calculation formula is as follows: ; In the formula, For the first The microarray site and the first Euclidean distance between microarray sites; , For respectively the first The physical spatial horizontal and vertical coordinates of each microarray site; , For respectively the first The physical spatial coordinates of each of the 400 microarray sites were used to calculate the Euclidean distances between each pair of sites, resulting in a 400x400 neighborhood distance matrix. The elements in the corresponding rows and columns represent the physical spatial distance between two sites. The Euclidean distance between sites in adjacent rows or columns is approximately 400 micrometers, and the Euclidean distance between diagonally adjacent sites is approximately 565.7 micrometers. The Euclidean distance between sites separated by multiple rows or columns gradually increases with increasing spatial interval. A Gaussian weight matrix was constructed based on the neighborhood distance matrix as a spatial neighborhood regularization constraint. The construction rule is as follows: for each element in the distance matrix, its corresponding Gaussian weight value is calculated. This weight value decreases exponentially with the increase of the ratio of the squared distance to the squared distance decay parameter. The calculation formula is as follows: ; In the formula, For the first The microarray site and the first Gaussian weight values ​​between microarray sites; The distance attenuation parameter is set to 200 micrometers. Based on this value, the Gaussian weight between directly adjacent sites at a distance of 400 micrometers is approximately exp(-4), or about 0.018. The Gaussian weight between diagonally adjacent sites at a distance of about 565.7 micrometers is approximately exp(-8), or about 0.0003. The Gaussian weight between distant sites at a distance of more than 1000 micrometers approaches 0, and is considered 0 in actual calculations to save computation. For each site itself, its Gaussian weight is set to 1. All Gaussian weight values ​​form a 400-row, 400-column Gaussian weight matrix. This matrix has the following characteristics: for each site, sites within its finite neighborhood have non-zero weight values, while sites further away have zero weight; the weight values ​​decay rapidly with increasing distance, reflecting the physical assumption that the closer the spatial distance between sites, the more similar their dynamic parameters should be, based on the continuous spatial distribution of the fluid concentration field.

[0045] Step 62: Calculation of spatial smoothing correction gradient Using spatial neighborhood regularization constraints, the parameter difference penalty value of the initial dynamic parameters of each microarray site relative to its spatial neighborhood sites is calculated, i.e., the spatial smoothing correction gradient. The smoothing correction gradient for the binding rate constant is calculated as follows: For any microarray site, all 400 sites are traversed. For each neighboring site, the difference between the initial binding rate constant of that site and the initial binding rate constant of its neighbors is calculated. This difference is multiplied by the weight value between two points in the Gaussian weight matrix. The weighted differences of all neighboring sites are then summed. This sum is the spatial smoothing correction gradient of that site in the dimension of the binding rate constant. The smoothing correction gradient for the dissociation rate constant is obtained using the exact same calculation method, only the calculation object is replaced with the initial dissociation rate constant. The calculation formulas for both types of gradients are as follows: ; In the formula, , The first Spatial smooth correction gradients for binding rate constants and dissociation rate constants of individual microarray sites; , The first , No. The initial binding rate constant of each microarray site; , The first , No. The initial dissociation rate constant of each microarray site; the physical meaning and function of this gradient are as follows: When the initial binding rate constant of a site is significantly higher than that of its neighboring sites, the difference is positive for its neighboring sites. Multiplying by a positive weight yields a positive contribution term, and the sum of all neighboring sites is the positive gradient value. A positive gradient means that subsequent gradient descent updates will be adjusted in the direction of decreasing the binding rate constant of the site, making it move closer to the lower level of its surrounding neighborhood. When the initial binding rate constant of a site is significantly lower than that of its neighboring sites, the sum of all neighboring sites is the negative gradient value, and subsequent updates will increase the binding rate constant of the site. When the binding rate constant of the site is basically the same as that of all its neighboring sites, the positive and negative contributions cancel each other out, and the sum approaches 0.

[0046] Step 63: Spatial smoothing correction based on gradient descent update and acquisition of final target parameters Using spatial smoothing correction gradients, the initial dynamic parameters of all 400 microarray sites are updated sequentially via gradient descent to eliminate unreasonable spatial jumps in parameters between adjacent sites, obtaining the final target dynamic parameters with a smooth spatial transition. A regularization strength coefficient is introduced into the gradient descent update. This coefficient controls the strength of the spatial smoothing correction; the larger the coefficient, the stronger the smoothing effect, and the tighter the compression of parameter differences between adjacent sites; the smaller the coefficient, the weaker the smoothing effect, and the more original spatial differences are preserved. In this embodiment, the regularization strength coefficient is optimized to 0.1. The optimization basis for this value is: selecting chip test datasets with known spatial dynamic parameter gradients from historical data, comparing the spatial smoothing correction results under different regularization strength coefficients with the measured values ​​from independent verification experiments, and selecting the regularization strength coefficient that minimizes the root mean square deviation between the corrected parameter values ​​and the independent verification experiment values ​​as the optimal value. The specific gradient descent update operation for each microarray site is as follows: The updated binding rate constant is obtained by subtracting the product of the regularization intensity coefficient and the smoothing correction gradient of the binding rate constant from the initial binding rate constant. The updated dissociation rate constant is obtained by subtracting the product of the regularization intensity coefficient and the smoothing correction gradient of the dissociation rate constant from the initial dissociation rate constant, and the calculation formula is as follows: ; In the formula, , The first The final binding rate constant and dissociation rate constant after correction of each microarray site; As the regularization intensity coefficient, the gradient descent update operation is performed on all 400 microarray sites one by one to obtain the final target binding rate constant and target dissociation rate constant of each microarray site. The site numbers are taken from 1 to 400. After spatial smoothing correction, the final target dynamic parameters have the following physical rationality characteristics: In the adjacent site group distributed in strips along the fluid flow direction on the chip, since these sites are in a basically the same local fluid concentration environment, there should be no unreasonable abrupt jumps between their target binding rate constants. For example, there should not be a situation where the binding rate constant value of a certain row of sites is all about 50,000 per second per mole, while the corresponding value of the next adjacent row of sites suddenly jumps to 300 per second per mole. The spatial smoothing correction mechanism naturally eliminates such unreasonable spatial jumps caused by local data noise or fitting trapped in local extrema, so that the target binding rate constant and target dissociation rate constant on the chip surface exhibit a smooth spatial transition mode consistent with the continuous distribution characteristics of the fluid concentration field.

[0047] like Figure 2 As shown, embodiments of the present invention also provide a parallel dynamics parameter batch fitting system for high-density SPR microarrays, comprising: The acquisition module is used to acquire real-time sensing image data of each microarray site on the high-density SPR microarray chip, extract multi-dimensional time-series features from the real-time sensing image data to obtain a data complexity feature vector, and calculate the computational load weight of each microarray site through the data complexity feature vector. The partitioning module is used to cluster real-time sensor image data based on computational load weights and a dynamic task scheduling strategy to obtain multiple parallel fitting subtask clusters with balanced computational load. The module is used to extract the local residual distribution characteristics and spatial neighborhood heterogeneity characteristics of the sensing map data of each microarray site in the parallel fitting subtask cluster, and to construct a local heterogeneity perception model. The computation module is used to calculate the adaptive weight coefficients of each microarray site in the dynamic fitting process through the local heterogeneity sensing model; based on the adaptive weight coefficients, a weighted local dynamic fitting objective function corresponding to each microarray site is constructed. The iterative module is used to iteratively solve the weighted local dynamic fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; and to obtain the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip. The correction module is used to construct a spatial neighborhood regularization constraint term based on physical spatial coordinate information; the initial binding rate constant and the initial dissociation rate constant are spatially smoothed and corrected by the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site.

[0048] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0049] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0050] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0051] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A parallel batch fitting method for dynamic parameters of high-density SPR microarrays, characterized in that, The method includes: Step 1: Obtain real-time sensing image data of each microarray site on the high-density SPR microarray chip; extract multi-dimensional time-series features from the real-time sensing image data to obtain the data complexity feature vector; calculate the computational load weight of each microarray site using the data complexity feature vector. Step 2: Based on the computational load weight, a dynamic task scheduling strategy is used to cluster and divide the real-time sensor map data to obtain multiple parallel fitting sub-task clusters with balanced computational load. Step 3: Extract the local residual distribution features and spatial neighborhood heterogeneity features of the sensing map data of each microarray site in the parallel fitting subtask cluster, and construct a local heterogeneity perception model. Step 4: Calculate the adaptive weight coefficients of each microarray site during the dynamic fitting process using the local heterogeneity sensing model; based on the adaptive weight coefficients, construct the weighted local dynamic fitting objective function corresponding to each microarray site. Step 5: Iteratively solve the weighted local dynamics fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; obtain the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip. Step 6: Construct a spatial neighborhood regularization constraint term based on physical spatial coordinate information; perform spatial smoothing correction on the initial binding rate constant and the initial dissociation rate constant through the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site.

2. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 1, characterized in that, Real-time sensing image data of each microarray site on a high-density SPR microarray chip is acquired, and multi-dimensional time-series feature extraction is performed on the real-time sensing image data to obtain a data complexity feature vector. The computational load weights for each microarray site are calculated using the data complexity eigenvector, including: Based on the surface plasmon resonance effect excited by the metal sensitive film of the high-density SPR microarray chip by incident light at a specific angle, the initial optical response signal is obtained by capturing the intensity attenuation signal of reflected light in real time through a photodetector array. Based on the spatial mapping relationship between the initial optical response signal and the microarray sites on the high-density SPR microarray chip, pixel-level position analysis is performed to obtain the original timing response data of each microarray site. Based on the original time-series response data of each microarray site, background subtraction and baseline calibration are performed using a preset reference channel signal to obtain real-time sensing data of each microarray site. Based on real-time sensor data, the signal fluctuation variance, baseline drift slope and residual convergence steps are extracted to obtain the data complexity feature vector. The computational load weights of each microarray site are obtained by weighted summation based on the signal fluctuation variance, baseline drift slope, and residual convergence steps in the data complexity feature vector.

3. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 2, characterized in that, Based on computational load weighting, a dynamic task scheduling strategy is employed to cluster and partition real-time sensor image data, resulting in multiple parallel fitting sub-task clusters with balanced computational load, including: Based on the calculated load weight of each microarray site, the total load weight and average load weight of all microarray sites are calculated to obtain the cluster load baseline value. Based on the cluster load baseline and calculated load weight, a distribution strategy is adopted to allocate real-time sensor map data to the initial task queue, resulting in multiple initial fitting sub-task clusters. The load deviation value is calculated based on the sum of the actual loads of each cluster in the initial fitting sub-task cluster and the cluster load benchmark value. Real-time sensor data is then migrated between the initial fitting sub-task clusters based on the load deviation value to obtain multiple parallel fitting sub-task clusters with balanced computational load.

4. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 3, characterized in that, Local residual distribution features and spatial neighborhood heterogeneity features of the sensing map data of each microarray site in the parallel fitting subtask cluster are extracted to construct a local heterogeneity sensing model, including: The time-series difference signal is calculated based on the real-time sensor data of each microarray site in the parallel fitting sub-task cluster to obtain the local residual time-series sequence; Skewness and kurtosis are calculated based on the local residual time series to obtain the local residual distribution characteristics; spatial neighborhood sites are determined based on the array index of each microarray site to obtain the neighborhood site set; The response signal difference degree is calculated based on the real-time sensor map data of spatial neighborhood sites in the neighborhood site set to obtain the spatial neighborhood heterogeneity characteristics. Based on the characteristics of local residual distribution and spatial neighborhood heterogeneity, a multilayer perceptron network containing a residual distribution coding layer and a spatial heterogeneity coding layer is constructed to obtain a local heterogeneity perception model.

5. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 4, characterized in that, Adaptive weight coefficients for each microarray site during the dynamic fitting process are calculated using a local heterogeneity sensing model. Based on these adaptive weight coefficients, a weighted local dynamic fitting objective function for each microarray site is constructed, including: Based on the local heterogeneity perception model, forward propagation calculations are performed on the local residual distribution characteristics and spatial neighborhood heterogeneity characteristics to obtain the adaptive weight coefficients of each microarray site in the dynamic fitting process. Based on the adaptive weighting coefficients, the response signal deviations at each time point in the real-time sensor image data are weighted to obtain a weighted local response deviation sequence. Based on the weighted local response deviation sequence, an error sum-of-squares function with the binding rate constant and dissociation rate constant as optimization variables is constructed to obtain the weighted local dynamic fitting objective function for each microarray site.

6. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 5, characterized in that, The weighted local dynamics fitting objective function is iteratively solved to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip is obtained, including: Based on the weighted local dynamics fitting objective function, initial guesses of the binding rate constant and dissociation rate constant are set, and the Levenberg-Marquardt algorithm is used for iterative optimization to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site. Based on the design drawings of the dotted mask of the high-density SPR microarray chip, the center pixel coordinates of each microarray site are extracted to obtain the set of center pixel coordinates; Based on the mapping ratio between the central pixel coordinate set and the chip physical size, the physical spatial coordinate information of each microarray point on the surface of the high-density SPR microarray chip is obtained.

7. The parallel dynamic parameter batch fitting method for high-density SPR microarrays according to claim 6, characterized in that, A spatial neighborhood regularization constraint term is constructed based on physical spatial coordinate information; the initial binding rate constant and initial dissociation rate constant are spatially smoothed and corrected using the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant for each microarray site, including: Based on the physical space coordinate information, the Euclidean distance between each microarray point and its spatial neighboring points is calculated to obtain the neighborhood distance matrix; a Gaussian weight matrix is ​​constructed based on the neighborhood distance matrix to obtain the spatial neighborhood regularization constraint term; Based on the spatial neighborhood regularization constraint, the parameter difference penalty value of the initial binding rate constant and initial dissociation rate constant of each microarray site relative to the spatial neighborhood site is calculated to obtain the spatial smoothing correction gradient. Based on the spatial smoothing correction gradient, the initial binding rate constant and the initial dissociation rate constant are updated by gradient descent to obtain the final target binding rate constant and target dissociation rate constant for each microarray site.

8. A parallel batch fitting system for dynamic parameters of high-density SPR microarrays, the system implementing the method as described in any one of claims 1 to 7, characterized in that, include: The acquisition module is used to acquire real-time sensing image data of each microarray site on the high-density SPR microarray chip, extract multi-dimensional time-series features from the real-time sensing image data to obtain a data complexity feature vector, and calculate the computational load weight of each microarray site through the data complexity feature vector. The partitioning module is used to cluster real-time sensor image data based on computational load weights and a dynamic task scheduling strategy to obtain multiple parallel fitting subtask clusters with balanced computational load. The module is used to extract the local residual distribution characteristics and spatial neighborhood heterogeneity characteristics of the sensing map data of each microarray site in the parallel fitting subtask cluster, and to construct a local heterogeneity perception model. The computation module is used to calculate the adaptive weight coefficients of each microarray site in the dynamic fitting process through the local heterogeneity sensing model; based on the adaptive weight coefficients, a weighted local dynamic fitting objective function corresponding to each microarray site is constructed. The iterative module is used to iteratively solve the weighted local dynamic fitting objective function to obtain the initial binding rate constant and initial dissociation rate constant of each microarray site; and to obtain the physical spatial coordinate information of each microarray site on the surface of the high-density SPR microarray chip. The correction module is used to construct a spatial neighborhood regularization constraint term based on physical spatial coordinate information; the initial binding rate constant and the initial dissociation rate constant are spatially smoothed and corrected by the spatial neighborhood regularization constraint term to obtain the final target binding rate constant and target dissociation rate constant of each microarray site.

9. A computing device, characterized in that, include: One or more processors; A storage device for storing one or more programs, which, when executed by one or more processors, cause the one or more processors to implement the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 7.

Images