Gradient enhanced particle swarm solution-based prognosis-guided intensity modulated radiotherapy plan optimization method

Abstract

The invention provides a prognosis-guided intensity modulated radiation therapy plan optimization method based on gradient enhanced particle swarm solution. The method comprises the following steps: S10, constructing a prognosis-based radiation field flux intensity distribution optimization model; s20, based on the radiation field flux intensity distribution optimization model, randomly generating a group of flux distribution initialization particle swarm in a solution space; step S30, carrying out iterative updating on particles in the flux distribution initialization particle swarm by randomly comparing interaction information and gradient information so as to find a global optimal flux; and S40, judging whether a preset termination condition is met or not, if the preset termination condition is not met, turning to the step S30 to continue iteration, and if the preset termination condition is met, outputting the currently searched global optimal flux to obtain an intensity modulated radiotherapy plan. According to the method, the prognosis-guided personalized intensity modulated radiotherapy plan can be efficiently and quickly generated for the patient, and a technical foundation is laid for wide clinical application of prognosis-guided intensity modulated radiotherapy plan design.

Description

Technical Field

[0001] The invention relates to the technical field of radiotherapy, and in particular to a prognosis-guided intensity-modulated radiotherapy plan optimization method using a gradient-enhanced particle swarm solution. Background Art

[0002] The lung cancer prognosis-guided radiotherapy plan design method introduces the prognostic analytical model into the reverse optimization process of the plan, quantitatively represents the potential expected clinical control and risk, and provides a theoretical means to objectively balance the radiation dose of tumors and normal tissues and organs. It is expected to further improve the clinical benefits of lung cancer patients and is the frontier direction of the current development of radiotherapy technology. However, the inclusion of the prognostic model has brought challenges to the reverse optimization solution of the radiotherapy plan: the highly non-convex nature of the prognostic model makes the optimization objective have multi-extreme characteristics, and the conventional gradient algorithm solution under the physical dose optimal strategy cannot provide a global optimal plan, which is not conducive to the clinical treatment of patients.

[0003] In order to meet the clinical limited time window, efficiently finding the global optimal solution is the first problem to be solved in the automatic optimization of prognosis-guided plans. In this regard, in 2018, Modiri et al. first proposed the use of particle swarm optimization (PSO) to solve the non-convex problem of prognosis planning in the article "Individualized estimates of overall survival in radiation therapy plan optimization—A concept study". As a swarm intelligence algorithm, PSO can use group interaction to achieve global optimization of adaptive search, but due to its over-reliance on the best particles in the group, it has low diversity and is prone to premature convergence in high-dimensional space. Its feasibility has only been verified in three-dimensional conformal radiotherapy. In 2023, Yang et al. proposed the random contrastive interaction particle swarm optimization (RCIPSO) in the article "Random Contrastive Interaction for Particle Swarm Optimization in High-Dimensional Environment". As one of the new efficient large-scale swarm intelligence algorithms, it uses a random contrastive interaction strategy to enable the particle swarm to maintain a high diversity in high-dimensional space, thereby exploring the entire solution space to achieve global optimization, and provides a new theoretical method for the optimization of prognosis plans for intensity-modulated radiotherapy (IMRT). However, like other swarm intelligence algorithms, RCIPSO is also a method based on random search, which requires the algorithm to explore the solution space through a large number of function evaluations to ensure the stable convergence of the algorithm, and therefore cannot meet the time efficiency of clinical work.

[0004] In summary, the existing technology still needs to be improved in order to solve the optimization problem of IMRT prognostic guidance planning more quickly and accurately. Summary of the invention

[0005] The purpose of the present invention is to overcome the limitations of the above-mentioned prior art and provide a prognosis-guided intensity-modulated radiotherapy plan optimization method based on gradient enhanced particle swarm algorithm. The large-scale non-convex problem in the prognosis-guided plan optimization is solved by using a particle swarm algorithm based on gradient enhancement. It can efficiently and quickly generate personalized intensity-modulated radiotherapy plans for patients, and maximize the improvement of the expected patient prognosis of the optimized output plan.

[0006] According to a first aspect of the present invention, a prognosis-guided intensity-modulated radiotherapy plan optimization method using a gradient-enhanced particle swarm solution is provided, comprising the following steps: step S10: constructing a prognosis-based field flux intensity distribution optimization model; step S20: based on the field flux intensity distribution optimization model, randomly generating a group of flux distribution initialization particle groups in the solution space; step S30: iteratively updating particles in the flux distribution initialization particle group by randomly comparing interaction information and gradient information to find a global optimal flux; step S40: determining whether a preset termination condition is met; if not, proceeding to step S30 to continue iterating; if the preset termination condition is met, outputting the currently found global optimal flux to obtain an intensity-modulated radiotherapy plan.

[0007] In one implementation, a prognosis-based field flux intensity distribution optimization model is established through the following steps: determining the number of fields and field angles of the radiotherapy plan; based on the determined number of fields and field angles, generating and saving a dose deposition matrix through a dose calculation engine, taking a photon intensity flux map as the optimization solution object, and introducing the optimization parameters into the dose-based optimization model by calculating the field dose distribution; taking the best patient prognosis as the optimization guide, calculating the probability of each prognostic event caused by the radiotherapy dose through the prognostic model, and weighting the probability of each prognostic event to form a prognostic target item; setting a photon flux intensity smoothing target item to improve the smoothness of the field flux intensity distribution; setting dose constraints and dose-volume constraints to ensure coverage and uniformity of the target area dose, while limiting the dose borne by the endangered organs; weighting the prognostic target item and the photon flux intensity smoothing target item to form a total objective function, and combining the dose constraints and dose-volume constraints to form a prognosis-based field flux intensity distribution optimization model.

[0008] In another implementation, in step S10, the field flux intensity distribution optimization model is converted into an unconstrained model using a penalty function method.

[0009] In another implementation, in step S20, based on the field flux intensity distribution optimization model, a group of flux distribution initialization particle groups are randomly generated in the solution space, including: initializing the position of the particle group by randomly generating N different flux vectors in the solution space, and setting the initial velocity of each particle; calculating the fitness of each particle in the initialized particle group and comparing them, and taking the particle with the smallest fitness value as the initial global optimal flux.

[0010] In another implementation, the fitness of each particle in the initialized particle group is calculated by: bringing the position of each particle in the initialized particle group into the field flux intensity distribution optimization model for calculation; and using the calculated function value as the fitness of each particle.

[0011] In another implementation, in step S30, particles in the flux distribution initialized particle group are iteratively updated by randomly comparing interaction information and gradient information to find the global optimal flux, including: each particle in the flux distribution initialized particle group obtains a new particle position by learning the position of the high-quality flux distribution in the entire solution space and the fitness descent direction represented by the local gradient information, wherein the search direction of the particle at time t+1 is jointly determined by the particle inertia, interaction information and gradient enhancement information, and the particle inertia is characterized by the speed of the particle at the last moment, reflecting the ability of the particle to inherit the previous speed; the particle interaction information enables the particle to learn the search direction of excellent individuals in the group and find the area where the optimal flux distribution is located in the global scope, and the gradient enhancement information is the vector of the fastest descent direction of the particle at the current position, which is conducive to strengthening the particle's development of the local area; after all particles are updated, the fitness value of the updated particle at the new position is evaluated; the global optimal flux is updated according to the fitness value of the updated particle at the new position, and the updated global optimal flux is used as the currently found global optimal flux.

[0012] In another implementation, the particle position is updated as follows:

[0013] a. Initialize the flux distribution for each particle X in the particle population i Randomly select TS particles (X j , j≠i) to form an interactive topology and calculate the fitness value of each particle in the topology. j )≤f(X i ), particle X j For X i If there are two or more dominant particles in the current topology, particle X i Will be updated according to the following steps, otherwise, particle X i The individual advantage will be retained and jump directly to the next particle update;

[0014] b. Particle X i Learn the two dominant particles X with the largest fitness difference in the current topology best and X worst The search direction is used as the interactive information, and its expression is as follows:

[0015]

[0016] in, represents the particle interaction information, Used to balance learning intensity, r1 and r2 are random vectors;

[0017] c. Calculate particle X iThe direction of the steepest descent at the current position, that is, the vector along the negative gradient direction at the current position:

[0018]

[0019] in, is the derivative of the objective function with respect to flux intensity, and r is a random vector;

[0020] d. According to particle inertia Interaction Information and gradient boosting information Update the particle's velocity and position:

[0021]

[0022] Among them, c1, c2 and c3 are the weights of the three components.

[0023] According to the second aspect of the present invention, there is provided an application of a prognosis-guided intensity modulated radiotherapy plan optimization method using a gradient enhanced particle swarm solver, which uses the method provided by the first aspect of the present invention to obtain an intensity modulated radiotherapy plan and perform quality control of the intensity modulated radiotherapy plan.

[0024] Compared with the prior art, the present invention has the following beneficial effects:

[0025] (1) The optimization algorithm of the present invention allows particles to simultaneously learn group interaction information and the direction of steepest descent, combining the global optimization capability of the heuristic algorithm and the local search capability of the gradient algorithm, thereby further improving the convergence speed of the algorithm.

[0026] (2) The optimization algorithm of the present invention can effectively solve large-scale non-convex problems and expand the prognosis-guided plan design from the current three-dimensional conformal radiotherapy to intensity-modulated radiotherapy.

[0027] (3) The scheme of the present invention can efficiently and quickly generate a prognosis-guided personalized IMRT plan for patients, maximize the expected patient prognosis of the optimized output plan, and lay a technical foundation for the widespread clinical application of prognosis-guided IMRT plan design. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 This is a flowchart of the steps of the prognosis-guided intensity-modulated radiotherapy plan optimization method using gradient-enhanced particle swarm solving according to an embodiment of the present invention.

[0029] Figure 2 This is a flowchart of the steps of a prognosis-guided intensity-modulated radiotherapy plan optimization method using gradient-enhanced particle swarm solution according to another embodiment of the present invention.

[0030] Figure 3 For Figure 2Correspondingly, a general flow chart of a prognosis-guided intensity-modulated radiotherapy plan optimization method using gradient-enhanced particle swarm optimization according to another embodiment of the present invention.

[0031] Figure 4 This is an optimization flow chart based on a gradient enhanced particle swarm algorithm according to another embodiment of the present invention.

[0032] Figure 5 The iterative convergence curves of 10 plans are optimized for lung cancer patients according to another embodiment of the present invention.

[0033] Figure 6 It is a schematic diagram of the convergence function values ​​of six different algorithm optimization plans and the expected total radiotherapy risks of the obtained plans according to another embodiment of the present invention. DETAILED DESCRIPTION

[0034] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0035] See also Figure 1 A prognosis-guided intensity-modulated radiotherapy plan optimization method using a gradient-enhanced particle swarm solution according to an embodiment of the present invention comprises the following steps:

[0036] Step S10: constructing a prognosis-based radiation field flux intensity distribution optimization model;

[0037] Step S20: Based on the field flux intensity distribution optimization model, a group of flux distribution initialization particle groups are randomly generated in the solution space;

[0038] Step S30: iteratively updating particles in the flux distribution initialization particle group by randomly comparing interaction information and gradient information to find the global optimal flux;

[0039] Step S40: Determine whether the preset termination condition is met. If not, go to step S30 to continue iteration. If the preset termination condition is met, output the currently found global optimal flux to obtain an intensity modulated radiotherapy plan.

[0040] See also Figure 2 The prognosis-guided intensity-modulated radiotherapy plan optimization method using gradient-enhanced particle swarm optimization in an embodiment of the present invention further includes step S50: performing plan quality verification and evaluation.

[0041] Optionally, a prognosis-based field flux intensity distribution optimization model is established through the following steps: determining the number of fields and field angles of the radiotherapy plan; based on the determined number of fields and field angles, generating and saving a dose deposition matrix through a dose calculation engine, taking the photon intensity flux map as the optimization solution object, and introducing the optimization parameters into the dose-based optimization model by calculating the field dose distribution; taking the best patient prognosis as the optimization guide, calculating the probability of each prognostic event caused by the radiotherapy dose through the prognostic model, and weighting the probability of each prognostic event to form a prognostic target item; setting a photon flux intensity smoothing target item to improve the smoothness of the field flux intensity distribution; setting dose constraints and dose-volume constraints to ensure the coverage and uniformity of the target area dose, while limiting the dose borne by the endangered organs; weighting the prognostic target item and the photon flux intensity smoothing target item to form a total objective function, and combining the dose constraints and dose-volume constraints to form a prognosis-based field flux intensity distribution optimization model.

[0042] Optionally, in step S10, the field flux intensity distribution optimization model is converted into an unconstrained model using a penalty function method.

[0043] Optionally, in step S20, based on the field flux intensity distribution optimization model, a group of flux distribution initialization particle groups are randomly generated in the solution space, including: initializing the position of the particle group by randomly generating N different flux vectors in the solution space, and setting the initial velocity of each particle; calculating the fitness of each particle in the initialized particle group and comparing them, and taking the particle with the smallest fitness value as the initial global optimal flux.

[0044] Optionally, the fitness of each particle in the initialized particle group is calculated in the following manner: the position of each particle in the initialized particle group is brought into the optimization model based on the field flux intensity distribution for calculation; and the calculated function value is used as the fitness of each particle.

[0045] Optionally, in step S30, particles in the flux distribution initialized particle group are iteratively updated by randomly comparing interaction information and gradient information to find the global optimal flux, including: each particle in the flux distribution initialized particle group obtains a new particle position by learning the position of the high-quality flux distribution in the entire solution space and the fitness descent direction represented by the local gradient information, wherein the search direction of the particle at time t+1 is jointly determined by the particle inertia, interaction information and gradient enhancement information, and the particle inertia is characterized by the speed of the particle at the previous moment, reflecting the ability of the particle to inherit the previous speed; the particle interaction information enables the particle to learn the search direction of excellent individuals in the group and find the area where the optimal flux distribution is located in the global scope, and the gradient enhancement information is the vector of the fastest descent direction of the particle at the current position, which is conducive to strengthening the particle's development of the local area; after all particles are updated, the fitness value of the updated particle at the new position is evaluated; the global optimal flux is updated according to the fitness value of the updated particle at the new position, and the updated global optimal flux is used as the currently found global optimal flux.

[0046] Optionally, the particle position is updated as follows:

[0047] a. Initialize the flux distribution for each particle X in the particle population i Randomly select TS particles (X j , j≠i) to form an interactive topology and calculate the fitness value of each particle in the topology. j )≤f(X i ) is called particle X j For X i If there are two or more dominant particles in the current topology, particle X i Will be updated according to the following steps, otherwise, particle X i The individual advantage will be retained and jump directly to the next particle update;

[0048] b. Particle X i Learn the two dominant particles X with the largest fitness difference in the current topology best and X worst The search direction is used as the interactive information, and its expression is as follows:

[0049]

[0050] in, represents the particle interaction information, Used to balance learning intensity, r1 and r2 are random vectors;

[0051] c. Calculate particle X i The direction of the steepest descent at the current position, that is, the vector along the negative gradient direction at the current position:

[0052]

[0053] in, is the derivative of the objective function with respect to flux intensity, and r is a random vector;

[0054] d. According to particle inertia Interaction Information and gradient boosting information Update the particle's velocity and position:

[0055]

[0056] Among them, c1, c2 and c3 are the weights of the three components.

[0057] Specifically, the scheme of the present invention is further described according to the following embodiments:

[0058] Example 1

[0059] According to one embodiment of the present invention, a prognosis-guided intensity modulated radiotherapy plan optimization method using a gradient-enhanced particle swarm algorithm is provided. In this embodiment, the patient's best prognosis is used to guide the plan optimization, and a gradient-enhanced particle swarm algorithm is used to solve the large-scale non-convex problem in the optimization, so as to achieve efficient and rapid generation of personalized intensity modulated radiotherapy plans. Figure 3 As shown, the method specifically comprises the following steps:

[0060] S10: Inputting the prognosis-based field flux intensity distribution optimization model into the algorithm process.

[0061] Optionally, a prognosis-based radiation field flux intensity distribution optimization model is established by the following steps:

[0062] S11: Determine the number of radiation fields and their angles in the radiotherapy plan, use the dose calculation engine to generate and save the dose deposition matrix, use the photon intensity flux map as the optimization solution object, and use (to represent the calculated dose distribution) to introduce the optimization parameters into the dose-based optimization model.

[0063] S12: With the best patient prognosis as the optimization guide, the prognostic model is used to calculate the probability of each prognostic event caused by radiotherapy dose, and these prognostic probabilities are weighted to form the prognostic target item.

[0064] S13: Set the photon flux intensity smoothing target item to maximize the smoothness of the optimized flux intensity distribution, which is beneficial to the implementation of radiotherapy equipment and improves the efficiency of plan execution.

[0065] S14: Set dose constraints and dose-volume constraints to ensure coverage and uniformity of target doses, while limiting the dose received by organs at risk to within a clinically acceptable range.

[0066] S15: Weighting each objective function to form a total objective function, and combining it with the constraint function to form an intensity modulated radiotherapy plan optimization model, that is, a radiation field flux intensity distribution optimization model based on prognosis.

[0067] It should be understood that in this embodiment, the overall objective function includes the above-mentioned prognostic objective term and smoothing objective term, the constraint function includes a dose constraint term, a dose-volume constraint term and may further include a flux non-negativity constraint term, etc.; the setting of relevant parameters in the optimization model includes the weight of each prognostic event in the prognostic objective term, the weight of the smoothing objective term, etc., and the basic data of the patient is substituted; the constraint term is processed using a penalty function, and the optimization model is converted into an unconstrained model.

[0068] The optimization model established through steps S11-S15, i.e., the prognosis-based field flux intensity distribution optimization model, can achieve the best prognosis for patients while meeting existing clinical requirements for target area dose coverage and organ-at-risk dose limitation.

[0069] S20: Establish a gradient-enhanced particle swarm algorithm as a plan optimization engine to obtain intensity-modulated radiotherapy plans.

[0070] In this step, a particle swarm algorithm based on gradient enhancement achieves rapid convergence on a global scale by allowing particles to learn both group interaction information and gradient information simultaneously.

[0071] Specifically, the optimization steps of a particle swarm algorithm based on gradient enhancement are as follows:

[0072] S21: Initialize the position of the particle swarm and speed

[0073] Initialize the position of the particle swarm by randomly generating N different flux vectors in the solution space and set the initial velocity of each particle

[0074] S22: Calculate and compare the fitness of the initial particle swarm, and record the particle with the smallest fitness value as the initial global optimal flux G.

[0075] The fitness is the function value calculated by substituting the particle position (i.e., flux vector) into the optimization model. The smaller the fitness, the better the individual prognosis and / or the more favorable the radiotherapy plan executed by the machine can be achieved by the flux intensity corresponding to the particle position.

[0076] S23: Start iterative update, each particle The new particle position is obtained by learning the location of the high-quality flux distribution in the entire solution space and the fitness decrease direction represented by the local gradient information

[0077] In step S23, the search direction of the particle at time t+1 is determined by the particle inertia, interaction information and gradient enhancement information.

[0078] Specifically, particle inertia is represented by the particle's velocity at the previous moment, reflecting the particle's ability to inherit the previous velocity; particle interaction information enables particles to learn the search direction of excellent individuals in the group and find the area with the optimal flux distribution on a global scale; and gradient enhancement information is the vector of the particle's fastest descent direction at its current position, which is conducive to strengthening the particle's development of local areas.

[0079] S24: After all particles are updated, the fitness values ​​of these particles in the new positions are evaluated, and the global optimal flux G is updated.

[0080] S25: Determine whether the preset termination condition is met. If not, go to step S23 and enter the next round of iteration; otherwise, the algorithm ends and the current optimal solution G is used as the final flux distribution output to obtain the intensity-modulated radiotherapy plan.

[0081] The preset termination condition is used to determine whether the optimization algorithm stops. The preset termination condition may be, for example, that the difference between the optimal solutions G between iterations is small enough or that the maximum number of iterations is reached.

[0082] S30: Conduct plan quality verification and evaluation.

[0083] In this step, the plan can be evaluated jointly by physical dose indicators and prognostic indicators.

[0084] Example 2

[0085] In a lung cancer application example, the prognosis-guided intensity-modulated radiotherapy plan optimization method provided by the present invention using a gradient-enhanced particle swarm solution includes the following process: Figure 3 and Figure 4 As shown:

[0086] (1) The prognosis-guided field flux intensity distribution optimization model was input and converted into an unconstrained model using the penalty function method.

[0087] Among them, the input field flux intensity distribution optimization model can be established through the following steps:

[0088] (1.1) The dose deposition matrix A is obtained by using the radiation field information in the original clinical plan and the quadrant infinite beam algorithm (QIB) in the built-in intensity modulated radiotherapy plan design module IMRTP (IMRT Planning) of the radiotherapy research software CERR (Computational Environment for Radiotherapy Research) developed and published by the University of Washington. The photon intensity flux map x is used as the optimization solution object, and d = A x represents the calculated dose distribution.

[0089] (1.2) Optimize the best patient outcome: local control failure (LRF) of lung cancer caused by radiotherapy radiation dose, ≥ grade 3 radiation-induced cardiac events (CE 3+ ), radiation pneumonitis events of grade ≥3 (RILT 3+ ) and grade ≥2 radiation esophagitis events (ARIE 2+ ) are weighted to form a prognostic target item, and the probability of occurrence of each prognostic event is calculated by the corresponding prognostic model. These prognostic models can use trained models or other existing models. Specifically, the function expression corresponding to the prognostic target item is:

[0090]

[0091] Among them, θ s，m is the probability that organ s will induce risk event m after receiving dose d The corresponding weight.

[0092] Specifically, since each weight in the prognostic target item represents the relative importance of each prognosis, it can be jointly determined by the clinician and the patient. Except for 2, the rest of θ R Both are set to 1.

[0093] (1.3) Set the photon flux intensity smoothing objective: Optionally, the smoothness is measured by the second-order difference of the flux in the X and Y directions, and its expression is:

[0094]

[0095] Among them, B represents the field set, h and k represent the resolution of flux in the X direction and Y direction respectively.

[0096] (1.4) Setting dose constraints and dose-volume constraints C: Dose constraints and dose-volume constraints can be given by the NCCN clinical practice guidelines and clinical practice experience to impose dose statistical level constraints on the target area and organs at risk.

[0097] In this application example, the dose constraints and dose-volume constraints mainly set are D 95％ and D max , the V of the heart 50Gy and D mean , V of lung 20Gy and D mean , D of esophagus mean and D max , spinal cord D max ,The specific settings of these constraints can be given based on clinical ,guidelines and physician experience.

[0098] (1.5) Each objective function is weighted to form the total objective function F, and combined with the constraint function to form the intensity modulated radiotherapy plan optimization model, that is, its mathematical expression is:

[0099]

[0100] St

[0101] (1)x≥0.

[0102] (2) d = Ax.

[0103] (3)C(d)≤0.

[0104] Among them, f outcome and f smooth are the prognostic objective term and smooth objective term constructed above, c is the dose constraint and dose-volume constraint function, and γ is the weight corresponding to the smooth term. Optionally, the smooth term weight γ is set to 0.1.

[0105] (2) Establish a particle swarm algorithm based on gradient enhancement.

[0106] (2.1) Initialize the positions of N particles and speed Use uniform sampling to randomly generate N different flux vectors in the solution space. As the position of the initial particle group, set the initial speed of each particle is 0.

[0107] Optionally, for a 1000-dimensional optimization problem, N is set to 900.

[0108] (2.2) Calculate the fitness of the initial particle group and compare them. Record the particle with the smallest fitness value as the initial global optimal flux G: Substitute the particle position into the optimization model and calculate the function value as the fitness. The particle with the smallest fitness indicates that its position can achieve the best individual prognosis in the current group and (or) is more conducive to the radiotherapy plan executed by the machine. It will be recorded and saved.

[0109] (2.3) Start iterative update and update the speed for each particle and location

[0110] (2.3.1) For each particle X i Randomly select TS particles (X j , j≠i) to form an interactive topology and calculate the fitness value of each particle in the topology. j )≤f(X i ) is called particle X j For X i If there are two or more dominant particles in the current topology, particle X i Will be updated according to the following steps; otherwise, particle X i The individual advantage will be retained and jump directly to the next particle update.

[0111] Optionally, TS is set to 25.

[0112] (2.3.2) Particle X i Learn the two dominant particles X with the largest fitness difference in the current topology best and X worst The search direction is used as the interactive information, and its expression is as follows:

[0113]

[0114] In the formula, Used to balance learning intensity, r1 and r2 are random vectors.

[0115] Optionally, It is set to 0.3, and r1 and r2 are randomly sampled from a uniform distribution in [0, 1].

[0116] (2.3.3) Calculate particle X i The direction of the steepest descent at the current position, that is, the vector along the negative gradient direction at the current position:

[0117]

[0118] In the formula, is the derivative of the objective function with respect to flux intensity, and r is a random vector.

[0119] Optionally, the derivative of the objective function with respect to the flux intensity can be obtained by the chain rule, and r is obtained by random sampling from a uniform distribution in [0, 1].

[0120] (2.3.4) According to particle inertia Interaction Information and gradient boosting information Update particle velocity and position:

[0121]

[0122] In the formula, c1, c2 and c3 are the weights of the three components.

[0123] Optionally, c1, c2 and c3 are set to 0.6, 0.6 and 0.8 respectively.

[0124] (2.4) Update the global optimal flux G: After all particles are updated, evaluate the fitness values ​​of these particles at their new positions, and record the position of the particle with the lowest current fitness as the optimal flux G.

[0125] (2.5) Determine whether the preset termination condition is met: If not, go to step (2.3) and enter the next round of iteration; otherwise, the algorithm ends and the current optimal solution G is output as the final flux distribution to obtain the intensity-modulated radiotherapy plan.

[0126] Optionally, the preset termination condition may be preset to satisfy f(x t )-f(x t+5 )≤10 -6 .

[0127] (3) Conduct plan quality verification and evaluation.

[0128] Optionally, the plans are evaluated jointly by physical dose metrics and prognostic metrics.

[0129] In order to further verify the technical effect, the prognosis-guided intensity-modulated radiotherapy plan automatic optimization method based on the gradient-enhanced particle swarm algorithm according to the embodiment of the present invention regenerates the radiotherapy plan for 10 lung cancer IMRT plans and compares them with the clinical plan. The convergence stability of the algorithm is verified by repeatedly optimizing the plan and observing the change range of its convergence curve, see Figure 5 The iterative convergence curves of 10 optimized plans for one patient are shown.

[0130] In addition, the target convergence state, expected total risk of radiotherapy and optimization time of the optimization plan were respectively optimized by comparing the GradRCIPSO algorithm with the active set, sequential quadratic programming (SQP), interior point method, gradient descent method and RCIPSO algorithms, and the synchronization advantage between the solution accuracy and efficiency of the algorithm was jointly evaluated. Figure 6 The values ​​of the convergence functions of the optimizations of six different algorithms and the expected total radiotherapy risk of the resulting plans are presented.

[0131] Furthermore, Table 1 shows the results of the Wilcoxon rank sum test of the dosimetric statistical indicators and expected radiotherapy risks of the two plans for 10 patients. The total risk in the figure is represented by LRF2+ , RILT 3+ , C.E. 3+ and ARIE 2+ The sum of the probabilities of an event occurring.

[0132] Table 1 Comparison of dosimetry and prognosis between the optimized plan and the clinical plan in 10 patients with lung cancer

[0133]

[0134] Analysis of the above results shows that the method of the present invention can successfully generate a prognosis-guided radiotherapy plan for patients. The obtained plan has a comparable dose statistical level to the original clinical plan and the expected total radiotherapy risk is reduced by 13.82%.

[0135] In addition, the algorithm GradRCIPSO of the present invention has stable solution quality and significantly better solution accuracy than the four gradient algorithms (P<0.01). Compared with the swarm intelligence algorithm RCIPSO without gradient enhancement, GradRCIPSO can reduce the average optimization time from 149 minutes to 29 minutes while maintaining similar solution quality (P>0.05), increasing the optimization speed by about 5 times.

[0136] Another aspect of the present invention provides an application of a prognosis-guided intensity modulated radiotherapy plan optimization method using a gradient enhanced particle swarm solver, which uses the method provided in the first aspect of the present invention to obtain an intensity modulated radiotherapy plan and perform quality control of the intensity modulated radiotherapy plan.

[0137] In summary, compared with the prior art, the beneficial effects of the present invention are:

[0138] (1) The optimization algorithm of the present invention allows particles to simultaneously learn group interaction information and the direction of steepest descent, combining the global optimization capability of the heuristic algorithm and the local search capability of the gradient algorithm, thereby further improving the convergence speed of the algorithm.

[0139] (2) The optimization algorithm of the present invention can effectively solve large-scale non-convex problems and expand the prognosis-guided plan design from the current three-dimensional conformal radiotherapy to intensity-modulated radiotherapy.

[0140] (3) The scheme of the present invention can efficiently and quickly generate a prognosis-guided personalized IMRT plan for patients, maximize the expected patient prognosis of the optimized output plan, and lay a technical foundation for the widespread clinical application of prognosis-guided IMRT plan design.

[0141] It should be noted that the present invention may be a method, a system, an apparatus and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for causing a processor to implement various aspects of the present invention.

[0142] Computer readable storage medium can be a tangible device that can hold and store instructions used by an instruction execution device. Computer readable storage medium can be, for example, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination thereof. More specific examples (non-exhaustive list) of computer readable storage medium include: a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a static random access memory (SRAM), a portable compact disk read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanical encoding device, for example, a punch card or a convex structure in a groove on which instructions are stored, and any suitable combination thereof. The computer readable storage medium used here is not interpreted as a transient signal itself, such as a radio wave or other freely propagating electromagnetic wave, an electromagnetic wave propagated by a waveguide or other transmission medium (for example, a light pulse by an optical fiber cable), or an electrical signal transmitted by a wire.

[0143] The computer-readable program instructions described herein can be downloaded from a computer-readable storage medium to each computing / processing device, or downloaded to an external computer or external storage device via a network, such as the Internet, a local area network, a wide area network, and / or a wireless network. The network can include copper transmission cables, optical fiber transmissions, wireless transmissions, routers, firewalls, switches, gateway computers, and / or edge servers. The network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and forwards the computer-readable program instructions for storage in the computer-readable storage medium in each computing / processing device.

[0144] The computer program instructions for performing the operation of the present invention may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages, such as Smalltalk, C++, Python, etc., and conventional procedural programming languages, such as "C" language or similar programming languages. Computer-readable program instructions may be executed entirely on a user's computer, partially on a user's computer, as an independent software package, partially on a user's computer, partially on a remote computer, or entirely on a remote computer or server. In the case of a remote computer, the remote computer may be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or may be connected to an external computer (e.g., using an Internet service provider to connect via the Internet). In some embodiments, an electronic circuit, such as a programmable logic circuit, a field programmable gate array (FPGA), or a programmable logic array (PLA), may be personalized by utilizing the state information of the computer-readable program instructions, and the electronic circuit may execute the computer-readable program instructions, thereby realizing various aspects of the present invention.

[0145] Various aspects of the present invention are described herein with reference to the flow charts and / or block diagrams of the methods, devices (systems) and computer program products according to embodiments of the present invention. It should be understood that each box of the flow chart and / or block diagram and the combination of each box in the flow chart and / or block diagram can be implemented by computer-readable program instructions.

[0146] These computer-readable program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device, thereby producing a machine, so that when these instructions are executed by the processor of the computer or other programmable data processing device, a device that implements the functions / actions specified in one or more boxes in the flowchart and / or block diagram is generated. These computer-readable program instructions can also be stored in a computer-readable storage medium, and these instructions cause the computer, programmable data processing device, and / or other equipment to work in a specific manner, so that the computer-readable medium storing the instructions includes a manufactured product, which includes instructions for implementing various aspects of the functions / actions specified in one or more boxes in the flowchart and / or block diagram.

[0147] Computer-readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device so that a series of operating steps are performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the instructions executed on the computer, other programmable data processing apparatus, or other device to implement the functions / actions specified in one or more boxes in the flowchart and / or block diagram.

[0148] The flow charts and block diagrams in the accompanying drawings show the possible architecture, functions and operations of the systems, methods and computer program products according to multiple embodiments of the present invention. In this regard, each box in the flow chart or block diagram can represent a part of a module, a program segment or an instruction, and a part of the module, a program segment or an instruction contains one or more executable instructions for realizing the specified logical function. In some alternative implementations, the functions marked in the box can also occur in a different order from the order marked in the accompanying drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the block diagram and / or flow chart, and the combination of the boxes in the block diagram and / or flow chart can be implemented by a dedicated hardware-based system that performs the specified function or action, or can be implemented by a combination of dedicated hardware and computer instructions. It is well known to those skilled in the art that it is equivalent to implement it by hardware, implement it by software, and implement it by combining software and hardware.

[0149] Embodiments of the present invention have been described above, and the above description is exemplary, not exhaustive, and is not limited to the disclosed embodiments. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The selection of terms used herein is intended to best explain the principles of the embodiments, practical applications, or technical improvements in the marketplace, or to enable other persons of ordinary skill in the art to understand the embodiments disclosed herein. The scope of the present invention is defined by the appended claims.

Claims

1. A prognosis-guided intensity-modulated radiotherapy plan optimization method using gradient-enhanced particle swarm optimization, characterized in that: The following steps are involved: step S10: Construct a prognosis-based radiation field flux intensity distribution optimization model; Step S20: Based on the field flux intensity distribution optimization model, a group of flux distribution initialization particle groups are randomly generated in the solution space; Step S30: iteratively updating particles in the flux distribution initialization particle group by randomly comparing interaction information and gradient information to find the global optimal flux; Step S40: Determine whether the preset termination condition is met. If not, go to step S30 to continue iteration. If the preset termination condition is met, output the currently found global optimal flux to obtain an intensity modulated radiotherapy plan.

2. The method according to claim 1, characterized in that The following steps are used to establish the prognosis-based radiation field flux intensity distribution optimization model: Determine the number and angle of radiation fields in the radiotherapy plan; Based on the determined number of radiation fields and radiation field angles, the dose deposition matrix is ​​generated and saved through the dose calculation engine. The photon intensity flux map is used as the optimization solution object. The optimization parameters are introduced into the dose-based optimization model by calculating the radiation field dose distribution. Taking the best patient prognosis as the optimization guide, the probability of each prognostic event caused by radiotherapy dose is calculated through the prognostic model, and the probability of each prognostic event is weighted to form the prognostic target item; Set the photon flux intensity smoothing objective item to improve the smoothness of the field flux intensity distribution; Set dose constraints and dose-volume constraints to ensure target dose coverage and uniformity while limiting the dose to organs at risk; The prognosis objective item and the photon flux intensity smoothing objective item are weighted to form the total objective function, and the dose constraint and dose-volume constraint items are combined to form a prognosis-based field flux intensity distribution optimization model.

3. The method according to claim 1, characterized in that In step S10, the field flux intensity distribution optimization model is converted into an unconstrained model using a penalty function method.

4. The method according to claim 1, characterized in that In step S20, based on the field flux intensity distribution optimization model, a group of flux distribution initialization particle groups are randomly generated in the solution space, including: Initialize the position of the particle swarm by randomly generating N different flux vectors in the solution space, and set the initial velocity of each particle; The fitness of each particle in the initialized particle swarm is calculated and compared, and the particle with the smallest fitness value is taken as the initial global optimal flux.

5. The method according to claim 4, characterized in that The fitness of each particle in the initialized particle swarm is calculated in the following way: The position of each particle in the initialized particle group is brought into the optimization model based on the field flux intensity distribution for calculation; The calculated function value is used as the fitness of each particle.

6. The method according to claim 1, characterized in that In step S30, the particles in the flux distribution initialization particle group are iteratively updated by randomly comparing the interaction information and the gradient information to find the global optimal flux, including: Each particle in the flux distribution initialized particle group obtains a new particle position by learning the position of the high-quality flux distribution in the entire solution space and the fitness decline direction represented by the local gradient information. The search direction of the particle at time t+1 is jointly determined by the particle inertia, interaction information and gradient enhancement information. The particle inertia is characterized by the particle's speed at the last moment, reflecting the particle's ability to inherit the previous speed; the particle interaction information enables the particle to learn the search direction of the excellent individuals in the group and find the area with the optimal flux distribution in the global range. The gradient enhancement information is the vector of the fastest descent direction of the particle at the current position, which is conducive to strengthening the particle's development of the local area; After all particles have been updated, evaluate the fitness values ​​of the updated particles in their new positions; The global optimal flux is updated according to the fitness value of the updated particle in the new position, and the updated global optimal flux is used as the currently found global optimal flux.

7. The method according to claim 6, characterized in that The method for updating particle positions is as follows: a. Initialize the flux distribution for each particle X in the particle population i Randomly select TS particles (X j , j≠i) to form an interactive topology and calculate the fitness value of each particle in the topology. j )≤f(X i ) is called particle X j For X i If there are two or more dominant particles in the current topology, particle X i Will be updated according to the following steps, otherwise, particle X i The individual advantage will be retained and jump directly to the next particle update; b. Particle X i Learn the two dominant particles X with the largest fitness difference in the current topology best and X worst The search direction is used as the interactive information, and its expression is as follows: in, represents the particle interaction information, Used to balance learning intensity, r1 and r2 are random vectors; c. Calculate particle X i The direction of the steepest descent at the current position, that is, the vector along the negative gradient direction at the current position: in, is the derivative of the objective function with respect to flux intensity, and r is a random vector; d. According to particle inertia Interaction Information and gradient boosting information Update the particle's velocity and position: Among them, c1, c2 and c3 are the weights of the three components.

8. An application of a prognosis-guided intensity-modulated radiotherapy plan optimization method using a gradient-enhanced particle swarm solution, characterized in that: The method according to any one of claims 1 to 7 is used to obtain an intensity modulated radiotherapy plan and to perform quality control of the intensity modulated radiotherapy plan.

Images