美国数学建模优秀论文03年B题

The Genetic Algorithm-Based Optimization Approach for

Gamma Unit Treatment Planning System

1 Restatement of the Problem

Stereotactic radiosurgery delivers a single high dose of ionizing radiation to a radiographically well-defined, small intracranial 3D brain tumor without delivering any significant fraction of the prescribed dose to the surrounding brain tissue. One of the three often-used modalities is the gamma knife unit. The gamma knife unit delivers a single high dose of ionizing radiation emanating from 201 cobalt-60 unit sources through a heavy helmet. All 201 beams simultaneously intersect at the isocenter, resulting in a spherical (approximately) dose distribution at the effective dose levels. Irradiating the isocenter to deliver dose is termed a “shot”. Shots can be represented as different spheres. Four interchangeable outer collimator helmets with beam channel diameters of 4, 8, 14, and 18 mm are available for irradiating different size volumes. For a target volume larger than one shot, multiple shots can be used to cover the entire target. By combining multiple shots of radiation, the treatment plan can be customized to treat lesions of varying sizes and shapes. In practice, most target volumes are treated with 1 to 15 shots.

Gamma knife treatment plans are conventionally produced using a manual iterative approach. In each iteration, the planner attempts to determine: (1) the number of shots, (2) the shot sizes, (3) the shot locations, and (4) the shot exposed times (weights) that would adequately cover the target and spare critical structures. For large or irregularly shaped treatment volumes, this process becomes rather tedious and time consuming. Also, the quality of the plan produced often depends upon both the patience and the experience of the user. Consequently, a number of researchers have studied techniques for automating the gamma knife treatment planning process, e.g.,

[1, 2]. The algorithms that have been tested include simulated annealing [3, 4], and mixed integer programming, and a nonlinear programming [5-8], etc.

In treatment planning problems, the objective is to deliver a homogeneous (uniform) dose of radiation to the tumor (typically called the target) area while avoiding unnecessary damage to the surrounding tissue and organs. One approach approximates each radiation shot as a sphere [9], thus reducing the problem to one of geometric coverage. Then, a sphere-packing approach [10] can be used to determine the shot locations and sizes. In this paper, considering the physical limitations and biological uncertainties involved in the gamma knife therapy process, we will also formulate the optimal treatment planning for a gamma knife unit as a sphere-packing problem, and propose a reasonable algorithm to find a solution.

美国数学建模优秀论文03年B题

2 Assumptions

To account for all physical limitations and biological uncertainties involved in the gamma knife therapy process, we make several assumptions as follows:

(A1) The shape of the target is not too irregular, and the target volume is a bounded. As a rule of thumb, the target to be treated should be less than 3.5 cm in all dimensions. Its three-dimensional (3D) digital image, usually consisting of millions of points, can be obtained from a CT or MRI.

(A2) Considering the target volume as a 3D grid of points, we divide this grid into two subsets, the subset of points in and out of the target, is denoted as T and N, respectively.

(A3) Four interchangeable outer collimator helmets with beam channel diameters w={4, 8, 14, 18} mm are available for irradiating different size volumes. We use (xs,ys,zs) to denote the coordinates of the center location of the shot, ts,w to denote the time (weight) that each shot is exposed. The total dose delivered is a linear function of the ts,w. For a target volume larger than one shot, multiple shots can be used to cover the entire target. In the gamma knife case, there is a bound of the number, n, of the shots. The typical range of values for n is:

1≤n≤15. (1)

(A4) Neurosurgeons commonly use isodose curves as a means of judging the quality of a treatment plan, and may wish to impose a requirement that the entire target is surrounded by an isodose line of x%, i.e., 30~70%. In our model, we will use an isodose line of 50%, which means that a constraint that the 50% line must surround the target.

(A5) The dose cloud was approximated as a spherically symmetric distribution by averaging the profiles along the x, y, and z axes. Other effects will be ignored in the following optimization formulations.

(A6) The total dose deposited in the target and critical organ should be more than a fraction, P, of the total dose delivered by a plan. Typically, we choose the values for P:

25%≤P≤40%.

3 Optimization Models

3.1 Analysis of the Problem

美国数学建模优秀论文03年B题

The goal of radiosurgery is to deplete tumor cells while preserving normal structures. In general, an optimal treatment plan is designed to meet the following requirements:

(R1) Match specified isodose contours to the target volumes.

(R2) Match specified dose-volume constraints of the target and critical organ. (R3) Constrain dose to specified normal tissue points below tolerance doses. (R4) Minimize the integral dose to the entire volume of normal tissues or organs. (R5) Minimize the dose gradient across the target volume.

(R6) Minimize the maximum dose to critical volumes.

Also, we have the following constraints:

(C1) Prohibit shots from protruding outside the target.

(C2) Prohibit shots from overlapping (to avoid hot spots).

(C3) Cover the target volume with effective dosage as much as possible. But at

least 90% of the target volume must be covered by shots.

(C4) Use as few shots as possible.

In order to satisfy the above requirements and constraints, we can formulate the optimal treatment planning for a gamma knife unit as a sphere-packing problem. The optimization is performed in two steps. First, for the sphere-packing problem, a geometry-based heuristic is used to produce a reasonable configuration of shot number, sizes and locations. Next, a dose-based optimization is used to produce the final treatment plan.

3.2 Geometry-Based Heuristic for the Sphere-Packing Problem

The geometry-based heuristic is designed to quickly produce a reasonable configuration of shot number, sizes and locations. In our geometry-based heuristic, each shot of radiation is modeled as a sphere, and the medial axis transform, known as the “skeleton”, of the target volume is used to guide the placement of the shots. The skeleton is frequently used in shape analysis and other related areas [11-13]. In our case, we will just use the skeleton to quickly find good locations of the shots.

Our geometry-based heuristic is in three stages. First, we generate the skeleton using a 3D skeleton algorithm. Then, we place shots and choose their sizes along the skeleton to maximize a measure of our objective. This process is done by a genetic algorithm (GA)-based shot placement approach. After the number of focusing helmets to be included in the treatment plan, the planning produces a list of the possible helmet combinations and a suggested number of shots to use for the dose-based optimization outlined in Section 3.3.

3.2.1 Skeleton Generation

In our geometry-based heuristic, a 3D skeleton algorithm that follows similar procedures to that in [5] is adopted. Here, we use a morphologic thinning approach

美国数学建模优秀论文03年B题

[14] to create the skeleton as opposed to the Euclidean distance technique. The first step in the skeleton generation is to compute the contour map containing distance information from the point to a nearest target boundary. Then, based on the contour map, several known skeleton extraction methods in the literature [5, 11-14] can be used. Since the method in [5] is simple and fast, we use it for our skeleton generation.

3.2.2 Genetic Algorithm-Based Shot Placement

In the shot placement algorithm, we restrict our attention to points on the skeleton. Our approach is to search a good location of the point to place a shot from all points of the skeleton. We start from a special type of skeleton points, an end point, that help determining the shot size and the location, as shown in Fig. 1, while in [5] uses two special types of skeleton points, i.e., an end point and a cross point. Then, based on the end points of the skeleton, we randomly find a best location to place a shot by using the GA-based shot placement algorithm.

Fig. 1: Examples of the end-points.

We should first give the definitions of an end point [5].

Definition 1. A point is an end-point if: (1) It is in the skeleton; (2) It has only one V-neighbor in the skeleton.

Starting from an end-point, we then look for the best point to place a shot and determine the shot size by using GA [15, 16]. In the GA-based shot placement algorithm, we must solve the following problems:

(1) The encoding method. In general, the bit string (0-1’s) encoding is the most common method adopted by GA researchers because of its simplicity and tractability. However, in this case, if we directly encode the point coordinates (xs,ys,zs) into a bit string, after the crossover and mutation operation, it will result in generating some points that are not in the skeleton. In order to solve this problem, we build a corresponding table to the relationship of the point number and the point coordinates (xs,ys,zs) in the skeleton, as shown in Table 1. As such, we just encode the point number into a bit string. We just select m points from all points of the skeleton to form a population. One point is a chromosome.

Table 1. The relationship of the point number and the point coordinates

美国数学建模优秀论文03年B题

(xs,ys,zs) in the skeleton.

Point number

1

.

.

.

M Point coordinates (xs,ys,zs) (x1,y1,z1) . . . (xM,yM,zM)

(2) Performance evaluation. The key to the GA-based approach is the evaluation function. Ideally, we would like to place shots that cover the entire region, without overdosing within (or outside) of the target. Overdosing occurs outside the target if we choose a shot size that is too large for the current location, and hence the shot protrudes from the target, which is the constraint (C1). Overdosing occurs within the target if we place two shots too close together for their chosen sizes, which is the constraint (C2).

Before defining a fitness function, we give some definitions need to be used as follows:

Definition 2.

(1) Fraction: A target part that is not large enough to be destroyed by the smallest

shot without any harm to the surrounding normal tissue.

(2) Span: The minimum distance between the current location and the end-point at

which we started.

(3) Radium: The approximate Euclidean distance to the target boundary.

We would like to ensure that Span, Radium and w (the shot size) are as close as possible. Therefore, we choose a fitness function that is a sum of squared differences between these three quantities. The fitness function can ensure that the generating Fraction is the smallest after every shot is placed on the target [5].

Fit= s,r(x,y,z)+ s,w(x,y,z,w)+ r,w(x,y,z,w), (2) where,

s,r(x,y,z)=(span(x,y,z) radium(x,y,z))2

s,w(x,y,z,w)=(span(x,y,z) w)2

r,w(x,y,z,w)=(radium(x,y,z) w)2 (3) (4) (5)

The (3) ensures that we pack the target volume as well as possible, that is the current span between shots should be close to the distance to the closest target boundary. The (4) is used to choose a helmet size that fits the skeleton best for the

美国数学建模优秀论文03年B题

current location. The (5) favors a location that is the appropriate distance from the target boundary for the current shot size.

(3) Genetic operators. Based on the encoding method, we develop the genetic operators in the GA. They are crossover and mutation operation.

Crossover / Recombination operation. Crossover / Recombination is a process of exchanging genetic information, which is important for the entire search process. We adopt one-point crossover operation in the GA. The crossover points are randomly set. Mutation operation. Any change in a gene is called a mutation. Here, we used the point mutation.

Then, we proposed the following GA-based shot placement Algorithm:

S1) Find a skeleton using the method in Section 3.2.1, and all the end-points according to the Definition 1. Considering one of the end-points as a start point; S2) Randomly searching all the points in the skeleton using the GA to find the location and size of the best shots as follows:

S2a) Generate randomly m (m=100) points from all the points (e.g., M=1000)

in the skeleton. The m points are the chromosomes. Set the crossover rate

pc=0.95, the mutation rate pm=0.05, the desired Fit, and the number

nof generations g;

S2b) Calculate the Fit (2) of all the m points in the skeleton;

nS2c) If the algorithm has run g steps, or one of the m points satisfies the

desired Fit, the GA stops (at this time, the best shot is chosen);

Else encode m points into the bit strings;

S2d) To do crossover and mutation operation to the m bit strings, go to S2b). S3) Considering the rest of the target in whole as a new target, repeat S1-S2) until the

rest of the target are fractions (At this time, all best shots are found).

3.3 Dose-Based Optimization

After we obtained the number, sizes and the locations of shots, we develop a dose-based optimization method for Gamma Knife treatment plans.

To set up the dose engine, we needed to determine a functional form for the dose delivered at a point (i,j,k), from the shot centered at (xs,ys,zs). The complete dose distribution can be calculated as a sum of contributions from all of the individual shots of radiation. The dose for all (i,j,k) can be computed using [8]:

美国数学建模优秀论文03年B题

D(i,j,k)=

(s,w)∈S×W∑ts,wDs,w(xs,ys,zs,i,j,k)

, (6)

where Ds,w(xs,ys,zs,i,j,k) is the dose delivered to the point (i,j,k) by the shot of size w centered at (xs,ys,zs) with a delivery time of ts,w. Note that Ds,w(xs,ys,zs,i,j,k) is a complicated (non-convex) function, we should find a nonlinear function to approximate it. A better nonlinear function has been noted in the literature to approximate this dose distribution [8]. We therefore used the functional form as follows:

Ds,w(xs,ys,zs,i,j,k)=∑λi(1 ∫i=12x1 ∞e x2dx)

(7) where,

x=t ri

σi,

and λi,γi, and σi are coefficients [8].

In order to meet the requirement (R1), based on assumption (A4), the optimization formulation should impose a constraint on the 50% isodose line that must surround the target. As such, we can impose strict lower and upper bounds on the dose allowed in the target, namely for all (i,j,k)∈T, the dose D(i,j,k) satisfies:

0.5≤D(i,j,k)≤2, (i,j,k)∈T (8)

In order to meet the requirement (R2), based on assumption (A3), during a certain treatment plan, no more than n shots are to be used, so in each card (any section of the target), card({(s,w)∈S×Wts,w>0})≤n. (9)

Note that the value of tolerance doses of the normal tissue points is q,

D(i,j,k)<q, (i,j,k)∈N.

For a treatment plan, the number of shots is no more than 15, so the tolerance doses of a specified normal tissue point should be:

q=15/201=7.46%,

so, we have

0≤D(i,j,k)≤7.46%, (i,j,k)∈N. (10)

美国数学建模优秀论文03年B题

In order to meet the requirement (R3), based on assumption (A6), the tolerance dose radio of the total dose deposited in the target and critical organ to the total dose delivered by a plan is:

∑D(i,j,k)

(i,j,k)∈T

(i,j,k)∈N+T∑D(i,j,k)≥P, P∈[0.25,0.4]. (11)

Also, in order to satisfy the constrain (C3), the gamma knife treatment unit should cover the target volume with effective dosage as much as possible. But at least 90% of the target volume must be covered by shots. We denote:

V: the total volume of target,

Vs: the total effective dosage volume of target, whose dose value at the point is more than 0.5,

f: the effective dosage rate, which satisfies the following inequation:

V90%≤f=s≤100%V. (12)

The exposed time of each shot ts,w should meet:

ts,w≥0. (13)

If we introduce a binary variable δs,w that indicates whether shot s uses width w or not, i.e.,

1, if in one of n shotsδs,w= 0, otherwise .

Moreover, we have the constraints (C1), (C2) and (C4).

Given these constraints (8)-(13), based on the requirement (R4), its goal is to minimize the dose outside of the target. Then, the first optimization model can be formulated:

min

Objective:

s.t.(i,j,k)∈N∑D(i,j,k) (14) 0.5≤D(i,j,k)≤2, (i,j,k)∈T

card({(s,w)∈S×Wts,w>0})≤n

0≤D(i,j,k)≤7.46%, (i,j,k)∈N

美国数学建模优秀论文03年B题

(i,j,k)∈T∑D(i,j,k)

∑D(i,j,k)≥P

, P∈[0.25,0.4] (i,j,k)∈N+T

90%≤f≤100%

ts,w≥0

Also, in order to meet the requirement (R5): minimize the dose gradient across the target volume, the treatment plan needs to be both conformal and homogeneous. It is easy to specify homogeneity in the models simply by imposing lower and upper bounds on the dose delivered to points in the target T. Typically, however, the imposition of rigid bounds leads to plans that are overly homogeneous and not conformal enough, that is, they provide too much dose outside the target. To overcome this problem, the notion of “underdose” (UD) was suggested in [8]. UD measures how much the delivered dose is below the prescribed dose on the target points. In our models, we either constrain the UD to be less than a prespecified value, or attempt to minimize the total UD.

In practical application, rather than calculating the dose at every point, it is easy to accurately estimate the total dose delivered by a plan based solely on the ts,w variables and other pre-calculated constants. An upper bound is also placed on the dose to the target. Given these constraints, the optimizer seeks to minimize the total underdosage in the target. A point is considered to be underdosed if it receives less than the prescribed isodose θ, which for the example formulation is assumed to be 1. We actually use the optimization process to model UD. UD is constrained to be:

UD(i,j,k)=max(0,1 D(i,j,k))

at every point in the target. Provided minimizing UD, we can implement this construct using linear constraints:

θ≤UD(i,j,k)+D(i,j,k), (i,j,k)∈T

0≤UD(i,j,k), (i,j,k)∈T

So, in order to meet the requirement (R5), we have the second complete formulation:

min

Objective: (i,j,k)∈T∑UD(i,j,k) (15)

美国数学建模优秀论文03年B题

s.t.0.5≤D(i,j,k)≤2, (i,j,k)∈T

card({(s,w)∈S×Wts,w>0})≤n

0≤D(i,j,k)≤7.46%, (i,j,k)∈N

(i,j,k)∈T∑D(i,j,k)

∑D(i,j,k)≥P

, P∈[0.25,0.4] (i,j,k)∈N+T

90%≤f≤100%

ts,w≥0

θ≤UD(i,j,k)+D(i,j,k), (i,j,k)∈T

0≤UD(i,j,k), (i,j,k)∈T

And, in order to meet the requirement (R6), we have the third formulation:

Objective: min

s.t.D(i,j,k), (i,j,k)∈T, when δs,w=0 (16) 0.5≤D(i,j,k)≤2, (i,j,k)∈T

card({(s,w)∈S×Wts,w>0})≤n

0≤D(i,j,k)≤7.46%, (i,j,k)∈N

(i,j,k)∈T∑D(i,j,k)

∑D(i,j,k)≥P

, P∈[0.25,0.4] (i,j,k)∈N+T

90%≤f≤100%

ts,w≥0

All the formulations are based on the assumptions that the neurosurgeon can use the volume of the target in order to a priori determine a realistic upper bound n on the number of shots required for treatment.

Several issues need to be resolved to create models that are practical, implementable, and solvable (in a reasonable time frame). Two main approaches are proposed in the literatures [5-8], namely using mixed integer programming and non-linear programming, to simultaneously optimize all of the variables.

4 Simulation Results and Model Testing

美国数学建模优秀论文03年B题

We developed an Optimization Package to implement the algorithms of our models in MATLAB. Also, a numerous of computer simulations have been done using different shape or size targets.

Before simulation, we generate the 3D target volume with the Painting Software Tool and export it to the optimization program. Then, the optimization program can provide various computed results and 2D/3D figures.

We plot a dose volume histogram for the four different helmets using (7) to examine its correctness, as shown in Fig. 2. The histogram depicts the fraction of the volume that receives a particular dose for the target volumes. The fit is best for the small shots and decreases slightly in accuracy for the larger ones. The lines show the fraction of the target and critical organ that receives a particular dosage.

Fig. 2. A dose volume histogram for four different helmets.

Generally speaking, the shape of the target is not too irregular, so we choose five typical shapes of the targets in different sizes. In Fig. 3, we illustrate the maximum section of a typically bean-shaped target, whose maximum size along all dimensions is 35 mm. Using the skeleton generation algorithm in Section 3.2.1, we get the corresponding skeleton as shown in Fig. 4. Then, we apply the GA-based shot placement algorithm, resulting in three shots for the target, one 14mm helmet and two 8mm helmets. The locations and sizes of the helmets in 2D are indicated in Fig. 5, while 3D shot placements are shown in Fig. 6.

美国数学建模优秀论文03年B题

Fig.3. The maximum section of the target. Fig. 4. The skeleton.

Fig.5. The locations and sizes of Fig. 6. The 3D shot placements in the target.

the helmets in 2D.

For this target, we also plot six different isodose lines: 30%, 40%, 50%, 60%, 70%, and 100%, as shown in Fig. 7(a)-(f). In Fig. 7, the thick (red) line is the target outline, and the thin (black) line is the isodose line. In Fig. 7(c), the 50% isodose line covers all the points of the target, while in Fig. 7(f), for the 100% isodose line, no point of the shots exceeds the boundary of the target.

We also have the same simulations to other four shapes of the target, and the 3D shot placements in other four targets are as shown in Fig. 8(a)-(d).

The optimized plans for all of the five shapes of the targets are shown in Table 2. We also obtain the minimum target doses and the percent coverage (by 50%, 80% and 100% isodose) for the five shapes of the targets.

美国数学建模优秀论文03年B题

(a) The 30% isodose. (b) The 40% isodose. (c) The 50% isodose.

(d) The 60% isodose. (e) The 70% isodose. (f) The 100% isodose.

Fig. 7. The specified isodose lines of different values.

Table 2. Optimized plans – Five targets.

Target

Number

1 (Fig. 6)

2 (Fig. 8a)

3 (Fig. 8b) Maximum Section Width 35 mm 26 mm 20 mm Helmet Sizes 18 mm4 mm 8mm 4mm 4mm

8mm

14mm

4mm

4mm Number Minimum of ShotsTarget Dose1 5 4 3 1 1 1 6 2 Percent Coverage (by Isodose ) 9 0%4 (Fig. 8c) 5 (Fig. 8d) 10 mm 8 mm 0.45 0.67 98.9% 100% 82.3% 97.3% 69.4% 56.9%

美国数学建模优秀论文03年B题

(a)

(b)

(c) (d)

Fig. 8. The 3D shot placements in other four targets.

From all of the above simulation results, we know that the geometry-based heuristic with the GA optimization approach is a useful tool for assisting in the selection of the appropriate number of shots and helmet sizes. Also, they indicate that our model exceeds the predefined quality of the treatment planning.

5 Sensitivity Analysis

We are wondering how our model be applied to other conditions.

(1) Does the model have the ability to the sensitive structures ?

In this case, we should apply more dose constraints, such as an upper bound on either the mean dose or the maximum dose to the sensitive structures closing to the target.

(2) Can we treat the tumor at the 30%, 40 %, 60% or 70% isodose level ?

In Fig. 7, with less isodose line, the dose outside of the target volume decreases rapidly, resulting in a reduction in the integral dose to normal brain. While, with

美国数学建模优秀论文03年B题

higher isodose line, the isodose line can not cover all the points of the target. For the nodular regions or the sensitive organs, the higher isodose coverage level should be specified.

(3) Apart from the time considerations, the outstanding question from an optimization viewpoint is the issue of global versus local optimality. How about our model ?

First, the GA-based shot placement algorithm can ensure that the generating Fraction is the smallest one after every shot is placed on the target. For the whole target, it also makes the sum of all Fractions to be the smallest after all the shots are placed on the target.

(4) For the shot placement, when we have several comparably optimization schemes, how to choose the best one from them

?

For example, for a 10 mm diameter sphere target, we have two comparably optimization schemes, as shown in Table 3. The first one is to place one 8 mm shot, and the second one is to place six 4 mm shots. If we only simply consider the treatment as a sphere-packing problem, the better choice is the first one. However, in practical treatment, we should consider the diffuse regions where no shot is irradiated. In these regions, the summing dose value of the point is more than 1, resulting increasing of the total effective dosage. In this case, we will adopt the second one.

Table 3. Two comparably optimization schemes for a 10 mm diameter sphere target.

Scheme Maximum Helmet Number Minimum Percent Coverage (by Number Section Sizes of ShotsTarget Isodose )

Width Dose 100%

100%4.3%69.4%

(5) Will a point outside the tumor whose dose value is more than 1 appear ?

In Fig. 9, though the shots have not protruded outside the target, which is constraint (C1), some points outside the tumor overdose. This occurs due to too irregular shape of the target, which is not avoidable. In this case, we should make a choice how to choose an optimization planning under some constraints.

Fig. 9. The 100% isodose of the target in Fig. 8(a).

美国数学建模优秀论文03年B题

(6) How about the robustness of our model ?

The robustness of the technique is of critical importance. For many cases optimized thus far, despite this somewhat arbitrary starting point, high-quality dose distributions have been obtained in all cases.

6 Strengths and Limitations

6.1 Strengths

(1) Our optimization based automated approach will generate not only more uniform, better treatment plans, but also produce them in less time than is currently used.

(2) The geometry-based approach is based on skeletonization ideas from computational graphics, which can speed the process of shot placement.

(3) The GA-based shot placement algorithm can guide the planner in selecting the number of shots of radiation and the appropriate collimator helmet sizes. Also it can quickly place a shot which provides a significant reduction in the time required for each optimization. Moreover, the GA-based shot placement algorithm can ensure the global and local optimality simultaneously.

(4) The model parameters can be tuned to improve solution speed and robustness. Given the practical patient data, the algorithm can quickly provide the best strategies for practical usage.

(5) The graphical interface is a very intuitive way to demonstrate the isodose curve and the treatment effects of the planning.

6.2 Limitations

There are several problems that our model must deal with.

(1) Our GA-based shot placement approach depends on the skeleton. That is to say, the skeleton is a key factor for the effectiveness of the algorithm. Furthermore, other better methods to find the skeleton should be discussed.

(2) Though our model can be used to different shapes and sizes of targets, whether it fits to too irregular target needs to be further examined.

(3) In our model, we use the function in [8] to approximate the dose calculation. Some other different methods of dose calculation should be examined to improve its accuracy, because the dose calculation will affect our optimization method.

(4) While the mixed integer programming and non-linear programming approaches have proven very effective in practice, there is no guarantee that there is not a better treatment plan. Some more intelligent algorithms, such as neural network-based dynamic programming algorithm could be considered.

(5) We have not tested our model on the practical patient data. Some practical issues should be considered. Further improvement needs to be done so that the model