Real-time monitoring of tumor microenvironment hypoxia level and radiotherapy regimen optimization system

Abstract

The present application relates to a kind of tumor microenvironment oxygen deficiency degree real-time monitoring and radiotherapy scheme optimization system, specifically, medical radiation therapy technical field, the scheme is by continuously collecting and processing low-resolution physiological image, real-time reconstruction high-definition three-dimensional oxygen deficiency dynamic distribution map in tumor, and then, utilize graph neural network technology to predict future short-term oxygen deficiency evolution trend and its uncertainty, based on this prediction and quantitative confidence, future oxygen deficiency target area is accurately mapped to treatment coordinate system, and the constraint strength of deformation registration is adaptively adjusted, finally, according to biological effect model with spatial confidence weight, real-time solution and output forward-looking, risk controllable beam adjustment instruction, drive radiotherapy equipment dynamic optimization irradiation dose projection, to realize the adaptive precision strike to the tumor oxygen deficiency area that changes in a twinkling, aims at improving the killing effect of radiotherapy to stubborn oxygen deficiency cell, while reducing the risk of mis-illumination to surrounding normal tissue.

Description

Technical Field

[0001] This invention relates to the field of medical radiotherapy technology, and more specifically, to a system for real-time monitoring of the hypoxia level of the tumor microenvironment and optimization of radiotherapy protocols. Background Technology

[0002] In clinical oncology radiotherapy practice, image-guided radiotherapy has become a core technology for precision radiotherapy. During treatment, the tumor target area is not static but undergoes complex deformations and displacements in response to the patient's physiological activities such as breathing, heartbeat, and gastrointestinal peristalsis. More importantly, the tumor microenvironment, especially the distribution and degree of hypoxic regions, exhibits significant spatiotemporal dynamic heterogeneity. This hypoxic state is not only affected by macroscopic physiological movements but is also closely related to the rapidly changing microscopic physiological processes within the tumor, such as blood perfusion and metabolic consumption. However, current mainstream image-guided techniques, such as cone-beam CT, primarily provide anatomical structural location information. While functional imaging techniques, such as positron emission tomography (PET) hypoxia imaging, can provide valuable hypoxia distribution maps before treatment, their image acquisition time is several minutes long, and there is an inherent time delay in the uptake and clearance of radioactive tracers. The obtained images are essentially averaged static snapshots of biological processes over tens of minutes, and cannot capture the dynamic evolution of hypoxic regions over seconds. Therefore, during each single treatment irradiation session that lasts several minutes, the clinic is actually irradiating a hypothesized, constant biological target area, which is fundamentally contradictory to the objective reality of real-time changes in hypoxic regions.

[0003] Existing technologies face significant bottlenecks in addressing the dynamic nature of tumor hypoxia. The core issue lies in the lack of a non-invasive, real-time monitoring and spatial mapping method for hypoxia levels that can be synchronized with radiotherapy equipment, achieving high spatiotemporal resolution. Firstly, at the monitoring level, existing technologies such as oxygen electrode measurements are invasive and can only provide single-point information. While techniques like functional magnetic resonance imaging (fMRI) based on blood oxygenation levels can non-invasively reflect changes in blood oxygenation, their temporal resolution is typically on the order of seconds, with limited spatial resolution and susceptibility to motion artifacts. Furthermore, they struggle to operate stably in the complex electromagnetic environment of treatment rooms and cannot achieve millisecond-level synchronous acquisition and communication with radiotherapy equipment such as medical linear accelerators at the hardware level. Secondly, at the information fusion and mapping level, even if dynamic hypoxia signals can be acquired, there is still the challenge of real-time and accurate registration of multi-source, asynchronous, and multimodal data with the treatment bed coordinate system. Existing image registration algorithms are mostly designed for rigid or deformable registration between fractions, which takes a long time to compute and cannot meet the needs of real-time tracking and mapping within fractions. This means that dynamic hypoxia information cannot be directly and accurately mapped into the dose calculation grid of the treatment plan, thus failing to provide reliable spatial guidance for real-time beam adjustment. Due to the lack of this core technology, radiotherapy protocol optimization loses the ability to target the most critical radiation resistance target area, namely the transient hypoxia region, in real time and with precision, which seriously restricts the further improvement of radiotherapy efficacy. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a system for real-time monitoring of hypoxia in the tumor microenvironment and optimization of radiotherapy plans. The system utilizes a hypoxia map reconstruction module, a dynamic prediction module, a coordinate mapping module, and a beam optimization module to solve the problems mentioned in the background.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: specifically including: a hypoxia map reconstruction module, a dynamic prediction module, a coordinate mapping module, and a beam optimization module, wherein;

[0006] Hypoxia Atlas Reconstruction Module: Receives low-resolution real-time hypoxia image sequences from medical imaging equipment (the medical imaging equipment is a fast plane wave excited blood oxygen level dependent magnetic resonance imaging sequence acquisition device integrated into the radiotherapy system. This device continuously acquires raw K-space data, which is reconstructed by fast Fourier transform and outputs a low-resolution real-time hypoxia image sequence). The low-resolution real-time hypoxia image sequence is then input into a preset spatiotemporal constrained generative adversarial network model. The spatiotemporal constrained generative adversarial network model is trained and inferred by calculating a comprehensive loss function including spatiotemporal gradient consistency loss, generating and outputting a high spatiotemporal resolution three-dimensional hypoxia atlas sequence.

[0007] Dynamic prediction module: Connected to the output of hypoxia map reconstruction module, it receives high spatiotemporal resolution 3D hypoxia map sequence, converts each time frame data of the high spatiotemporal resolution 3D hypoxia map sequence into spatiotemporal map structure data and inputs it into gated spatiotemporal map convolutional network for processing. The gated spatiotemporal map convolutional network performs forward propagation to output the predicted hypoxia mean map for a future set time point, and calculates the corresponding prediction uncertainty variance map based on multiple forward propagation sampling. Finally, it outputs a prediction result data packet containing the predicted hypoxia mean map and the prediction uncertainty variance map.

[0008] Coordinate mapping module: Connected to the output of the dynamic prediction module, it obtains the prediction result data packet, uses the predicted hypoxia mean map in the prediction result data packet as the floating image to be registered, and uses the pre-stored treatment plan CT image as the reference image to perform non-rigid deformation registration calculation. In the process of optimizing the deformation field, the non-rigid deformation registration calculation uses the prediction uncertainty variance map to weight the smoothness constraint of the deformation field, and calculates and generates the target coordinate mapping, final deformation field and registration confidence heatmap of the hypoxia area in the radiotherapy equipment coordinate system at future time.

[0009] Beam optimization module: Connected to the output of coordinate mapping module, it receives target coordinate mapping, final deformation field and registration confidence heatmap, constructs beam parameter optimization objective function based on confidence weighted bioequivalent dose model, solves and generates adjustment amounts for the weight of each beam subfield and the position of the multi-leaf grating blades in real time according to the optimization objective function, and outputs beam adjustment control commands corresponding to the adjustment amounts to the actuator of radiotherapy equipment.

[0010] In a preferred embodiment, the hypoxia map reconstruction module includes a spatiotemporally constrained generative adversarial network model comprising a generator network. The generator network receives a low-resolution real-time hypoxia image sequence as input and performs spatiotemporal feature encoding and fusion processing on the input low-resolution real-time hypoxia image sequence, specifically:

[0011] The generator network first uses a three-dimensional convolutional coding unit to perform convolution operations on a low-resolution real-time hypoxia image sequence that contains multiple consecutive time frames within a sliding time window centered on the current time point to be reconstructed, simultaneously in the spatial and slice dimensions, to extract the deep spatial features corresponding to each time frame within the sliding time window.

[0012] Subsequently, the generator network uses a one-dimensional temporal convolutional unit to convolve the sequence composed of deep spatial features of each time frame within the sliding time window along the time dimension, thereby extracting the temporal evolution features of the sliding time window.

[0013] Finally, the generator network fuses the deep spatial features and temporal evolution features of the current time point to be reconstructed to generate a latent feature map that integrates spatiotemporal information corresponding to the current time point to be reconstructed.

[0014] In a preferred embodiment, the spatiotemporal gradient consistency loss included in the comprehensive loss function is composed of a weighted sum of spatial gradient consistency components, temporal gradient consistency components, and cyclic consistency components; the weighted sum is obtained by multiplying the spatial gradient consistency component by a first weighting coefficient, the temporal gradient consistency component by a second weighting coefficient, and the cyclic consistency component by a third weighting coefficient, and then adding them together.

[0015] The calculation process of the spatial gradient consistency component is as follows: the generator network calculates the gradient magnitude distribution along the three-dimensional spatial direction of any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the latent feature map, and calculates the divergence between the gradient magnitude distribution and a prior spatial gradient magnitude distribution that conforms to the characteristics of biological tissue and is learned from a pre-acquired high-resolution gold standard image.

[0016] The calculation process of the temporal gradient consistency component is as follows: The difference between the maps corresponding to two adjacent time points in the high spatiotemporal resolution 3D hypoxia map sequence generated by the generator network is calculated to obtain the first difference image; the difference between the images corresponding to two time points in the low-resolution real-time hypoxia image sequence is calculated, and this difference is upsampled to the same size as the high spatiotemporal resolution 3D hypoxia map to obtain the second difference image; the gradient magnitudes of the first and second difference images along the time dimension are calculated respectively; finally, the first norm distance between the gradient magnitudes of the first and second difference images is calculated.

[0017] The calculation process of the cyclic consistency component is as follows: a preset three-dimensional Gaussian blur downsampling operator, used to simulate the imaging process of medical imaging equipment, is used to process any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the generator network to obtain a simulated low-resolution image; the first norm distance between the simulated low-resolution image and the image at the corresponding time point in the input low-resolution real-time hypoxia image sequence is calculated.

[0018] In a preferred embodiment, the spatiotemporal constrained generative adversarial network model is trained by optimizing the comprehensive loss function. After training, forward inference is performed on the input low-resolution real-time hypoxia image sequence. The forward inference process independently performs spatiotemporal feature encoding and fusion processing on each time point in the input low-resolution real-time hypoxia image sequence, and the decoder of the generator network upsamples and decodes the latent feature map to generate a high spatiotemporal resolution three-dimensional hypoxia map corresponding to the current time point, and finally outputs a high spatiotemporal resolution three-dimensional hypoxia map sequence.

[0019] In a preferred embodiment, the specific operation of converting each time frame data of the high spatiotemporal resolution three-dimensional hypoxia map sequence into spatiotemporal map structure data in the dynamic prediction module is as follows:

[0020] A1. Define each voxel in the three-dimensional hypoxia map of each time frame as an independent graph node, and construct a feature vector containing multi-dimensional information for each node. The feature vector consists of five components in sequence: the hypoxia value of the voxel at the current time point, the gradient component of the voxel along the X-axis in the three-dimensional coordinate system, the gradient component along the Y-axis, the gradient component along the Z-axis, and the approximate instantaneous rate of change obtained by calculating the difference between the hypoxia value of the voxel at the current time point and the immediately preceding time point.

[0021] A2. Based on the dual criteria, construct the connecting edges between graph nodes to form a spatiotemporal graph structure; the dual criteria include the spatial proximity criterion and the feature similarity criterion; according to the spatial proximity criterion, establish connecting edges between each node and a predetermined number of voxel nodes that are closest to its spatial coordinates in three-dimensional space; according to the feature similarity criterion, calculate the cosine similarity between the feature vectors of any two nodes; if the cosine similarity calculation result exceeds the preset first similarity threshold, then establish a connecting edge between the two nodes.

[0022] A3. In the time dimension, directed connection edges are established between nodes corresponding to voxels at different consecutive time points to represent the state evolution path of the voxel at the same spatial coordinate position.

[0023] By independently performing the above A1-A3 operations on each time frame in the high spatiotemporal resolution three-dimensional hypoxia map sequence, a corresponding spatiotemporal graph structure data containing the feature set of all nodes in that time frame and the set of connecting edges between nodes is generated for each time frame. These spatiotemporal graph structure data arranged in chronological order together constitute a dynamically evolving spatiotemporal graph data sequence that serves as the input to a gated spatiotemporal graph convolutional network.

[0024] In a preferred embodiment, the specific process by which the gated spatiotemporal graph convolutional network processes and calculates the prediction uncertainty variance map is as follows:

[0025] The gated spatiotemporal graph convolutional network consists of alternating stacked spatial graph convolutional layers and gated temporal convolutional layers. In the spatial graph convolutional layer, for the current central node to be processed, the processing procedure is as follows:

[0026] B1. Apply a shared learnable weight matrix to the feature vector of the central node and the feature vector of each of its neighboring nodes to obtain the transformation features of the central node and the transformation features of each neighboring node.

[0027] B2. Concatenate the transformation features of the center node with the transformation features of a neighbor node to form a concatenated feature vector; perform a dot product operation on the concatenated feature vector and a learnable attention vector, and then process the result through a leaky linear rectified activation function to obtain the original attention score of the neighbor node relative to the center node; repeat this process for all neighbor nodes to obtain a set of original attention scores.

[0028] B3. Normalize the original attention scores using an exponential function to convert each original score into a value between zero and one. These values ​​are the attention coefficients corresponding to each neighbor node.

[0029] B4. Multiply the transformation feature of each neighboring node by its corresponding attention coefficient, and sum the product of all neighboring nodes to obtain the spatial aggregation feature of the central node.

[0030] In the gated temporal convolutional layer, for each node, its spatial aggregated features output by the spatial graph convolutional layer at multiple consecutive historical time steps are arranged in chronological order to form the spatial aggregated feature sequence of that node; the processing procedure of this layer is as follows:

[0031] Two sets of independent learnable one-dimensional convolutional kernels are used to perform one-dimensional convolution operations on the input spatial aggregated feature sequence. The convolution result of the first set of convolutional kernels is applied with a hyperbolic tangent activation function to obtain the feature extraction component. The convolution result of the second set of convolutional kernels is applied with a sigmoid activation function to obtain the gated signal component. Finally, the feature extraction component and the gated signal component are multiplied element-wise to obtain the hidden state feature of the node at the current time step.

[0032] The multi-step prediction process is as follows: The gated spatiotemporal graph convolutional network unfolds in the time dimension, utilizing spatiotemporal graph data sequences from multiple consecutive time steps in the past. This data is processed through alternating stacked spatial graph convolutional layers and gated temporal convolutional layers to obtain the hidden state features of each node at the last input time step. Then, a linear transformation layer consisting of a learnable weight matrix and bias terms maps the hidden state features of each node to a predicted hypoxia value for a future set time step. This linear transformation layer constitutes the final prediction output layer of the gated spatiotemporal graph convolutional network. Finally, the predicted hypoxia values ​​of all nodes are rearranged according to their original spatial locations, thus forming a predicted hypoxia mean map for a future set time point.

[0033] The specific process of calculating the prediction uncertainty variance map based on multiple forward propagation sampling is as follows:

[0034] After training the gated spatiotemporal graph convolutional network, a random deactivation operation is kept active on the linear transformation weights of its spatial graph convolutional layers and the final prediction output layer. For the same set of historical spatiotemporal graph data sequences as input, multiple independent forward propagation calculations are performed. Due to the randomness of the random deactivation operation, each forward propagation yields a slightly different prediction result sample. Next, a discretization interval partitioning scheme covering the possible hypoxia value range is set, mapping continuous hypoxia prediction values ​​to corresponding intervals. For each voxel position in the predicted hypoxia mean map, the frequency distribution of its predicted value falling into each interval in multiple samplings is statistically analyzed. Finally, the prediction uncertainty at that position is quantified by calculating the statistical information entropy of this frequency distribution. The information entropy is specifically calculated as the sum of the following two parts:

[0035] The first part is to calculate the average frequency of each interval across all samples for all discretized intervals, then calculate the product of the average frequency and the natural logarithm of the average frequency, and finally sum the product of all intervals and take the negative value.

[0036] The second part involves calculating the information entropy of the frequency distribution in each independent sample, and then averaging the information entropy results of all multiple samples. The uncertainty values ​​of all voxel positions calculated in this way together constitute the prediction uncertainty variance map.

[0037] The final output of the dynamic prediction module is a prediction result data package containing a prediction mean hypoxia map and a prediction uncertainty variance map.

[0038] In a preferred embodiment, the specific process of performing non-rigid deformation registration calculation in the coordinate mapping module is as follows:

[0039] First, a dense spatial deformation mapping is defined from the coordinate system of the reference image to the coordinate system of the floating image. This mapping is described by a three-dimensional displacement field, and the three-dimensional displacement field is parameterized using a free deformation model based on cubic B-splines. Specifically, a three-dimensional control point grid with uniform spacing is defined on the reference image, and each control point in the three-dimensional control point grid is associated with a three-dimensional displacement vector. The displacement of any point in space is obtained by weighted interpolation of the displacements of its neighboring control points by a predetermined number of control points using cubic B-spline basis functions. The optimization variable is the three-dimensional displacement vector of all control points in the three-dimensional control point grid.

[0040] Secondly, with the goal of maximizing the local structural correlation between the deformed floating image and the reference image, the image similarity metric loss is calculated. The image similarity metric uses local cross-correlation, and its calculation process is as follows: In the reference image, a local neighborhood is defined with each voxel as the center, and the mean intensity of all voxels in the reference image and the mean intensity of all voxels in the deformed floating image within this neighborhood are calculated. Then, the deviation between the intensity of the reference image and its mean, and the deviation between the intensity of the deformed floating image and its mean, are calculated at each position within this neighborhood. Next, the square of the product of these two deviations is calculated, and then the sum of the squares of these two deviations is calculated. Finally, the square of the deviation product is divided by the product of the sums of the squares of the two deviations to obtain the local correlation ratio at the central voxel. The local correlation ratios of all voxels in the entire reference image are summed, and the sum is negative, which is taken as the total image similarity metric loss.

[0041] Next, a deformation field regularization constraint term weighted by uncertainty information is constructed. This term is the integral result of the deformation field over the spatial region covered by the reference image. Its integrand consists of a weighted sum of three parts, each of which is related to the local geometric properties of the deformation field:

[0042] The first part is used to penalize the degree to which the local Jacobian matrix of the deformation field deviates from the orthogonal transformation; the local Jacobian matrix represents the spatial gradient of the deformation field at a point, and is composed of the identity matrix plus the gradient matrix of the displacement field at that point;

[0043] The second part is used to penalize local volume changes caused by deformation; local volume changes are characterized by the determinant value of the local Jacobian matrix.

[0044] The third part, the uncertainty weighted smoothing term, is calculated as follows: First, the predicted uncertainty variance map is resampled to the same coordinate system as the reference image through the inverse mapping of the current deformation field, resulting in a resampled uncertainty map. Then, an uncertainty gating function is defined, which maps the uncertainty value at each position in the resampled uncertainty map to a weight between zero and one. Specifically, a positive steepness coefficient and a predefined uncertainty threshold are set. The uncertainty value at the current position is subtracted from the threshold and then multiplied by the steepness coefficient. The resulting product is input into the negative exponent of the natural exponential function for calculation. The result is then incremented by one and the reciprocal is taken to obtain the weight at that position. Finally, the Frobenius norm of the displacement field gradient matrix is ​​calculated, and the square of this Frobenius norm is weighted using the weights.

[0045] The integral result of the deformation field regularization constraint term is the cumulative amount of the weighted sum of the three parts over the spatial region;

[0046] Finally, the optimization objective function for non-rigid deformation registration calculation is composed of the image similarity measurement loss and the deformation field regularization constraint term weighted by uncertainty. The optimal control point displacement vector is solved by minimizing this objective function, and then the final deformation field is obtained.

[0047] In a preferred embodiment, the specific process of calculating and generating the target coordinate mapping, final deformation field, and registration confidence heatmap of the hypoxic region in the radiotherapy equipment coordinate system at future times is as follows:

[0048] The final deformation field is solved using an optimization strategy combining multi-resolution hierarchical optimization and the Gauss-Newton trust-region method. First, the deformation field is initialized and optimized on a coarser-resolution 3D control point mesh to obtain a low-resolution deformation field estimate. Then, the resolution of the control point mesh is progressively refined, and the deformation field obtained at the current coarser resolution is used as the initial value to continue optimization at a finer resolution until the preset maximum resolution is reached, resulting in a high-precision final deformation field. In each single optimization iteration at each resolution level, the Gauss-Newton method is used to approximate the Hessian matrix of the objective function with respect to the optimization variables, and a local quadratic model is constructed. The optimal update step size for the current iteration is determined by solving a trust-region subproblem.

[0049] After obtaining the final deformation field, its inverse mapping is calculated using numerical methods. The spatial coordinates of all voxels in the hypoxic region defined by the hypoxic mean map predicted at future time in the floating image are transformed to the coordinate system of the reference image, i.e. the coordinate system of the radiotherapy equipment, through inverse mapping, thereby obtaining the target coordinate mapping.

[0050] The registration confidence heatmap is generated by comprehensively evaluating the reliability of the registration results at each spatial location in the reference image; specifically, for each location in the reference coordinate system, its confidence is the product of the following three components:

[0051] The first component is the prediction confidence, which is calculated based on the value at that position in the resampled uncertainty map, specifically by using a negative natural exponential function with the uncertainty value as the exponent.

[0052] The second component is the deformation rationality, which is calculated based on the determinant of the local Jacobian matrix of the final deformation field at that position. Specifically, it is calculated by a negative natural exponential function with the absolute value of the determinant minus one as the exponent.

[0053] The third component is the local matching degree, calculated based on the image similarity metric of the local neighborhood at that location. The calculation process for this component is as follows: First, following the calculation process of local cross-correlation, calculate the image similarity metric value of the local neighborhood at that location, i.e., the local correlation ratio; then, take the negative of this similarity metric value to obtain the local similarity loss value at that location; next, find the maximum value of the local similarity loss value among all locations in the reference image, called the global maximum local loss value; finally, divide the local similarity loss value at that location by the global maximum local loss value to obtain a ratio, and then subtract this ratio from one to get the result, which is the local matching degree component.

[0054] By multiplying these three components, a registration confidence heatmap with values ​​ranging from zero to one is finally generated.

[0055] In a preferred embodiment, the specific process of constructing the objective function for beam parameter optimization based on the confidence-weighted bioequivalent dose model in the beam optimization module is as follows:

[0056] First, the time-integrated physical dose that will be received at each spatial location in the coordinate system of the radiotherapy equipment within the future optimized time window is calculated. The calculation of this time-integrated physical dose is based on the subfield dose rate weights of all subfields of all beams within the optimized time window, and the unit dose rate distribution determined by the change of the position of the multi-leaf grating blades over time. The contribution of all subfields of all beams within the optimized time window is integrated and summed.

[0057] Secondly, a spatiotemporally dependent biological effect model is constructed to calculate the bioequivalent dose received by each voxel in the target hypoxic region within a future optimized time window; the target hypoxic region is defined by a target coordinate mapping; specifically, the model includes two effective radiosensitivity parameters related to the hypoxia level, namely the alpha effective parameter and the beta effective parameter; the alpha effective parameter is obtained by averaging the reciprocal of the oxygen enhancement ratio of the voxel at each moment within the optimized time window over time;

[0058] The beta effective parameter is obtained by averaging the time-inverse of the square of the oxygen enhancement ratio of the voxel at each moment within the optimized time window; the specific formula for calculating the bioequivalent dose is: the alpha effective parameter multiplied by the time-integrated physical dose, plus the beta effective parameter multiplied by a treatment time-related intra-session repair factor, and then multiplied by the square of the time-integrated physical dose.

[0059] Finally, a confidence-weighted prescription gap function is constructed, defining the prescription biological effect as the clinically required target biological effect value; for each voxel in the target hypoxic region, its prescription gap function consists of two terms:

[0060] The first item is the penalty for failing to meet the high confidence region standard. Its value is the value of the registration confidence heatmap corresponding to the voxel, multiplied by the square of the difference between the prescription biological effect of the voxel and its actual calculated bioequivalent dose.

[0061] The second term is the over-irradiation penalty term for low-confidence areas, which is a preset risk avoidance coefficient multiplied by a factor minus the value of the registration confidence heatmap, and then multiplied by the square of the actual calculated bioequivalent dose of the voxel; the beam parameter optimization objective function is the sum of the prescription difference functions of all voxels in the target hypoxic region.

[0062] In a preferred embodiment, the specific process of solving and generating the adjustment amounts for the weights of each beam subfield and the positions of the multi-leaf grating blades in real time according to the optimization objective function is as follows:

[0063] The real-time solution problem of the beam parameter optimization objective function is modeled as a stochastic optimal control problem in the finite-time domain. The state variables of this problem include the cumulative machine hops of each beam, the position and velocity of the multi-leaf grating blades, and the gantry angle. The control variables include the dose rate change rate of each beam subfield and the target velocity command of the multi-leaf grating blades. The optimization objective is to minimize an expected cumulative loss while satisfying the physical and kinematic constraints of the radiotherapy equipment. This expected cumulative loss consists of the sum of the integrals of three parts within the optimization time window:

[0064] The first part is the sum of confidence-weighted prescription gap function values ​​for all voxels within the target hypoxic region;

[0065] The second part is the sum of dose penalty function values ​​for all voxels in the regions at risk of organ failure;

[0066] The third part is a control variable smoothness regularization function value multiplied by a regularization weight coefficient;

[0067] A model predictive control framework is employed, combining differential dynamic programming algorithms and unscented transform techniques to solve stochastic optimal control problems in real-time rolling time domains. During each control cycle, the following operations are performed:

[0068] Based on the current status of the radiotherapy equipment and the latest received predicted mean hypoxia map, predicted uncertainty variance map, and registration confidence heatmap, the stochastic optimal control problem is reconstructed and solved in a shortened prediction time domain; the finite-time deterministic optimal control problem is solved using a differential dynamic programming algorithm to obtain the optimal control command sequence for the current control cycle.

[0069] From the optimal control command sequence obtained from the solution, the command of the first control cycle is extracted as the adjustment amount for the weight of the beam subfield and the position of the multileaf grating blades in the current control cycle. The adjustment amount is input into a forward simulator for safety verification. If the verification is successful, the beam adjustment control command packet containing the timestamp, beam identifier, subfield weight adjustment amount, and multileaf grating blade target position sequence is sent to the actuator of the radiotherapy equipment. If the verification fails, a safety fallback strategy is triggered. The specific content of the safety fallback strategy is as follows: while keeping the current beam parameters unchanged, a pre-calculated low-risk static optimization plan based on historical average predicted hypoxia distribution is immediately invoked to generate corresponding backup control commands and issue a warning. Finally, the command packet that has passed the safety verification is output to the actuator of the radiotherapy equipment through a low-latency communication interface. The actuator includes a dose rate controller, a multileaf grating driver, and a gantry controller.

[0070] The beneficial effects of this invention are as follows: This scheme continuously acquires and processes low-resolution physiological images to reconstruct a high-definition three-dimensional dynamic distribution map of hypoxia within the tumor in real time. Then, it uses graph neural network technology to predict the hypoxia evolution trend and its uncertainty in the near future. Based on this prediction and quantified confidence, the future hypoxic target area is accurately mapped to the treatment coordinate system, and the constraint strength of deformation registration is adaptively adjusted. Finally, based on the biological effect model with spatial confidence weights, a forward-looking and risk-controllable beam adjustment command is solved and output in real time, driving the radiotherapy equipment to dynamically optimize the irradiation dose, thereby achieving adaptive and precise targeting of the rapidly changing hypoxic tumor area. This aims to improve the killing effect of radiotherapy on refractory hypoxic cells while reducing the risk of false irradiation of surrounding normal tissues. Attached Figure Description

[0071] Figure 1 This is a flowchart of the method of the present invention;

[0072] Figure 2 This is a block diagram of the system structure of the present invention. Detailed Implementation

[0073] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0074] In the description of this application, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more features. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0075] In the description of this application, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this application is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed in this application.

[0076] Example 1

[0077] This embodiment provides, for example Figure 1 The method for real-time monitoring of hypoxia in the tumor microenvironment and optimization of radiotherapy protocols, as shown, includes the following steps:

[0078] Step S1: Receive a low-resolution real-time hypoxia image sequence from a medical imaging device and input the low-resolution real-time hypoxia image sequence into a preset spatiotemporal constraint generative adversarial network model. The spatiotemporal constraint generative adversarial network model is trained and inferred by calculating a comprehensive loss function including spatiotemporal gradient consistency loss, and generates and outputs a high spatiotemporal resolution three-dimensional hypoxia map sequence.

[0079] Step S2: Connect to the output of the hypoxia map reconstruction module, receive the high spatiotemporal resolution three-dimensional hypoxia map sequence, convert each time frame data of the high spatiotemporal resolution three-dimensional hypoxia map sequence into spatiotemporal map structure data and input it into the gated spatiotemporal map convolutional network for processing. The gated spatiotemporal map convolutional network performs forward propagation to output the predicted hypoxia mean map for a future set time point, and calculates the corresponding prediction uncertainty variance map based on multiple forward propagation sampling. Finally, it outputs a prediction result data packet containing the predicted hypoxia mean map and the prediction uncertainty variance map.

[0080] Step S3: Connect to the output of the dynamic prediction module, obtain the prediction result data packet, use the predicted hypoxia mean map in the prediction result data packet as the floating image to be registered, use the pre-stored treatment plan CT image as the reference image, and perform non-rigid deformation registration calculation. In the process of optimizing the deformation field, the non-rigid deformation registration calculation uses the prediction uncertainty variance map to weight the smoothness constraint of the deformation field, and calculates and generates the target coordinate mapping of the hypoxia area in the radiotherapy equipment coordinate system at future time, the final deformation field and the registration confidence heatmap.

[0081] Step S4: Connect to the output of the coordinate mapping module, receive the target coordinate mapping, the final deformation field and the registration confidence heatmap, construct the beam parameter optimization objective function based on the confidence weighted biological equivalent dose model, solve and generate the weight of each beam subfield and the adjustment amount of the multi-leaf grating blade position in real time according to the optimization objective function, and output the beam adjustment control command corresponding to the adjustment amount to the actuator of the radiotherapy equipment.

[0082] It also provides, for example Figure 2 The system shown is a real-time monitoring system for the degree of hypoxia in the tumor microenvironment and an optimization system for radiotherapy plans. Specifically, it includes: a hypoxia atlas reconstruction module, a dynamic prediction module, a coordinate mapping module, and a beam optimization module.

[0083] Hypoxia Atlas Reconstruction Module: Receives low-resolution real-time hypoxia image sequences from medical imaging equipment (the medical imaging equipment is a fast plane wave excited blood oxygen level dependent magnetic resonance imaging sequence acquisition device integrated into the radiotherapy system. This device continuously acquires raw K-space data, which is reconstructed by fast Fourier transform and outputs a low-resolution real-time hypoxia image sequence). The low-resolution real-time hypoxia image sequence is then input into a preset spatiotemporal constrained generative adversarial network model. The spatiotemporal constrained generative adversarial network model is trained and inferred by calculating a comprehensive loss function including spatiotemporal gradient consistency loss, generating and outputting a high spatiotemporal resolution three-dimensional hypoxia atlas sequence.

[0084] Dynamic prediction module: Connected to the output of hypoxia map reconstruction module, it receives high spatiotemporal resolution 3D hypoxia map sequence, converts each time frame data of the high spatiotemporal resolution 3D hypoxia map sequence into spatiotemporal map structure data and inputs it into gated spatiotemporal map convolutional network for processing. The gated spatiotemporal map convolutional network performs forward propagation to output the predicted hypoxia mean map for a future set time point, and calculates the corresponding prediction uncertainty variance map based on multiple forward propagation sampling. Finally, it outputs a prediction result data packet containing the predicted hypoxia mean map and the prediction uncertainty variance map.

[0085] Coordinate mapping module: Connected to the output of the dynamic prediction module, it obtains the prediction result data packet, uses the predicted hypoxia mean map in the prediction result data packet as the floating image to be registered, and uses the pre-stored treatment plan CT image as the reference image to perform non-rigid deformation registration calculation. In the process of optimizing the deformation field, the non-rigid deformation registration calculation uses the prediction uncertainty variance map to weight the smoothness constraint of the deformation field, and calculates and generates the target coordinate mapping, final deformation field and registration confidence heatmap of the hypoxia area in the radiotherapy equipment coordinate system at future time.

[0086] Beam optimization module: Connected to the output of coordinate mapping module, it receives target coordinate mapping, final deformation field and registration confidence heatmap, constructs beam parameter optimization objective function based on confidence weighted bioequivalent dose model, solves and generates adjustment amounts for the weights of each beam subfield and the position of the multi-leaf grating blades in real time according to the optimization objective function, and outputs beam adjustment control commands corresponding to the adjustment amounts to the actuator of radiotherapy equipment.

[0087] In this embodiment, it should be specifically noted that in the hypoxia map reconstruction module, the spatiotemporal constrained generative adversarial network model includes a generator network. The generator network receives a low-resolution real-time hypoxia image sequence as input and performs spatiotemporal feature encoding and fusion processing on the input low-resolution real-time hypoxia image sequence, specifically as follows:

[0088] The generator network first uses a 3D convolutional coding unit to perform convolution operations on a low-resolution real-time hypoxia image sequence within a sliding time window containing multiple consecutive time frames centered on the current time point to be reconstructed, simultaneously in the spatial and slice dimensions, to extract the deep spatial features corresponding to each time frame within the sliding time window; the size of the sliding time window can be set according to the image acquisition frequency and physiological movement cycle, for example, a window containing 5 to 7 consecutive time frames, to ensure that it can cover a complete respiratory or blood flow fluctuation cycle information;

[0089] Subsequently, the generator network uses a one-dimensional temporal convolutional unit to convolve the sequence composed of deep spatial features of each time frame within the sliding time window along the time dimension, and extracts the temporal evolution features of the sliding time window; the kernel size and stride settings of the one-dimensional temporal convolutional unit must be set to ensure that they can effectively capture the dynamic change patterns of low-resolution real-time hypoxia image sequences in the time dimension.

[0090] Finally, the generator network fuses the deep spatial features and temporal evolution features of the current time point to be reconstructed to generate a latent feature map corresponding to the current time point to be reconstructed, which incorporates spatiotemporal information. The fusion operation includes, but is not limited to, the concatenation operation of feature tensors or the element-wise addition operation. The fused latent feature map is input to the decoder part of the generator network for upsampling and reconstruction.

[0091] The spatiotemporal gradient consistency loss included in the comprehensive loss function is composed of a weighted sum of the spatial gradient consistency component, the temporal gradient consistency component, and the cyclic consistency component. The weighted sum is obtained by multiplying the spatial gradient consistency component by a first weight coefficient, the temporal gradient consistency component by a second weight coefficient, and the cyclic consistency component by a third weight coefficient, and then adding them together.

[0092] The calculation process of the spatial gradient consistency component is as follows: The generator network calculates the gradient magnitude distribution along the three-dimensional spatial direction of any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the latent feature map, and calculates the divergence between this gradient magnitude distribution and a prior spatial gradient magnitude distribution that conforms to the characteristics of biological tissue and is learned from a pre-acquired high-resolution gold standard image; the divergence is preferably KL divergence, and the high-resolution gold standard image is obtained by registering high-resolution hypoxia PET-CT images and treatment plan CT images before treatment;

[0093] The calculation process of the temporal gradient consistency component is as follows: The difference between the maps corresponding to two adjacent time points in the high spatiotemporal resolution 3D hypoxia map sequence generated by the generator network is calculated to obtain the first difference image; the difference between the images corresponding to two time points in the low-resolution real-time hypoxia image sequence is calculated, and this difference is upsampled to the same size as the high spatiotemporal resolution 3D hypoxia map to obtain the second difference image; the gradient magnitudes of the first and second difference images along the time dimension are calculated respectively; finally, the first norm distance between the gradient magnitudes of the first and second difference images is calculated; the first norm distance, i.e., the L1 norm, is used to constrain the reconstructed high-resolution sequence to maintain consistency with the physiological changes reflected by the original low-resolution signal in terms of the degree of temporal variation, avoiding the introduction of unreasonable temporal jitter;

[0094] The calculation process of the cycle consistency component is as follows: a preset three-dimensional Gaussian blur downsampling operator, used to simulate the imaging process of medical imaging equipment, is used to process any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the generator network to obtain a simulated low-resolution image; the first norm distance between the simulated low-resolution image and the image at the corresponding time point in the input low-resolution real-time hypoxia image sequence is calculated; this process constitutes a self-consistent cycle, which forces the high-resolution image generated by the network to be as close as possible to the original observation data after the simulated imaging degradation, thereby ensuring the fidelity of the reconstruction result;

[0095] The spatiotemporal constrained generative adversarial network (GAN) model is trained by optimizing the comprehensive loss function. After training, it performs forward inference on the input low-resolution real-time hypoxia image sequence. In the forward inference process, spatiotemporal feature encoding and fusion are performed independently for each time point in the input low-resolution real-time hypoxia image sequence. The decoder of the generator network upsamples and decodes the latent feature map to generate a high spatiotemporal resolution 3D hypoxia map corresponding to the current time point, and finally outputs a high spatiotemporal resolution 3D hypoxia map sequence. Each map in the output high spatiotemporal resolution 3D hypoxia map sequence corresponds one-to-one with each frame of the input low-resolution real-time hypoxia image sequence at the time point. Furthermore, this high spatiotemporal resolution 3D hypoxia map sequence serves as the sole input data source for the subsequent dynamic prediction module, providing unified and continuous high-quality hypoxia dynamic distribution information for the subsequent dynamic prediction, coordinate mapping, and beam optimization of the entire system.

[0096] In this embodiment, it is specifically necessary to explain the specific operation in the dynamic prediction module of converting each time frame data of the high spatiotemporal resolution three-dimensional hypoxia map sequence into spatiotemporal map structure data as follows:

[0097] A1. Define each voxel in the 3D hypoxia map of each time frame as an independent graph node, and construct a feature vector containing multidimensional information for each node. The feature vector consists of five components in sequence: the hypoxia value of the voxel at the current time point, the gradient component of the voxel along the X-axis, the gradient component along the Y-axis, the gradient component along the Z-axis in the 3D spatial coordinate system, and the approximate instantaneous rate of change obtained by calculating the difference between the hypoxia value of the voxel at the current time point and the immediately preceding time point. The spatial gradient component is approximated by calculating the difference between the hypoxia value of the voxel and the adjacent voxels in the corresponding direction. The instantaneous rate of change approximates the short-term dynamics of the hypoxia state at this location, together providing key information for subsequent graph learning to describe local spatial heterogeneity and temporal evolution.

[0098] A2. Based on a dual criterion, connection edges are constructed between graph nodes to form a spatiotemporal graph structure. The dual criterion includes a spatial proximity criterion and a feature similarity criterion. According to the spatial proximity criterion, a connection edge is established between each node and a predetermined number of voxel nodes that are closest to it in three-dimensional space. The predetermined number is, for example, twenty-six corresponding to the twenty-six neighborhood connection method. This method ensures that each node and its directly physically adjacent voxels are structurally connected, encoding the basic anatomical adjacency relationship. According to the feature similarity criterion, the cosine similarity between the feature vectors of any two nodes is calculated. If the calculated cosine similarity exceeds a preset first similarity threshold (for example, the threshold can be set between 0.7 and 0.9, preferably 0.8), then a connection edge is established between the two nodes. This criterion allows functionally similar but spatially separable hypoxic regions (such as cell clusters located in different parts of the tumor but with similar degrees of hypoxia) to be connected in the graph, thereby modeling non-local, functional long-range associations, which is crucial for capturing the overall evolution of hypoxic patterns.

[0099] A3. In the time dimension, directed connections are established between nodes corresponding to voxels at different consecutive time points to represent the same spatial coordinate position, so as to represent the state evolution path of that position itself; this edge clearly describes the causal or continuous change trajectory of the hypoxia state of each specific spatial point over time, which is the basis for learning temporal dynamics.

[0100] By independently performing the above A1-A3 operations on each time frame in the high spatiotemporal resolution three-dimensional hypoxia map sequence, a spatiotemporal graph structure data corresponding to each time frame is generated, which includes the feature set of all nodes in that time frame and the set of connecting edges between nodes. These spatiotemporal graph structure data arranged in chronological order together constitute a dynamically evolving spatiotemporal graph data sequence that serves as the input to a gated spatiotemporal graph convolutional network.

[0101] The specific process of using a gated spatiotemporal graph convolutional network to process and calculate the uncertainty variance map is as follows:

[0102] The gated spatiotemporal graph convolutional network consists of alternating stacked spatial graph convolutional layers and gated temporal convolutional layers. In the spatial graph convolutional layer, for the current central node to be processed, the processing procedure is as follows:

[0103] B1. Apply a shared learnable weight matrix to the feature vector of the central node and the feature vector of each of its neighboring nodes to obtain the transformation features of the central node and the transformation features of each neighboring node.

[0104] B2. Concatenate the transformation features of the center node with the transformation features of a neighbor node to form a concatenated feature vector; perform a dot product operation on the concatenated feature vector and a learnable attention vector, and then process the result through a leaky linear rectified activation function to obtain the original attention score of the neighbor node relative to the center node; repeat this process for all neighbor nodes to obtain a set of original attention scores; the attention vector is a parameter learned during network training, and its function is to evaluate the importance of different neighbor node features to the center node;

[0105] B3. Normalize the original attention scores using an exponential function, converting each original score into a value between zero and one, and ensuring that the sum of the values ​​of all neighboring nodes after the conversion is one. These values ​​are the attention coefficients corresponding to each neighboring node. This step realizes the normalization and probabilization of the weights, so that the aggregation operation has a clear weighted summation interpretation.

[0106] B4. Multiply the transformation feature of each neighboring node by its corresponding attention coefficient, and sum the product of all neighboring nodes to obtain the spatial aggregation feature of the central node, thereby capturing the non-local and non-uniform spatial interaction patterns in the graph. Through this adaptive attention mechanism, the network can dynamically and discriminatively integrate information from spatially neighboring nodes and feature-similar nodes, rather than averaging with fixed weights. This significantly improves the ability to model the interaction relationships of complex hypoxic spatial patterns.

[0107] In the gated temporal convolutional layer, for each node, its spatial aggregated features output by the spatial graph convolutional layer at multiple consecutive historical time steps are arranged in chronological order to form the spatial aggregated feature sequence of that node; the processing procedure of this layer is as follows:

[0108] Two independent sets of learnable one-dimensional convolutional kernels are used to perform one-dimensional convolution operations on the input spatial aggregated feature sequence. The convolution result of the first set of kernels is subjected to a hyperbolic tangent activation function to obtain the feature extraction component. The convolution result of the second set of kernels is subjected to a sigmoid activation function to obtain the gated signal component. Finally, the feature extraction component and the gated signal component are multiplied element-wise to obtain the hidden state feature of the node at the current time step. The value range of the gated signal component is between 0 and 1. Its function is similar to a filter, controlling how much historical feature information can flow into the current state, thereby effectively learning long-term temporal dependencies and alleviating the gradient vanishing problem.

[0109] The multi-step prediction process is as follows: The gated spatiotemporal graph convolutional network unfolds in the time dimension, utilizing spatiotemporal graph data sequences from multiple consecutive time steps in the past. This data is processed through alternating stacked spatial graph convolutional layers and gated temporal convolutional layers to obtain the hidden state features of each node at the last input time step. Then, a linear transformation layer composed of a learnable weight matrix and bias terms maps the hidden state features of each node to a predicted hypoxia value for a future set time step. This linear transformation layer constitutes the final prediction output layer of the gated spatiotemporal graph convolutional network. The predicted hypoxia values ​​of all nodes are rearranged according to their original spatial positions, thus forming a predicted hypoxia mean map for the future set time point. The prediction time span for the future set time step can be set according to the system response time and physiological change rate of the treatment device, for example, predicting the state within the next 0.5 to 2 seconds to provide a prospective measure for beam adjustment.

[0110] The specific process of calculating the prediction uncertainty variance map based on multiple forward propagation sampling is as follows:

[0111] After training the gated spatiotemporal graph convolutional network, the random deactivation operation is kept active on the linear transformation weights of its spatial graph convolutional layers and the final prediction output layer. "Keeping active" here means that during model inference (prediction), the random deactivation function of these layers is not turned off, but continues to randomly discard some neuron connections with the same deactivation probability as during training. This is crucial for achieving Monte Carlo sampling. For the same set of historical spatiotemporal graph data sequences as input, multiple independent forward propagation calculations are performed (e.g., 50 forward propagations). Due to the randomness of the random deactivation operation, each forward propagation yields a slightly different prediction result sample. Next, a discretization interval division scheme covering the possible hypoxia value range is set (e.g., dividing the hypoxia value range into 10 equal intervals). Continuous hypoxia prediction values ​​are mapped to the corresponding intervals. For each voxel position in the predicted hypoxia mean map, the frequency distribution of its predicted value falling into each interval in multiple samplings is statistically analyzed. Finally, the prediction uncertainty at that position is quantified by calculating the statistical information entropy of this frequency distribution. The information entropy is specifically calculated as the sum of the following two parts:

[0112] The first part is to calculate the average frequency of each interval across all samples for all discretized intervals, then calculate the product of the average frequency and the natural logarithm of the average frequency, and finally sum the product of all intervals and take the negative value. This part reflects the uncertainty of the average distribution of multiple predictions, i.e., "cognitive uncertainty", which mainly stems from the uncertainty of the model itself regarding the weights.

[0113] The second part involves calculating the information entropy of the frequency distribution in each independent sample (i.e., the sum of the products of the frequencies of each interval and the natural logarithms of those frequencies, then taking the negative value), and then calculating the average of this information entropy result for all multiple samples. This part reflects the average of the uncertainty of the distribution of a single prediction, i.e., "random uncertainty," which mainly originates from data noise. The uncertainty values ​​calculated for all voxel positions using this method together constitute the prediction uncertainty variance map. This prediction uncertainty variance not only indicates the reliability of the prediction, but its spatial distribution information can also be used to guide downstream modules in making confidence-weighted decisions.

[0114] The final output of the dynamic prediction module is a prediction result data package containing a prediction mean hypoxia map and a prediction uncertainty variance map.

[0115] In this embodiment, it is necessary to specifically explain the process of performing non-rigid deformation registration calculation in the coordinate mapping module as follows:

[0116] First, a dense spatial deformation mapping is defined from the coordinate system of the reference image to the coordinate system of the floating image. This mapping is described by a three-dimensional displacement field, which is parameterized using a free deformation model based on cubic B-splines. Specifically, a three-dimensional control point grid with uniform spacing is defined on the reference image, and each control point in the grid is associated with a three-dimensional displacement vector. The displacement of any point in space is obtained by weighted interpolation of the displacements of its neighboring control points using a predetermined number of cubic B-spline basis functions. The optimization variable is the three-dimensional displacement vector of all control points in the three-dimensional control point grid. The predetermined number of neighboring control points specifically refers to taking four consecutive control points along each of the three spatial dimensions in the smallest grid cell containing that point, for a total of sixty-four control points. The cubic B-spline basis functions are cubic polynomial spline functions. This ensures that the generated displacement field has second-order continuity, thus producing smooth deformation. Second-order continuity means that the displacement field and its first derivative (slope) and second derivative (curvature) are continuous in space, which effectively avoids unnatural folds or cracks in the deformation field and ensures the physical rationality of the mapping. Weighted interpolation refers to, for a point in space, first determining the grid cell in which it is located and its normalized coordinates within the cell, then using cubic B-spline basis functions to calculate the weight of each of the above sixty-four control points under the normalized coordinates, and finally multiplying the displacement vectors of all control points by their corresponding weights and summing them to obtain the displacement of that point. The basis functions have local support, that is, each control point only affects the deformation in its finite neighborhood, which makes the deformation adjustment local and facilitates fine optimization of specific regions.

[0117] Secondly, driven by maximizing the local structural correlation between the deformed floating image and the reference image, an image similarity metric loss is calculated. The image similarity metric employs local cross-correlation, and its calculation process is as follows: In the reference image, a local neighborhood is defined centered on each voxel point. The mean intensity of all voxels in the reference image and the mean intensity of all voxels in the deformed floating image within this neighborhood are calculated. Then, the deviations between the intensity of the reference image and its mean, and the deviations between the intensity of the deformed floating image and its mean, are calculated at each location within this neighborhood. Next, the square of the product of these two deviations is calculated, followed by the sum of the squares of each deviation. Finally, the square of the deviation product is divided by the product of the sums of the squares of the two deviations. The local correlation ratio at the central voxel point is obtained. The size of the local neighborhood can be set according to the image resolution and the desired registration locality. For example, a cube neighborhood with a side length of five voxels centered at the voxel point can be used. This neighborhood size balances the sensitivity to local structural changes and the robustness to image noise. If the neighborhood is too small, noise interference will be introduced. If the neighborhood is too large, local details will be blurred, affecting the registration accuracy. The local correlation ratios of all voxel points on the entire reference image are summed, and the sum is negative, which is used as the total image similarity loss metric. Local cross-correlation is not sensitive to linear changes in image intensity, so it is suitable for registering functional hypoxia maps with CT images of anatomical structures.

[0118] Next, a deformation field regularization constraint term weighted by uncertainty information is constructed. This term is the integral result of the deformation field over the spatial region covered by the reference image. Its integrand consists of a weighted sum of three parts, each of which is related to the local geometric properties of the deformation field:

[0119] The first part penalizes the degree to which the local Jacobian matrix of the deformation field deviates from the orthogonal transformation. The local Jacobian matrix represents the spatial gradient of the deformation field at a point and is composed of the identity matrix plus the gradient matrix of the displacement field at that point. The calculation of this part involves weighting the square of the trace of the matrix obtained by multiplying the transpose of the local Jacobian matrix by itself and then subtracting the identity matrix, in order to constrain excessive shear and stretch deformation. This penalty term encourages the deformation field to be approximated locally as a rigid rotation, which helps to maintain the relative shape of the structure.

[0120] The second part penalizes local volume changes caused by deformation; these local volume changes are characterized by the determinant value of the local Jacobian matrix, which, when close to one, indicates approximately isochoric deformation. The calculation of this part involves squaring the result after subtracting one from the determinant value and then weighting it to encourage approximately isochoric deformation that conforms to the physical properties of biological tissue. This is based on the reasonable assumption that biological soft tissue is incompressible over short periods of time, preventing physically unrealistic extreme compression or expansion.

[0121] The third part, the uncertainty-weighted smoothing term, is calculated as follows: First, the predicted uncertainty variance map is resampled to the same coordinate system as the reference image through the inverse mapping of the current deformation field, resulting in a resampled uncertainty map. Then, an uncertainty gating function is defined, which maps the uncertainty value at each position in the resampled uncertainty map to a weight between zero and one. Specifically, a positive steepness coefficient and a predefined uncertainty threshold are set. The uncertainty value at the current position is subtracted from the threshold and multiplied by the steepness coefficient. The resulting product is input into the negative exponent of the natural exponential function for calculation. The result is then incremented by one and the reciprocal is taken to obtain the weight at that position. The uncertainty threshold can be set based on the statistical distribution of historical predicted uncertainties, for example, by taking the median (e.g., the range of values ​​in the predicted uncertainty variance map). The median (e.g., 0.3) and the steepness coefficient control the smoothness of the transition of the weights from zero to one (e.g., the steepness coefficient can be set to 5.0 to produce a relatively sharp but still smooth transition near the threshold); finally, the Frobenius norm of the displacement field gradient matrix is ​​calculated, which is the square root of the sum of the squares of all elements of the matrix, and is used to measure the severity of the deformation. The square of this Frobenius norm is weighted by the weights; thus, in the high uncertainty region, since the weights approach one, a strong smoothing constraint is imposed on the displacement field gradient at this location to prevent unbelievable severe deformation; in the low uncertainty region, since the weights approach zero, the smoothing constraint is relaxed to allow for more flexible deformation adjustment based on high confidence image information, realizing adaptive regularization based on prediction confidence, making the deformation field flexible in "confident" areas and conservative in "doubtful" areas;

[0122] The integral result of the deformation field regularization constraint term is the cumulative amount of the three-part weighted sum over the spatial region; the weighting coefficients corresponding to the three parts need to be determined in model training or empirical parameter tuning to balance the relationship between similarity driving and various physical constraints.

[0123] Finally, the optimization objective function for non-rigid deformation registration calculation is composed of the image similarity measurement loss and the deformation field regularization constraint term weighted by uncertainty. The optimal control point displacement vector is solved by minimizing this objective function, and then the final deformation field is obtained.

[0124] The specific process for calculating and generating the target coordinate mapping, final deformation field, and registration confidence heatmap of the hypoxic region in the radiotherapy equipment coordinate system at future time points is as follows:

[0125] The final deformation field is solved using an optimization strategy combining multi-resolution hierarchical optimization and the Gauss-Newton trust region method. First, the deformation field is initialized and optimized on a coarse-resolution 3D control point mesh to obtain a low-resolution deformation field estimate. Then, the resolution of the control point mesh is progressively refined, using the deformation field obtained at the current coarse resolution as the initial value, and optimization continues at finer resolutions until the preset maximum resolution is reached, resulting in a high-precision final deformation field. Multi-resolution optimization typically sets three to four resolution levels, with the control point mesh spacing halved at each level to balance optimization efficiency and avoid local minima. The coarse-resolution levels capture large-scale deformation. The strategy of optimizing local details at finer resolution levels significantly improves the stability and convergence speed of the optimization process. In each iteration at each resolution level, the Gauss-Newton method is used to approximate the Hessian matrix of the objective function with respect to the optimization variables, and a local quadratic model is constructed. The optimal update step size for the current iteration is determined by solving a trust region subproblem. The mathematical description of the trust region subproblem is: finding an update vector composed of the sequentially arranged adjustments of all control point displacements that minimizes the following objective function value, while simultaneously satisfying that the Euclidean norm of the update vector does not exceed a dynamically set positive number, i.e., the trust region radius. The objective function consists of two terms... Composition: The first term is a scalar obtained by multiplying half of the update vector by its transpose, then by the Hessian matrix approximated by the Gauss-Newton method, and finally by the update vector. The second term is a scalar obtained by multiplying the transpose of the gradient vector of the objective function at the current solution with respect to the control point displacement by the update vector. Here, the update vector represents the search direction and step size; the gradient vector contains the partial derivatives of the objective function with respect to each control point displacement component; the Hessian matrix describes the curvature information of the objective function near the current solution; and the trust region radius limits the magnitude of the update in a single iteration to ensure the stability of the optimization process. The initial trust region radius can be set to 10 times the mean norm of the current control point displacement vector. %, as a reasonable starting point; by solving this subproblem, an update vector is obtained, which is used to update the control point displacement and evaluate the descent of the actual objective function. At the same time, the trust region radius of the next iteration is dynamically adjusted according to the degree of matching between the model prediction and the actual descent. The adjustment strategy of the trust region radius is as follows: if the descent predicted by the model matches the actual situation well (for example, the actual descent reaches 70% or more of the predicted descent), the trust region radius is expanded to accelerate convergence; if the prediction is poor, the trust region radius is reduced to improve stability. This dynamic adjustment mechanism enables the optimization process to advance quickly with large steps in flat areas and to explore cautiously with small steps in complex terrain.

[0126] After obtaining the final deformation field, its inverse mapping is calculated using numerical methods. The spatial coordinates of all voxels in the hypoxic region defined by the hypoxic mean map of the future time in the floating image are transformed to the coordinate system of the reference image, i.e., the coordinate system of the radiotherapy equipment, through inverse mapping, thereby obtaining the target coordinate mapping. The hypoxic region can be obtained by segmenting it by setting a hypoxic threshold on the predicted hypoxic mean map. This hypoxic threshold can be set according to the tumor type and clinical experience. For example, the region with a hypoxic value higher than the mean of the whole image plus one standard deviation is defined as the hypoxic target area.

[0127] The registration confidence heatmap is generated by comprehensively evaluating the reliability of the registration results at each spatial location in the reference image; specifically, for each location in the reference coordinate system, its confidence is the product of the following three components:

[0128] The first component is the prediction confidence level, which is calculated based on the value at that position in the uncertainty map after resampling. The lower the uncertainty value, the higher the value of this component. Specifically, it is calculated using a negative natural exponential function with the uncertainty value as the exponent. The adjustment parameter in the exponential function is used to control the strength of the influence of uncertainty on the confidence level. The larger the adjustment parameter, the stronger the attenuation effect of uncertainty on the confidence level.

[0129] The second component is the deformation rationality, which is calculated based on the determinant of the local Jacobian matrix of the final deformation field at this position. The closer the determinant value is to one, the higher the value of this component. Specifically, it is calculated by a negative natural exponential function with the absolute value of the determinant value minus one as the exponent. This component is used to exclude physically unreasonable regions of extreme compression or expansion deformation.

[0130] The third component is the local matching degree, calculated based on the image similarity metric of the local neighborhood at that location. A higher value indicates better alignment quality. The calculation process for this component is as follows: First, following the calculation process of local cross-correlation, calculate the image similarity metric of the local neighborhood at that location, i.e., the local correlation ratio. Then, take the negative of this similarity metric to obtain the local similarity loss value at that location. Next, find the maximum value of the local similarity loss value among all locations in the reference image, called the global maximum local loss value. Finally, divide the local similarity loss value at that location by the global maximum local loss value to obtain a ratio. Subtract this ratio from one to get the result, which is the local matching degree component. This construction ensures that the component is close to 1 at the best matching location (where the local loss is minimal) and close to 0 at the worst matching location. This component directly reflects the local accuracy of image registration; the smaller the loss, the higher the matching degree.

[0131] By multiplying these three components, a registration confidence heatmap with values ​​ranging from zero to one is generated. The closer the value is to one, the more reliable the mapping result at that location. These three components evaluate reliability from three independent dimensions: predicted source, deformation physics, and image matching. Their product result integrates multiple aspects of information and is more robust than a single indicator. This heatmap can intuitively display the confidence level of different spatial locations in the registration result, providing key input for risk trade-offs in subsequent beam optimization. For example, during beam optimization, regions with high confidence can be assigned higher dose irradiation priority, while regions with low confidence can adopt a more conservative strategy.

[0132] In this embodiment, it is specifically necessary to explain the process of constructing the objective function for beam parameter optimization based on the confidence-weighted biological equivalent dose model in the beam optimization module as follows:

[0133] First, the time-integrated physical dose to be received at each spatial location in the radiotherapy equipment coordinate system within the future optimized time window is calculated. This time-integrated physical dose is calculated by integrating and summing the contributions of all subfields of all beams within the optimized time window, based on the subfield dose rate weights of all subfields of all beams within the optimized time window and the unit dose rate distribution determined by the position of the multi-leaf grating blades over time. The length of the optimized time window can be set according to the system latency and prediction effectiveness, for example, from 0.5 seconds to 2 seconds, and the integration step size is consistent with the subfield weight update cycle, such as 10 milliseconds.

[0134] Secondly, a spatiotemporally dependent biological effect model is constructed to calculate the bioequivalent dose received by each voxel in the target hypoxic region within a future optimized time window. The target hypoxic region is defined by a target coordinate mapping. The calculation of this bioequivalent dose depends not only on the time-integrated physical dose but also on the real-time hypoxia level of the voxel throughout the entire optimized time window. The real-time hypoxia level is represented by a sequence of hypoxia values ​​at the voxel's location that evolves over time, obtained by mapping the predicted hypoxia mean map sequence at multiple time points within the future optimized time window using the final deformation field onto the radiotherapy equipment coordinate system, as output by the dynamic prediction module. The temporal resolution of the hypoxia value sequence is consistent with the output frame rate of the predicted hypoxia mean map sequence, for example, 10 frames per second, thus characterizing the hypoxia dynamics within the optimized window. Specifically, the model includes two effective radiosensitivity parameters related to the hypoxia level: the alpha effective parameter and the beta effective parameter. The alpha effective parameter is obtained by averaging the reciprocal of the oxygen enhancement ratio of the voxel at each moment within the optimized time window over time.

[0135] The beta effective parameter is obtained by averaging the reciprocal of the square of the oxygen enhancement ratio at each moment within the optimized time window for this voxel. The oxygen enhancement ratio is a parameter characterizing the tissue's resistance to radiation under hypoxic conditions relative to normoxic conditions; its value is a function of the hypoxia level, and the specific functional relationship can be obtained by fitting in vitro cell experimental data, for example, increasing exponentially with increasing hypoxia level. The specific formula for calculating the bioequivalent dose is: alpha effective parameter multiplied by the time-integrated physical dose, plus beta effective parameter multiplied by a treatment time-related intra-fractional repair factor, and then multiplied by the square of the time-integrated physical dose. The intra-fractional repair factor is a parameter used to correct changes in biological effects caused by prolonged irradiation time in a single treatment; its value can be calculated using the time factor formula in the linear square model.

[0136] Finally, a confidence-weighted prescription gap function is constructed as a core component of the beam parameter optimization objective function, defining the prescription biological effect as the clinically required target biological effect value; for each voxel in the target hypoxic region, its prescription gap function consists of two terms:

[0137] The first item is a penalty for non-compliance in high-confidence areas. Its value is the value of the registration confidence heatmap corresponding to the voxel, multiplied by the square of the difference between the prescribed biological effect of the voxel and the actual calculated bioequivalent dose (when the difference is positive). This item ensures that in areas with high confidence of location, the actual delivered bioequivalent dose is as close as possible to the target value of the prescribed biological effect.

[0138] The second item is the over-irradiation penalty for low-confidence areas. Its value is a preset risk aversion coefficient multiplied by a factor subtracted from the registration confidence heatmap value, and then multiplied by the square of the actual calculated bioequivalent dose for that voxel. The risk aversion coefficient is used to balance the irradiation risk of uncertain areas, and its value can be set in the range of 0.1 to 0.5 according to clinical risk preferences, for example, set to 0.2. The beam parameter optimization objective function is the sum of the prescription difference functions of all voxels in the target hypoxic area.

[0139] The specific process of solving the objective function in real time and generating the weights of each beam subfield and the adjustment amounts for the positions of the multi-leaf grating blades is as follows:

[0140] The real-time solution problem of the beam parameter optimization objective function is modeled as a stochastic optimal control problem in the finite-time domain. The state variables of this problem include the cumulative machine hops of each beam, the position and velocity of the multi-leaf grating blades, and the gantry angle. The control variables include the dose rate change rate of each beam subfield and the target velocity command of the multi-leaf grating blades. The optimization objective is to minimize an expected cumulative loss while satisfying the physical and kinematic constraints of the radiotherapy equipment. This expected cumulative loss consists of the sum of the integrals of three parts within the optimization time window:

[0141] The first part is the sum of confidence-weighted prescription gap function values ​​for all voxels within the target hypoxic region;

[0142] The second part is the sum of dose penalty function values ​​for all voxels in the organ at risk regions, for example, using a dose-volume penalty function;

[0143] The third part is a control variable smoothness regularization function value multiplied by a regularization weight coefficient, which is used to suppress high-frequency jitter of control commands and ensure stable operation of the equipment.

[0144] The expectation operation in the expected cumulative loss is performed on the prediction uncertainty distribution provided by the dynamic prediction module. It calculates the cumulative loss for each of the multiple possible sample hypoxia mean spectrum obtained from the prediction uncertainty distribution, and then approximates the cumulative loss by taking the average of these sample losses. The number of sampling groups can be set according to the real-time requirements of the calculation, for example, 20 samplings. The optimization problem also needs to comply with system dynamic constraints, which describe that the rate of change of the state variable with time is a function of the current state, control variables, and time. At the same time, it also needs to satisfy machine physics constraints, including but not limited to the dose rate limit, the maximum speed and acceleration limit of the multi-leaf grating blades, and the frame angular velocity limit.

[0145] A model predictive control framework is employed, combining differential dynamic programming algorithms and unscented transform techniques to solve stochastic optimal control problems in real-time rolling time domains. During each control cycle, the following operations are performed:

[0146] Based on the current radiotherapy equipment status and the latest received predicted mean hypoxia map, predicted uncertainty variance map, and registration confidence heatmap, a stochastic optimal control problem is reconstructed and solved in a shortened prediction time domain. The control period is synchronized with the system data update period, for example, 50 milliseconds. The shortened prediction time domain is typically less than or equal to the optimization time window. During the solution process, the unscented transformation technique is used to approximate the propagation of prediction uncertainty through system dynamics to efficiently calculate the expected cumulative loss. The unscented transformation approximates the propagation of the probability distribution by selecting a carefully designed set of sampling points; for example, for an n-dimensional state, 2n+1 sampling points are used. The differential dynamic programming algorithm is used to solve this finite-time deterministic optimal control problem to obtain the optimal control command sequence for the current control period. The optimal control command sequence contains the control variable values ​​for multiple future control steps.

[0147] From the optimal control command sequence obtained from the solution, the command of the first control cycle is extracted as the adjustment amount for the weight of the beam subfield and the position of the multi-leaf grating blades in the current control cycle. The adjustment amount is input into a forward simulator for safety verification to ensure that it does not violate the dose rate mutation limit, the maximum speed limit of the multi-leaf grating blades, and the anti-collision rules between blades. The specific process of safety verification is as follows: the subfield weight change sequence and the multi-leaf grating blade position trajectory implied by the adjustment amount over a future period are substituted into a fast-calculated dose delivery and equipment motion simulation model to check whether it exceeds any preset machine safety hard constraints. The fast simulator uses a simplified pencil beam or convolution algorithm for dose calculation and simulates the blade position based on the kinematic model to complete the verification in milliseconds. If the verification passes, a timestamp and beam identifier will be included. The beam adjustment control command package containing the subfield weight adjustment amount and the multi-leaf grating target position sequence is sent to the actuator of the radiotherapy equipment. If the verification fails, a safety fallback strategy is triggered. The specific content of the safety fallback strategy is as follows: while keeping the current beam parameters unchanged, a pre-calculated low-risk static optimization plan based on the historical average predicted hypoxia distribution is immediately invoked to generate corresponding backup control commands and issue a warning. The goal of the low-risk static optimization plan is to strictly meet all machine constraints while covering the historical average hypoxia area. Its calculation does not depend on real-time prediction, so it serves as a reliable backup plan. Finally, the command package or backup command that has passed the safety verification is output to the actuator of the radiotherapy equipment through a low-latency communication interface. The actuator includes a dose rate controller, a multi-leaf grating driver, and a gantry controller.

[0148] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0149] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0150] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0151] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0152] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0153] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

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

Claims

1. A system for real-time monitoring of hypoxia in the tumor microenvironment and optimization of radiotherapy protocols, characterized in that, Specifically, it includes: The system includes a hypoxia map reconstruction module, a dynamic prediction module, a coordinate mapping module, and a beam optimization module. Hypoxia atlas reconstruction module: Receives low-resolution real-time hypoxia image sequences from medical imaging equipment and inputs them into a preset spatiotemporal constraint generative adversarial network model. The spatiotemporal constraint generative adversarial network model is trained and inferred by calculating a comprehensive loss function including spatiotemporal gradient consistency loss, and generates and outputs a high spatiotemporal resolution three-dimensional hypoxia atlas sequence. Dynamic prediction module: Connected to the output of hypoxia map reconstruction module, it receives high spatiotemporal resolution 3D hypoxia map sequence, converts each time frame data of the high spatiotemporal resolution 3D hypoxia map sequence into spatiotemporal map structure data and inputs it into gated spatiotemporal map convolutional network for processing. The gated spatiotemporal map convolutional network performs forward propagation to output the predicted hypoxia mean map for a future set time point, and calculates the corresponding prediction uncertainty variance map based on multiple forward propagation sampling. Finally, it outputs a prediction result data packet containing the predicted hypoxia mean map and the prediction uncertainty variance map. Coordinate mapping module: Connected to the output of the dynamic prediction module, it obtains the prediction result data packet, uses the predicted hypoxia mean map in the prediction result data packet as the floating image to be registered, and uses the pre-stored treatment plan CT image as the reference image to perform non-rigid deformation registration calculation. In the process of optimizing the deformation field, the non-rigid deformation registration calculation uses the prediction uncertainty variance map to weight the smoothness constraint of the deformation field, and calculates and generates the target coordinate mapping, final deformation field and registration confidence heatmap of the hypoxia area in the radiotherapy equipment coordinate system at future time. Beam optimization module: Connected to the output of coordinate mapping module, it receives target coordinate mapping, final deformation field and registration confidence heatmap, constructs beam parameter optimization objective function based on confidence weighted bioequivalent dose model, solves and generates adjustment amounts for the weights of each beam subfield and the position of the multi-leaf grating blades in real time according to the optimization objective function, and outputs beam adjustment control commands corresponding to the adjustment amounts to the actuator of radiotherapy equipment.

2. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 1, characterized in that: In the hypoxia map reconstruction module, the spatiotemporal constrained generative adversarial network model includes a generator network. The generator network receives a low-resolution real-time hypoxia image sequence as input and performs spatiotemporal feature encoding and fusion processing on the input low-resolution real-time hypoxia image sequence, specifically: The generator network first uses a three-dimensional convolutional coding unit to perform convolution operations on a low-resolution real-time hypoxia image sequence that contains multiple consecutive time frames within a sliding time window centered on the current time point to be reconstructed, simultaneously in the spatial and slice dimensions, to extract the deep spatial features corresponding to each time frame within the sliding time window. Subsequently, the generator network uses a one-dimensional temporal convolutional unit to convolve the sequence composed of deep spatial features of each time frame within the sliding time window along the time dimension, thereby extracting the temporal evolution features of the sliding time window. Finally, the generator network fuses the deep spatial features and temporal evolution features of the current time point to be reconstructed to generate a latent feature map that integrates spatiotemporal information corresponding to the current time point to be reconstructed.

3. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 2, characterized in that: The spatiotemporal gradient consistency loss included in the comprehensive loss function is composed of a weighted sum of spatial gradient consistency components, temporal gradient consistency components, and cyclic consistency components; the weighted sum is obtained by multiplying the spatial gradient consistency component by a first weight coefficient, the temporal gradient consistency component by a second weight coefficient, and the cyclic consistency component by a third weight coefficient, and then adding them together. The calculation process of the spatial gradient consistency component is as follows: the generator network calculates the gradient magnitude distribution along the three-dimensional spatial direction of any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the latent feature map, and calculates the divergence between the gradient magnitude distribution and a prior spatial gradient magnitude distribution that conforms to the characteristics of biological tissue and is learned from a pre-acquired high-resolution gold standard image. The calculation process of the temporal gradient consistency component is as follows: The difference between the maps corresponding to two adjacent time points in the high spatiotemporal resolution 3D hypoxia map sequence generated by the generator network is calculated to obtain the first difference image; the difference between the images corresponding to two time points in the low-resolution real-time hypoxia image sequence is calculated, and this difference is upsampled to the same size as the high spatiotemporal resolution 3D hypoxia map to obtain the second difference image; the gradient magnitudes of the first and second difference images along the time dimension are calculated respectively; finally, the first norm distance between the gradient magnitudes of the first and second difference images is calculated. The calculation process of the cyclic consistency component is as follows: a preset three-dimensional Gaussian blur downsampling operator, used to simulate the imaging process of medical imaging equipment, is used to process any map in the high spatiotemporal resolution three-dimensional hypoxia map sequence generated by the generator network to obtain a simulated low-resolution image; the first norm distance between the simulated low-resolution image and the image at the corresponding time point in the input low-resolution real-time hypoxia image sequence is calculated.

4. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 3, characterized in that: The spatiotemporal constrained generative adversarial network model is trained by optimizing the comprehensive loss function. After training, it performs forward inference on the input low-resolution real-time hypoxia image sequence. In the forward inference process, spatiotemporal feature encoding and fusion processing are performed independently on each time point in the input low-resolution real-time hypoxia image sequence. The decoder of the generator network upsamples and decodes the latent feature map to generate a high spatiotemporal resolution three-dimensional hypoxia map corresponding to the current time point. Finally, a high spatiotemporal resolution three-dimensional hypoxia map sequence is output.

5. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 4, characterized in that: In the dynamic prediction module, the specific operation of converting each time frame data of the high spatiotemporal resolution three-dimensional hypoxia map sequence into spatiotemporal map structure data is as follows: A1. Define each voxel in the 3D hypoxia map of each time frame as an independent graph node, and construct a feature vector containing multidimensional information for each node; The vector consists of five components in sequence: the hypoxia value of the voxel at the current time point, the gradient component of the voxel along the X-axis in the three-dimensional coordinate system, the gradient component along the Y-axis, the gradient component along the Z-axis, and the approximate instantaneous rate of change obtained by calculating the difference between the hypoxia value of the voxel at the current time point and the immediately preceding time point. A2. Construct connecting edges between graph nodes based on dual criteria to form a spatiotemporal graph structure; the dual criteria include spatial proximity criterion and feature similarity criterion; Based on the spatial proximity criterion, a connection edge is established between each node and a predetermined number of voxel nodes whose spatial coordinates are closest to each other in three-dimensional space. According to the feature similarity criterion, the cosine similarity between the feature vectors of any two nodes is calculated. If the cosine similarity calculation result exceeds the preset first similarity threshold, a connection edge is established between the two nodes. A3. In the time dimension, directed connection edges are established between nodes corresponding to voxels at different consecutive time points to represent the state evolution path of the voxel at the same spatial coordinate position. By independently performing the above A1-A3 operations on each time frame in the high spatiotemporal resolution three-dimensional hypoxia map sequence, a corresponding spatiotemporal graph structure data containing the feature set of all nodes in that time frame and the set of connecting edges between nodes is generated for each time frame. These spatiotemporal graph structure data arranged in chronological order together constitute a dynamically evolving spatiotemporal graph data sequence that serves as the input to a gated spatiotemporal graph convolutional network.

6. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 5, characterized in that: The specific process by which the gated spatiotemporal graph convolutional network processes and calculates the prediction uncertainty variance map is as follows: The gated spatiotemporal graph convolutional network consists of alternating stacked spatial graph convolutional layers and gated temporal convolutional layers. In the spatial graph convolutional layer, for the current central node to be processed, the processing procedure is as follows: B1. Apply a shared learnable weight matrix to the feature vector of the central node and the feature vector of each of its neighboring nodes to obtain the transformation features of the central node and the transformation features of each neighboring node. B2. Concatenate the transformation features of the center node with the transformation features of a neighbor node to form a concatenated feature vector; perform a dot product operation on the concatenated feature vector and a learnable attention vector, and then process the result through a leaky linear rectified activation function to obtain the original attention score of the neighbor node relative to the center node. Repeat this process for all neighboring nodes to obtain a set of raw attention scores; B3. Normalize the original attention scores using an exponential function to convert each original score into a value between zero and one. These values ​​are the attention coefficients corresponding to each neighbor node. B4. Multiply the transformation feature of each neighboring node by its corresponding attention coefficient, and sum the product of all neighboring nodes to obtain the spatial aggregation feature of the central node. In the gated temporal convolutional layer, for each node, its spatial aggregated features output by the spatial graph convolutional layer at multiple consecutive historical time steps are arranged in chronological order to form the spatial aggregated feature sequence of that node; the processing procedure of this layer is as follows: Two sets of independent learnable one-dimensional convolution kernels are used to perform one-dimensional convolution operations on the input spatial aggregated feature sequence respectively; the hyperbolic tangent activation function is applied to the convolution result of the first set of convolution kernels to obtain the feature extraction components; Apply a sigmoid activation function to the convolution result of the second set of convolution kernels to obtain the gated signal component; finally, multiply the feature extraction component and the gated signal component element by element to obtain the hidden state feature of the node at the current time step. The multi-step prediction process is as follows: The gated spatiotemporal graph convolutional network is unfolded in the time dimension. It uses the spatiotemporal graph data sequence of multiple consecutive time steps in the past and processes it through alternating stacked spatial graph convolutional layers and gated temporal convolutional layers to obtain the hidden state features of each node at the last input time step. Then, through a linear transformation layer consisting of a learnable weight matrix and bias terms, the hidden state features of each node are mapped to the hypoxia value prediction for a future set time step; the linear transformation layer constitutes the final prediction output layer of the gated spatiotemporal graph convolutional network; the predicted hypoxia values ​​of all nodes are rearranged according to their original spatial positions, thus forming the predicted hypoxia mean map for a future set time point. The specific process of calculating the prediction uncertainty variance map based on multiple forward propagation sampling is as follows: After training the gated spatiotemporal graph convolutional network, a random deactivation operation is kept active on the linear transformation weights of its spatial graph convolutional layers and the final prediction output layer. For the same set of historical spatiotemporal graph data sequences as input, multiple independent forward propagation calculations are performed. Due to the randomness of the random deactivation operation, each forward propagation yields a slightly different prediction result sample. Next, a discretization interval partitioning scheme covering the possible hypoxia value range is set, mapping continuous hypoxia prediction values ​​to corresponding intervals. For each voxel position in the predicted hypoxia mean map, the frequency distribution of its predicted value falling into each interval in multiple samplings is statistically analyzed. Finally, the prediction uncertainty at that position is quantified by calculating the statistical information entropy of this frequency distribution. The information entropy is specifically calculated as the sum of the following two parts: The first part is to calculate the average frequency of each interval across all samples for all discretized intervals, then calculate the product of the average frequency and the natural logarithm of the average frequency, and finally sum the product of all intervals and take the negative value. The second part involves calculating the information entropy of the frequency distribution in each independent sample, and then averaging the information entropy results of all multiple samples. The uncertainty values ​​of all voxel positions calculated in this way together constitute the prediction uncertainty variance map. The final output of the dynamic prediction module is a prediction result data package containing a prediction mean hypoxia map and a prediction uncertainty variance map.

7. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 6, characterized in that: The specific process of performing non-rigid deformation registration calculation in the coordinate mapping module is as follows: First, a dense spatial deformation mapping is defined from the coordinate system of the reference image to the coordinate system of the floating image. This mapping is described by a three-dimensional displacement field, and the three-dimensional displacement field is parameterized using a free deformation model based on cubic B-splines. Specifically, a three-dimensional control point grid with uniform spacing is defined on the reference image, and each control point in the three-dimensional control point grid is associated with a three-dimensional displacement vector. The displacement of any point in space is obtained by weighted interpolation of the displacements of its neighboring control points by a predetermined number of control points using cubic B-spline basis functions. The optimization variable is the three-dimensional displacement vector of all control points in the three-dimensional control point grid. Secondly, with the goal of maximizing the local structural correlation between the deformed floating image and the reference image, the image similarity metric loss is calculated. The image similarity metric uses local cross-correlation, and its calculation process is as follows: In the reference image, a local neighborhood is defined with each voxel as the center, and the mean intensity of all voxels in the reference image and the mean intensity of all voxels in the deformed floating image within this neighborhood are calculated. Then, the deviation between the intensity of the reference image and its mean, and the deviation between the intensity of the deformed floating image and its mean, are calculated at each position within this neighborhood. Next, the square of the product of these two deviations is calculated, and then the sum of the squares of these two deviations is calculated. Finally, the square of the deviation product is divided by the product of the sums of the squares of the two deviations to obtain the local correlation ratio at the central voxel. The local correlation ratios of all voxels in the entire reference image are summed, and the sum is negative, which is taken as the total image similarity metric loss. Next, a deformation field regularization constraint term weighted by uncertainty information is constructed. This term is the integral result of the deformation field over the spatial region covered by the reference image. Its integrand consists of a weighted sum of three parts, each of which is related to the local geometric properties of the deformation field: The first part is used to penalize the degree to which the local Jacobian matrix of the deformation field deviates from the orthogonal transformation; the local Jacobian matrix represents the spatial gradient of the deformation field at a point, and is composed of the identity matrix plus the gradient matrix of the displacement field at that point; The second part is used to penalize local volume changes caused by deformation; local volume changes are characterized by the determinant value of the local Jacobian matrix. The third part, the uncertainty weighted smoothing term, is calculated as follows: First, the predicted uncertainty variance map is resampled to the same coordinate system as the reference image through the inverse mapping of the current deformation field, resulting in a resampled uncertainty map. Then, an uncertainty gating function is defined, which maps the uncertainty value at each position in the resampled uncertainty map to a weight between zero and one. Specifically, a positive steepness coefficient and a predefined uncertainty threshold are set. The uncertainty value at the current position is subtracted from the threshold and then multiplied by the steepness coefficient. The resulting product is input into the negative exponent of the natural exponential function for calculation. The result is then incremented by one and the reciprocal is taken to obtain the weight at that position. Finally, the Frobenius norm of the displacement field gradient matrix is ​​calculated, and the square of this Frobenius norm is weighted using the weights. The integral result of the deformation field regularization constraint term is the cumulative amount of the weighted sum of the three parts over the spatial region; Finally, the optimization objective function for non-rigid deformation registration calculation is composed of the image similarity measurement loss and the deformation field regularization constraint term weighted by uncertainty. The optimal control point displacement vector is solved by minimizing this objective function, and then the final deformation field is obtained.

8. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 7, characterized in that: The specific process for calculating and generating the target coordinate mapping, final deformation field, and registration confidence heatmap of the hypoxic region in the radiotherapy equipment coordinate system at future time points is as follows: The final deformation field is solved using an optimization strategy that combines multi-resolution hierarchical optimization with the Gauss-Newton trust region method. First, the deformation field is initialized and optimized on a coarser resolution 3D control point mesh to obtain a low-resolution deformation field estimate. Then, the resolution of the control point mesh is gradually refined, and the deformation field obtained at the current coarser resolution is used as the initial value to continue optimization at a finer resolution until the preset maximum resolution is reached, resulting in a high-precision final deformation field. In each single optimization iteration at each resolution level, the Gauss-Newton method is used to approximate the Hessian matrix of the objective function with respect to the optimization variables, and a local quadratic model is constructed. The optimal update step size for the current iteration is determined by solving a trust region subproblem. After obtaining the final deformation field, its inverse mapping is calculated using numerical methods. The spatial coordinates of all voxels in the hypoxic region defined by the hypoxic mean map predicted at future time in the floating image are transformed to the coordinate system of the reference image, i.e. the coordinate system of the radiotherapy equipment, through inverse mapping, thereby obtaining the target coordinate mapping. The registration confidence heatmap is generated by comprehensively evaluating the reliability of the registration results at each spatial location in the reference image; specifically, for each location in the reference coordinate system, its confidence is the product of the following three components: The first component is the prediction confidence, which is calculated based on the value at that position in the resampled uncertainty map, specifically by using a negative natural exponential function with the uncertainty value as the exponent. The second component is the deformation rationality, which is calculated based on the determinant of the local Jacobian matrix of the final deformation field at that position. Specifically, it is calculated by a negative natural exponential function with the absolute value of the determinant minus one as the exponent. The third component is the local matching degree, calculated based on the image similarity metric of the local neighborhood at that location. The calculation process for this component is as follows: First, following the calculation process of local cross-correlation, calculate the image similarity metric value of the local neighborhood at that location, i.e., the local correlation ratio; then, take the negative of this similarity metric value to obtain the local similarity loss value at that location; next, find the maximum value of the local similarity loss value among all locations in the reference image, called the global maximum local loss value; finally, divide the local similarity loss value at that location by the global maximum local loss value to obtain a ratio, and then subtract this ratio from one to get the result, which is the local matching degree component. By multiplying these three components, a registration confidence heatmap with values ​​ranging from zero to one is finally generated.

9. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 8, characterized in that: In the beam optimization module, the specific process of constructing the objective function for beam parameter optimization based on the confidence-weighted bioequivalent dose model is as follows: First, the time-integrated physical dose that will be received at each spatial location in the coordinate system of the radiotherapy equipment within the future optimized time window is calculated. The calculation of this time-integrated physical dose is based on the subfield dose rate weights of all subfields of all beams within the optimized time window, and the unit dose rate distribution determined by the change of the position of the multi-leaf grating blades over time. The contribution of all subfields of all beams within the optimized time window is integrated and summed. Secondly, a spatiotemporally dependent biological effect model is constructed to calculate the biologically equivalent dose received by each voxel in the target hypoxic region within a future optimized time window; the target hypoxic region is defined by a target coordinate mapping. Specifically, the model includes two effective radiosensitivity parameters related to hypoxia level: the alpha effective parameter and the beta effective parameter. The alpha effective parameter is obtained by averaging the reciprocal of the oxygen enhancement ratio of the voxel at each moment within the optimization time window. The beta effective parameter is obtained by averaging the time-inverse of the square of the oxygen enhancement ratio of the voxel at each moment within the optimized time window; the specific formula for calculating the bioequivalent dose is: the alpha effective parameter multiplied by the time-integrated physical dose, plus the beta effective parameter multiplied by a treatment time-related intra-session repair factor, and then multiplied by the square of the time-integrated physical dose. Finally, a confidence-weighted prescription gap function is constructed, defining the prescription biological effect as the clinically required target biological effect value; for each voxel in the target hypoxic region, its prescription gap function consists of two terms: The first item is the penalty for failing to meet the high confidence region standard. Its value is the value of the registration confidence heatmap corresponding to the voxel, multiplied by the square of the difference between the prescription biological effect of the voxel and its actual calculated bioequivalent dose. The second term is the over-irradiation penalty term for low-confidence areas, which is a preset risk avoidance coefficient multiplied by a factor minus the value of the registration confidence heatmap, and then multiplied by the square of the actual calculated bioequivalent dose of the voxel; the beam parameter optimization objective function is the sum of the prescription difference functions of all voxels in the target hypoxic region.

10. The system for real-time monitoring of tumor microenvironment hypoxia and optimization of radiotherapy protocols according to claim 9, characterized in that: The specific process of solving the objective function in real time and generating the weights of each beam subfield and the adjustment amounts for the positions of the multi-leaf grating blades is as follows: The real-time solution problem of the beam parameter optimization objective function is modeled as a stochastic optimal control problem in the finite-time domain. The state variables of this problem include the cumulative machine hops of each beam, the position and velocity of the multi-leaf grating blades, and the gantry angle. The control variables include the dose rate change rate of each beam subfield and the target velocity command of the multi-leaf grating blades. The optimization objective is to minimize an expected cumulative loss while satisfying the physical and kinematic constraints of the radiotherapy equipment. This expected cumulative loss consists of the sum of the integrals of three parts within the optimization time window: The first part is the sum of confidence-weighted prescription gap function values ​​for all voxels within the target hypoxic region; The second part is the sum of dose penalty function values ​​for all voxels in the regions at risk of organ failure; The third part is a control variable smoothness regularization function value multiplied by a regularization weight coefficient; A model predictive control framework is employed, combining differential dynamic programming algorithms and unscented transform techniques to solve stochastic optimal control problems in real-time rolling time domains. During each control cycle, the following operations are performed: Based on the current status of the radiotherapy equipment and the latest received predicted mean hypoxia map, predicted uncertainty variance map, and registration confidence heatmap, the stochastic optimal control problem is reconstructed and solved in a shortened prediction time domain; the finite-time deterministic optimal control problem is solved using a differential dynamic programming algorithm to obtain the optimal control command sequence for the current control cycle. From the optimal control command sequence obtained from the solution, the command of the first control cycle is extracted as the adjustment amount for the weight of the beam subfield and the position of the multileaf grating blades in the current control cycle. The adjustment amount is input into a forward simulator for safety verification. If the verification is successful, the beam adjustment control command packet containing the timestamp, beam identifier, subfield weight adjustment amount, and multileaf grating blade target position sequence is sent to the actuator of the radiotherapy equipment. If the verification fails, a safety fallback strategy is triggered. The specific content of the safety fallback strategy is as follows: while keeping the current beam parameters unchanged, a pre-calculated low-risk static optimization plan based on historical average predicted hypoxia distribution is immediately invoked to generate corresponding backup control commands and issue a warning. Finally, the command packet that has passed the safety verification is output to the actuator of the radiotherapy equipment through a low-latency communication interface. The actuator includes a dose rate controller, a multileaf grating driver, and a gantry controller.

Images