Estimation method of dynamic parameters of hepatocellular carcinoma PETCT imaging based on improved MCMC algorithm and its application

Abstract

The present invention relates to a method for estimating the kinetic parameters of hepatocellular carcinoma PET / CT imaging based on an improved MCMC algorithm and its application, belonging to the technical field of liver cancer auxiliary diagnosis. The present invention first obtains prior information and sample information. Secondly, based on Bayesian theory, the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer is derived, and the Markov Chain Monte Carlo (MCMC) method is used to sample the parameter posterior distribution. The sampling algorithm adopted is an improved MCMC algorithm. The global optimization and local optimization of the flower pollination algorithm and the adaptive substitution probability are introduced into the standard MH algorithm to improve the convergence efficiency and accuracy of the MH algorithm, and estimate the kinetic parameters. Finally, the differences in the tracer kinetic parameters between the hepatocellular carcinoma tissue images and the surrounding healthy liver tissue images are analyzed and compared, and the Akaike Information Criterion (AIC) and Root Mean Square Error (RMSE) are used to compare the performance of the algorithm before and after improvement.

Description

Technical Field

[0001] The present invention relates to a method for estimating dynamic parameters of hepatocellular carcinoma PETCT imaging based on an improved MCMC algorithm and an application thereof, and belongs to the technical field of liver cancer auxiliary diagnosis. Background Art

[0002] Primary liver cancer is a highly prevalent and extremely harmful malignant tumor. Its mortality rate ranks third among digestive system malignancies, with approximately 110,000 deaths annually, accounting for 45% of all liver cancer deaths worldwide. Primary liver cancer primarily originates in the epithelial or mesenchymal tissues of the liver. The vast majority of primary liver cancers are hepatocellular carcinomas (HCC). With the development of positron emission tomography (PET) technology, tracer kinetics have become an important tool for studying human functional information. Using tracer kinetic models, it is possible to quantitatively analyze and calculate physiologically significant indicators such as glucose metabolism rate, which plays a crucial role in the study of diseases such as hepatocellular carcinoma. Currently, methods such as the least squares method, commonly used for quantitatively calculating liver tracer kinetic parameters, require the user to pre-set initial values ​​for each parameter. This is difficult for users to accurately set initial values, which directly impacts the reliability of the kinetic modeling results. Therefore, the present invention utilizes an improved MCMC algorithm based on an optimization algorithm to calculate the kinetic model of hepatocellular carcinoma tracers to improve the reliability of parameter estimation. Summary of the Invention

[0003] The present invention provides a method for estimating the kinetic parameters of hepatocellular carcinoma PETCT imaging based on an improved MCMC algorithm and its application, so as to solve the problem of inaccurate estimation of the kinetic parameters of PETCT imaging. The present invention improves the reliability of the estimation of the kinetic parameters of PETCT imaging.

[0004] The technical solution of the present invention is: in a first aspect, a method for estimating dynamic parameters of hepatocellular carcinoma PETCT imaging based on an improved MCMC algorithm, wherein the specific steps of the method are as follows:

[0005] Step 1, obtain prior information and sample information;

[0006] Step 2: Based on Bayesian theory, the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer was derived, and the Markov chain Monte Carlo (MCMC) method was used to sample the posterior distribution of the parameters. The sampling algorithm adopted was the improved MH algorithm to estimate the kinetic parameters.

[0007] As a further solution of the present invention, the improved MH algorithm in Step 2 is specifically as follows:

[0008] (1) Select an initial value θ0 that has physical significance;

[0009] (2) According to the current value θt+1 ,

[0022] , t , , generate a candidate value θ from the proposed distribution q(·|θ t ); * ;

[0010] (3) Calculate the acceptance rate α according to the formula;

[0011]

[0012] (4) Generate a random variable u from the uniform distribution u(0, 1). When α ≥ u, accept the updated θ of the candidate value t+1 = θ * , otherwise θ t+1 = θ t ; p(θ t+1 |y) represents the posterior probability distribution of the kinetic model parameter θ t+1 ; t represents the number of iterations;

[0013] (5) Repeat steps (2) to (4) multiple times and select the optimal value θ best ;

[0014] (6) Introduce the probability conversion formula p(t) of the flower pollination algorithm;

[0015]

[0016] where, R is a random number in the uniform distribution (0, 1), p min = 0.2, p max = 0.8, t is the current number of iterations, T is the maximum number of iterations, Fitness max,t is the maximum fitness value in the t-th iteration process, Fitness min,t is the minimum fitness value in the t-th iteration process; Fitness t represents the fitness value of the current t-th iteration;

[0017] (7) If rand < p(t), update the position using the global optimization formula;

[0018]

[0019] where, take β = 1, s0 = 0.1, Γ is the Gamma function, rand is a random number in the uniform distribution (0, 1);​​​​​​​​​​​​localbest -θ′ t )

[0023] θ localbest is a local optimal value near the current iteration number t;

[0024] (8) Determine whether the updated value is out of bounds. If so, return a random position within the parameter range.

[0025] (9) Determine whether to accept the candidate location update based on the original MH algorithm acceptance probability;

[0026] (10) Determine whether the current termination condition is met. If the condition is met, the algorithm ends; otherwise, a new round of iteration is performed.

[0027] (11) Repeat steps (6) to (10) until convergence, then stop the iteration and output the results.

[0028] In the second aspect, the PETCT imaging kinetic parameters of hepatocellular carcinoma estimated by the method described in the first aspect are used to assist in the diagnosis of liver cancer. When used, the differences in the tracer kinetic parameters between the hepatocellular carcinoma tissue images and the surrounding healthy liver tissue images are analyzed and compared to assist in the diagnosis of liver cancer.

[0029] The beneficial effects of the present invention are as follows: based on Bayesian theory, the present invention derives the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer, uses the Markov Chain Monte Carlo (MCMC) method to sample the parameter posterior distribution, adopts an improved MCMC algorithm as the sampling algorithm, introduces the global optimization and local optimization of the flower pollination algorithm and the adaptive substitution probability into the standard MH algorithm, thereby improving the convergence efficiency and accuracy of the MH algorithm, estimating the kinetic parameters, and improving the reliability of the estimation of the kinetic parameters of PETCT imaging. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 The figure is a flow chart of the overall method of the present invention. DETAILED DESCRIPTION

[0031] Example 1: The patient data used in this invention are as follows: Patient data were provided by the PET / CT Center of the First People's Hospital. They included 31 case data from 18 HCC patients. Written informed consent was obtained from all patients.

[0032] The PET / CT scanning protocol used in this study was as follows: A whole-body CT scan (120 kV, 100 mA) was performed 60 minutes after 18F-fluorodeoxyglucose (FDG) injection, covering 6 to 8 sites, with each site undergoing a 1.5-minute PET scan. All CT images were reconstructed onto a 512 × 512 matrix using a filtered back-projection algorithm. Attenuation correction was performed based on the CT data, and all PET images were reconstructed onto a 200 × 200 matrix using the ordered subset expectation maximization principle. Dynamic PET image acquisition took 5 minutes, and the PET imaging data were reconstructed into 15-second × 4 frames and 60-second × 4 frames using dynamic software.

[0033] like Figure 1 As shown, in the first aspect, a method for estimating dynamic parameters of hepatocellular carcinoma PETCT imaging based on an improved MCMC algorithm is provided, wherein the specific steps of the method are as follows:

[0034] Step 1, obtain prior information and sample information;

[0035] The data processing method used for the prior information and sample information in the present invention is as follows: the time-activity curves of hepatocellular carcinoma, surrounding healthy liver tissue, aorta and portal vein are extracted from 8-frame serial PET / CT images.

[0036] Step 2: Based on Bayesian theory, the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer was derived, and the Markov chain Monte Carlo (MCMC) method was used to sample the posterior distribution of the parameters. The sampling algorithm adopted was the improved MH algorithm to estimate the kinetic parameters.

[0037] The tracer kinetic model used in the derivation of the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer based on the Bayesian theory is as follows:

[0038] 18F-FDG is a drug tracer that has similar properties to glucose in the first two metabolic stages, and the concentration of the radioactive isotope 18F in its molecule can be detected by PET instruments. The three-compartment, four-parameter model, as a kinetic model describing the metabolic process of 18F-FDG, is widely used in the study of tissues such as myocardial and skeletal muscle.

[0039] C B (t) = f A *C a (t)+(1-f A )*C p (t) (1)

[0040]

[0041]

[0042] Solving the differential equation yields the following relationship:

[0043]

[0044] Among them: symbol represents the convolution operation, and

[0045]

[0046]

[0047] Among them, C B (t) represents the total blood 18 The function of F-FDG concentration as a function of time, C a (t) in aortic blood 18 The function of F-FDG concentration as a function of time, C p (t) is the amount of blood in the portal vein 18 The concentration of F-FDG changes as a function of time. C1(t) represents the unphosphorylated F-FDG in the interstitial space or in cells. 18 The function of F-FDG concentration changing with time, C2(t) represents 18 The concentration of FDG-6-P (FDG-6-Phosphate), the product of F-FDG phosphorylation in tissue, changes with time. k1-k4 are metabolic rate constants, and the hepatic arterial blood flow ratio parameter fa is used. Since it is impossible to measure the FDG concentration C1(t) and FDG-6-P concentration in tissue separately, C T (t) represents the function of the change of FDG concentration in the total tissue (liver cancer tissue or normal liver tissue) over time obtained from the PET image.

[0048] Bayesian method: Bayesian theory was first proposed by British scholar Bayes in the mid-18th century. After the 1950s, Bayesian theory gained attention and began to be widely used in statistical decision-making and other problems. Using Bayes' theorem, we can estimate the posterior probability distribution of the parameter based on the prior information of the parameter θ and the likelihood function of the observed data y:

[0049]

[0050] where θ = (k1, k2, k3, k4, fa) are the parameters of the dynamic model, p(y|θ) is the likelihood function of the observed data y, and p(θ) is the prior probability distribution of the parameters θ = (k1, k2, k3, k4, fa).

[0051] Since the problem is complex and the Bayesian formula cannot be directly applied to solve the parameters, the Markov Chain Monte Carlo method (MCMC) is used to couple the Bayesian method and the tracer dynamics three-compartment model. Currently, the MCMC method has become an effective tool for processing complex high-dimensional integral operations. However, as the parameter dimension increases, the posterior distribution becomes more complex, and the sample acceptance rate of the standard MH algorithm will decrease, causing sampling to stagnate. The present invention adopts the MH algorithm improved by the flower pollination optimization algorithm. The overall algorithm process is as follows: Figure 1 shown.

[0052] The improved MH algorithm in Step 2 is specifically as follows:

[0053] (1) Select an initial value θ0 that has physical significance;

[0054] (2) According to the current value θ t , from the proposed distribution q(·|θ t ) generates candidate values ​​θ * ;

[0055] (3) Calculate the acceptance rate α according to the formula;

[0056]

[0057] (4) Generate a random variable u from the uniform distribution u(0,1), and when α≥u, accept the candidate value update θ t+1 =θ * , otherwise θ t+1 =θ t ;p(θ t+1 |y) represents the dynamic model parameter θ t+1 The posterior probability distribution of , t represents the number of iterations;

[0058] (5) Repeat steps (2) to (4) multiple times to select the optimal value θ best ;

[0059] (6) Introducing the probability conversion formula p(t) of the flower pollination algorithm;

[0060]

[0061] Among them, R is a random number in uniform distribution (0,1), p min =0.2,p max =0.8, t is the current number of iterations, T is the maximum number of iterations, Fitness max,t is the maximum fitness value during the tth iteration, min,t is the minimum fitness value in the tth iteration process; Fitnesst denotes the fitness value at the current \(t\)-th iteration;

[0062] (7) If \(rand < p(t)\), update the position using the global optimization formula;

[0063]

[0064] where \(\beta = 1\), \(s_0 = 0.1\), \(\Gamma\) is the Gamma function, and \(rand\) is a random number uniformly distributed on \((0, 1)\);

[0065]

[0066] If \(rand\geq p(t)\), update the position using the local optimization formula;

[0067] \(\theta'\) t+1 =\(\theta'\) t +rand×(\(\theta\) localbest -\(\theta'\) t )

[0068] \(\theta\) localbest is a local optimal value near the current iteration number \(t\);

[0069] (8) Determine whether the updated value is out of bounds. If it is out of bounds, return a random position within the parameter range;

[0070] (9) Determine whether to accept the candidate position update according to the acceptance probability of the original MH algorithm;

[0071] (10) Determine whether the current situation meets the termination condition. If it meets the condition, the algorithm ends; otherwise, a new round of iteration is performed;

[0072] (11) Repeat steps (6) - (10) until convergence, then stop the iteration and output the result.

[0073] In the second aspect, the hepatocellular carcinoma PETCT imaging kinetic parameters estimated by the method described in the first aspect above are applied to assist in the diagnosis of liver cancer. When in use, the differences in the tracer kinetic parameters of the hepatocellular carcinoma tissue image and the surrounding healthy liver tissue image are analyzed and compared to assist in the diagnosis of liver cancer.

[0074] Result analysis: Independent samples t-tests were performed using SPSS, and ROC analysis was performed using Medcalc version 13.0.0.0 software. If \(P < 0.05\), it was considered to have statistical differences.

[0075] Evaluation index:

[0076] AIC = 2K - 2ln(L)

[0077]

[0078] Among them, K is the number of parameters, L is the likelihood function, m is the number of measurements, y i is the measured value, f i is the model value.

[0079] Table 1 shows the results of kinetic parameters for PET / CT imaging of HCC and surrounding normal liver tissue using the unmodified MH algorithm. In HCC images, k1 (P = 0.022), k2 (P < 0.001), and fa (P < 0.01) were significantly higher than those in surrounding normal liver tissue. Compared with surrounding normal liver tissue, k3 (P = 0.059) and k4 (0.218) were increased in HCC images, but this was not statistically significant. HCC is primarily supplied by the hepatic artery rather than the portal vein. This study demonstrated that fa (0.712 ± 0.196 vs. 0.361 ± 0.128) was significantly higher in HCC images than in normal liver tissue. These values ​​are consistent with previous studies of kinetic parameters in liver PET / CT imaging.

[0080] Table 2 shows the results of kinetic parameters for PET / CT imaging of hepatocellular carcinoma and surrounding normal liver tissue using the improved MH algorithm. Kinetic parameters k1 (P = 0.017), k2 (P < 0.001), k3 (P < 0.001), k4 (P = 0.044), and fa (P < 0.001) for hepatocellular carcinoma images were higher than those for normal liver tissue images, consistent with previous studies. Furthermore, differences were observed in all five kinetic parameters. Table 3 shows that the improved MH algorithm achieved lower AIC and RMSE values. Based on the P values ​​and evaluation metrics, the improved MH algorithm outperformed the original MH algorithm.

[0081] Table 1. Kinetic parameter results of PET / CT imaging of hepatocellular carcinoma and normal liver tissue without the improved MH algorithm

[0082]

[0083] Table 2. Kinetic parameter results of PET / CT imaging of hepatocellular carcinoma and normal liver tissue using the improved MH algorithm

[0084]

[0085] Table 3. Comparison of evaluation indicators between the improved MH algorithm and the original MH algorithm

[0086]

[0087] The specific embodiments of the present invention are described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Various changes can be made within the knowledge of ordinary technicians in this field without departing from the scope of the present invention.

Claims

1. A method for estimating dynamic parameters of hepatocellular carcinoma PETCT imaging based on an improved MCMC algorithm, characterized by: The specific steps of the method are as follows: Step1. Obtain prior information and sample information; Step2. Based on Bayesian theory, the posterior distribution of the kinetic parameters of the hepatocellular carcinoma PET / CT imaging tracer is derived. The Markov Chain Monte Carlo (MCMC) method is used to sample the posterior distribution of the parameters. The sampling algorithm adopted is the improved Metropolis-Hastings (MH) algorithm to estimate the kinetic parameters; The specific improved MH algorithm in Step2 is as follows: (1) Select an initial value θ0 with physical meaning; (2) According to the current value θ t , from the proposed distribution q(·|θ t ) generates candidate values ​​θ * ; (3) Calculate the acceptance rate α according to the formula; (4) Generate a random variable u from the uniform distribution u((0,1)). When α≥u, accept the candidate value update θ t+1 =θ * , otherwise θ t+1 =θ t ;p((θ t+1 |y) represents the dynamic model parameter θ t+1 The posterior probability distribution of , t represents the number of iterations; (5) Repeat steps (2) to (4) multiple times to select the optimal value θ best ; (6) Introduce the probability transformation formula p(t) of the flower pollination algorithm; Where R is a random number in uniform distribution (0,1), p min =0.2,p max =0.8, t is the current number of iterations, T is the maximum number of iterations, Fitness max,t is the maximum fitness value during the tth iteration, min,t is the minimum fitness value in the tth iteration process; Fintess t Indicates the fitness value of the current t-th iteration; (7) If rand < p((t), update the position using the global optimization formula; where β = 1, s0 = 0.1, Γ is the Gamma function, and rand is a uniformly distributed random number in (0, 1); If rand ≥ p(t), update the position using the local optimization formula; θ′ t+1 =θ′ t +rand×((θ localbest -θ′ t ) θ localbest is a local optimal value near the current iteration number t; (8) Determine whether the updated value is out of bounds. If it is out of bounds, return a random position within the parameter range; (9) Determine whether to accept the candidate position update according to the acceptance probability of the original MH algorithm; (10) Determine whether the current situation meets the termination condition. If it meets the condition, the algorithm ends; otherwise, a new round of iteration is performed; [[ID= ​

Images