Patients
This study was approved by the Institutional Review Committee (IRB) of the First People’s Hospital of Yun-nan Province (No. 2017YYLH035), informed consent was obtained from all the patients, and all the methods were performed in accordance with the Declaration of Helsinki.
We recruited 28 patients with clinically suspected HCC, and a 5-min dynamic PET/CT scan was added before conventional PET/CT. Ten patients were excluded because they had a non-HCC pathological diagnosis (n = 5), lacked a pathological diagnosis (n = 3), and had suboptimal imaging quality (n = 2).
Eighteen patients (17 males and 1 female) who had pathologically confirmed HCC were finally included in this study. Sixteen patients had a single lesion, and two patients had two lesions. A total of 20 HCC tumors that were confirmed by surgery (n = 14) or biopsy (n = 6) were used in this study, and the long axis of these tumors was 1.9–15.0 cm (average 6.5 ± 3.6).
Dynamic PET/CT
All examinations were performed using a Philips Ingenuity TF PET/CT scanner (Cleveland, OH, USA), and a Philips IntelliSpace Portal v7.0.4.20175 was used for postprocessing. In summary, after the patients had fasted for at least 6 h, blood glucose was verified. A low-dose liver CT scan (120 kV, 100 mAs) was performed for attenuation correction and image fusion. A 5-min dynamic PET scan was performed over the liver region after intravenous administration of 5.5 MBq/kg 18F-FDG. Dynamic PET data were divided into 16 frames using the following sampling schedule: 12 frames of 5 s and 4 frames of 60 s each. Dynamic PET images were then reconstructed using the standard ordered subsets expectation maximization (OSEM) algorithm.
Regions of interest (ROIs) were drawn manually in the CT images of each patient, including the HCC, background liver tissue, aorta and portal vein (the extrahepatic portal vein rather than the intrahepatic portal vein), with manual slice-by-slice adjustment. An ROI was copied to the PET/CT images after image fusion, and time activity curves (TACs) consisting of the maximum SUV (SUVmax) extracted from each frame were generated. Figure 1 shows the ROIs drawn on a transaxial dynamic PET/CT image of a patient with HCC.
Kinetic modeling
A dual-input three-compartment model was used to assess the steady-state hepatic metabolism of 18F-FDG, as shown in Fig. 2 [17, 18].
In Fig. 2, k1 (ml/min/ml) represents the rate constant of 18F-FDG from the blood to the liver tissue, and k2 represents the clearance rate back to the blood. k3 is the rate constant of further phosphorylation of 18F-FDG to 18F-FDG-6-phosphate and k4 is the dephosphorylation rate of phosphatase. CB(t) represents the 18F-FDG concentration in the blood:
$${C}_{B}\left(t\right)=HPI\times A\left(t\right)+\left(1-HPI\right)\times P\left(\mathrm{t}\right),$$
(1)
where A(t) represents the 18F-FDG concentration in the hepatic artery, and P(t) represents the 18F-FDG concentration in the portal vein. The HPI represents the hepatic artery perfusion index (the ratio of arterial blood volume to total blood volume). CM(t) represents the free-state 18F-FDG concentration and the metabolized 18F-FDG-6-phosphate concentration in the liver tissue compartment. CT(t) represents the curve of the tracer concentration in the tissue measured from the PET image over time and is the output function of the kinetic model:
$${C}_{T}\left(\mathrm{t}\right)=\frac{{K}_{1}}{{\alpha }_{2}-{\alpha }_{1}}\times \left[\left({k}_{3}+{k}_{4}-{\alpha }_{1}\right){e}^{-{\alpha }_{1}t}+\left({\alpha }_{2}-{k}_{3}-{k}_{4}\right){e}^{-{\alpha }_{2}t}\right]\otimes{C}_{B}\left(t\right)$$
(2)
where \({\alpha}_{2} \, {\text{and}} \, {\alpha}_{1}\) can be described as follows:
$${\alpha }_{1}=\frac{{k}_{2}+{k}_{3}+{k}_{4}-\sqrt{{({k}_{2}+{k}_{3}+{k}_{4})}^{2}-4{k}_{2}\times {k}_{4}}}{2}$$
(3)
$${\alpha }_{2}=\frac{{k}_{2}+{k}_{3}+{k}_{4}+\sqrt{{({k}_{2}+{k}_{3}+{k}_{4})}^{2}-4{k}_{2}\times {k}_{4}}}{2}$$
(4)
Estimation of the kinetic parameters
This study proposes an improved algorithm, the DCGSA, for the estimation of liver kinetic parameters. In the GSA, a solution space about the objective function is initialized randomly, with each individual in the space as one feasible solution in the objective function. With constant movement, individuals will move toward the most mass, which is the optimal solution of the search space [16]. The convergence speed of the GSA is faster, which makes it fall into a local optimum without global search [19].
In the DCGSA, the dynamic adjustment strategy is introduced for the gravitational constant of GSA to improve the algorithm exploration capacity and mining capacity. At the same time, inertia weights and chaotic sequences are added to the particle speed update process to avoid falling into a local optimum and to improve the global search ability. The overall flow of the DCGSA is shown in Fig. 3. The DCGSA can be simply divided into four parts: the initialization phase, the evaluation phase, the acceleration phase, and the update position phase.
Initialization phase
To start the search process for the DCGSA, the initial population of N individuals is randomly generated in the search space, which represents a set of model parameters (k1, k2, k3, k4 and the HPI). The position of each particle is represented by the following:
$${X}_{i} =\left({x}_{i}^{1}, {x}_{i}^{2},\cdots ,{x}_{i}^{d},\cdots ,{x}_{i}^{D}\right)\quad \text{for}\,\, i=1,\cdots ,N$$
(5)
where \({x}_{i}^{D}\) is the position of the ith individual in the D dimension.
Evaluation phase
The individual is evaluated according to the fitness function, defined as the sum of the squared errors squared errors of the experimental data and the fitted data:
$${fit}_{i}\left(t\right)= {\sum }_{i=1}^{N}{\left( {{C}_{T}}^{^{\prime}}\left(t\right)- {C}_{i}^{m}\left(t\right) \right)}^{2},$$
(6)
where N is the number of individuals, i is the index, \({{C}}{(}{{t}}{)}\) is the actual concentration of 18F-FDG in the obtained tissue, and \({{C}}_{{i}}^{{m}}{(}{{t}}{)}\) is the estimated concentration of 18F-FDG in the obtained tissue. After that, the individual continuously updates its inertial mass \({M}_{i}\left(t\right)\) during movement by the following equations, that is, to find the best model parameters in the search space:
$${{m}}_{{i}}\left({{t}}\right)= \frac{{{fi}}{{t}}_{{i}}\left({{t}}\right)- \text{worst}\left({{t}}\right)}{{\text{best}}\left({{t}}\right)-\text{worst}\left({{t}}\right)}$$
(7)
$${{M}}_{{i}}\left({{t}}\right)= \frac{{{m}}_{{i}}{(t)}}{{\sum }_{{j=1}}^{{N}}{{{m}}}_{{j}}{(t)}},$$
(8)
where best(t) and worst(t) are the best and worst fitness, respectively:
$$\mathrm{best}\left(t\right)= \underset{\mathrm{j}\in 1,\cdots ,\mathrm{N}}{\mathrm{min}}{\mathrm{fit}}_{\mathrm{j}}\left(\mathrm{t}\right)$$
(9)
$$\mathrm{worst}\left(t\right)= \underset{\mathrm{j}\in \{1,\cdots ,\mathrm{N}\}}{\mathrm{max}}{\mathrm{fit}}_{\mathrm{j}}(\mathrm{t})$$
(10)
Acceleration phase
The individual will be subjected to the force of other individuals in the search space. Based on the law of universal gravitation, the force \({\text{F}}_{\text{ij}}^{\text{d}}\left({\text{t}}\right)\) of individual i by individual j is:
$${F}_{ij}^{d}\left(t\right)=G\left(t\right)\frac{{M}_{i}\left(t\right)\times {M}_{j}\left(t\right)}{{R}_{ij}\left(t\right)+ \varepsilon }\left({x}_{i}^{d}\left(t\right) - {x}_{j}^{d}\left(t\right)\right)$$
(11)
where Mi(t) and Mj(t) are the inertial masses of individuals i and j at time t, respectively, that is, the TACs calculated by a set of kinetic parameters. \(\epsilon\) is a small value to prevent errors. Rij(t) is the Euclidean distance between individuals i and j. G(t) is a gravitational constant that decreases with time t and is described by Eq. (12); it can affect the force and acceleration of the individual:
$$G\left(t\right)={G}_{0} {e}^{-\alpha \frac{t}{T}}$$
(12)
where \({\text{G}}_{0}\) is an initial value, α is a constant, t is the current number of iterations, and T is the maximum number of iterations. The traditional gravitational constant was not fully explored in the early stage of iteration, and fell into a local optimum [20]. Lei et al. [21] introduced a self-adaptive gravitational constant, and Seyedali et al. [22] used chaotic mapping to adjust the gravitational constant.
Therefore, this paper proposes a new gravitational constant expression, introducing an improved dynamic adjustment strategy, which contains a random variable as follows:
$${G}{{^{\prime}}\left(t\right)}={G}_{0} {e}^{-\alpha \frac{t}{{T}^{1.5}} \times \left({rand}_{t}+\frac{t}{T}\right)},$$
(13)
where \({\text{rand}}_{\text{t}}\) is a random number in (0,1). \(G{^{\prime}}\left(t\right)\) adopts a large step and long movement in the early stage of iteration to increase the particle exploration ability and enough time for optimization. Random variables can abruptly change the gravitational constant during iterations, improving the ability of particles to jump out of the local optimum. The process of moving in small steps is adopted in the later stage of iteration, which effectively avoids the premature convergence of particles and improves the mining capacity.
The resultant force of each particle is calculated as follows:
$${F}_{i}^{d}\left(t\right)=\sum_{j\in Kbest,j\ne i}{rand}_{j}{F}_{ij}^{d}\left(t\right)$$
(14)
where \({\text{Kbest}}\) decreases with the number of iterations, the initial value is N, and \({\text{rand}}_{\text{j}}\) is a random number between the interval (0, 1). The acceleration of particle i at time t is calculated as follows:
$${a}_{i}^{d}\left(t\right)= \frac{{F}_{i}^{d}(t)}{{M}_{i}(t)}$$
(15)
Update position phase
In each iteration, the particles update their velocity and position as follows:
$${v}_{i}^{d}\left(t+1\right)= {rand}_{i} \times {v}_{i}^{d}\left(t\right)+ {a}_{i}^{d}\left(t\right)$$
(16)
$${x}_{i}^{d}\left(t+1\right)= {x}_{i}^{d}\left(t\right)+{v}_{i}^{d}(t+1)$$
(17)
where \({\text{rand}}_{\text{i}}\) is a uniformly distributed random number between the interval (0, 1). However, some particles move too fast in the moving process, thus flying out of the solution space. Li et al. [23] introduced inertia weight instead of random variables to restrict the particle velocity. Gao et al. [24] replaced random variables with chaotic sequences. Therefore, this paper adds inertia weight and chaotic sequence to the particle speed updating process to further limit the particle speed:
$${v}_{i}^{d}\left(t+1\right)={W}^{d}(t)\times {v}_{i}^{d}\left(t\right)+ {a}_{i}^{d}\left(t\right)$$
(18)
$${x}_{i}^{d}\left(t+1\right)={x}_{i}^{d}\left(t\right)+{v}_{i}^{d}(t+1)\times c(i)$$
(19)
$${W}^{d}\left(t\right)={\omega }_{max}-\frac{{\omega }_{max}-{\omega }_{min}}{T}\times t,$$
(20)
where \({\omega}_{\text{max}}\) and \({\omega}_{\text{min}}\) are the inertia weights (\({\omega}_{\text{max}}\)=0.7, \({\omega}_{\text{min}}\)=0.1, which is mainly based on experience), and c(i) is a chaotic sequence. A larger inertia weight can improve the exploration ability, and a smaller inertia weight can improve the mining ability. Moreover, the ergodicity and dynamics of the chaotic sequence have the ability to jump out of the local optimum.
Statistical analysis
Statistical analysis was performed using MedCalc version 13.0.0.0 (MedCalc software, Ostend, Belgium). The derived parameters are expressed as the mean ± standard deviation. The Student’s t-test was used to compare the estimated parameters between HCCs and background liver tissues. The box plot was used to assess the consistency of the estimated k1 and k3 for the different methods. The diagnostic performance of k1 and k3 among the three methods was compared using receiver operating characteristic (ROC) curve analysis. P < 0.05 indicated significant differences. The fitting quality of the TACs among the three methods was compared using the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), and a smaller value represents better curve fitting.
