Gradient enhanced particle swarm optimization for inverse planning of intensity modulated radiation therapy

Abstract

The application provides a gradient-enhanced particle swarm solution-based prognosis-guided intensity modulated radiotherapy plan optimization method, which comprises the following steps: S10, constructing a fluence intensity distribution optimization model based on prognosis; S20, randomly generating a group of fluence distribution initialization particle groups in a solution space based on the fluence intensity distribution optimization model; S30, iteratively updating the particles in the fluence distribution initialization particle groups by randomly comparing interaction information and gradient information, so as to find a global optimal fluence; S40, judging whether a preset termination condition is met, if the preset termination condition is not met, the iteration is continued in step S30, and if the preset termination condition is met, the global optimal fluence currently found is outputted to obtain an intensity modulated radiotherapy plan. The application can efficiently and quickly generate a prognosis-guided individualized intensity modulated radiotherapy plan for a patient, and lays a technical foundation for the wide clinical application of the prognosis-guided intensity modulated plan design.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of radiotherapy, in particular to a gradient-enhanced particle swarm optimization method for prognosis-guided intensity modulated radiotherapy plan optimization. BACKGROUND

[0002] The lung cancer prognosis-guided radiotherapy plan design method introduces a prognosis analysis model into the plan inverse optimization process, quantitatively represents the potential expected clinical control and risk, and provides a theoretical means for objectively balancing the irradiation dose of tumor and normal tissue organs, which is expected to further improve the clinical benefit of lung cancer patients, and is a frontier direction of current radiotherapy technology development. However, the inclusion of the prognosis model brings challenges to the inverse optimization solution of the radiotherapy plan: the highly non-convexity of the prognosis model causes the optimization objective to have multiple extreme values, and the conventional gradient algorithm under the physical dose optimal strategy cannot provide a globally optimal plan, which is not conducive to the clinical treatment of patients.

[0003] Under the condition of meeting the clinical limited time window, efficiently finding the global optimal solution is the first problem of the prognosis guided plan automatic optimization. In this regard, in 2018, Modiri et al. first proposed to use the particle swarm optimization (PSO) algorithm 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 kind of swarm intelligence algorithm, PSO can realize global optimization of adaptive search by using group interaction, but due to the over-reliance on the best particle of the group, it has low diversity and is prone to premature convergence in high-dimensional space. Its feasibility can only be verified in three-dimensional conformal radiotherapy. In 2023, Yang et al. proposed a random contrastive interaction particle swarm optimization (RCIPSO) algorithm 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 can maintain high diversity in high-dimensional space through random contrastive interaction strategy, so as to explore the entire solution space and realize global optimization, providing a new theoretical method for prognosis planning optimization of intensity-modulated radiotherapy (IMRT). However, like other swarm intelligence algorithms, RCIPSO is also a random search-based method, which requires a large number of function evaluations to explore the solution space to ensure the stable convergence of the algorithm, so it cannot meet the time efficiency of clinical work.

[0004] In summary, the existing technology still needs to be improved to solve the optimization problem of IMRT prognosis guided plan more quickly and accurately. SUMMARY

[0005] The purpose of the present application is to overcome the limitations of the above-mentioned prior art and provide a gradient-enhanced particle swarm optimization method for prognosis guided intensity-modulated radiotherapy plan optimization. The gradient-enhanced particle swarm optimization algorithm is used to solve the large-scale non-convex problem in prognosis guided plan optimization, which can efficiently and quickly generate individualized 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 application, a prognosis-guided intensity modulated radiotherapy plan optimization method based on gradient enhanced particle swarm optimization is provided, comprising the following steps: step S10: constructing a fluence intensity distribution optimization model based on prognosis; step S20: generating a group of fluence distribution initialization particle populations in the solution space based on the fluence intensity distribution optimization model; step S30: iteratively updating the particles in the fluence distribution initialization particle populations by randomly comparing the interaction information and the gradient information to find the global optimal fluence; step S40: determining whether the preset termination condition is met, if the preset termination condition is not met, going to step S30 to continue iteration, and if the preset termination condition is met, outputting the current global optimal fluence to obtain the intensity modulated radiotherapy plan.

[0007] In an implementation manner, the fluence intensity distribution optimization model based on prognosis is established by the following steps: determining the number and angle of the fields of the radiotherapy plan; generating a dose deposition matrix by a dose calculation engine based on the determined number and angle of the fields and saving the dose deposition matrix, taking a photon intensity fluence map as an optimization solving object, introducing optimization parameters into a dose-based optimization model by calculating the field dose distribution; taking the best patient prognosis as an optimization guide, calculating the probabilities of each prognosis event caused by the radiotherapy dose by a prognosis model, and weighting the probabilities of each prognosis event to form a prognosis objective term; setting a photon fluence intensity smoothing objective term to improve the smoothness of the fluence intensity distribution; setting a dose constraint and a dose-volume constraint term to ensure the coverage and uniformity of the target dose and limit the dose borne by the organs at risk; weighting the prognosis objective term and the photon fluence intensity smoothing objective term to form a total objective function, and combining the dose constraint and the dose-volume constraint term to form the fluence intensity distribution optimization model based on prognosis.

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

[0009] In another implementation manner, in step S20, generating a group of fluence distribution initialization particle populations in the solution space based on the fluence intensity distribution optimization model comprises: initializing the position of the particle population by randomly generating N different fluence vectors in the solution space, and setting the initial speed of each particle; calculating the fitness of each particle in the initialization particle population and comparing the fitness, and taking the particle with the minimum fitness value as the initial global optimal fluence.

[0010] In another implementation manner, the fitness of each particle in the initialization particle population is calculated by: bringing the position of each particle in the initialization particle population into the fluence intensity distribution optimization model based on prognosis to calculate; and taking the calculated function value as the fitness of each particle.

[0011] In another implementation, in step S30, the particles in the flux distribution initialization particle swarm are iteratively updated by randomly comparing interaction information and gradient information to find the global optimal flux. This includes: each particle in the flux distribution initialization particle swarm learns the position of the best flux distribution in the entire solution space and the fitness descent direction represented by the local gradient information to obtain a new particle position. 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 represented by the particle's velocity at the previous moment, reflecting the particle's ability to inherit the previous velocity. The particle interaction information enables the particle to learn the search direction of the best individuals in the swarm and find the region where the optimal flux distribution is located globally. The gradient enhancement information is the vector of the fastest descent direction of the particle at the current position, which is beneficial to enhance the particle's development of the local region. After all particles have been updated, the fitness value of the updated particles at the new position is evaluated. The global optimal flux is updated according to the evaluated fitness value of the updated particles at the new position, and the updated global optimal flux is taken as the currently found global optimal flux.

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

[0013] a. Initialize each particle X in the particle swarm for the flux distribution. i Randomly select TS particles (X) j (j≠i) form an interactive topology, and calculate the fitness value of each particle in the topology. When f(X) j )≤f(X i When ), particle X is called j For X i The dominant particle, if there are two or more dominant particles in the current topology, particle X i The update will proceed according to the following steps; otherwise, particle X... i Individual advantages will be preserved, and the process will jump directly to the next particle update.

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

[0015]

[0016] in, Represents particle interaction information. To balance the learning intensity, r1 and r2 are random vectors;

[0017] c. Calculate particle X iThe steepest descent direction at the current position, i.e. the vector in the negative gradient direction at the current position:

[0018]

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

[0020] d. According to the particle inertia Interaction information and gradient-enhanced information The velocity and position of the particle are updated as:

[0021]

[0022] wherein, c1, c2 and c3 are weights for measuring the three components.

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

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

[0025] (1) The optimization algorithm of the present application can further improve the convergence speed of the algorithm by allowing the particle to learn the group interaction information and the steepest descent direction simultaneously, combining the global optimization ability of the heuristic algorithm and the local search ability of the gradient algorithm.

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

[0027] (3) The present application can efficiently and quickly generate a prognosis-guided individualized intensity modulated radiotherapy plan for a patient, and maximize the improvement of the expected patient prognosis of the optimized output plan, thereby laying a technical foundation for the wide clinical application of prognosis-guided intensity modulated plan design. BRIEF DESCRIPTION OF DRAWINGS

[0028] Figure 1 The step flowchart of the gradient-enhanced particle swarm solution-based prognosis-guided intensity modulated radiotherapy plan optimization method of the embodiments of the present application.

[0029] Figure 2 The step flowchart of the gradient-enhanced particle swarm solution-based prognosis-guided intensity modulated radiotherapy plan optimization method of another embodiment of the present application.

[0030] Figure 3 The step flowchart of the gradient-enhanced particle swarm solution-based prognosis-guided intensity modulated radiotherapy plan optimization method of another embodiment of the present application. Figure 2The overall flow chart of the gradient enhanced particle swarm solution-based prognosis guided intensity modulated radiotherapy plan optimization method of another embodiment of the application.

[0031] Figure 4 The optimization flow chart based on the gradient enhanced particle swarm algorithm of another embodiment of the application.

[0032] Figure 5 The iteration convergence curve of the 10 optimized plans of the lung cancer patient of another embodiment of the application.

[0033] Figure 6 The convergence function value of the six different algorithm optimized plans and the expected radiotherapy total risk of the obtained plans of another embodiment of the application. DETAILED DESCRIPTION

[0034] The technical solutions in the embodiments of the application will be clearly and completely described below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the protection scope of the application.

[0035] Referring to Figure 1 The gradient enhanced particle swarm solution-based prognosis guided intensity modulated radiotherapy plan optimization method of the embodiment of the application comprises the following steps:

[0036] Step S10: constructing a fluence intensity distribution optimization model based on prognosis;

[0037] Step S20: based on the fluence intensity distribution optimization model, randomly generating a group of fluence distribution initialization particle swarm in the solution space;

[0038] Step S30: iteratively updating the particles in the fluence distribution initialization particle swarm by randomly comparing the interaction information and the gradient information, so as to find the global optimal fluence;

[0039] Step S40: judging whether the preset termination condition is met, if the preset termination condition is not met, going to step S30 to continue iteration, and if the preset termination condition is met, outputting the currently found global optimal fluence to obtain the intensity modulated radiotherapy plan.

[0040] Referring to Figure 2 The gradient enhanced particle swarm solution-based prognosis guided intensity modulated radiotherapy plan optimization method of the embodiment of the application further comprises step S50: performing plan quality verification and evaluation.

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

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

[0043] Optionally, in step S20, based on the fluence intensity distribution optimization model, a group of fluence distribution initialization particle populations are randomly generated in the solution space, including: initializing the position of the particle population by randomly generating N different fluence vectors in the solution space, while setting the initial speed of each particle; calculating the fitness of each particle in the initialization particle population and comparing, taking the particle with the minimum fitness value as the initial global optimal fluence.

[0044] Optionally, the fitness of each particle in the initialization particle population is calculated by: bringing the position of each particle in the initialization particle population into the fluence intensity distribution optimization model based on the fluence intensity distribution optimization model for calculation; taking the calculated function value as the fitness of each particle.

[0045] Optionally, in step S30, the particles in the flux distribution initialization particle swarm are iteratively updated by randomly comparing interaction information and gradient information to find the global optimal flux. This includes: each particle in the flux distribution initialization particle swarm learns the position of the high-quality flux distribution in the entire solution space and the fitness descent direction represented by the local gradient information to obtain a new particle position. 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 represented by the particle's velocity at the previous moment, reflecting the particle's ability to inherit the previous velocity. The particle interaction information enables the particle to learn the search direction of the excellent individuals in the swarm and find the region where the optimal flux distribution is located globally. The gradient enhancement information is the vector of the particle's steepest descent direction at the current position, which is beneficial to enhance the particle's development of the local region. After all particles have been updated, the fitness value of the updated particles at the new position is evaluated. The global optimal flux is updated according to the evaluated fitness value of the updated particles at the new position, and the updated global optimal flux is taken as the currently found global optimal flux.

[0046] Alternatively, the particle position can be updated as follows:

[0047] a. Initialize each particle X in the particle swarm for the flux distribution. i Randomly select TS particles (X) j (j≠i) form an interactive topology, and calculate the fitness value of each particle in the topology. When f(X) j )≤f(X i When ), particle X is called j For X i The dominant particle, if there are two or more dominant particles in the current topology, particle X i The update will proceed according to the following steps; otherwise, particle X... i Individual advantages will be preserved, and the process will jump directly to the next particle update.

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

[0049]

[0050] in, Represents particle interaction information. To balance the learning intensity, r1 and r2 are random vectors;

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

[0052]

[0053] wherein, is the derivative of the objective function with respect to the fluence intensity, r is a random vector;

[0054] d. According to the particle inertia interaction information and gradient enhanced information update the velocity and position of the particle:

[0055]

[0056] wherein, c1, c2 and c3 are weights measuring the three components.

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

[0058] Embodiment 1

[0059] A gradient-enhanced particle swarm solving prognosis-guided intensity modulated radiotherapy plan optimization method is provided according to an embodiment of the application. In this embodiment, the best prognosis-guided plan of the patient is optimized, and a gradient-enhanced particle swarm algorithm is used to solve the large-scale non-convex problem in optimization, so as to realize efficient and fast generation of individualized intensity modulated radiotherapy plan. Referring to the method shown in the figure, the method specifically comprises the following steps: Figure 3

[0060] S10: input the fluence intensity distribution optimization model based on prognosis in the algorithm flow.

[0061] Optionally, the fluence intensity distribution optimization model based on prognosis is established by the following steps:

[0062] S11: determine the number and angle of the radiotherapy plan field, generate the dose deposition matrix by using the dose calculation engine and save it, take the photon intensity fluence map as the optimization solving object, and introduce the optimization parameter into the dose-based optimization model by (which represents the calculation of dose distribution).

[0063] S12: take the best patient prognosis as the optimization guide, calculate the probability of each prognosis event caused by radiotherapy dose by using the prognosis model, and weight these prognosis probabilities to form the prognosis objective term.

[0064] S13: set the photon fluence intensity smoothing objective term to maximize the smoothing degree of the optimized fluence intensity distribution, which is conducive to the implementation of the radiotherapy equipment and improves the plan execution efficiency.

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

[0066] S15: Weight each objective function to form a total objective function, and combine the constraint functions to form an optimization model for IMRT plan, i.e. a fluence intensity distribution optimization model based on prognosis.

[0067] It should be understood that in this embodiment, the total objective function includes the prognosis objective term and the smoothing objective term described above, the constraint function includes the dose constraints and the dose-volume constraints, and can further include a fluence non-negativity constraint term, etc.; the setting of the related parameters in the optimization model includes the weights of each prognosis event in the prognosis objective term, the weight of the smoothing objective term, etc., and are substituted into the basic data of the patient; the constraint terms are processed using a penalty function to convert the optimization model into an unconstrained model.

[0068] The optimization model, i.e. the fluence intensity distribution optimization model based on prognosis, established by the steps S11-S15 can achieve the best prognosis of the patient while meeting the existing clinical requirements for target dose coverage and dose limitation of organs at risk.

[0069] S20: Establish a gradient-enhanced particle swarm algorithm as a plan optimization engine to obtain an IMRT plan.

[0070] In this step, the gradient-enhanced particle swarm algorithm realizes fast convergence in the global range by allowing particles to learn both group interaction information and gradient information.

[0071] Specifically, the optimization steps of the gradient-enhanced particle swarm algorithm are as follows:

[0072] S21: Initialize the position and the speed

[0073] The position of the particle swarm is initialized by randomly generating N different fluence vectors in the solution space, and the initial speed of each particle is set

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

[0075] The fitness is the function value calculated after substituting the particle position (i.e. the fluence vector) into the optimization model. The smaller the fitness, the better the prognosis of the fluence intensity corresponding to the position of the particle and / or the more beneficial the plan for machine execution of radiotherapy.

[0076] S23: Start iteration and update each particle The new particle position is obtained by learning the position of the excellent flux distribution in the entire solution space and the fitness descending direction represented by the local gradient information

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

[0078] Specifically, the particle inertia is represented by the velocity of the particle at the previous moment, which reflects the ability of the particle to inherit the previous velocity; the particle interaction information enables the particle to learn the search direction of the excellent individual in the group and find the region where the optimal flux distribution is located in the global range; and the gradient enhancement information is a vector in the steepest descent direction of the particle at the current position, which is beneficial to strengthen the development of the particle in the local region.

[0079] S24: After all the particles are updated, the fitness values of the particles at 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 iteration; otherwise, the algorithm ends, and the current optimal solution G is output as the flux distribution obtained by solving the intensity modulated radiotherapy plan.

[0081] The preset termination condition is used to determine whether the optimization algorithm is stopped. For example, the preset termination condition is that the optimal solution G between iterations is small enough or reaches the maximum number of iterations.

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

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

[0084] Example 2

[0085] In a lung cancer application example, the prognosis-guided intensity modulated radiotherapy plan optimization method solved by the gradient-enhanced particle swarm provided by the application includes the following processes, combined with Figure 3 and Figure 4 as shown:

[0086] (1) Input the field flux intensity distribution optimization model based on prognosis guidance and convert it into an unconstrained model by the penalty function method.

[0087] The input field flux intensity distribution optimization model can be established by the following steps:

[0088] (1.1) Obtain dose deposition matrix A by using the field information in the original clinical plan and using the quadrant infinite beam algorithm (QIB) in the built-in IMRT planning module IMRTP (IMRT Planning) of the radiotherapy research software CERR (Computational Environment for Radiotherapy Research) developed and disclosed by the University of Washington, USA, and then performing dose calculation. Photon fluence map x is taken as the optimization solution object, and d = A x represents the calculated dose distribution.

[0089] (1.2) Take the best patient prognosis as the optimization guide: the occurrence probabilities of lung cancer local control failure events (LRF) caused by radiotherapy radiation dose, ≥3-grade radiation heart events (CE 3+ ), ≥3-grade radiation pneumonitis events (RILT 3+ ), and ≥2-grade radiation esophagitis events (ARIE 2+ ) are weighted to form a prognosis objective term, and the occurrence probabilities of each prognosis event are calculated by a corresponding prognosis model. These prognosis models can use already trained models or other existing models. Specifically, the function expression of the prognosis objective term is:

[0090]

[0091] Wherein, θ s，m is the weight corresponding to the probability of organ s inducing risk event m after receiving dose d irradiation.

[0092] Specifically, since each weight in the prognosis objective term represents the relative importance between each prognosis, it can be determined by the clinician and the patient together. Alternatively, in addition to θ R being set to 2, the rest θ are set to 1.

[0093] (1.3) Set the photon fluence intensity smoothing objective term: optionally, the smoothing degree is measured by the second-order difference of the fluence in the X and Y directions, and the expression is:

[0094]

[0095] Wherein, B represents the field set, and h and k represent the resolution of the fluence in the X and Y directions, respectively.

[0096] (1.4) Set the dose constraint and dose-volume constraint term C: the dose constraint and dose-volume constraint can be given by the NCCN clinical practice guidelines and clinical practice experience, and the dose statistical level constraints are performed on the target and organs at risk.

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

[0098] (1.5) The weighted sum of each objective function forms the overall objective function F, which, combined with the constraint function, forms the intensity-modulated radiotherapy (IMRT) planning optimization model. 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 Let c represent the prognostic target term and the smoothing target term constructed above, respectively, where c is the dose constraint and dose-volume constraint function, and γ is the weight of the smoothing term. Optionally, the weight γ of the smoothing term 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 N distinct flux vectors are randomly generated in the solution space using a uniform sampling method. The initial position of the particle swarm is used to set the initial velocity of each particle. It is 0.

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

[0108] (2.2) 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: 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 swarm and / or is more conducive to the radiotherapy plan executed by the machine. Record and save it.

[0109] (2.3) Start iteration update, update the velocity of each particle and position

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

[0111] Optionally, TS is set to 25.

[0112] (2.3.2) The particle X i learns the search direction of the two dominant particles X best and X worst with the largest fitness difference in the current topology as interaction information, which is expressed as follows:

[0113]

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

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

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

[0117]

[0118] In the formula, is the derivative of the objective function with respect to the 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 chain rule, and r is randomly sampled in a uniform distribution of [0, 1].

[0120] (2.3.4) Update the particle velocity and position according to the particle inertia interaction information and gradient enhancement information :

[0121]

[0122] wherein c1, c2 and c3 are 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) updating the global optimal flux G: after all particles are updated, the fitness values of the particles at the new positions are evaluated, and the position of the particle with the minimum current fitness value is recorded as the optimal flux G.

[0125] (2.5) judging whether a preset termination condition is met: if not, going to step (2.3) to enter the next round of iteration; otherwise, the algorithm ends, and the current optimal solution G is output as the flux distribution obtained by the final solving, so that an intensity modulated radiotherapy plan is obtained.

[0126] Optionally, the preset termination condition can be preset as meeting f(x t )-f(x t+5 )≤10 -6 .

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

[0128] Optionally, the plan is evaluated by physical dose indicators and prognosis indicators.

[0129] In order to further verify the technical effect, the gradient enhanced particle swarm algorithm-based prognosis guided intensity modulated radiotherapy plan automatic optimization method according to the embodiment of the application is used to regenerate radiotherapy plans for 10 lung cancer IMRT plans and compare them with clinical plans. The convergence stability of the algorithm is verified by repeatedly optimizing the plans and observing the change range of the convergence curves, see Figure 5 An iteration convergence curve of a plan optimized for a patient for 10 times is shown.

[0130] In addition, the target convergence state, expected total radiotherapy risk and optimization time of the plans optimized by the six algorithms, i.e., the gradient enhanced particle swarm algorithm (GradRCIPSO), active set, sequential quadratic programming (SQP), interior point method, gradient descent method and RCIPSO, are compared, so as to evaluate the synchronous advantages of the solving precision and efficiency of the algorithm, see Figure 6 The optimization convergence function values and expected total radiotherapy risks of the plans obtained by the six different algorithms are shown.

[0131] Further, Table 1 shows the results of Wilcoxon rank sum test of the dosimetric statistical indicators and expected radiotherapy risks of the two plans of the 10 patients, and the total risk is represented by LRF2+ RILT 3+ CE 3+ and ARIE 2+ the sum of the probabilities of the events.

[0132] Table 1 Dose statistics and prognosis comparison of optimized plans and clinical plans of 10 lung cancer patients

[0133]

[0134] According to the analysis of the above results, the method can successfully generate a prognosis-guided radiotherapy plan for the patient, and the dose statistics of the obtained plan are equivalent to those of the original clinical plan, and the total risk of radiotherapy is reduced by 13.82%.

[0135] In addition, the algorithm GradRCIPSO has stable solving quality, and the solving precision is significantly better than that of the four gradient algorithms (P<0.01). Compared with the population intelligence algorithm RCIPSO without gradient enhancement, GradRCIPSO can reduce the average optimization time from 149 min to 29 min while having similar solving quality (P>0.05), and the optimization speed is improved by about 5 times.

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

[0137] Compared with the prior art, the application has the following beneficial effects:

[0138] (1) The optimization algorithm of the application scheme combines the global optimization ability of the heuristic algorithm and the local search ability of the gradient algorithm by allowing particles to learn group interaction information and steepest descent direction simultaneously, thereby further improving the convergence speed of the algorithm.

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

[0140] (3) The application scheme can efficiently and quickly generate a prognosis-guided individualized intensity-modulated radiotherapy plan for a patient, and maximally improve the expected prognosis of the optimized output plan, thereby laying a technical foundation for the wide clinical application of prognosis-guided intensity-modulated plan design.

[0141] It should be noted that the application can be a method, system, device and / or computer program product. The computer program product can include a computer readable storage medium having computer readable program instructions loaded thereon for causing a processor to implement various aspects of the application.

[0142] Computer readable storage media can be tangible storage media which can retain and store instructions for use by an instruction execution device. Computer readable storage media can be, for example, but is not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of computer readable storage media include the following: a portable computer diskette, 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 disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission media (e.g., light pulses passing through a fiber-optic cable), or electrical signals transmitted through a wire.

[0143] Computer readable program instructions described herein can be downloaded to respective computing / processing devices from a computer readable storage medium or to an external computer or external storage device via a network, for example, the Internet, a local area network, a wide area network and / or a wireless network. The network can comprise copper transmission cables, optical transmission fibers, wireless transmission, routers, firewalls, switches, gateway computers and / or edge servers. A network adapter card or network interface in each computing / processing device receives computer readable program instructions from the network and forwards the computer readable program instructions for storage in a computer readable storage medium within the respective computing / processing device.

[0144] Computer readable program instructions for carrying out operations of the present application can be assembler instructions, instruction-set-architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state-setting data, or either source code or object code written in any combination of one or more programming languages, including an object oriented programming language such as Smalltalk, C++ or the like and conventional procedural programming languages, such as the "C" programming language or similar programming languages. The computer readable program instructions can execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer can be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection can be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, programmable logic circuitry, field-programmable gate array (FPGA), or programmable logic array (PLA) can execute the computer readable program instructions by utilizing state information of the computer readable program instructions to personalize the electronic circuitry, in order to perform aspects of the present application.

[0145] The computer readable program instructions can also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions / acts specified in the flowchart and / or block diagram block or blocks.

[0146] The computer readable program instructions can also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions / acts specified in the flowchart and / or block diagram block or blocks.

[0147] The computer readable program instructions can also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus or other device to produce a computer implemented process such that the instructions which execute on the computer, other programmable data processing apparatus, or other device implement the functions / acts specified in the flowchart and / or block diagram block or blocks.

[0148] The flow diagrams and the block diagrams in the drawings are presented to illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present application. In this regard, each block in the flow diagrams and the block diagrams can represent a module, segment, or portion of instructions, which comprises one or more executable instructions for implementing the specified logic functions. In some alternative implementations, the functions noted in the blocks can occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks can sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and / or flow diagrams, and combinations thereof, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and

[0149] Embodiments of the application have been described above. The descriptions are intended to be illustrative, and not restrictive, of possible embodiments of the application. Many modifications and variations of the described embodiments are possible, given the benefit of the present disclosure, without departing from the scope and spirit of the described embodiments. The scope of the application is defined by the appended claims.

Claims

1. A gradient-enhanced particle swarm solution guided optimization method for predicting intensity modulated radiotherapy plan optimization, characterized in that, The method comprises the following steps: Step S10: Construct a fluence intensity distribution optimization model based on prognosis; Step S20: Randomly generate a group of fluence distribution initialization particle populations in the solution space based on the fluence intensity distribution optimization model; Step S30: Iteratively update the particles in the fluence distribution initialization particle population by randomly comparing the interaction information and the gradient information to find the global optimal fluence; Step S40: Determine whether the preset termination condition is met, if the preset termination condition is not met, go to step S30 to continue iteration, and if the preset termination condition is met, output the current global optimal fluence to obtain the intensity modulated radiotherapy plan; The fluence intensity distribution optimization model based on prognosis is established by the following steps: Determine the number and angle of the radiation fields of the radiotherapy plan; Based on the determined number and angle of the radiation fields, generate a dose deposition matrix through a dose calculation engine and save it, take a photon intensity fluence map as an optimization solving object, and introduce optimization parameters into a dose-based optimization model by calculating the radiation field dose distribution; Take the best patient prognosis as the optimization guide, calculate the probability of each prognosis event caused by the radiotherapy dose through the prognosis model, and weight the probability of each prognosis event to form a prognosis objective term; Set a photon fluence intensity smoothing objective term to improve the smoothness of the fluence intensity distribution of the radiation field; Set a dose constraint and a dose-volume constraint term to ensure the coverage and uniformity of the target dose and limit the dose borne by the organs at risk; Weight the prognosis objective term and the photon fluence intensity smoothing objective term to form a total objective function, and combine the dose constraint and the dose-volume constraint term to form the fluence intensity distribution optimization model based on prognosis; The expression of the smoothness is: Wherein, B represents a set of radiation fields, h and k represent the resolution of the fluence in the X direction and the Y direction respectively.

2. The method of claim 1, wherein, In step S10, the fluence intensity distribution optimization model is converted into an unconstrained model by using the penalty function method.

3. The method of claim 1, wherein, In step S20, a group of fluence distribution initialization particle populations are randomly generated in the solution space based on the fluence intensity distribution optimization model, which comprises: Initialize the position of the particle population by randomly generating N different fluence vectors in the solution space, and set the initial speed of each particle; Calculate the fitness of each particle in the initialization particle population and compare it, and take the particle with the minimum fitness value as the initial global optimal fluence.

4. The method of claim 3, wherein, The fitness of each particle in the initialization particle population is calculated in the following way: Bring the position of each particle in the initialization particle population into the fluence intensity distribution optimization model based on the radiation field to calculate it; Take the calculated function value as the fitness of each particle.

5. The method of claim 1, wherein, In step S30, the particles in the fluence distribution initialization particle population are iteratively updated by randomly comparing the interaction information and the gradient information to find the global optimal fluence, which comprises: The flux distribution initialization each particle in the particle swarm learns the position of the excellent flux distribution in the whole solution space and the fitness reduction direction represented by the local gradient information, to obtain a new particle position, wherein the search direction of the particle at t+1 moment is determined by the particle inertia, the interaction information and the gradient enhancement information, the particle inertia is represented by the velocity of the particle at the last moment, which embodies the ability of the particle to inherit the previous velocity; the particle interaction information enables the particle to learn the search direction of the excellent individual in the group, to find the area where the optimal flux distribution is located in the global range, and the gradient enhancement information is the vector of the steepest descent direction of the particle at the current position, which is conducive to strengthening the development of the particle in the local area; After all the 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 taken as the current global optimal flux.

6. The method of claim 5, wherein, The method for updating the particle position is as follows: a. Initialize each particle X in the population of flux distribution i Randomly select TS particles X j Construct the interaction topology, j≠i, and calculate the fitness value of each particle in the topology, when f(X j )≤f(X i ), the particle X j is called the dominant particle of X i , if there are two or more dominant particles in the current topology, the particle X i will be updated according to the following steps, otherwise, the particle X i will keep the individual advantage and jump directly to the next particle update; b. particles X i Learning the two dominant particles X with the largest fitness difference in the current topology best and X worst The search direction of X as the interactive information, whose expression is as follows: wherein, represents particle interaction information, for balancing learning intensity, r1 and r2 are random vectors; c. Compute particle X i In the steepest descent direction at the current position, i.e. the vector in the negative gradient direction at the current position: wherein is the derivative of the objective function with respect to the flux intensity, r is a random vector; d. According to the particle inertia Interaction information And gradient enhanced information Update the velocity and position of the particle: Wherein, c1, c2 and c3 are the weights of the three components.

7. Use of a gradient-enhanced particle swarm solution for a prognosis-guided optimization method of intensity modulated radiotherapy planning, characterized in that, The method of any one of claims 1 to 6 is used to obtain the intensity modulated radiotherapy plan, and the quality control of the intensity modulated radiotherapy plan is performed.

Images