# Integrated Prediction Framework for Clinical Scores of Cognitive Functions in ESRD Patients

Authors: Yutao Zhang; Quan Sheng; Xidong Fu; Haifeng Shi; Zhuqing Jiao
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

The clinical scores are applied to determine the stage of cognitive function in patients with end-stage renal disease (ESRD). However, accurate clinical scores are hard to come by. This paper proposed an integrated prediction framework with GPLWLSV to predict clinical scores of cognitive functions in ESRD patients. GPLWLSV incorporated three parts, graph theoretic algorithm (GTA) and principal component analysis (PCA), whale optimization algorithm with Levy flight (LWOA), and least squares support vector regression machine (LSSVRM). GTA was adopted to extract features from the brain functional networks in ESRD patients, while PCA was used to select features. LSSVRM was built to explore the relationship between the selected features and the clinical scores of ESRD patients. Whale optimization algorithm (WOA) was introduced to select better parameters of the kernel function in LSSVRM; it aims to improve the exploration competence of LSSVRM. Levy flight was used to optimize the ability to jump out of local optima in WOA and improve the convergence of coefficient vectors in WOA, which lead to an increase in the generalization ability and convergence speed of WOA. The results validated that the prediction accuracy of GPLWLSV was higher than that of several comparable frameworks, such as GPSV, GPLSV, and GPWLSV. In particular, the average of root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) between the predicted scores and the actual scores of ESRD patients was 2.40, 2.06, and 9.83%, respectively. The proposed framework not only can predict the clinical scores more accurately but also can capture imaging markers associated with decline of cognitive function. It helps to understand the potential relationship between structural changes in the brain and cognitive function of ESRD patients.


---

## Body

## 1. Introduction

In recent years, the mortality rate of patients with end-stage renal disease (ESRD) has been increasing year by year. ESRD is the second highest increase of any disease, which is one of the top risks to human health. Although ESRD patients represent 0.1% of the global population, World Public Health Organization spends 2%-3% of its total expenditure on the treatment of ESRD. ESRD not only places a severe economic burden on society and families but also places a significant mental burden on patients [1]. Cognitive impairment is common in ESRD patients, especially in the aspects of orientation, attention, and executive ability [2]. Bugnicourt et al. [3] found that the incidence of cognitive impairment in patients undergoing hemodialysis was as high as 30% ∼ 60%. Cognitive impairment may affect late treatment of ESRD patients, such as dietary adjustment and medication compliance. Cognitive impairment in ESRD patients is often ignored in clinical practice, and its pathophysiological mechanism has not been fully elucidated [4]. In consequence, it is of great significance to study the cognitive impairment mode of patients with ESRD and clarify the exact pathophysiological mechanism.As neuroimaging technology develops by leaps and bounds, the neuropathological mechanism of ESRD can be learned from the perspective of central nervous systems. For example, Liang [5] et al. observed the brain neural activity of ESRD patients with the help of resting-state fMRI technology. Compared with normal controls, the functional activity of ESRD patients significantly decreased in bilateral frontal parietal temporal lobe, suggesting that the abnormal connection of brain functional networks in ESRD patients could result in cognitive disorders. In virtue of Diffusion Tensor Imaging (DTI) technology, Chou et al. [6] pointed out that long-term hemodialysis would aggravate cerebral interstitial edema in ESRD patients and cause demyelination of pontine axons, indicating that hemodialysis may lead to extensive white matter damage in the brain. Chai et al. [7] found that the volume of gray matter in the left lobe and bilateral putamen of ESRD patients was significantly lower than that of healthy subjects through voxel-based morphometry. This suggests that changes in the volume of gray matter in the left lobe bring about cognitive impairment probably. These researches relying on neuroimaging technology have greatly deepened our understanding of the cognitive function of ESRD patients, but it is still impossible to determine the stage of cognitive function of ESRD patients.Clinically, neurologists often judge the stage of cognitive function of ESRD patients according to the scores of the Montreal Cognitive Assessment Scale (MoCA), also known as clinical scores. Accurate prediction of clinical scores is beneficial to estimating the stage of cognitive functions of ESRD patients. Jiang et al. [8] calculated the correlation coefficient between topological attribute parameters of brain functional networks and clinical scores of cognitive functions in ESRD patients by Pearson correlation analysis. They mainly focused on the imaging markers affecting the cognitive function of ESRD patients and could not judge the current status of cognitive function well. Lu et al. [9] proposed a prediction method of clinical scores based on the brain functional networks. By virtue of simple definition and small computation, clustering coefficients in the brain functional networks were extracted as features to predict clinical scores, which were referred to judge the current stage of cognitive function. However, this method is subjective in the process of extracting features; at the same time, the influence of other features on the cognitive function of patients is ignored. Yang et al. [10] predicted clinical scores by support vector regression machine (SVRM). Regrettably, SVRM has high volatility and low prediction accuracy in the process of predicting clinical scores.On the above considerations, we propose an integrated framework to predict clinical scores of ESRD patients. The main work is as follows: firstly, fMRI data were preprocessed to construct the brain functional networks, and graph theoretic algorithm (GTA) was adopted to extract features. Secondly, principal component analysis (PCA) was used to select features. Then, the least squares support vector regression machine (LSSVRM) was built to explore the relationship between the selected features and the clinical scores of ESRD patients. Meanwhile, the whale optimization algorithm (WOA) was introduced to select better parameters of kernel function in LSSVRM so as to enhance the exploration competence of LSSVRM. Finally, Levy flight replaced the traditional selection of WOA and in the meantime optimized the ability to jump out of local optimum in WOA. The framework called GPLWLSV was constructed to predict the clinical scores of ESRD patients and then determine their current stage of cognitive function.


## 2. Data and Methods

### 2.1. Research Framework

Figure1 shows our prediction framework, which mainly includes the following steps: (1) The original resting-state fMRI images were preprocessed by Data Processing Assistant for Resting-State fMRI (DPARSF). (2) Time series from preprocessed images were abstracted to construct the brain functional networks. (3) The area under the curve (AUC) of topological attribute parameters in the brain functional networks was extracted as features through GTA. (4) PCA was used to filter redundant features and to retain important features. (5) LSSVRM was built to retain features. (6) WOA was improved by Levy flight to optimize the selection strategy of kernel function parameters in LSSVRM. (7) GPLWLSV was constructed to predict clinical scores of cognitive functions in ESRD patients.Figure 1
A flowchart of prediction framework.


### 2.2. Experimental Data and Pretreatment

A total of 50 patients with ESRD were admitted to Changzhou Second People’s Hospital Affiliated to Nanjing Medical University from May 2021 to March 2022, including 27 male and 23 female individuals aged 49.12 ± 8.23 years. Synchronously, a total of 40 normal controls were also admitted to the same hospital, including 22 male and 18 female individuals aged 47.26 ± 7.01 years. There were no significant differences (P>0.05) in the gender ratio, age, and education level between them. One hour before fMRI examination, the stage of cognitive function of all subjects was assessed by trained neurologists who did not know the data of subjects via clinical scores of cognitive functions. Table 1 gives the demographic information of these two groups of subjects.Table 1
Demographic information of subjects.
Gender (male/female)Age (years,x¯±s)Education years (years,x¯±s)Clinical scores (points,x¯±s)ESRD patients (n = 50)27/2349.12 ± 8.2311.13 ± 2.0521.47 ± 2.75Normal controls (n = 40)22/1847.26 ± 7.0111.24 ± 2.1327.38 ± 1.35t/χ20.0091.1330.387−13.412P>0.050.2600.7780.000All subjects were scanned with GE Discovery MR 750W 3.0T superconducting MR scanner with 32-channel head and neck joint coil. The head of subject was fixed with a cushion to reduce the artifacts produced by head movement. Gradient echo plane echo imaging (GRE-EPI) sequence was used to collect fMRI images. During the collection process, subjects were required to keep quiet and awake and try not to think. Machine scanning parameters are as follows: repetition time (TR) = 2000 ms, echo time (TE) = 40 ms, field of view (FOV) = 24 cm, flip angle (FA) = 90°, matrix size = 64 × 64, and layer thickness = 6 mm.After achieving the fMRI images of all subjects, DPARSF (available athttp://rfmri.org/dpabi) participated in preprocessing. The specific steps are as follows: (a) Image format was transformed: the correct format can only be opened by DPARSF. (b) The first 10 time points were deleted: errors may occur when the instrument and subjects were not in a stable state. (c) Time points were corrected: data of time points at different levels were corrected to the same time point to maintain consistency. (d) Brain positions were corrected: the brain must remain in the same position throughout the examination. (e) Space was normalized: put the brain of subject into a standard space for the convenience of follow-up statistics and reports. (f) Spatial structure differences were reduced through smooth: better image results were beneficial to the validity of statistical tests. (g) Linear drift was removed: changes in time and temperature caused errors in the inspection instrument. (h) Filtering was performed: the frequency range is 0.01∼0.08 Hz. The resting fMRI signal after low-frequency filtering has important physiological significance and may reflect spontaneous neural activity. (i) Final time series was obtained: the mean BOLD time series of head motion parameters, white matter, and cerebrospinal fluid were regression and removed.


### 2.3. Principle of GP

#### 2.3.1. Features Were Extracted by GTA

The construction of brain functional networks is as follows. Firstly, the brain of each subject was divided into 90 brain regions on the basis of Automated Anatomical Labeling (AAL). Secondly, the Pearson correlation coefficient was calculated between the time series of two brain regions;K-matrix was born, which is a 90 × 90 symmetric matrix with all 1 s on the diagonal. Finally, the Fisher-Z transformation is performed on the elements in the K-matrix divided by the diagonal, which is converted to the Z-value close to the normal distribution to generate the Z-matrix:(1)K=1p1,2⋯p1,90p2,1⋱⋮⋮⋱⋮p90,1⋯⋯1,(2)Z=1z1,2⋯z1,90z2,1⋱⋮⋮⋱⋮z90,1⋯⋯1,(3)zij=12ln1+pij−ln1−piji≠j,where pij represents the Pearson correlation coefficient between the time series of the i-th brain region and the j-th brain region.Taking matrix sparsity as the threshold, the preprocessing forZ-matrix was binary. The matrix sparsity was set to 0.1–0.4 with a span of 0.01. Within the threshold range of the matrix sparsity, the AUC of topological attribute parameters was calculated by GTA as features, which reflected the overall change of brain functional networks and reduced the influence of individual differences. These features include the AUC of global efficiency (Eglobal), local efficiency (Elocal), clustering coefficient (Cp), characteristic path length (Lp), standardized clustering coefficient (γ), standardized characteristic path length (λ), and small-world properties (σ).


#### 2.3.2. Features Were Selected by PCA

In the process of machine learning, a large number of features need to be extracted and converted into various data that can be processed by computers [11]. Increasing the number of features can improve the effect of machine learning. The trouble is that as the number of features increases, the dimension of feature vectors will also add. This not only makes machine learning more difficult but can also lead to overfitting, affecting the final accuracy. In view of this situation, it is necessary to filter out some unimportant features or combine some related features [12]. Features were selected by PCA that keep the information contained in the original data as much as possible while reducing the dimension. The specific steps are as follows.The first step is to standardize the matrix. The AUC of Eglobal, Elocal, Cp, Lp,γ, λ, and σ of ESRD patients was marked as α1∼α7. Constructing a matrix called A, A = [α1, α2, α3, α4, α5, α6, α7], and αj = [x1j, x2j, x3j…… xnj]T.B-matrix is the normalization of A-matrix,B = [β1, β2, β3, β4, β5, β6, β7], and βj = [ y1j, y2j, y3j…… ynj] T [13]:(4)yij=xij−x¯xij−x¯2/n−1,where x¯ is the mean value of samples in αj and n is the number of samples.The second step is to construct the covariance matrix ofB-matrix. Computing the covariance between βi and βj in β1 ∼ β7,(5)rij=covX,Y=∑i=1nXi−X¯Yi−Y¯n−1,where Xi is the value of the i-th sample; X¯ represents the mean value of samples; n is the number of samples.Constructing the covariance matrix calledR,(6)R=r11r12⋯r1mr21⋱r2m⋮⋱⋮rm1rm2⋯rmm,where m represents the number of features and its value is 7.The third step is to extract the original features. Computing the eigenvalues and eigenvectors of theR-matrix, the principal component is the eigenvector of the R-matrix, and the magnitude of the corresponding eigenvalue indicates the importance of the eigenvector. Calculating the weight of the features, λ1∼λ7 are the eigenvalues of the R-matrix, λj represents the dominance of αj in the A-matrix, and the weight of αj in the whole feature set called α1∼α7 can be expressed as [14](7)wi=λi∑j=17λj.The weight of extracted features is higher than 0.6 set in this paper.Finally, the extracted features are transformed to new features. Finding the eigenvectorvj in the R-matrix corresponding to αj, the new feature cj is transformed from vj [15]:(8)cj=B⋅vj.


### 2.4. Principle of LSSVRM

The principle of SVRM is to map data vectors from low-dimensional space to high-dimensional space. SVRM constructs decision functions with the aid of the principle of risk minimization. SVRM is suited to solve nonlinear problems of small samples and high dimension of feature vectors. LSSVRM is simplified on the basis of SVRM. There are two main points: (1) changing constraints in SVRM to improve computing efficiency; (2) selecting different decision functions to reduce the operation time. LSSVRM maps nonlinear vector Φ(x) to high-dimensional space and transforms it into a linear regression problem, as shown in formula (9) [16]:(9)y=ωTΦx+b,where ω is the n-dimensional weight vector; b is the deviation.Following the principle of risk minimization, linear regression can be transformed into an optimal problem [17]:(10)minJω,ξ=12ωTω+12J∑i=1mξi2,(11)yi=ωTΦxi+b+ξi,where J represents the normalization function; ξi represents the error; i represents the i-the dimension of the space vector; ω represents the weight vector. Since ω belongs to the high-dimensional space, it cannot be solved directly, so the kernel function is introduced [18]:(12)y=∑i=1naiKx,xi+b.In order to ensure the efficiency of operation, RBF is chosen as the kernel function, and its function is expressed as follows [19]:(13)Kxi,xj=exp−|xi−xj|2σ2.It can be seen from the above process thatJ and σ2 have the greatest influence on the LSSVRM. In searching for optimal J and σ2, WOA is replaced by LWOA to seek the optimal solution.


### 2.5. Principle of LWOA

WOA is a new optimization algorithm that simulates the predatory behavior of whales. Whales locate and surround prey through their scent. We can define a certain number of virtual whales as search agents, assuming that the location of the scent of prey is at or near the current optimal location. We look for the optimal solution as the next location of the whale by comparing the feasible solutions of various search agents. Other search agents update their location to complete the strategy of finding the optimal solution.Whales constantly adjust their position according to the position of prey during hunting. To describe this strategy of hunting, the following mathematical model is presented [20]:(14)M1=CX∗t−Xt,Xt+1=X∗t−HM1,where t is the number of current iterations; X(t) is the coordinate vector of the current whale; X(t+1) is the target coordinate vector after the next iteration; X∗t is the coordinate vector of the best solution at present. If there is a better feasible solution, X∗t should be updated immediately. C and H are coefficients, which are acquired by formulas (22) and (23) [21]:(15)H=2hr1−h,C=2r2,where r1and r2 are random numbers between 0 and 1; h is computed from formula (24) [21]:(16)h=2−2tTmax.For purpose of finding a better target position to approach the prey, the whale will randomly use any whale coordinate vector to replace the whale coordinate vector of the next iteration so as to achieve the purpose of deviating from the prey. This avoids falling into local optimality. The following mathematical model is shown [22]:(17)Xt+1=X∗t−HM1,H<1,p<0.5,Xrandt−HM2,H⩾1,p<0.5,X∗t+M3eblcos2πl,p⩾0.5,M1=CX∗t−Xt,M2=CXrand−Xt,M3=X∗t−Xt,where t is the current iteration times; X∗t is the optimal position vector so far; Xrand is the random position vector of whale; X(t) is the current position vector of whale; b is the constant, and default is 1; the role of b is to control the hunting path shape; l is acquired by the following formulas:(18)l=h2−1r3+1,h2=−1−tTmax,where r3 is a random number between 0 and 1; t represents the current iteration number; Tmax represents the maximum iteration number. When t ≥ 0.5Tmax, H is always less than 1. At this point, the whale enters the attack mode and no longer deviates from the prey through search agents.Due to the few adjustment parameters of WOA, it has a fast rate of convergence and a certain ability to jump out of the local optimum. It is worth noting that WOA can be optimized further. The reasons are as follows: Firstly, WOA searches by means of random system. Excessive reliance on random system limits the convergence speed of WOA. Secondly, WOA is subject to coefficient vector. WOA will lose the ability to jump out of the local optimum when the number of iterations reaches half of the maximum number of iterations set earlier. Consequently, WOA is accompanied by a certain risk of falling into the local optimum, leading to inaccurate results of prediction [23].The defects of WOA can be solved by an improved WOA with Levy flight named LWOA. Levy flight is a kind of random search that relies on Levy distribution, which has been applied many times in the optimization field in recent years. Levy flight is able to improve cuckoo and particle swarm optimization algorithms [23–25] and so on. LWOA owns a faster convergence speed and higher convergence accuracy; LWOA has a better ability to jump out of the local optimum. The specific steps of improvement are as follows.WOA is improved by Levy flight, and formula (22) is replaced by the following formula [16, 26]:(19)H=2hLevyλ−h,where Levy(λ) means that it obeys the Levy distribution with parameter λ [16, 27]:(20)Levy∼μ=t−λ.Due to the complexity of Levy flight, the Mantegna algorithm is adopted to simulate it, and its mathematical expression is as follows [16, 28]:(21)s=μν1/β,where μ and v obey the normal distribution with parameters σμ and σv [16, 29]:(22)μ∼N0,σμ2,ν∼N0,σν2,(23)σμ=Γ1+βsinπβ/2Γ1+β/2β2β−1/2,(24)σν=1.For higher operation efficiency of the algorithm,β is a constant 1.5, and σμ is a constant 0.7.The coefficient vectorH in WOA converges linearly with certain limitations. WOA should be promoted to jump out of local optimum; formula (24) changed by the following formula [16]:(25)h=2e0.15−log10t/Tmax4,where t represents the current iteration number and Tmax represents the maximum iteration number.In the early stage of iteration, the value ofh will decrease slowly as the increase of iteration times, which is conducive to the global search out of local optimum. At the end of iteration, the value of h will decrease exponentially to improve the ability of rapid local search.


## 2.1. Research Framework

Figure1 shows our prediction framework, which mainly includes the following steps: (1) The original resting-state fMRI images were preprocessed by Data Processing Assistant for Resting-State fMRI (DPARSF). (2) Time series from preprocessed images were abstracted to construct the brain functional networks. (3) The area under the curve (AUC) of topological attribute parameters in the brain functional networks was extracted as features through GTA. (4) PCA was used to filter redundant features and to retain important features. (5) LSSVRM was built to retain features. (6) WOA was improved by Levy flight to optimize the selection strategy of kernel function parameters in LSSVRM. (7) GPLWLSV was constructed to predict clinical scores of cognitive functions in ESRD patients.Figure 1
A flowchart of prediction framework.


## 2.2. Experimental Data and Pretreatment

A total of 50 patients with ESRD were admitted to Changzhou Second People’s Hospital Affiliated to Nanjing Medical University from May 2021 to March 2022, including 27 male and 23 female individuals aged 49.12 ± 8.23 years. Synchronously, a total of 40 normal controls were also admitted to the same hospital, including 22 male and 18 female individuals aged 47.26 ± 7.01 years. There were no significant differences (P>0.05) in the gender ratio, age, and education level between them. One hour before fMRI examination, the stage of cognitive function of all subjects was assessed by trained neurologists who did not know the data of subjects via clinical scores of cognitive functions. Table 1 gives the demographic information of these two groups of subjects.Table 1
Demographic information of subjects.
Gender (male/female)Age (years,x¯±s)Education years (years,x¯±s)Clinical scores (points,x¯±s)ESRD patients (n = 50)27/2349.12 ± 8.2311.13 ± 2.0521.47 ± 2.75Normal controls (n = 40)22/1847.26 ± 7.0111.24 ± 2.1327.38 ± 1.35t/χ20.0091.1330.387−13.412P>0.050.2600.7780.000All subjects were scanned with GE Discovery MR 750W 3.0T superconducting MR scanner with 32-channel head and neck joint coil. The head of subject was fixed with a cushion to reduce the artifacts produced by head movement. Gradient echo plane echo imaging (GRE-EPI) sequence was used to collect fMRI images. During the collection process, subjects were required to keep quiet and awake and try not to think. Machine scanning parameters are as follows: repetition time (TR) = 2000 ms, echo time (TE) = 40 ms, field of view (FOV) = 24 cm, flip angle (FA) = 90°, matrix size = 64 × 64, and layer thickness = 6 mm.After achieving the fMRI images of all subjects, DPARSF (available athttp://rfmri.org/dpabi) participated in preprocessing. The specific steps are as follows: (a) Image format was transformed: the correct format can only be opened by DPARSF. (b) The first 10 time points were deleted: errors may occur when the instrument and subjects were not in a stable state. (c) Time points were corrected: data of time points at different levels were corrected to the same time point to maintain consistency. (d) Brain positions were corrected: the brain must remain in the same position throughout the examination. (e) Space was normalized: put the brain of subject into a standard space for the convenience of follow-up statistics and reports. (f) Spatial structure differences were reduced through smooth: better image results were beneficial to the validity of statistical tests. (g) Linear drift was removed: changes in time and temperature caused errors in the inspection instrument. (h) Filtering was performed: the frequency range is 0.01∼0.08 Hz. The resting fMRI signal after low-frequency filtering has important physiological significance and may reflect spontaneous neural activity. (i) Final time series was obtained: the mean BOLD time series of head motion parameters, white matter, and cerebrospinal fluid were regression and removed.


## 2.3. Principle of GP

### 2.3.1. Features Were Extracted by GTA

The construction of brain functional networks is as follows. Firstly, the brain of each subject was divided into 90 brain regions on the basis of Automated Anatomical Labeling (AAL). Secondly, the Pearson correlation coefficient was calculated between the time series of two brain regions;K-matrix was born, which is a 90 × 90 symmetric matrix with all 1 s on the diagonal. Finally, the Fisher-Z transformation is performed on the elements in the K-matrix divided by the diagonal, which is converted to the Z-value close to the normal distribution to generate the Z-matrix:(1)K=1p1,2⋯p1,90p2,1⋱⋮⋮⋱⋮p90,1⋯⋯1,(2)Z=1z1,2⋯z1,90z2,1⋱⋮⋮⋱⋮z90,1⋯⋯1,(3)zij=12ln1+pij−ln1−piji≠j,where pij represents the Pearson correlation coefficient between the time series of the i-th brain region and the j-th brain region.Taking matrix sparsity as the threshold, the preprocessing forZ-matrix was binary. The matrix sparsity was set to 0.1–0.4 with a span of 0.01. Within the threshold range of the matrix sparsity, the AUC of topological attribute parameters was calculated by GTA as features, which reflected the overall change of brain functional networks and reduced the influence of individual differences. These features include the AUC of global efficiency (Eglobal), local efficiency (Elocal), clustering coefficient (Cp), characteristic path length (Lp), standardized clustering coefficient (γ), standardized characteristic path length (λ), and small-world properties (σ).


### 2.3.2. Features Were Selected by PCA

In the process of machine learning, a large number of features need to be extracted and converted into various data that can be processed by computers [11]. Increasing the number of features can improve the effect of machine learning. The trouble is that as the number of features increases, the dimension of feature vectors will also add. This not only makes machine learning more difficult but can also lead to overfitting, affecting the final accuracy. In view of this situation, it is necessary to filter out some unimportant features or combine some related features [12]. Features were selected by PCA that keep the information contained in the original data as much as possible while reducing the dimension. The specific steps are as follows.The first step is to standardize the matrix. The AUC of Eglobal, Elocal, Cp, Lp,γ, λ, and σ of ESRD patients was marked as α1∼α7. Constructing a matrix called A, A = [α1, α2, α3, α4, α5, α6, α7], and αj = [x1j, x2j, x3j…… xnj]T.B-matrix is the normalization of A-matrix,B = [β1, β2, β3, β4, β5, β6, β7], and βj = [ y1j, y2j, y3j…… ynj] T [13]:(4)yij=xij−x¯xij−x¯2/n−1,where x¯ is the mean value of samples in αj and n is the number of samples.The second step is to construct the covariance matrix ofB-matrix. Computing the covariance between βi and βj in β1 ∼ β7,(5)rij=covX,Y=∑i=1nXi−X¯Yi−Y¯n−1,where Xi is the value of the i-th sample; X¯ represents the mean value of samples; n is the number of samples.Constructing the covariance matrix calledR,(6)R=r11r12⋯r1mr21⋱r2m⋮⋱⋮rm1rm2⋯rmm,where m represents the number of features and its value is 7.The third step is to extract the original features. Computing the eigenvalues and eigenvectors of theR-matrix, the principal component is the eigenvector of the R-matrix, and the magnitude of the corresponding eigenvalue indicates the importance of the eigenvector. Calculating the weight of the features, λ1∼λ7 are the eigenvalues of the R-matrix, λj represents the dominance of αj in the A-matrix, and the weight of αj in the whole feature set called α1∼α7 can be expressed as [14](7)wi=λi∑j=17λj.The weight of extracted features is higher than 0.6 set in this paper.Finally, the extracted features are transformed to new features. Finding the eigenvectorvj in the R-matrix corresponding to αj, the new feature cj is transformed from vj [15]:(8)cj=B⋅vj.


## 2.3.1. Features Were Extracted by GTA

The construction of brain functional networks is as follows. Firstly, the brain of each subject was divided into 90 brain regions on the basis of Automated Anatomical Labeling (AAL). Secondly, the Pearson correlation coefficient was calculated between the time series of two brain regions;K-matrix was born, which is a 90 × 90 symmetric matrix with all 1 s on the diagonal. Finally, the Fisher-Z transformation is performed on the elements in the K-matrix divided by the diagonal, which is converted to the Z-value close to the normal distribution to generate the Z-matrix:(1)K=1p1,2⋯p1,90p2,1⋱⋮⋮⋱⋮p90,1⋯⋯1,(2)Z=1z1,2⋯z1,90z2,1⋱⋮⋮⋱⋮z90,1⋯⋯1,(3)zij=12ln1+pij−ln1−piji≠j,where pij represents the Pearson correlation coefficient between the time series of the i-th brain region and the j-th brain region.Taking matrix sparsity as the threshold, the preprocessing forZ-matrix was binary. The matrix sparsity was set to 0.1–0.4 with a span of 0.01. Within the threshold range of the matrix sparsity, the AUC of topological attribute parameters was calculated by GTA as features, which reflected the overall change of brain functional networks and reduced the influence of individual differences. These features include the AUC of global efficiency (Eglobal), local efficiency (Elocal), clustering coefficient (Cp), characteristic path length (Lp), standardized clustering coefficient (γ), standardized characteristic path length (λ), and small-world properties (σ).


## 2.3.2. Features Were Selected by PCA

In the process of machine learning, a large number of features need to be extracted and converted into various data that can be processed by computers [11]. Increasing the number of features can improve the effect of machine learning. The trouble is that as the number of features increases, the dimension of feature vectors will also add. This not only makes machine learning more difficult but can also lead to overfitting, affecting the final accuracy. In view of this situation, it is necessary to filter out some unimportant features or combine some related features [12]. Features were selected by PCA that keep the information contained in the original data as much as possible while reducing the dimension. The specific steps are as follows.The first step is to standardize the matrix. The AUC of Eglobal, Elocal, Cp, Lp,γ, λ, and σ of ESRD patients was marked as α1∼α7. Constructing a matrix called A, A = [α1, α2, α3, α4, α5, α6, α7], and αj = [x1j, x2j, x3j…… xnj]T.B-matrix is the normalization of A-matrix,B = [β1, β2, β3, β4, β5, β6, β7], and βj = [ y1j, y2j, y3j…… ynj] T [13]:(4)yij=xij−x¯xij−x¯2/n−1,where x¯ is the mean value of samples in αj and n is the number of samples.The second step is to construct the covariance matrix ofB-matrix. Computing the covariance between βi and βj in β1 ∼ β7,(5)rij=covX,Y=∑i=1nXi−X¯Yi−Y¯n−1,where Xi is the value of the i-th sample; X¯ represents the mean value of samples; n is the number of samples.Constructing the covariance matrix calledR,(6)R=r11r12⋯r1mr21⋱r2m⋮⋱⋮rm1rm2⋯rmm,where m represents the number of features and its value is 7.The third step is to extract the original features. Computing the eigenvalues and eigenvectors of theR-matrix, the principal component is the eigenvector of the R-matrix, and the magnitude of the corresponding eigenvalue indicates the importance of the eigenvector. Calculating the weight of the features, λ1∼λ7 are the eigenvalues of the R-matrix, λj represents the dominance of αj in the A-matrix, and the weight of αj in the whole feature set called α1∼α7 can be expressed as [14](7)wi=λi∑j=17λj.The weight of extracted features is higher than 0.6 set in this paper.Finally, the extracted features are transformed to new features. Finding the eigenvectorvj in the R-matrix corresponding to αj, the new feature cj is transformed from vj [15]:(8)cj=B⋅vj.


## 2.4. Principle of LSSVRM

The principle of SVRM is to map data vectors from low-dimensional space to high-dimensional space. SVRM constructs decision functions with the aid of the principle of risk minimization. SVRM is suited to solve nonlinear problems of small samples and high dimension of feature vectors. LSSVRM is simplified on the basis of SVRM. There are two main points: (1) changing constraints in SVRM to improve computing efficiency; (2) selecting different decision functions to reduce the operation time. LSSVRM maps nonlinear vector Φ(x) to high-dimensional space and transforms it into a linear regression problem, as shown in formula (9) [16]:(9)y=ωTΦx+b,where ω is the n-dimensional weight vector; b is the deviation.Following the principle of risk minimization, linear regression can be transformed into an optimal problem [17]:(10)minJω,ξ=12ωTω+12J∑i=1mξi2,(11)yi=ωTΦxi+b+ξi,where J represents the normalization function; ξi represents the error; i represents the i-the dimension of the space vector; ω represents the weight vector. Since ω belongs to the high-dimensional space, it cannot be solved directly, so the kernel function is introduced [18]:(12)y=∑i=1naiKx,xi+b.In order to ensure the efficiency of operation, RBF is chosen as the kernel function, and its function is expressed as follows [19]:(13)Kxi,xj=exp−|xi−xj|2σ2.It can be seen from the above process thatJ and σ2 have the greatest influence on the LSSVRM. In searching for optimal J and σ2, WOA is replaced by LWOA to seek the optimal solution.


## 2.5. Principle of LWOA

WOA is a new optimization algorithm that simulates the predatory behavior of whales. Whales locate and surround prey through their scent. We can define a certain number of virtual whales as search agents, assuming that the location of the scent of prey is at or near the current optimal location. We look for the optimal solution as the next location of the whale by comparing the feasible solutions of various search agents. Other search agents update their location to complete the strategy of finding the optimal solution.Whales constantly adjust their position according to the position of prey during hunting. To describe this strategy of hunting, the following mathematical model is presented [20]:(14)M1=CX∗t−Xt,Xt+1=X∗t−HM1,where t is the number of current iterations; X(t) is the coordinate vector of the current whale; X(t+1) is the target coordinate vector after the next iteration; X∗t is the coordinate vector of the best solution at present. If there is a better feasible solution, X∗t should be updated immediately. C and H are coefficients, which are acquired by formulas (22) and (23) [21]:(15)H=2hr1−h,C=2r2,where r1and r2 are random numbers between 0 and 1; h is computed from formula (24) [21]:(16)h=2−2tTmax.For purpose of finding a better target position to approach the prey, the whale will randomly use any whale coordinate vector to replace the whale coordinate vector of the next iteration so as to achieve the purpose of deviating from the prey. This avoids falling into local optimality. The following mathematical model is shown [22]:(17)Xt+1=X∗t−HM1,H<1,p<0.5,Xrandt−HM2,H⩾1,p<0.5,X∗t+M3eblcos2πl,p⩾0.5,M1=CX∗t−Xt,M2=CXrand−Xt,M3=X∗t−Xt,where t is the current iteration times; X∗t is the optimal position vector so far; Xrand is the random position vector of whale; X(t) is the current position vector of whale; b is the constant, and default is 1; the role of b is to control the hunting path shape; l is acquired by the following formulas:(18)l=h2−1r3+1,h2=−1−tTmax,where r3 is a random number between 0 and 1; t represents the current iteration number; Tmax represents the maximum iteration number. When t ≥ 0.5Tmax, H is always less than 1. At this point, the whale enters the attack mode and no longer deviates from the prey through search agents.Due to the few adjustment parameters of WOA, it has a fast rate of convergence and a certain ability to jump out of the local optimum. It is worth noting that WOA can be optimized further. The reasons are as follows: Firstly, WOA searches by means of random system. Excessive reliance on random system limits the convergence speed of WOA. Secondly, WOA is subject to coefficient vector. WOA will lose the ability to jump out of the local optimum when the number of iterations reaches half of the maximum number of iterations set earlier. Consequently, WOA is accompanied by a certain risk of falling into the local optimum, leading to inaccurate results of prediction [23].The defects of WOA can be solved by an improved WOA with Levy flight named LWOA. Levy flight is a kind of random search that relies on Levy distribution, which has been applied many times in the optimization field in recent years. Levy flight is able to improve cuckoo and particle swarm optimization algorithms [23–25] and so on. LWOA owns a faster convergence speed and higher convergence accuracy; LWOA has a better ability to jump out of the local optimum. The specific steps of improvement are as follows.WOA is improved by Levy flight, and formula (22) is replaced by the following formula [16, 26]:(19)H=2hLevyλ−h,where Levy(λ) means that it obeys the Levy distribution with parameter λ [16, 27]:(20)Levy∼μ=t−λ.Due to the complexity of Levy flight, the Mantegna algorithm is adopted to simulate it, and its mathematical expression is as follows [16, 28]:(21)s=μν1/β,where μ and v obey the normal distribution with parameters σμ and σv [16, 29]:(22)μ∼N0,σμ2,ν∼N0,σν2,(23)σμ=Γ1+βsinπβ/2Γ1+β/2β2β−1/2,(24)σν=1.For higher operation efficiency of the algorithm,β is a constant 1.5, and σμ is a constant 0.7.The coefficient vectorH in WOA converges linearly with certain limitations. WOA should be promoted to jump out of local optimum; formula (24) changed by the following formula [16]:(25)h=2e0.15−log10t/Tmax4,where t represents the current iteration number and Tmax represents the maximum iteration number.In the early stage of iteration, the value ofh will decrease slowly as the increase of iteration times, which is conducive to the global search out of local optimum. At the end of iteration, the value of h will decrease exponentially to improve the ability of rapid local search.


## 3. Results

### 3.1. Experimental Settings

Table2 shows the AUC of topology attribute parameters of the brain functional networks of ESRD patients and normal controls calculated by GTA. Within the whole matrix sparsity threshold range, the AUC of γ, σ, and Elocal in ESRD patients was significantly lower than that in normal controls, with statistical significance (P<0.05). Nevertheless, there were no significant differences (P>0.05) in the AUC of λ, Cp, Lp, and Eglobal.Table 2
Comparison of AUC of global topological parameters between two groups (mean ± SD).
ParameterESRD patients (n = 50)Normal controls (n = 40)tPγ0.646 ± 0.0710.669 ± 0.056−1.7140.004λ0.326 ± 0.1230.329 ± 0.0240.7100.480σ0.589 ± 0.0650.607 ± 0.065−1.3840.008Cp0.175 ± 0.0140.176 ± 0.016−0.2960.768Lp0.543 ± 0.0270.553 ± 0.069−0.9660.337Eglobal0.171 ± 0.0060.170 ± 0.011−0.4770.635Elocal0.230 ± 0.0070.241 ± 0.006−0.2530.001Table3 shows the corresponding weight of each feature in the feature set. The proportion of AUC of Elocal was the highest, up to 65.31%. The proportion of AUC of Eglobal was next to Elocal, accounting for 30.89%. The proportion of AUC of γ, λ, σ, Cp, and Lp was less than 5%. The proportion of γ was the lowest, only 0.03%.Table 3
Weight of AUC of global topological parameters of ESRD patients in feature set.
Parameterγ (%)λ (%)σ (%)Cp (%)Lp (%)Eglobal (%)Elocal (%)wi × 100%0.030.070.170.283.2530.8965.31The experimental results can be seen from Tables2 and 3. The AUC of Elocal in ESRD patients was significantly lower than that in the normal controls, with statistically significant differences (P<0.05). Meanwhile, the AUC of Elocal accounted for the highest proportion in the feature set. Obviously, the AUC of Elocal was selected as the feature to construct GPSV [30], GPLSV [31], GPWLSV [32], and GPLWLSV for predicting the clinical scores of cognitive functions in ESRD patients [33].The prediction accuracy of various frameworks was evaluated by a tenfold cross-validation method. Firstly, the AUC of Elocal of 50 ESRD patients and the corresponding clinical scores were collected as data set calledD. Secondly, D was divided into 10 mutually exclusive subsets of the same size, D = D1∪D2∪ D3∪D4∪D5∪D6∪D7∪D8∪D9∪D10, and Di∪Dj = ∅. Each subset Di was separated from D through stratified sampling, which ensures consistency of data distribution. Each subset Di contains five samples. Taking the union of 9 subsets as the training set and the remaining subset as the test set, 10 groups of training sets and test sets are formed. Different frameworks were trained and tested for 10 times, and the average of 10 groups of test results was calculated.


### 3.2. Experimental Results

The root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) were selected as the testing standards of prediction accuracies. Table4 shows the prediction accuracies of some representative frameworks. All the comparable predictive frameworks were evaluated by tenfold cross-validation method. The test results show that the prediction accuracy of GPLWLSV is higher than that of GPSV, GPLSV, and GPWLSV. The average of RMSE of GPLWLSV was 2.40, 1.04, 0.93, and 0.46 points lower than that of GPSV, GPLSV, and GPWLSV, respectively. The average of MAE of GPLWLSV was 2.06, which was 0.82, 0.74, and 0.42 points lower than that of GPSV, GPLSV, and GPWLSV, respectively. MAPE can reflect the relative errors of frameworks better than MAE. The average of MAPE of GPLWLSV was 9.83%, lower than 10%, which was 4.10%, 3.75%, and 1.93% lower than that of GPSV, GPLSV, and GPWLSV, respectively. The bar chart in Figure 2 shows that the prediction accuracy of GPLWLSV is better than that of GPSV, GPLSV, and GPWLSV intuitively.Table 4
Prediction accuracies of comparable frameworks.
Prediction frameworkThe test groupRMSEMAEMAPEGPSV13.50642.38350.1364GPSV23.4823.09510.1454GPSV33.26482.84990.1267GPSV43.61543.17920.1542GPSV52.55082.29650.1109GPSV62.27241.79190.079GPSV74.1153.68250.1788GPSV82.82042.6590.1277GPSV95.72474.05230.1999GPSV103.07072.78520.1334Average3.44232.87750.1392GPLSV15.26344.1970.2149GPLSV24.66283.18990.1305GPLSV33.48192.81210.1491GPLSV43.16112.78370.1385GPLSV52.10461.88840.0898GPLSV62.1622.12460.0994GPLSV73.35453.01990.1524GPLSV83.98913.7250.18GPLSV92.28051.77330.0812GPLSV102.88792.48380.1214Average3.33482.79980.1357GPWLSV14.24083.47820.1695GPWLSV23.32932.81780.124GPWLSV33.06572.41780.1188GPWLSV43.30272.98390.1466GPWLSV51.31251.13880.0529GPWLSV61.92721.41140.0618GPWLSV72.94122.71960.1311GPWLSV82.78992.68650.1284GPWLSV92.67582.32820.1092GPWLSV103.03132.80140.1335Average2.86162.47840.1176GPLWLSV13.91323.10260.1561GPLWLSV23.25062.83340.1258GPLWLSV32.73072.43340.1075GPLWLSV42.99762.37380.1216GPLWLSV51.21921.11110.053GPLWLSV60.7180.64970.0313GPLWLSV72.51672.11560.1036GPLWLSV82.07091.96360.0914GPLWLSV92.17571.94670.0916GPLWLSV102.41632.03020.1007Average2.40092.0560.0983Figure 2
Average of prediction accuracies of various frameworks. (a) The RMSE and MAE of various frameworks. (b) The MAPE of various frameworks.
(a)(b)Figure3 shows the comparison between the predicted scores of various frameworks and the actual scores. The solid blue line represents the actual scores, and the solid purple line represents the predicted scores. As can be seen from the figure, different from GPSV and GPLSV, GPWLSV and GPLWLSV can fit most of the scores of tests set well, and the predicted scores are closer to the actual scores in the case of large fluctuation. The strong fluctuation of actual scores will lead to a large error between the actual scores and the predicted scores of GPSV and GPLSV, while the predicted scores of GPWLSV and GPLWLSV are stable relatively. This is due to the powerful optimization ability of WOA.Figure 3
Actual scores and predicted scores of various frameworks. (a) GPSV, (b) GPLSV, (c) GPWLSV, and (d) GPLWLSV.
(a)(b)(c)(d)Specifically, WOA optimizes the penalty factor and kernel parameters of LSSVRM, and we only need to adjust two parameters to improve the generalization ability and prediction accuracy of LSSVRM. The prediction accuracy of GPLWLSV is higher than that of GPWLSV for two reasons. Firstly, WOA adjusts parameters through a random system, and excessive reliance on a random system limits GPWLSV. Secondly, WOA is subject to coefficient vector. When the number of iterations reaches half of the maximum set earlier, WOA will lose its ability to jump out of the local optimum. At this time, WOA may fall into local optimum to some extent, which will lead to inaccurate prediction results. LWOA replaces the random adjustment of WOA with Levy flight; LWOA improves the convergence mode of coefficient vector of WOA; LWOA optimizes the ability of WOA to jump out of local optimum and thus improves the generalization ability and convergence speed of WOA. As a result, the prediction accuracy of GPLWLSV is higher than GPWLSV.


### 3.3. Discriminative Brain Regions

Node efficiency is one of the measures of node centrality. The higher the node efficiency is, the stronger the capacity of information transmission of the node is and the more important its position in the network is. The higher the efficiency of a node is, the stronger the information transmission capacity of the node is and the more important the node is in the network [34]. For the sake of finding out the key brain regions which affect the cognitive function of patients with ESRD [35], the multiple regression method was taken to analyze the relationship between the node efficiency and clinical scores of ESRD patients [36]. The results showed that the node efficiency of the left amygdala (AMYG.L) was negatively correlated with the clinical scores (r = −0.562, P=0.014) significantly, as shown in Figure 4(a), while the node efficiency of the right parahippocampal gyrus (PHG.R) was positively correlated with the clinical scores (r = 0.551, P=0.035) significantly, as shown in Figure 4(b). In addition, no correlation was found between node efficiency and clinical scores in other brain regions (P>0.05). The amygdala is mainly involved in emotional processing [37], and previous studies have also shown that depression of varying degrees is common in patients with ESRD [38].Figure 4
Relationship between node efficiency and clinical scores.
(a)(b)Given the above, it can be speculated that the amygdala node efficiency may be related to the depressed mood of ESRD patients, and the emotional instability and intense reaction of patients may lead to the increase of amygdala node efficiency. Qin et al. [39] found that patients with type 2 diabetes had multiple brain regions with increased node efficiency. This phenomenon is explained as the compensatory mechanism of the network. The right parahippocampal gyrus is associated with learning and memory functions [40], and its reduction of node efficiency may be connected with cognitive impairment in ESRD patients.To find out the relationship between each brain region and the cognitive function of ESRD patients, the node efficiency of 90 brain regions was calculated and contrastedd between the patients and the normal controls. The brain regions with significant differences were found; the result is shown in Figure 5. Blue nodes represent the brain regions with reduced node efficiency in ESRD patients compared to the normal controls. These included left insula (INS.L), right insula (INS.R), left median cingulate and paracingulate gyrus (DCG.L), right median and paracingulate gyrus (DCG.R), left hippocampus (HIP.L), right hippocampus (HIP.R), left parahippocampal gyrus (PHG.L), right parahippocampal gyrus (PHG.R), right transverse temporal gyrus (HES.R), left superior temporal gyrus (STG.L), left temporal pole: superior temporal gyrus (TPOsup.L), right temporal pole: superior temporal gyrus (TPOsup.R), left temporal pole: middle temporal gyrus (TPOmid.L), and right temporal pole: middle temporal gyrus (TPOmid.R). The red nodes represent the brain regions with increased node efficiency in ESRD patients compared to the normal controls. These include left amygdala (AMYG.L), right amygdala (AMYR.R), right calcarine fissure and surrounding cortex (CAL.R), left cuneus (CUN.L), right cuneus (CUN.R), left superior occipital gyrus (SOG.L), and left middle occipital gyrus (MOG.L).Figure 5
Brain regions with statistically significant differences in node efficiency between the two groups. (a) Left lateral view, (b) right lateral view, (c) left lateral view, (d) right lateral view, and (e) dorsal view of the whole brain.The brain is a highly optimized complex system capable of integrating information to deal with changes in cognitive needs. According to research findings, the brain regions with reduced node efficiency in ESRD patients are mainly located in the paralimbic network. The hippocampus (HIP) and parahippocampal gyrus (PHG) are related to learning and memory functions [40]; meanwhile, the median and paracingulate gyrus (DCG) are involved in cognitive control function [41]. Chou et al. [42] applied DTI technology to view the changes in brain functional networks in ESRD patients, and the results also found that the node efficiency in the paralimbic network was significantly decreased compared with the normal controls. Notably, ESRD patients had higher node efficiency in brain regions associated with visual networks, which are responsible for processing visual information. These regions are about attention, visual memory, and other neurocognitive functions [43]. Accordingly, it can be speculated that the increased node efficiency in these brain regions may also be a compensatory mechanism [39]. In summary, node efficiency can effectively distinguish ESRD patients from the normal controls, thus achieving accurate classification.


## 3.1. Experimental Settings

Table2 shows the AUC of topology attribute parameters of the brain functional networks of ESRD patients and normal controls calculated by GTA. Within the whole matrix sparsity threshold range, the AUC of γ, σ, and Elocal in ESRD patients was significantly lower than that in normal controls, with statistical significance (P<0.05). Nevertheless, there were no significant differences (P>0.05) in the AUC of λ, Cp, Lp, and Eglobal.Table 2
Comparison of AUC of global topological parameters between two groups (mean ± SD).
ParameterESRD patients (n = 50)Normal controls (n = 40)tPγ0.646 ± 0.0710.669 ± 0.056−1.7140.004λ0.326 ± 0.1230.329 ± 0.0240.7100.480σ0.589 ± 0.0650.607 ± 0.065−1.3840.008Cp0.175 ± 0.0140.176 ± 0.016−0.2960.768Lp0.543 ± 0.0270.553 ± 0.069−0.9660.337Eglobal0.171 ± 0.0060.170 ± 0.011−0.4770.635Elocal0.230 ± 0.0070.241 ± 0.006−0.2530.001Table3 shows the corresponding weight of each feature in the feature set. The proportion of AUC of Elocal was the highest, up to 65.31%. The proportion of AUC of Eglobal was next to Elocal, accounting for 30.89%. The proportion of AUC of γ, λ, σ, Cp, and Lp was less than 5%. The proportion of γ was the lowest, only 0.03%.Table 3
Weight of AUC of global topological parameters of ESRD patients in feature set.
Parameterγ (%)λ (%)σ (%)Cp (%)Lp (%)Eglobal (%)Elocal (%)wi × 100%0.030.070.170.283.2530.8965.31The experimental results can be seen from Tables2 and 3. The AUC of Elocal in ESRD patients was significantly lower than that in the normal controls, with statistically significant differences (P<0.05). Meanwhile, the AUC of Elocal accounted for the highest proportion in the feature set. Obviously, the AUC of Elocal was selected as the feature to construct GPSV [30], GPLSV [31], GPWLSV [32], and GPLWLSV for predicting the clinical scores of cognitive functions in ESRD patients [33].The prediction accuracy of various frameworks was evaluated by a tenfold cross-validation method. Firstly, the AUC of Elocal of 50 ESRD patients and the corresponding clinical scores were collected as data set calledD. Secondly, D was divided into 10 mutually exclusive subsets of the same size, D = D1∪D2∪ D3∪D4∪D5∪D6∪D7∪D8∪D9∪D10, and Di∪Dj = ∅. Each subset Di was separated from D through stratified sampling, which ensures consistency of data distribution. Each subset Di contains five samples. Taking the union of 9 subsets as the training set and the remaining subset as the test set, 10 groups of training sets and test sets are formed. Different frameworks were trained and tested for 10 times, and the average of 10 groups of test results was calculated.


## 3.2. Experimental Results

The root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) were selected as the testing standards of prediction accuracies. Table4 shows the prediction accuracies of some representative frameworks. All the comparable predictive frameworks were evaluated by tenfold cross-validation method. The test results show that the prediction accuracy of GPLWLSV is higher than that of GPSV, GPLSV, and GPWLSV. The average of RMSE of GPLWLSV was 2.40, 1.04, 0.93, and 0.46 points lower than that of GPSV, GPLSV, and GPWLSV, respectively. The average of MAE of GPLWLSV was 2.06, which was 0.82, 0.74, and 0.42 points lower than that of GPSV, GPLSV, and GPWLSV, respectively. MAPE can reflect the relative errors of frameworks better than MAE. The average of MAPE of GPLWLSV was 9.83%, lower than 10%, which was 4.10%, 3.75%, and 1.93% lower than that of GPSV, GPLSV, and GPWLSV, respectively. The bar chart in Figure 2 shows that the prediction accuracy of GPLWLSV is better than that of GPSV, GPLSV, and GPWLSV intuitively.Table 4
Prediction accuracies of comparable frameworks.
Prediction frameworkThe test groupRMSEMAEMAPEGPSV13.50642.38350.1364GPSV23.4823.09510.1454GPSV33.26482.84990.1267GPSV43.61543.17920.1542GPSV52.55082.29650.1109GPSV62.27241.79190.079GPSV74.1153.68250.1788GPSV82.82042.6590.1277GPSV95.72474.05230.1999GPSV103.07072.78520.1334Average3.44232.87750.1392GPLSV15.26344.1970.2149GPLSV24.66283.18990.1305GPLSV33.48192.81210.1491GPLSV43.16112.78370.1385GPLSV52.10461.88840.0898GPLSV62.1622.12460.0994GPLSV73.35453.01990.1524GPLSV83.98913.7250.18GPLSV92.28051.77330.0812GPLSV102.88792.48380.1214Average3.33482.79980.1357GPWLSV14.24083.47820.1695GPWLSV23.32932.81780.124GPWLSV33.06572.41780.1188GPWLSV43.30272.98390.1466GPWLSV51.31251.13880.0529GPWLSV61.92721.41140.0618GPWLSV72.94122.71960.1311GPWLSV82.78992.68650.1284GPWLSV92.67582.32820.1092GPWLSV103.03132.80140.1335Average2.86162.47840.1176GPLWLSV13.91323.10260.1561GPLWLSV23.25062.83340.1258GPLWLSV32.73072.43340.1075GPLWLSV42.99762.37380.1216GPLWLSV51.21921.11110.053GPLWLSV60.7180.64970.0313GPLWLSV72.51672.11560.1036GPLWLSV82.07091.96360.0914GPLWLSV92.17571.94670.0916GPLWLSV102.41632.03020.1007Average2.40092.0560.0983Figure 2
Average of prediction accuracies of various frameworks. (a) The RMSE and MAE of various frameworks. (b) The MAPE of various frameworks.
(a)(b)Figure3 shows the comparison between the predicted scores of various frameworks and the actual scores. The solid blue line represents the actual scores, and the solid purple line represents the predicted scores. As can be seen from the figure, different from GPSV and GPLSV, GPWLSV and GPLWLSV can fit most of the scores of tests set well, and the predicted scores are closer to the actual scores in the case of large fluctuation. The strong fluctuation of actual scores will lead to a large error between the actual scores and the predicted scores of GPSV and GPLSV, while the predicted scores of GPWLSV and GPLWLSV are stable relatively. This is due to the powerful optimization ability of WOA.Figure 3
Actual scores and predicted scores of various frameworks. (a) GPSV, (b) GPLSV, (c) GPWLSV, and (d) GPLWLSV.
(a)(b)(c)(d)Specifically, WOA optimizes the penalty factor and kernel parameters of LSSVRM, and we only need to adjust two parameters to improve the generalization ability and prediction accuracy of LSSVRM. The prediction accuracy of GPLWLSV is higher than that of GPWLSV for two reasons. Firstly, WOA adjusts parameters through a random system, and excessive reliance on a random system limits GPWLSV. Secondly, WOA is subject to coefficient vector. When the number of iterations reaches half of the maximum set earlier, WOA will lose its ability to jump out of the local optimum. At this time, WOA may fall into local optimum to some extent, which will lead to inaccurate prediction results. LWOA replaces the random adjustment of WOA with Levy flight; LWOA improves the convergence mode of coefficient vector of WOA; LWOA optimizes the ability of WOA to jump out of local optimum and thus improves the generalization ability and convergence speed of WOA. As a result, the prediction accuracy of GPLWLSV is higher than GPWLSV.


## 3.3. Discriminative Brain Regions

Node efficiency is one of the measures of node centrality. The higher the node efficiency is, the stronger the capacity of information transmission of the node is and the more important its position in the network is. The higher the efficiency of a node is, the stronger the information transmission capacity of the node is and the more important the node is in the network [34]. For the sake of finding out the key brain regions which affect the cognitive function of patients with ESRD [35], the multiple regression method was taken to analyze the relationship between the node efficiency and clinical scores of ESRD patients [36]. The results showed that the node efficiency of the left amygdala (AMYG.L) was negatively correlated with the clinical scores (r = −0.562, P=0.014) significantly, as shown in Figure 4(a), while the node efficiency of the right parahippocampal gyrus (PHG.R) was positively correlated with the clinical scores (r = 0.551, P=0.035) significantly, as shown in Figure 4(b). In addition, no correlation was found between node efficiency and clinical scores in other brain regions (P>0.05). The amygdala is mainly involved in emotional processing [37], and previous studies have also shown that depression of varying degrees is common in patients with ESRD [38].Figure 4
Relationship between node efficiency and clinical scores.
(a)(b)Given the above, it can be speculated that the amygdala node efficiency may be related to the depressed mood of ESRD patients, and the emotional instability and intense reaction of patients may lead to the increase of amygdala node efficiency. Qin et al. [39] found that patients with type 2 diabetes had multiple brain regions with increased node efficiency. This phenomenon is explained as the compensatory mechanism of the network. The right parahippocampal gyrus is associated with learning and memory functions [40], and its reduction of node efficiency may be connected with cognitive impairment in ESRD patients.To find out the relationship between each brain region and the cognitive function of ESRD patients, the node efficiency of 90 brain regions was calculated and contrastedd between the patients and the normal controls. The brain regions with significant differences were found; the result is shown in Figure 5. Blue nodes represent the brain regions with reduced node efficiency in ESRD patients compared to the normal controls. These included left insula (INS.L), right insula (INS.R), left median cingulate and paracingulate gyrus (DCG.L), right median and paracingulate gyrus (DCG.R), left hippocampus (HIP.L), right hippocampus (HIP.R), left parahippocampal gyrus (PHG.L), right parahippocampal gyrus (PHG.R), right transverse temporal gyrus (HES.R), left superior temporal gyrus (STG.L), left temporal pole: superior temporal gyrus (TPOsup.L), right temporal pole: superior temporal gyrus (TPOsup.R), left temporal pole: middle temporal gyrus (TPOmid.L), and right temporal pole: middle temporal gyrus (TPOmid.R). The red nodes represent the brain regions with increased node efficiency in ESRD patients compared to the normal controls. These include left amygdala (AMYG.L), right amygdala (AMYR.R), right calcarine fissure and surrounding cortex (CAL.R), left cuneus (CUN.L), right cuneus (CUN.R), left superior occipital gyrus (SOG.L), and left middle occipital gyrus (MOG.L).Figure 5
Brain regions with statistically significant differences in node efficiency between the two groups. (a) Left lateral view, (b) right lateral view, (c) left lateral view, (d) right lateral view, and (e) dorsal view of the whole brain.The brain is a highly optimized complex system capable of integrating information to deal with changes in cognitive needs. According to research findings, the brain regions with reduced node efficiency in ESRD patients are mainly located in the paralimbic network. The hippocampus (HIP) and parahippocampal gyrus (PHG) are related to learning and memory functions [40]; meanwhile, the median and paracingulate gyrus (DCG) are involved in cognitive control function [41]. Chou et al. [42] applied DTI technology to view the changes in brain functional networks in ESRD patients, and the results also found that the node efficiency in the paralimbic network was significantly decreased compared with the normal controls. Notably, ESRD patients had higher node efficiency in brain regions associated with visual networks, which are responsible for processing visual information. These regions are about attention, visual memory, and other neurocognitive functions [43]. Accordingly, it can be speculated that the increased node efficiency in these brain regions may also be a compensatory mechanism [39]. In summary, node efficiency can effectively distinguish ESRD patients from the normal controls, thus achieving accurate classification.


## 4. Discussion

The purpose of this study was to look for imaging markers related to the decline of cognitive function in ESRD patients, which could be modeled to predict the clinical scores of the patients, and the clinical scores could determine the current stage of cognitive function. As a consequence, the GPLWLSV is constructed. The experimental results show that the GPLWLSV achieves the optimal effect of prediction than comparable frameworks, indicating its effectiveness. The GPLWLSV has three main advantages.First of all, the imaging markers are accurate. By right of the principle of GP, GPLWLSV found the imaging markers related to the decline of cognitive function in ESRD patients. Compared with the normal controls, ESRD patients showed a significant decrease in the AUC ofγ, σ, and Elocal at multiple sparsity. This result is consistent with a previous study based on resting-state fMRI networks [44]. γ and Elocal mainly affect specific information processing or fault tolerance rate of the network [45]. The decrease of γ and Elocal indicates that information transmission efficiency of different brain regions is reduced in ESRD patients, which impairs the ability to manage the brain potentially [46]. This provides a new perspective and potential imaging markers for understanding the underlying pathophysiological mechanisms of cognitive impairment in ESRD patients.Secondly, the operation speed is faster. GPLSV has a shorter operation time than GPSV, although the prediction accuracy is similar. This is because LSSVRM is an improvement on SVRM. It changes inequality constraints into equality constraints in SVRM and transforms solving quadratic programming problems into solving linear equations, speeding up the operation speed greatly. In a word, the LSSVRM enables GPLWLSV to run more efficiently.Finally, the prediction accuracy is higher. The prediction accuracy of GPWLSV and GPLWLSV is higher than that of GPSV and GPLSV significantly. This is because WOA possesses powerful optimization ability; it optimized the strategy of selecting two parameters in the kernel function of SVRM and LSSVRM, thus improving the prediction accuracy of GPWLSV and GPLWLSV. However, WOA still has deficiencies. Firstly, WOA adopts a random system to adjust parameters and relies on random system excessively. Secondly, WOA is limited by the coefficient vector. When the number of iterations reaches half of the maximum iterations set earlier, WOA will lose its ability to jump out of local optimum. In brief, WOA may fall into local optimum to some extent, leading to low prediction accuracy. GPLWLSV chose Levy flight instead of WOA to adjust parameters. LWOA improved the convergence mode of coefficient vector of WOA and hoisted the capacity of WOA to jump out of local optimum, and the convergence speed of WOA was improved. To sum up, the prediction accuracy of GPLWLSV was higher than comparable frameworks.The prediction framework still has room for further improvement. The main reason is that the size of the data set in this paper is small, which affects the training effectiveness of GPLWLSV to a certain extent. Researches show that DTI, DWI, DKI, and other biomarkers of different modes are complementary to each other to some extent. Compared with single-mode analysis, multimode analysis can often achieve better results of prediction by integrating complementary information of different modes [47]. In future studies, more clinical data of different modes will be collected, leading to a more in-depth and comprehensive study on a larger data set.


## 5. Conclusion

In this paper, an integrated prediction framework with GPLWLSV was constructed to predict the clinical scores of cognitive functions in ESRD patients. According to the principle of GP, the framework found the imaging markers related to the decline of cognitive function in ESRD patients; meanwhile, GP determined that the AUC of Elocal accounted for the highest proportion of all the features. Elocal had a great influence on the cognitive function of ESRD patients. The inequality constraints in SVRM are simplified, and the operation speed is improved. LWOA was introduced to optimize the parameter selection strategy in LSSVRM, and the prediction accuracy was improved. In clinical diagnosis, it is often necessary to analyze the clinical scores of a large number of patients with ESRD in order to judge the stage of cognitive function. In general, our framework can obtain imaging markers related to the decline of cognitive function accurately and give consideration to work efficiency and accuracy simultaneously.



---

*Source: 8124053-2022-08-09.xml*