Proton therapy beam spot weight end-to-end prediction method based on heterogeneous graph neural network driven by physical prior

Abstract

This invention discloses an end-to-end prediction method for beam spot weights in proton therapy based on a physical prior-driven heterogeneous graph neural network, relating to the field of radiotherapy technology. This method acquires the patient's daily CT scan, target area and organ-at-risk contours, and radiation field configuration; converts the HU value to relative stopping power; constructs a beam spot-voxel heterogeneous bipartite graph and pre-calculates dose contributions for edge construction; constructs a hierarchical compressed graph to perform coarse-grained message passing; performs fine-grained message passing after embedding and unfolding; and finally outputs the beam spot weights and generates a standard DICOM file. This invention achieves end-to-end prediction of beam spot weights, eliminates the inverse optimization step, meets the online treatment time requirements, incorporates physical priors to ensure the rationality of the results, and improves the efficiency and clinical applicability of treatment plan generation.

Description

Technical Field

[0001] This invention relates to the field of radiotherapy technology, and in particular to an end-to-end prediction method for proton therapy beam spot weights based on a physical prior-driven heterogeneous graph neural network. Background Technology

[0002] Proton therapy, with its unique Bragg peak physics, can concentrate energy release in the tumor area while significantly reducing radiation damage to surrounding normal tissues, making it an important direction for the development of tumor radiotherapy. Online adaptive proton therapy, as a next-generation treatment technology, requires the treatment plan to be regenerated based on the CT images taken that day within minutes of the patient maintaining the treatment position. This adapts to changes in patient position, tumor regression, and other anatomical alterations, further improving treatment precision. Traditional proton therapy plan generation relies on iterative optimization algorithms, requiring repeated adjustments to beam spot weights to meet the dual requirements of target dose coverage and organ protection. This process typically takes tens of minutes or even hours, failing to meet the real-time requirements of online adaptive therapy.

[0003] Existing deep learning-based proton therapy planning optimization methods mainly fall into two categories. One category uses a 3D convolutional neural network to predict the 3D dose distribution within the patient's body, then uses a traditional inverse optimization engine to deduce executable beam spot weights. While this method can quickly generate dose distributions, the inverse optimization step is essentially an iterative process, becoming a bottleneck for the speed of online adaptive therapy. Furthermore, the separation of dose prediction and weight derivation easily introduces errors. The other category is beam spot weight fine-tuning methods, which use the dose residuals of the reference plan as input to predict weight adjustments. This type of method heavily relies on a high-quality reference plan, lacks the ability to generate treatment plans from scratch, and cannot handle situations with significant changes in anatomical structures. In addition, flux mapping prediction methods for photon intensity-modulated radiotherapy (IMRT) are based on completely different physical models. Photons do not have Bragg peaks or energy layers, and their longitudinal dose distribution characteristics differ fundamentally from those of protons, making them unsuitable for direct transfer to proton therapy scenarios.

[0004] Existing deep learning methods generally lack effective utilization of the physical laws governing proton therapy. Three-dimensional convolutional neural networks encode beam information into fixed-size feature maps, losing the physical coupling relationships between beams and failing to accurately model the dose superposition effect between adjacent beams. Most methods treat treatment plan generation as a purely data-driven regression problem, without explicitly incorporating prior physical knowledge of proton therapy. This results in some predictions not conforming to the physical laws of proton dose distribution, requiring extensive manual corrections by clinicians before they can be used in actual treatment. Furthermore, performing neural network inference directly on the original beam and voxel-scale graphs incurs enormous computational and memory overhead, making it difficult for existing methods to meet the real-time requirements of clinical practice while maintaining accuracy. Summary of the Invention

[0005] The present invention proposes an end-to-end prediction method for proton therapy beam spot weights based on a physical prior-driven heterogeneous graph neural network to solve the problems mentioned in the prior art.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: an end-to-end prediction method for proton therapy beam spot weights based on a physical prior-driven heterogeneous graph neural network, comprising the following steps: The system acquires the patient's CT images from the same day, the target area and the outline of organs at risk obtained by automatic segmentation, and the field configuration file exported by the treatment planning system. The Schneider piecewise linear model is used to convert the HU values ​​of the CT images into a three-dimensional distribution of relative stopping power. Construct a speckle-voxel isomer bipartite map, extract all candidate speckles from the field configuration file as speckle nodes, sample voxels from the target area, organs at risk, and normal tissues within 20mm of the target area from CT images as voxel nodes, if the total number of voxels in the target area exceeds 50,000, perform isotropic downsampling at equal intervals of 2 along the three axes, and allocate sampling quotas for organs at risk according to volume ratio. The pen beam algorithm based on the single Gaussian approximation pre-calculates the dose contribution value of each beam spot to each voxel. The dose contribution value is the product of the depth dose distribution function and the transverse Gaussian scattering kernel. A directed edge from the beam spot to the voxel is established only when the dose contribution value is not less than 0.001 times the maximum dose contribution of the beam spot. A hierarchical compressed graph structure is constructed. Voxel nodes are grouped according to their anatomical structure and super voxel nodes are generated through k-means clustering. Bead nodes are grouped according to their radiation field and energy layer and bead cluster nodes are generated through DBSCAN clustering. Aggregation edges are established between bead clusters and super voxels. Four rounds of coarse-grained message passing are performed on the hierarchical compressed graph. Bidirectional message passing is carried out between super voxel nodes and cluster nodes through a heterogeneous graph attention network. At the same time, horizontal and vertical physical interaction information is passed between cluster nodes to update node embeddings. Perform an embedding unrolling operation, broadcast the embedding of the speckle cluster node to all its member speckle nodes, broadcast the embedding of the super voxel node to all its member voxel nodes, and concatenate it with their original features to generate the initial node embedding for the fine-grained stage through a linear transformation. Two rounds of fine-grained message passing are performed on the original speckle-voxel heterogeneous bipartite graph. Independent reasoning is performed by shooting field grouping. Bidirectional message passing is performed between speckle nodes and voxel nodes to update speckle node embeddings. The final embedding of each beamspot node is processed by the output head of a 3-layer multilayer perceptron. The activation function uses ReLU to ensure that the output is non-negative. The predicted weight value of each beamspot is output, and a beamspot weight file in standard DICOMRT-Plan format is generated.

[0007] Furthermore, it also includes the step of constructing physical interaction edges between beam clusters, establishing lateral interaction edges for spatially adjacent beam clusters within the same energy layer, and establishing vertical interaction edges between beam clusters in adjacent energy layers with an energy difference of no more than 5 MeV. ; in, Cluster of spots with clusters of spots Horizontal interaction edge weights between them Cluster of spots The coordinates of the centroid in the plane of equal centers, Cluster of spots The coordinates of the centroid in the plane of equal centers, This is the half-width of the beam spot lateral scattering corresponding to this energy layer.

[0008] Furthermore, it also includes a multi-loss function joint training step, which uses a weighted combination of weighted regression loss, physical consistency dose regularization loss and dose-volume histogram constraint loss as the total loss function. During the training process, a temperature annealing strategy is used to gradually approximate the accurate dose-volume histogram index constraint. The initial temperature value is set to 5.0, and it is multiplied by a decay factor of 0.5 every 50 training cycles until it drops to 0.1.

[0009] Furthermore, when constructing the speckle-voxel isomer bipartite map, the speckle node adopts a 7-dimensional feature vector, including the speckle nominal energy, the x-axis incident position of the isocenter plane, the y-axis incident position of the isocenter plane, the gantry angle of the field, the treatment bed angle, the half-width of the speckle lateral scattering, and the reference plan weight. All of the above features are directly extracted from the standard DICOMRT-Plan format file. The voxel node feature vector includes the relative stopping power value, the one-heat encoding of the anatomical structure label, the Euclidean distance to the geometric center of the target area, and the water equivalent depth in each field direction.

[0010] Furthermore, when grouping voxel nodes according to anatomical structure and generating super voxel nodes through clustering, voxel groups are divided according to anatomical structure labels based on the automatic segmentation results. Anatomical structures with voxel sizes reaching a set threshold are further subdivided using the k-means clustering algorithm combined with three-dimensional spatial coordinates. The number of clusters is adaptively set based on the number of structural voxels. Each cluster constitutes a super voxel node. The feature vector of the super voxel node integrates regional voxel statistical features, regional volume, anatomical labels, and geometric distance features. The aggregation edges between fasciculus and super voxels are filtered and retained based on edge weight threshold conditions.

[0011] Furthermore, when grouping beam spot nodes by radiation field and energy layer and generating beam spot cluster nodes through clustering, beam spots of the same energy layer within the same radiation field are grouped using the DBSCAN clustering algorithm. Each cluster forms a beam spot cluster node. The feature vector of the beam spot cluster node includes the arithmetic mean of the corresponding features of member beam spots, the number of member beam spots, and the width and height of the bounding rectangle of the spatial distribution of member beam spots. The vertical interaction edge weight between beam spot clusters of adjacent energy layers is calculated as follows: ; in, Cluster of spots with clusters of spots Vertical interaction edge weights between them For energy The beam pattern at the water equivalent depth The dose value at the central axis depth, Cluster of spots The average energy, Cluster of spots The average energy, It is 1.2 times the larger range of the two depth dose distribution curves.

[0012] Furthermore, when performing coarse-grained message passing on the hierarchical compressed graph, the original features of the super voxel nodes and cluster nodes are mapped to a unified latent space through independent linear transformations, and four rounds of bidirectional message passing are performed. In each round, the cluster nodes first transmit dose information to the super voxel nodes, and the edge weights are calculated by combining node embedding and edge weights with an attention mechanism, and the super voxel nodes are updated by weighted aggregation. Then, the super voxel nodes feed back dose constraint information to the cluster nodes, and the basic update message of the cluster nodes is obtained by using the same attention mechanism. Combined with the interaction information of adjacent cluster nodes, the cluster node embedding is updated by the feedforward network and layer normalization.

[0013] Furthermore, when performing fine-grained message passing on the original beam spot-voxel isomer bipartite graph, the entire planned beam spot is divided into 3 to 5 independent field groups according to the field affiliation. Each field group retains only voxel nodes that have effective edge associations with the beam spots within the group. Two rounds of fine-grained bidirectional message passing are performed independently on each field group. The original dose contribution value, water equivalent projection depth, and lateral distance are used as edge features and participate in attention calculation after logarithmic compression transformation. Each field group completes information interaction through coarse-grained global graph operations.

[0014] Furthermore, the multi-loss function joint training step adopts an adaptive batch processing strategy, dynamically controlling the batch size based on the differences in sample graph structure scale. Upper limits for the total number of nodes in the batch are set for both the coarse-grained and fine-grained stages, and samples are accumulated one by one until the corresponding upper limit is reached to complete batch construction. The graph structure data of all samples is pre-calculated and serialized for storage, and the data is directly read during the training stage to avoid online graph construction overhead. The AdamW optimizer is used to complete model training, and parameters such as the weight decay coefficient and learning rate are configured. The learning rate is dynamically adjusted using a cosine annealing restart scheduling strategy.

[0015] Furthermore, the online inference process first acquires the patient's CT images from the current day and obtains the target area and organ-at-risk contours through an automatic segmentation algorithm, converting the CT images into relative stopping power distributions; based on the preset field configuration, the dose contribution matrix of all beam spots is pre-calculated in parallel on the GPU, a hierarchical heterogeneous graph is constructed and input into the trained neural network model, and the prediction weights of each beam spot are obtained through forward inference; a standard DICOMRT-Plan file is generated according to the treatment plan system format requirements, and the overall inference process meets the clinical time window requirements.

[0016] Compared with existing technologies, the beneficial effects of this invention are: This invention achieves end-to-end prediction of beam spot weights in proton therapy, directly outputting executable beam spot weight values. It completely eliminates the time-consuming inverse optimization step in traditional methods, significantly shortening the treatment plan generation time and meeting the clinical time window requirements of online adaptive proton therapy. The entire inference process requires no manual intervention, reducing the workload of clinicians and contributing to the widespread clinical application of online adaptive proton therapy.

[0017] This invention employs a physics-prior-driven heterogeneous graph neural network architecture, constructing a graph topology based on the physical dose contribution relationship between beam spots and voxels, naturally encoding the physical causal relationships of proton therapy. It explicitly models the lateral scattering coupling of beam spots in the same energy layer and the longitudinal dose superposition effect of beam spots in adjacent energy layers, enabling the model to learn beam spot weight allocation patterns that conform to the physical laws of proton therapy, thus improving the physical consistency and clinical applicability of the prediction results.

[0018] This invention employs a hierarchical graph structure and a field-based grouping inference strategy. Through two-level compression of voxels and fasciculations, it significantly reduces the computational scale while preserving crucial spatial and physical information. The fine-grained inference method based on field grouping effectively controls memory usage, enabling the model to run stably on conventional medical computing devices. The combined coarse-grained and fine-grained message passing mechanism ensures both the satisfaction of global dose constraints and the capture of fine-grained dose distribution details.

[0019] This invention employs a multi-loss function joint training strategy, comprehensively considering weighted regression accuracy, physical dose consistency, and clinical dose-volume histogram constraints. The temperature annealing strategy enables the model to progressively approach precise clinical dose requirements during training, generating treatment plans that simultaneously meet clinical standards for uniform target dose coverage and dose limitations for organs at risk. All model inputs come from standard format files in the treatment planning system, requiring no additional data collection and facilitating integration into existing clinical workflows. Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating the overall prediction process of the end-to-end prediction method for proton therapy beam spot weights based on a physical prior-driven heterogeneous graph neural network proposed in this invention. Figure 2 Flowchart for constructing a hierarchical heterogeneous graph; Figure 3 Construct a mechanism graph for physical interaction edges; Figure 4 This is a flowchart of message passing in a two-stage graph neural network. Figure 5 The flowchart shows the multi-loss joint training and adaptive optimization process. Detailed Implementation

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

[0022] In the description of this invention, it should be understood that the terms "center," "longitudinal," "lateral," "length," "width," "thickness," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," "outer," "clockwise," and "counterclockwise," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0023] Furthermore, 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 indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include one or more of the stated features. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified. Furthermore, the terms "installed," "connected," and "linked" should be interpreted broadly; for example, they may refer to a fixed connection, a detachable connection, or an integral connection; they may refer to a mechanical connection or an electrical connection; they may refer to a direct connection or an indirect connection through an intermediate medium; and they may refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances. The invention will now be described in further detail with reference to the accompanying drawings.

[0024] Reference Figures 1 to 5 An end-to-end prediction method for proton therapy beam spot weights based on a physics-prior-driven heterogeneous graph neural network includes the following steps: The system acquires the patient's CT images from the same day, the target area and the outline of organs at risk obtained by automatic segmentation, and the field configuration file exported by the treatment planning system. The Schneider piecewise linear model is used to convert the HU values ​​of the CT images into a three-dimensional distribution of relative stopping power. A beam spot-voxel isomer bipartite map is constructed. All candidate beam spots are extracted from the field configuration file as beam spot nodes. Voxels of the target area, organs at risk, and normal tissue within 20 mm outside the target area are sampled from CT images as voxel nodes. If the total number of target area voxels exceeds 50,000, isotropic downsampling is performed at equal intervals of 2 along the three axes. The sampling quota for organs at risk is allocated according to the volume ratio and the total quota is capped at 100,000. The sampling quota for normal tissue is 10,000. The pencil beam algorithm based on single Gaussian approximation pre-calculates the dose contribution value of each beam spot to each voxel. The dose contribution value is the product of the depth dose distribution function and the lateral Gaussian scattering kernel. A directed edge from the beam spot to the voxel is established only when the dose contribution value is not less than 0.001 times the maximum dose contribution of the beam spot. The edge features include the dose contribution value, the water equivalent projection depth, and the lateral distance. A hierarchical compressed graph structure is constructed. Voxel nodes are grouped according to their anatomical structure and super voxel nodes are generated by k-means clustering. Bead nodes are grouped according to their radiation field and energy layer and bead cluster nodes are generated by DBSCAN clustering. Aggregation edges are established between bead clusters and super voxels. The weight of the aggregation edge is the sum of the dose contributions of all member voxels in the super voxel multiplied by the volume of a single voxel. Four rounds of coarse-grained message passing are performed on the hierarchical compressed graph. Bidirectional message passing is carried out between super voxel nodes and cluster nodes through a heterogeneous graph attention network. At the same time, horizontal and vertical physical interaction information is passed between cluster nodes to update node embeddings. Perform an embedding unrolling operation, broadcast the embedding of the speckle cluster node to all its member speckle nodes, broadcast the embedding of the super voxel node to all its member voxel nodes, and concatenate it with their original features to generate the initial node embedding for the fine-grained stage through a linear transformation. Two rounds of fine-grained message passing are performed on the original speckle-voxel heterogeneous bipartite graph. Independent reasoning is performed by shooting field grouping. Bidirectional message passing is performed between speckle nodes and voxel nodes to update speckle node embeddings. The final embedding of each fascicle node is processed by the output head of a 3-layer multilayer perceptron. The activation function uses ReLU to ensure that the output is non-negative. The predicted weight value of each fascicle is output, and a standard DICOMRT-Plan format fascicle weight file that can be directly imported into the treatment planning system is generated.

[0025] This invention also includes a step of constructing physical interaction edges between beam clusters: establishing lateral interaction edges for spatially adjacent beam clusters within the same energy layer, and establishing longitudinal interaction edges between beam clusters in adjacent energy layers with an energy difference not exceeding 5 MeV. ; in, Cluster of spots with clusters of spots The horizontal interaction edge weights between them are dimensionless. Cluster of spots The centroid coordinates in the isocentric plane, in mm; Cluster of spots The centroid coordinates in the isocentric plane, in mm; The lateral scattering half-width of the beam spot corresponding to the energy layer is measured in mm. It is used to quantify the degree of overlap in dose distribution between adjacent beam spot clusters, enabling the model to capture the lateral scattering coupling effect between beam spots within the same energy layer and achieve synergistic optimization of beam spot weights.

[0026] This invention also includes a multi-loss function joint training step, which uses a weighted combination of weighted regression loss, physical consistency dose regularization loss, and dose-volume histogram constraint loss as the total loss function. During training, a temperature annealing strategy is used to gradually approximate the accurate dose-volume histogram constraint. The initial temperature value is set to 5.0, and every 50 training cycles, it is multiplied by a decay factor of 0.5 until it drops to 0.1, ensuring that the predicted beam spot weights simultaneously meet clinical dose requirements and physical consistency, and avoiding weight distributions that do not conform to the physical laws of proton therapy.

[0027] In this invention, the construction of the speckle-voxel isomer bipartite graph specifically involves the following: the feature vector of each speckle node is 7-dimensional, including the speckle nominal energy, the x-axis incident position of the isocenter plane, the y-axis incident position of the isocenter plane, the gantry angle of the field, the treatment bed angle, the half-width of the speckle lateral scattering, and the reference plan weight. All features are directly extracted from the standard DICOMRT-Plan format file of the treatment planning system without additional acquisition or calculation. The feature vector of each voxel node includes the relative stopping power value, the one-heat encoding of the anatomical structure label, the Euclidean distance to the geometric center of the target area, and the water equivalent depth in each field direction. The water equivalent depth is obtained by line integral calculation of the relative stopping power values ​​of all voxels along the incident direction of the field.

[0028] In this invention, the process of grouping voxel nodes according to anatomical structure and generating super voxel nodes through clustering specifically involves grouping voxels according to anatomical structure labels based on the automatic segmentation results. For anatomical structures with more than 500 voxels, the k-means clustering algorithm is used to further subdivide them according to three-dimensional spatial coordinates. The number of clusters K is adaptively determined as the ratio of the number of voxels in the structure to 200, rounded down and not less than 1. Each cluster forms a super voxel node. The feature vector of the super voxel node includes the arithmetic mean and standard deviation of the relative stopping power of voxels in the region, the physical volume of the region, the anatomical structure label, and the Euclidean distance from the geometric center of the region to the center of the target area. The aggregation edge from the beam spot to the super voxel is retained only if the aggregation edge weight is not less than 0.01 times the maximum aggregation edge weight of the beam spot.

[0029] In this invention, the grouping of beam spot nodes by radiation field and energy layer and the generation of beam spot cluster nodes through clustering specifically involves: within the same radiation field, beam spots of the same energy layer are grouped using the DBSCAN clustering algorithm. The neighborhood radius is set to twice the lateral scattering half-width of the beam spot in that energy layer, and the minimum number of samples is set to 1 to ensure that isolated beam spots can also form independent cluster nodes. Each cluster forms a beam spot cluster node. The feature vector of the beam spot cluster node includes the arithmetic mean of the corresponding features of the member beam spots, the number of member beam spots, and the width and height of the bounding rectangle of the spatial distribution of the member beam spots. The vertical interaction edge weight between adjacent energy layer beam spot clusters is calculated as follows: ; in, Cluster of spots with clusters of spots The vertical interaction edge weights between them are dimensionless. For energy The beam pattern at the water equivalent depth The central axis depth dose value, in cGy / MU, is obtained from a lookup table established by offline physical measurements; Cluster of spots The average energy, expressed in MeV; Cluster of spots The average energy, expressed in MeV; It is 1.2 times the larger range of the two depth dose distribution curves, in mm, and is used to quantify the depth dose superposition effect of beam spots in different energy layers, so that the model can capture the coupling relationship of longitudinal dose distribution.

[0030] In this invention, the coarse-grained message passing on the hierarchical compressed graph specifically involves first mapping the original features of the super voxel nodes and cluster nodes to a unified 128-dimensional latent space using independent linear transformation matrices. Then, four rounds of bidirectional message passing are performed. In each round, the cluster nodes first send dose information to the super voxel nodes. An attention mechanism is used to calculate the weight of each edge. The attention score is calculated by concatenating the cluster embedding, the super voxel embedding, and the log-compressed aggregated edge weights using the LeakyReLU activation function. The weighted aggregation yields the update message of the super voxel nodes. Then, the super voxel nodes send dose constraint feedback information to the cluster nodes. The same attention mechanism is used for weighted aggregation to obtain the basic update message of the cluster nodes. Finally, the basic update message and the horizontal and vertical interaction information from adjacent clusters are combined, and the embedding of the cluster nodes is updated through a two-layer feedforward neural network and a layer normalization.

[0031] In this invention, the fine-grained message passing on the original beam spot-voxel isomer bipartite graph specifically involves dividing the beam spots of the entire plan into 3 to 5 independent field groups according to their field affiliation. For each field group, only voxel nodes with valid edges to at least one beam spot within that field are retained. Typically, each group involves 20,000 to 40,000 voxel nodes. Two rounds of fine-grained message passing are performed independently on each field group. The message passing process is consistent with the bidirectional message passing format in the coarse-grained stage but uses an independent set of learnable parameters. The edge features use the original dose contribution value, water equivalent projection depth, and lateral distance, which are logarithmically compressed and then used in attention calculation. The information between each field group has been fully exchanged through the global graph operations in the coarse-grained stage, so there will be no information fragmentation problem.

[0032] In this invention, the multi-loss function joint training step employs an adaptive batch processing strategy. Since the heterogeneous graphs of different patients vary by 3 to 5 times in the number of nodes and edges, the batch size is controlled by the total number of nodes. The upper limit for the total number of nodes in all samples within a batch is set at 50,000 in the coarse-grained stage and 200,000 in the fine-grained stage. Samples are added to the batch one by one until the upper limit is reached. Before training, the graph structure of all samples, including node features, edge indices, and edge features, is pre-calculated and serialized to an HDF5 format disk file. During training, this file is directly loaded to avoid the overhead of online graph construction. The optimizer used is AdamW, with a weight decay coefficient set to 1×10⁻⁶.-4 The learning rate adopts a cosine annealing restart scheduling strategy, with an initial learning rate of 1×10. -3 .

[0033] In this invention, the online inference process specifically involves: first, acquiring the patient's CT image from the current day and obtaining the target area and organ-at-risk contours using an automatic segmentation algorithm, which takes 1 to 3 minutes; converting the CT image into a relative stopping power distribution; and pre-calculating the dose contribution matrix of all beam spots in parallel on the GPU based on a preset field configuration, which takes 10 to 60 seconds; constructing a hierarchical heterogeneous graph and inputting it into a trained neural network model; and outputting the predicted weights of all beam spots through forward inference, which takes less than 1 second; and writing the predicted weights into a standard DICOMRT-Plan file according to the format required by the treatment planning system, which takes less than 1 second. The entire inference process is controlled within 2 minutes after automatic segmentation, meeting the clinical time window of less than 5 minutes required for online adaptive proton therapy.

[0034] This invention presents an end-to-end prediction method for proton therapy speckle weights based on a physics-prior-driven heterogeneous graph neural network. Centered on the physical laws of proton therapy, it constructs a speckle-voxel heterogeneous bipartite graph model. Through a hierarchical graph structure and a two-level message passing mechanism, it achieves a direct mapping from patient CT images to executable speckle weights. The method explicitly incorporates prior knowledge such as the dose contribution of the pen beam algorithm and the lateral and longitudinal physical coupling between specks, completely eliminating the traditional inverse optimization step and significantly shortening the treatment plan generation time. The model employs joint training with multiple loss functions, simultaneously satisfying physical consistency and clinical dose constraints, and is applicable to online adaptive proton therapy for tumors in different locations such as the head and neck and chest. The invention is further described in detail below through two specific embodiments.

[0035] Example 1: This embodiment is applied to an online adaptive proton therapy scenario for a patient with squamous cell carcinoma of the head and neck. The patient's tumor is located in the nasopharynx, invading the left parapharyngeal space, and surrounded by 12 organs at risk, including the brainstem, spinal cord, bilateral parotid glands, and optic nerve. The complex anatomical structure necessitates extremely high dose distribution accuracy. The treatment plan employs three coplanar radiation fields, each containing 30 energy layers, with an energy range of 70 MeV to 220 MeV, totaling approximately 8000 candidate beam spots. The system is deployed on a dedicated computing server in the hospital's radiotherapy department, equipped with a GPU acceleration unit, and interfaces with existing treatment planning systems.

[0036] In the online inference phase, the patient's CT images prior to treatment are first acquired. The image slice thickness is 3mm, the slice interval is 3mm, and the matrix size is 512×512. An automatic segmentation algorithm is used to extract the gross tumor target area, clinical target area, and the contours of various organs at risk. This step takes approximately 2 minutes. The Schneider piecewise linear model is used to convert the HU values ​​of the CT images into a three-dimensional distribution of relative stopping power, accurately characterizing the deceleration ability of different tissues to protons. A preset field configuration file is exported from the treatment planning system, and the parameter information of all candidate speckles is extracted to construct speckle node feature vectors, including the speckle nominal energy, isocentric plane incident position, gantry angle of the field, treatment table angle, speckle lateral scatter half-width, and reference plan weight. All parameters are directly read from the standard DICOM format file without additional calculation.

[0037] A beam-voxel isomer bipartite graph was constructed, and voxel nodes were sampled hierarchically. The original number of voxels in the target region was approximately 60,000. Isotropic downsampling was performed at equal intervals of 2 along the three axes, resulting in approximately 7,500 target region voxel nodes. Sampling quotas were allocated to organs at risk based on their volume as a percentage of the total volume of all organs at risk, with a total quota cap of 100,000. Sampling density was appropriately increased for critical organs at risk such as the brainstem and spinal cord. Normal tissue voxels within a 20mm radius extending from the target region were collected, with a sampling quota of 10,000, for a total of approximately 120,000 voxel nodes sampled throughout the entire plan. The feature vector of each voxel node includes the relative stopping power value, the one-heat encoding of the anatomical structure label, the Euclidean distance to the geometric center of the target region, and the water equivalent depth in each field direction. The water equivalent depth was calculated by line integral of the relative stopping power values ​​of all voxels along the incident direction of the field.

[0038] A single Gaussian approximation-based pen-beam algorithm pre-calculates the dose contribution of each beam spot to each voxel. The dose contribution is obtained by multiplying the depth dose distribution function by the lateral Gaussian scattering kernel. All calculations are performed in parallel on the GPU, with a dose calculation time of approximately 0.5 milliseconds for a single beam spot. The dose contribution matrix is ​​sparsified; directed edges from beam spots to voxels are established only when the dose contribution of a beam spot at a voxel is not less than 0.001 times the maximum dose contribution of that beam spot. The edge features include the dose contribution value, the water equivalent projection depth, and the lateral distance. After sparsification, the non-zero rate of the matrix is ​​reduced to approximately 2%, significantly reducing the subsequent computational load.

[0039] A hierarchical compressed graph structure was constructed, and voxel nodes were grouped according to anatomical structure labels. For anatomical structures with more than 500 voxels, k-means clustering was used to subdivide them according to three-dimensional spatial coordinates, with each cluster containing approximately 200 voxels, forming super voxel nodes. A total of approximately 250 super voxel nodes were generated throughout the project. Their feature vectors included the arithmetic mean and standard deviation of the relative stopping power of voxels within the region, the physical volume of the region, the anatomical structure label, and the Euclidean distance from the geometric center of the region to the target center. Bead nodes were grouped according to the field and energy layer, and the DBSCAN clustering algorithm was used to generate bead cluster nodes. The neighborhood radius was set to twice the half-width of the lateral scattering of the bead in that energy layer, and the minimum number of samples was set to 1 to ensure that isolated bead spots could also form independent cluster nodes. A total of approximately 400 bead cluster nodes were generated throughout the project. Their feature vectors included the arithmetic mean of the corresponding features of member bead spots, the number of member bead spots, and the size of the outer rectangle of the spatial distribution of member bead spots. Convergence edges are established between beam clusters and supervoxels. The convergence edge weight is the sum of the dose contributions of all member voxels within the supervoxel multiplied by the volume of a single voxel. Edges are retained only if the convergence edge weight is not less than 0.01 times the maximum convergence edge weight of the beam cluster. Simultaneously, physical interaction edges are constructed between beam clusters. Lateral interaction edges are established between spatially adjacent beam clusters within the same energy layer, and vertical interaction edges are established between beam clusters in adjacent energy layers with an energy difference not exceeding 5 MeV.

[0040] To perform coarse-grained message passing, the original features of the super-voxel nodes and cluster nodes are first mapped to a unified 128-dimensional latent space using independent linear transformation matrices. Four rounds of bidirectional message passing are then performed. In each round, the cluster node first sends dose information to the super-voxel node. An attention mechanism is used to calculate the weight of each edge. The attention score is obtained by concatenating the cluster embedding, the super-voxel embedding, and the log-compressed aggregated edge weights, and then applying the LeakyReLU activation function. After weighted aggregation of information from all adjacent clusters, the super-voxel node embedding is updated through layer normalization. Subsequently, the super-voxel node sends dose constraint feedback information to the cluster node, and the same attention mechanism is used for weighted aggregation to obtain the basic update message for the cluster node. Finally, the basic update message and the lateral and vertical interaction information from adjacent clusters are combined, processed through a two-layer feedforward neural network, and then the cluster node embedding is updated through layer normalization.

[0041] An embedding unrolling operation is performed, broadcasting the embedding of each speckle cluster node to all its member speckle nodes. After concatenation with the original features of the speckle itself, a linear transformation is performed to generate the initial speckle embedding for the fine-grained stage. Similarly, the embedding of each super voxel node is broadcast to all its member voxel nodes to generate the initial voxel embedding for the fine-grained stage.

[0042] Fine-grained message passing is performed, dividing the planned beamspots into three independent field groups based on their field affiliation, with approximately 2700 beamspot nodes in each group. For each field group, only voxel nodes with valid edges to at least one beamspot within that field are retained, involving approximately 30,000 voxel nodes per group. Two rounds of fine-grained message passing are performed independently for each field group. The message passing process is consistent with the bidirectional message passing in the coarse-grained stage, but uses an independent set of learnable parameters. Edge features use the original dose contribution values, water-equivalent projection depth, and lateral distance, which are logarithmically compressed before being used in attention calculations. Information between field groups has been fully exchanged through global graph operations in the coarse-grained stage, preventing information fragmentation.

[0043] The final embedding of each fascicle node is processed using a three-layer multilayer perceptron output head. The network dimensions are 128, 64, 32, and 1 respectively. The ReLU activation function is used to ensure non-negativity of the output, resulting in the predicted weight value for each fascicle. The predicted weights are written into a standard DICOMRT-Plan file according to the format required by the treatment planning system, which can be directly imported into the treatment machine for irradiation. The entire online inference process takes approximately 1.5 minutes after automatic segmentation, meeting the clinical time window requirements.

[0044] During model training, historical clinical treatment plans of 150 patients with head and neck tumors were used as the training dataset, with each case containing 30 fractionation plans, totaling 4500 training samples. A multi-loss function joint training strategy was adopted, with the total loss function being a weighted combination of weighted regression loss, physical consistency dose regularization loss, and dose-volume histogram constraint loss. A temperature annealing strategy was used during training, with an initial temperature of 5.0, which was multiplied by a decay factor of 0.5 every 50 training epochs, eventually decreasing to 0.1. An adaptive batch processing strategy was employed, with a maximum of 50,000 nodes per batch in the coarse-grained stage and 200,000 nodes in the fine-grained stage. The optimizer used was AdamW, with a weight decay factor set to 1×10⁻⁶. -4 The learning rate adopts a cosine annealing restart scheduling strategy, with an initial learning rate of 1×10. -3 The total training cycles are 400. Before training, the graph structure of all samples is pre-computed and serialized and stored on disk, and then loaded directly during training to improve efficiency.

[0045] This embodiment addresses the complex anatomy and numerous organs at risk associated with head and neck tumors by achieving rapid end-to-end prediction of beam spot weights. The physics-prior-driven graph structure naturally encodes the causal relationships of proton therapy, ensuring the physical consistency of the prediction results. Hierarchical compression and grouped reasoning strategies effectively control computational and memory overhead, while joint training with multiple loss functions ensures that the generated plan simultaneously meets the requirements of target dose coverage and organ protection at risk. The end-to-end design completely eliminates the inverse optimization step, significantly shortening the treatment plan generation time and providing a feasible technical solution for online adaptive proton therapy of head and neck tumors.

[0046] Example 2: This embodiment applies to an online adaptive proton therapy scenario for a patient with locally advanced non-small cell lung cancer. The patient's tumor is located in the upper lobe of the right lung, measuring approximately 4cm × 3cm, and is surrounded by eight organs at risk, including both lungs, heart, spinal cord, and esophagus, exhibiting anatomical changes due to respiratory motion. The treatment plan employs four non-coplanar radiation fields, each containing 25 energy layers, with an energy range of 80MeV to 230MeV, totaling approximately 10,000 candidate beam spots. The system is deployed on a computing cluster within the hospital's radiotherapy department, supporting parallel inference for multiple patients and data integration with a 4D-CT imaging system.

[0047] In the online inference phase, the patient's 4D-CT images from the current day are first acquired, and image data from 10 respiratory phases are collected. The end-expiratory phase is selected as the reference phase for treatment plan generation. An automatic segmentation algorithm is used to extract the contours of the gross tumor target area, clinical target area, and various organs at risk. The influence of respiratory motion on the target area location is considered, and the target area boundary is appropriately expanded. This step takes approximately 2.5 minutes. The Schneider piecewise linear model is used to convert the HU values ​​of the CT images into a three-dimensional distribution of relative stopping power, optimized for the low-density characteristics of lung tissue, and improving the accuracy of water equivalent depth calculation. A preset field configuration file is exported from the treatment planning system, and parameter information of all candidate speckles is extracted to construct speckle node feature vectors.

[0048] A beam-voxel isomer bipartite map was constructed, and voxel nodes were sampled hierarchically. The original number of voxels in the target area was approximately 80,000. Isotropic downsampling was performed at equal intervals of 2 along the three axes, resulting in approximately 10,000 target area voxel nodes. Sampling quotas were allocated to organs at risk based on their volume ratio, with a total quota cap of 100,000. Larger organs at risk, such as the lungs and heart, were allocated higher quotas. Normal tissue voxels within a 20mm radius extending from the target area were collected, with a sampling quota of 10,000, for a total of approximately 130,000 voxel nodes sampled throughout the entire program. The feature vector of each voxel node includes the relative stopping power value, the one-heat encoding of the anatomical structure label, the Euclidean distance to the geometric center of the target area, and the water equivalent depth in each field direction.

[0049] A single Gaussian approximation-based pen-beam algorithm pre-calculates the dose contribution of each beam spot to each voxel. The calculation considers the average effect of respiratory motion and uses the average relative stopping power distribution across multiple respiratory phases for dose calculation. All calculations are performed in parallel on the GPU. The dose contribution matrix is ​​sparsified, retaining only edges whose dose contribution is not less than 0.001 times the maximum dose contribution of that beam spot. Edge features include the dose contribution value, the water equivalent projection depth, and the lateral distance.

[0050] A hierarchical compressed graph structure is constructed, and voxel nodes are grouped according to anatomical structure labels. For anatomical structures with more than 500 voxels, k-means clustering is used to subdivide them according to three-dimensional spatial coordinates, with each cluster containing approximately 200 voxels, forming super voxel nodes. Approximately 300 super voxel nodes are generated throughout the entire project. Bead cluster nodes are grouped according to radiation field and energy layer, and DBSCAN clustering is used to generate bead cluster nodes. The neighborhood radius is set to twice the half-width of the lateral scattering of the bead cluster at that energy layer, and the minimum number of samples is set to 1. Approximately 500 bead cluster nodes are generated throughout the entire project. Aggregation edges are established between bead clusters and super voxels, and edges are retained only if the weight of the aggregation edge is not less than 0.01 times the maximum aggregation edge weight of the bead cluster. Simultaneously, lateral and longitudinal physical interaction edges between bead clusters are constructed.

[0051] To perform coarse-grained message passing, the original features of the supervoxel nodes and cluster nodes are first mapped to a unified 128-dimensional latent space using independent linear transformation matrices. Four rounds of bidirectional message passing are then performed. In each round, the cluster node first sends dose information to the supervoxel, and an attention mechanism is used to calculate edge weights. After weighted aggregation, the supervoxel embedding is updated. Next, the supervoxel sends dose constraint feedback information to the cluster node. Integrating lateral and vertical interaction information, the cluster embedding is updated using a feedforward neural network and layer normalization.

[0052] The embedding unrolling operation is performed, broadcasting the cluster embedding to the member cluster nodes and the super voxel embedding to the member voxel nodes. After concatenation with their respective original features, the initial embedding for the fine-grained stage is generated through linear transformation.

[0053] Fine-grained message passing is performed, dividing the entire plan's speckles into four independent field groups based on their field affiliation, with approximately 2500 speckle nodes in each group. For each field group, only voxel nodes with valid edges to at least one speckle within that field are retained, resulting in approximately 25,000 voxel nodes per group. Two rounds of fine-grained message passing are performed independently for each field group, using an independent set of learnable parameters. Edge features are logarithmically compressed before being used in attention calculations.

[0054] The final embedding of each fascicle node is processed by a three-layer multilayer perceptron output head. The ReLU activation function ensures non-negativity of the output, yielding the predicted weight value for each fascicle. The predicted weights are written into a standard DICOMRT-Plan file, which can be directly imported into the therapy machine for irradiation. The entire online inference process takes approximately 1.8 minutes after automatic segmentation, meeting clinical time window requirements.

[0055] During the model training phase, historical clinical treatment plans of 120 non-small cell lung cancer patients were used as the training dataset, with each case containing 25 fractionation plans, totaling 3000 training samples. The training strategy was consistent with that of the head and neck model, adjusting the weights of the dose-volume histogram constraint loss to suit the characteristics of thoracic tumors, prioritizing the satisfaction of key clinical indicators such as bilateral lung V20, average cardiac dose, and maximum spinal cord dose. Data augmentation techniques were employed during training, randomly translating, rotating, and scaling CT images to improve the model's generalization ability.

[0056] This embodiment addresses the characteristics of thoracic tumors, including respiratory motion, large organ size at risk, and dose sensitivity, achieving rapid and accurate beam weight prediction. A physics-prior-driven model architecture ensures that the dose distribution conforms to the physical laws of proton therapy, while a hierarchical graph structure and grouped inference strategy effectively improve inference speed. Joint training with multiple loss functions ensures that the generated plan minimizes the radiation dose to normal tissues while guaranteeing target area dose coverage. The end-to-end design meets the time requirements of online adaptive therapy, helps reduce the impact of respiratory motion on treatment accuracy, and improves the efficacy of proton therapy for thoracic tumors.

[0057] Reference Figure 1 This figure illustrates the global lifeline of end-to-end prediction of proton therapy beam weights. The entire process is divided into four core stages: data preprocessing, heterogeneous graph construction, two-stage neural network inference, and standard file output. High vigilance is crucial; artifacts or deviations in the acquisition and segmentation of the patient's CT images on the same day can lead to catastrophic consequences such as off-target effects or direct irradiation of organs. The conversion of HU values ​​to relative stopping power and the pre-calculation of the dose contribution matrix form the physical foundation of the model. Without the single Gaussian approximation of the physical prior, the subsequent output of the graph network will have no clinical value. The entire process is strictly constrained to within two minutes after automatic segmentation, which is the critical time window for online adaptive therapy. If insufficient system computing power leads to timeouts, changes in the patient's internal anatomy will render the generated treatment plan completely ineffective, potentially causing serious medical accidents. The generation of the standard DICOMRT-Plan file is the final and most critical output barrier, ensuring that the prediction results can be safely and accurately executed by the treatment machine.

[0058] Reference Figure 2This stage is the core defense mechanism for controlling the explosive increase in computational dimensionality. Without strict downsampling and clustering strategies, the scale of the original voxels and beam spots will instantly overwhelm the GPU memory, causing the online inference system to completely crash. The threshold for the total number of target voxels and the quota allocation for organs at risk must be rigidly enforced. Oversampling will cause computation timeouts, while undersampling will miss critical dose hotspots, leading to target recurrence or necrosis of normal tissue. Using the K-means algorithm for spatial subdivision of voxels based on anatomical structure dimensions, and using DBSCAN for beam spots based on scattering width, is the last line of defense for preserving local anatomical and physical features. If the clustering radius or isolated point parameters are set incorrectly, anatomical boundaries will become blurred, causing organs at risk that should be protected to suffer direct, unintended attacks from high-dose proton beams. Setting an edge preservation threshold severs ineffective weak physical connections, preventing background noise from drowning out the dominant dose contribution.

[0059] Reference Figure 3 This figure illustrates the modeling path of the physical coupling effect between beam spots within the same radiation field. This is a mandatory constraint network to prevent the model from outputting non-physical weight distributions. If the system lacks lateral interaction edges, it will be unable to perceive the lateral scattering overlap of adjacent beam spots within the same energy layer, inevitably causing extremely dangerous high-dose hotspots within the target area, leading to irreversible radiation necrosis of tissues. Vertical interaction edges are responsible for intercepting the risk of uncontrolled Bragg peak broadening and superposition between different energy layers. By extracting the depth-dose curves and range parameters of adjacent energy layers, the model can quantify the depth-dose coupling relationship between beam spots. Once this dimension of information interaction is blocked, the system will output highly fragmented and clinically unenforceable discrete weights, not only increasing the load on the treatment machine and causing hardware failures, but also directly threatening the patient's life safety because the dose distribution cannot be physically realized.

[0060] Reference Figure 4 The message passing framework is the neural center determining the accuracy of dose distribution, and it consists of coarse-grained and fine-grained deductions. If sufficient global information exchange is not achieved in the coarse-grained stage, the model will fall into a local optimum blind spot, leading to dose spillover to distal organs at risk. The attention mechanism established between the supervoxel and the fascicle cluster must accurately capture the main dose contribution paths; otherwise, the shift in attention weights will lead the system to incorrectly suppress the effective beam, resulting in severe underdose of the tumor target area. The subsequent embedding unfolding operation is a crucial action in projecting macroscopic physical constraints onto microscopic nodes. In the fine-grained stage, strict isolation reasoning according to the radiation field is essential, effectively preventing catastrophic interference from different incident direction characteristics. The activation function is forcibly locked to ReLU, an absolute red line in the physical world, to prevent the generation of negative weights; any attempt to bypass this restriction will cause the treatment planning system to crash and malfunction.

[0061] Reference Figure 5The training and optimization process has extremely low fault tolerance, directly determining whether the model can survive in a clinical setting. Simple regression loss is highly deceptive in such physical systems; the system forcibly introduces dose-volume histogram constraints and physical consistency regularization for joint intervention. Without this multi-dimensional suppression, the network will generate extreme treatment plans that are mathematically small but clinically absolutely unacceptable. Temperature annealing strategies must strictly decay with the training cycle; if the initial temperature is too low, gradient explosion will occur, preventing convergence; if it is not reduced to a minimum in the later stages, it will be impossible to approach the true clinical dose threshold. Due to the drastic fluctuations in patient data volume, an adaptive batch processing and disk pre-computation serialization strategy with controlled node count is necessary. Any unauthorized attempt to perform dynamic graph construction online will trigger memory overflow and training interruption, completely destroying the usability of the end-to-end architecture.

[0062] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. An end-to-end prediction method for proton therapy beam spot weights based on a physics-prior-driven heterogeneous graph neural network, characterized in that, Includes the following steps: The system acquires the patient's CT images from the same day, the target area and the outline of organs at risk obtained by automatic segmentation, and the field configuration file exported by the treatment planning system. The Schneider piecewise linear model is used to convert the HU values ​​of the CT images into a three-dimensional distribution of relative stopping power. Construct a speckle-voxel isomer bipartite map, extract all candidate speckles from the field configuration file as speckle nodes, sample voxels from the target area, organs at risk, and normal tissues within 20mm of the target area from CT images as voxel nodes, if the total number of voxels in the target area exceeds 50,000, perform isotropic downsampling at equal intervals of 2 along the three axes, and allocate sampling quotas for organs at risk according to volume ratio. The pen beam algorithm based on the single Gaussian approximation pre-calculates the dose contribution value of each beam spot to each voxel. The dose contribution value is the product of the depth dose distribution function and the transverse Gaussian scattering kernel. A directed edge from the beam spot to the voxel is established only when the dose contribution value is not less than 0.001 times the maximum dose contribution of the beam spot. A hierarchical compressed graph structure is constructed. Voxel nodes are grouped according to their anatomical structure and super voxel nodes are generated through k-means clustering. Bead nodes are grouped according to their radiation field and energy layer and bead cluster nodes are generated through DBSCAN clustering. Aggregation edges are established between bead clusters and super voxels. Four rounds of coarse-grained message passing are performed on the hierarchical compressed graph. Bidirectional message passing is carried out between super voxel nodes and cluster nodes through a heterogeneous graph attention network. At the same time, horizontal and vertical physical interaction information is passed between cluster nodes to update node embeddings. Perform an embedding unrolling operation, broadcast the embedding of the speckle cluster node to all its member speckle nodes, broadcast the embedding of the super voxel node to all its member voxel nodes, and concatenate it with their original features to generate the initial node embedding for the fine-grained stage through a linear transformation. Two rounds of fine-grained message passing are performed on the original speckle-voxel heterogeneous bipartite graph. Independent reasoning is performed by shooting field grouping. Bidirectional message passing is performed between speckle nodes and voxel nodes to update speckle node embeddings. The final embedding of each beamspot node is processed by the output head of a 3-layer multilayer perceptron. The activation function uses ReLU to ensure that the output is non-negative. The predicted weight value of each beamspot is output, and a beamspot weight file in standard DICOMRT-Plan format is generated.

2. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, It also includes the step of constructing physical interaction edges between beam clusters, establishing lateral interaction edges for spatially adjacent beam clusters within the same energy layer, and establishing vertical interaction edges between beam clusters in adjacent energy layers with an energy difference of no more than 5 MeV. ； in, Cluster of spots with clusters of spots Horizontal interaction edge weights between them Cluster of spots The coordinates of the centroid in the plane of equal centers, Cluster of spots The coordinates of the centroid in the plane of equal centers, This is the half-width of the beam spot lateral scattering corresponding to this energy layer.

3. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, It also includes a multi-loss function joint training step, which uses a weighted combination of weighted regression loss, physical consistency dose regularization loss and dose-volume histogram constraint loss as the total loss function. During the training process, a temperature annealing strategy is used to gradually approximate the accurate dose-volume histogram index constraint. The initial temperature value is set to 5.0, and the temperature is multiplied by a decay factor of 0.5 every 50 training cycles until it drops to 0.

1.

4. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, When constructing the beam spot-voxel isomer bipartite map, the beam spot nodes adopt a 7-dimensional feature vector, including the beam spot nominal energy, the x-axis incident position of the isocenter plane, the y-axis incident position of the isocenter plane, the gantry angle of the field, the treatment bed angle, the beam spot lateral scattering half-width, and the reference plan weight. All of the above features are directly extracted from the standard DICOMRT-Plan format file. The voxel node feature vector includes the relative stopping power value, the one-heat encoding of the anatomical structure label, the Euclidean distance to the geometric center of the target area, and the water equivalent depth in each field direction.

5. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, When grouping voxel nodes by anatomical structure and generating super voxel nodes by clustering, voxel groups are divided according to anatomical structure labels based on automatic segmentation results. Anatomical structures whose voxel size reaches a set threshold are subdivided using the k-means clustering algorithm combined with three-dimensional spatial coordinates. The number of clusters is adaptively set according to the number of structural voxels. Each cluster constitutes a super voxel node. The feature vector of the super voxel node integrates the statistical features of the region voxels, the region volume, the anatomical label and the geometric distance features. The aggregation edges between the sclereus and the super voxel are filtered and retained according to the edge weight threshold condition.

6. The end-to-end prediction method for proton therapy beam spot weights based on a physical prior-driven heterogeneous graph neural network according to claim 1, characterized in that, When grouping beam spot nodes by radiation field and energy layer and generating beam spot cluster nodes through clustering, beam spots of the same energy layer within the same radiation field are grouped using the DBSCAN clustering algorithm. Each cluster forms a beam spot cluster node. The feature vector of the beam spot cluster node includes the arithmetic mean of the corresponding features of the member beam spots, the number of member beam spots, and the width and height of the bounding rectangle of the spatial distribution of the member beam spots. The vertical interaction edge weight between beam spot clusters of adjacent energy layers is calculated as follows: ； in, Cluster of spots with clusters of spots Vertical interaction edge weights between them For energy The beam pattern at the water equivalent depth The dose value at the central axis depth, Cluster of spots The average energy, Cluster of spots The average energy, It is 1.2 times the larger range of the two depth dose distribution curves.

7. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, When performing coarse-grained message passing on the hierarchical compressed graph, the original features of super voxel nodes and cluster nodes are mapped to a unified latent space through independent linear transformations, and four rounds of bidirectional message passing are performed. In each round, the cluster nodes first transmit dose information to the super voxel nodes. The attention mechanism is used to calculate the edge weights by combining node embedding and edge weights, and the super voxel nodes are updated by weighted aggregation. Then, the super voxel nodes feed back dose constraint information to the cluster nodes. The same attention mechanism is used to obtain the basic update message of the cluster nodes. Combined with the interaction information of adjacent clusters, the cluster nodes are embedded and updated by feedforward network and layer normalization.

8. The end-to-end prediction method for proton therapy beam weights based on physical prior-driven heterogeneous graph neural networks according to claim 1, characterized in that, When performing fine-grained message passing on the original beam spot-voxel isomorphic bipartite graph, the entire planned beam spot is divided into 3 to 5 independent field groups according to the field affiliation. Each field group retains only voxel nodes that have effective edge associations with the beam spots within the group. Two rounds of fine-grained bidirectional message passing are performed independently on each field group. The original dose contribution value, water equivalent projection depth, and lateral distance are used as edge features and participate in attention calculation after logarithmic compression transformation. Each field group completes information interaction through coarse-grained global graph operations.

9. The end-to-end prediction method for proton therapy beam weights based on a physical prior-driven heterogeneous graph neural network according to claim 3, characterized in that, The multi-loss function joint training step adopts an adaptive batch processing strategy, dynamically controlling the batch size based on the differences in sample graph structure size, setting the upper limit of the total number of batch processing nodes for both the coarse-grained and fine-grained stages, and accumulating each sample to the corresponding node limit to complete batch construction. The graph structure data of all samples is pre-computed and serialized for storage. During the training phase, the data is directly read to avoid the overhead of online graph construction. The AdamW optimizer is used to complete the model training, and the weight decay coefficient and learning rate are configured. The learning rate is dynamically adjusted using a cosine annealing restart scheduling strategy.

10. The end-to-end prediction method for proton therapy beam weights based on a physical prior-driven heterogeneous graph neural network according to claim 1, characterized in that, The online inference process first acquires the patient's CT images from the current day and obtains the target area and organ-at-risk contours through an automatic segmentation algorithm, converting the CT images into relative stopping power distributions; based on the preset field configuration, the dose contribution matrix of all beam spots is pre-calculated in parallel on the GPU, a hierarchical heterogeneous graph is constructed and input into the trained neural network model, and the prediction weights of each beam spot are obtained through forward inference; a standard DICOMRT-Plan file is generated according to the treatment plan system format requirements, and the overall inference process meets the clinical time window requirements.

Images