Prediction and Optimization of Electrospun Polyacrylonitrile Fiber Diameter Based on Grey System Theory

This paper provides a new method for predicting the diameter of electrospun nanofibers. Based on the grey system theory, the effects of polyacrylonitrile mass fraction, voltage, flow rate, and receiving distance on fiber diameter were studied. The GM(1,1) (grey model) model and DNGM(1,1) (The DNGM (1,1) model is based on the whitening differential equation using parameters to Directly estimate the approximate Non-homogeneous sequence Grey prediction Model) model were established to predict fiber diameter by a single-factor change, and the results showed high prediction accuracy. The multi-variable grey model MGM(1,n) (MGM(1,n) is a Multivariate Grey prediction Model) was used for prediction of fiber diameter when multiple factors change simultaneously. The results showed that the average modeling fitting error is 8.62%. The background value coefficients of the MGM(1,n) model were optimized, the average modeling fitting error was reduced to 1.01%, and the average prediction error was reduced to 1.33%. Combining the fractional optimization with the background-value coefficient optimization, the optimal background-value coefficient α and the order r were selected. The results showed that the average modeling fitting error is 0.85%, and the average prediction error is 0.38%. The results demonstrate that the grey system theory can effectively predict the diameter of polyacrylonitrile electrospinning fibers with high prediction accuracy. This theory can increase the control of nanofiber diameters in production.


Introduction
Nanofiber membranes prepared via electrospinning have extremely high porosity, high surface-to-volume ratios [1], and excellent surface properties leading to broad application prospects in textiles, filtration, medicine, sensing, and other fields [2][3][4][5][6][7]. Its properties also depend on different electrospun parameters. The fiber diameter seriously affects the pore size per unit area and affects the specific area of the fiber membrane, thus affecting the filtration and gas permeability of the fiber membrane [8]. Therefore, the prediction of fiber diameter and its control is of great significance for the performance and improvement of nanofiber membranes [9].
At present, two methods are used for predicting and controlling the diameter of electrospun nanofibers. One is to study the relationship between the jet radius and its moving distance in the direction of stretching [10]. The other is to study the relationship between fiber diameter and process parameters [11][12][13].
Spivak [10] comprehensively studied the factors of inertial force, hydrostatic force, viscous resistance, electric-field force, and surface tension. The relationship scaled between the jet-flow model, the prediction accuracy of the model for the research object is continuously improved. Model accuracy is the key and difficult point in building a model. Many researchers have proposed methods to improve accuracy often using background value reconstruction [30], background value coefficient optimization [31], initial value optimization [32], and the new information optimization model [33]. When building a grey prediction model, in order to facilitate the model solution, an accumulation sequence is often represented by a sequence of mean generations. It is common to say that the background value α is taken as 0.5. However, according to the actual situation, selecting the appropriate background value coefficient α can effectively improve the prediction accuracy.
The models mentioned above are all integer-order models, but in practice many objects satisfy the characteristics of fractional order, and the methods of fractional-order accumulation can better reveal their characteristics. Wu and Liu proposed a fractional-order grey model FGM(1,1) (FGM(1,1) is a Fractional-order Gray prediction Model), the prediction error of which can be reduced by using fractional power in grey generation [34,35].
In this paper, the principle and control of electrospinning are analyzed. The effects of main factors such as mass fraction, voltage, flow rate, and receiving distance on fiber diameter were studied, and an effective prediction model is found for the prediction of electrospinning fiber diameter. The GM(1,1) model was used to predict the diameter of single-variable electrospinning fibers. The original system sequence prediction accuracy was explored in the multivariate prediction of electrospinning fiber diameter based on the MGM(1,n) model. In the modeling principle, the method of optimizing the background value is adopted to correct the prediction fitting error caused by the selection of the background value parameters. To further improve the prediction accuracy, the fractional order idea is introduced into the modeling mechanism. The fractional order and background value optimization are combined to establish an FMGM(1,n) (FMGM(1,n) is a Fractional order Multivariate Grey prediction Model) model that is more adaptive to the system sequence.

Preparation of PAN and Testing
The appropriate amount of DMF was measured in a beaker. The PAN powder was then dried in a low temperature oven for 2 h. The PAN powder was slowly added to the solvent and transferred to a three-necked flask. The three-necked flask was placed in a constant temperature water bath at 60 • C, and the sample was stirred with a magnetic stirrer for 4 h. The stirred solution exhibited a stable bright yellow solution. The solution was allowed to stand for a period of time for defoaming treatment, and then transferred to a sealed beaker. The beaker containing the solution was labeled.
The control-variable method was adopted to explore the effect of mass fraction factors on fiber diameter. The voltage was set to 18 kV, the distance between the nozzle and the collecting plate was set to 16 cm, the inner diameter of the nozzle was selected to be 0.57 mm, and the flow rate of the solution was set to 0.5 mL/h. PAN electrospinning has a good spinning mass fraction ranging from about 8% to 16%. When the mass fraction is low, the bead-like structure is likely to occur, and when the mass fraction is high, the spun fiber membrane is not easily formed. Therefore, the mass fraction of this experiment was selected to be in the middle of the range. The electrospinning experiment was carried out with a solution having a PAN mass fraction of 9%, 10%, 11%, 12%, 13%, 14%, 10.5%, 11.5%, and  12.5%. When studying the voltage factor, flow rate factor, and receiving distance factor, the control variable method was similarly used, and the fixed PAN mass fraction was 12%. When one factor was studied, the other three factors remained unchanged, and only the value of the factor being studied was changed. The parameter values of specific factors are detailed in Table 1. Table 1. Experimental data sheet for the effect of individual factors on fiber diameter. L: the distance between the nozzle and the collecting plate; M: the mass fraction of polyacrylonitrile (PAN) solution; S: the flow rate of the solution at the nozzle; V: the loading voltage at the nozzle; x (0) i : the corresponding value of nanofiber diameter.

Sample Detection
The fiber membrane prepared according to the above experimental conditions was placed in a petri dish and vacuum dried at a low temperature for 8 h. The nanofiber membrane was cut into a sample of appropriate size and fixed on the stage with a conductive paste. The electrospun fiber membrane was observed using an environmental scanning electron microscope (ESEM, FEI Corporation, Hillsborough, Oregon, USA). Ten fibers were randomly selected in the ESEM image, and the diameter of the fibers were measured by the image analysis function in Photoshop software (Version 13.0, Adobe Systems Software Ireland Ltd, San Jose, California, USA), and the average value was calculated as the fiber diameter of the experimental spinning. The data of the fiber diameter under specific factors are shown in Table 1.

Single-Factor Model Establishment and Experimental Verification
Under the premise of ensuring the quality of the electrospun fiber membrane, the parameter values of various factors were set within the maximum range. For example, when the loading voltage was lower than 10 kV, the fiber membrane on the collecting plate was clustered although the fiber could be spun. When the loading voltage was 10-12 kV, the ESEM showed more bead-like structures, which seriously affects the quality of the fiber membrane. The bead-like structure is shown in Figure 1. Electrospinning is unstable when the loading voltage is higher than 22 kV. From this, we determined that the normal range of the spinning voltage is 12 kV to 22 kV. Some fiber morphology ESEM images between 12 kV and 22 kV are shown in Figure 2.
The loading voltage at the nozzle was recorded as V unit (kV), the mass fraction of PAN solution was recorded as M unit (%), the flow rate of the solution at the nozzle was recorded as S unit (mL/h), the distance between the nozzle and the collecting plate was recorded as L unit (cm), and x (0) i indicates the corresponding value of nanofiber diameter (nm). In Table 1, No. 1-6 are the parameter values of each factor, which is the maximum range value under the premise of ensuring the quality of the fiber membrane; parameter value of each factor increased sequentially. This part of the data was used to build a grey model. In grey prediction theory, the theory uses at least two sets of data to test the accuracy of the model. In this paper, three sets of data were used to test the accuracy. The data in sequence numbers 7-9 are interpolation data within the respective factors, and this part of the data was used to test the model's prediction accuracy.  The loading voltage at the nozzle was recorded as V unit (kV), the mass fraction of PAN solution was recorded as M unit (%), the flow rate of the solution at the nozzle was recorded as S unit (mL/h), the distance between the nozzle and the collecting plate was recorded as L unit (cm), and (0) i x indicates the corresponding value of nanofiber diameter (nm). In Table 1, No. 1-6 are the parameter values of each factor, which is the maximum range value under the premise of ensuring the quality of the fiber membrane; parameter value of each factor increased sequentially. This part of the data was used to build a grey model. In grey prediction theory, the theory uses at least two sets of data to test the accuracy of the model. In this paper, three sets of data were used to test the accuracy. The data in sequence numbers 7-9 are interpolation data within the respective factors, and this part of the data was used to test the model's prediction accuracy. The loading voltage at the nozzle was recorded as V unit (kV), the mass fraction of PAN solution was recorded as M unit (%), the flow rate of the solution at the nozzle was recorded as S unit (mL/h), the distance between the nozzle and the collecting plate was recorded as L unit (cm), and (0) i x indicates the corresponding value of nanofiber diameter (nm). In Table 1, No. 1-6 are the parameter values of each factor, which is the maximum range value under the premise of ensuring the quality of the fiber membrane; parameter value of each factor increased sequentially. This part of the data was used to build a grey model. In grey prediction theory, the theory uses at least two sets of data to test the accuracy of the model. In this paper, three sets of data were used to test the accuracy. The data in sequence numbers 7-9 are interpolation data within the respective factors, and this part of the data was used to test the model's prediction accuracy.

Grey Relational Analysis
In the grey theory, if two factors change with the development of the system, then the trend of the two factors has a high similarity indicating that the correlation between the two is higher. On the contrary, if the two factors change with the development of the system, the trend has a low similarity the degree of correlation is low [36]. The four above groups of factors underwent relational analysis to explore the correlation between various factors and fiber diameter. Taking the relational analysis of voltage factor and fiber diameter as an example, the raw data are expressed as the following: Step 1: The physical units between the original data sequences are different. In order to make the data sequences comparable, it is necessary to preprocess the data sequence, eliminate the order of magnitude difference and dimension of the series, and retain the trend of data changes. The raw data are processed using the mean conversion equation as shown in Equations (1) and (2).
Step 2: As shown in Equation (3), the absolute value ∆ ij (t) of the difference corresponding to the parameter value of each factor is solved. We then selected the maximum value ∆max and the minimum value ∆min in the absolute value and solved the correlation coefficient L ij (t) using Equation (4): In Equation (4), k takes a value of 0.5, and the calculated absolute value has a maximum value of ∆max = 0.2426 and a minimum value of ∆min = 0.0492. Step 3: The final step solves the value of the degree of grey relation: From Equation (5), the grey relational degree between the voltage and the fiber diameter is Similarly, using the above formula, the value of the grey relational degree between the other three factors and the fiber diameter is obtained. The results are as follows: the grey relational degree value of the mass fraction R 2 = 0.5670, the grey relational degree of the flow rate R 3 = 0.6631, and the grey relational degree of the distance factor R 4 = 0.6752. Analysis and comparison of the above four factors and fiber-diameter grey relational degree indicated that the grey relational degree of each factor was greater than 0.5, indicating that these four factors have a greater impact on fiber diameter. The results show that it is reasonable to select these four factors to study the diameter of nanofibers and can be used for the subsequent grey prediction theory modeling.

Establish GM(1,1) Model and Verification
The GM(1,1) model was established by taking the voltage factor as an example when studying a single factor. The value of the fiber diameter corresponding to the uniform change of the voltage is recorded as a one-dimensional original data sequence. Its expression form is: In order to get a better system of growth, random interference from the original data sequence needs to be eliminated. The raw data is subjected to an accumulation process to obtain X (1) sequence. The accumulation process is as shown in Equation (6): The background value sequence Z (1) is obtained by processing X (1) with Equation (7). Here, α is the background value weight coefficient. The value of α is in the range (0, 1), usually for the convenience of calculation, the average generation series modeling calculation is generally adopted. That is, the value of α is taken as 0.5 [37].
In Equation (9), a is the main variable parameter, and b is the background value of the model. The parameter matrix of the a and b values is determined by the least squares method. The solution is as follows: In Equation (11), the matrices B and Y are written as follows: The solved parameters a and b are substituted into the whitenization Equation (8) to obtain the GM(1, 1) time response function Equation (13). On this basis, the discretization Equation (14) is obtained.
Substituting the initial values X (1) = X (0) and k = 0 into Equation (14), the value of the undetermined parameter C can be obtained. On this basis, an accumulated analog value sequence is obtained. The adjacent two values of the ∧ X (1) sequence are subtracted as shown in Equation (15).
The final reduction yields the simulated and predicted value sequence corresponding to the original data. The model error is evaluated using the average relative error MAPE (The mean absolute percentage error (MAPE) is a measure of prediction accuracy of a forecasting method in statistics) as shown in Equation (16).
In order to guarantee the quality of the product, the prediction error must be within 3%. If the average relative error exceeds the target error, then a corresponding model improvement is required. The Matlab program of the GM(1,1) model can be written such that the original data sequence under each factor is brought into the program. The fiber diameter fitting value and predicted value corresponding to each factor were obtained. The results are shown in Tables    It can be seen from Tables 2 and 3 that the fitting error and prediction error of the voltage factor are the smallest. It is shown that the fitting value and the measured value of the voltage factor have the best fitting effect. The average relative fitting error is within 2.5% for the four factors above except for the mass fraction factor; the average relative prediction error is less than 2%. This shows that the model has higher prediction accuracy. The forecasting requirements can be met in actual production activities. The average modeling fit accuracy is 3.35% for the mass fraction factor sequence, which exceeds the expected error target. The model needs research and improvement, and this further improvement is described in Section 3.3.

Model Improvement
The previous section showed that the prediction error of the mass fraction factor exceeds the target value. This might be because the GM(1,1) modeling mechanism and model structure are difficult to obtain satisfactory fitting and prediction accuracy for if the original data are an approximate non-homogeneous sequence. Therefore, it is necessary to improve the modeling mechanism of the grey model to adapt it to the current system sequence.
The study found that the DNGM(1,1) (The DNGM (1,1) model is based on the whitening differential equation using parameters to Directly estimate the approximate Non-homogeneous sequence Grey prediction Model) model is suitable for approximate non-homogeneous data sequences. The model is based on the solution of the whitenization equation applicable to the approximate non-homogeneous model. Kramer's law is then used to solve the corresponding coefficients, and the undetermined parameter values of the differential equations are obtained. A specific expression under the approximate non-homogeneous sequence-prediction model is obtained. Modeling avoids the error caused by directly using the difference equation to estimate the parameters, but it obtains a response function from the differential equation.
The DNGM(1,1) model is similar to the GM(1,1) modeling process except for the whitenization equation (shown in Equation (17)). The undetermined parameters in the equation are obtained from the intermediate parameters The intermediate parameters α, β and γ are obtained via Equations (18) and (19): On the basis of solving the intermediate parameters, the undetermined parameters a, b and c of the white differential equation are solved by Equation (20).
We can then solve the response expression of the white differential equation. The adjacent two values of the responsive sequence are subtracted; finally, the analog value Equation (21) under the DNGM(1,1) model is obtained.
The model error was evaluated using the average relative error MAPE. The corresponding Matlab program leads to the following in Table 4:  Table 4 shows that the fitting error of the mass fraction factor under the DNGM(1,1) model is 0.48%, and the prediction error is 1.62%. These conditions satisfy the expected target error requirement. The improved model is successful in predicting the diameter of the mass fraction factor.
In summary, under the GM(1,1) model, the average prediction error of the voltage factor is 0.81%, the average prediction error of the flow factor is 1.31%, and the average prediction error of the distance factor is 1.82%. Under the improved DNGM(1,1) model, the average prediction error of the mass fraction factor is 1.62%. The average prediction error of each factor is less than 3%. This indicates that the GM(1,1) model and the DNGM(1,1) model are successful in predicting the diameter of the electrospun nanofiber. Using any parameter within the allowable range of each factor, the corresponding nanofiber diameter can be accurately predicted.

Experimental Verification and Optimization of Multi-Variable Fiber Diameter Prediction
When multiple variables change simultaneously in actual production activities, the multi-variable grey model MGM(1,n) can be used to solve the fiber diameter prediction problem. To ensure the comparability of the electrospinning fiber diameter prediction experiment, the experiment needs to be carried out at the same formulation concentration. The concentration of the solution with a mass fraction of 12% is suitable, the surface of the fiber is smooth, and the fiber diameter distribution is uniform than other mass fractions. Therefore, the experimental conditions of a PAN mass fraction of 12% were selected to study the effect of the distance between the nozzle and the collector, the solution flow rate at the nozzle, and the loading voltage on the fiber diameter. Sequence numbers 1-6 in the following four-dimensional data were used for modeling, and sequence numbers 7-9 were used to test the prediction accuracy. The sequence of fiber diameter measurement values in Table 5 is recorded as x (0) 1 (nm), the voltage factor sequence is recorded as x (0) 2 (kV), and the distance between the head and the collector is recorded as x   Table 5 records the data on the changes in fiber diameter as these three different factors change simultaneously.  The multi-variable grey model MGM(1,n) is an extension of the GM(1,1) model under the n variables. Similar to the GM(1,1) modeling process, only the main modeling part of MGM(1,n) is presented here. The first-order ordinary differential equations of the MGM(1,n) model are shown in Equation (22).
In the MGM(1,n) model, n is the number of variables. The equation parameter used for solving and process transformation is similar to the GM(1,1) model and will not be described here. The sequence of analog values is expressed as Equation (23).
The analog value sequence is inversely accumulated to obtain a predicted value corresponding to the original data. The accuracy of the model was evaluated by taking the average relative error MAPE.
Based on the multi-variable grey theory, the MGM(1,4) model was established and recorded as Model 1. We then used Matlab for solving and analyzing the data in Table 5 to get the relative error of fiber diameter fitting and prediction (Figure 3).  In Figures 3, 4 and 5, term B represents the relative error of the fit, and C represents the predicted relative error. Figure 3 shows that the maximum error is 41.33%. Further calculation shows that the average fitting error of the MGM(1,4) model is 8.62%, which does not meet the required accuracy requirements. Model improvements are needed to accommodate the current multivariate prediction of nanofiber diameter problems.
The selection of the background value coefficientα affects the sequence of system background values, thereby affecting the sequence of predicted values. The automatic optimization and weighting method is adopted to find the best background value coefficient from the internal modeling principle; this will improve the accuracy of the model.
The automatic optimization-weight method is as follows: the initial value of α is 0, andα increases by a small amount of α Δ ( 0 α Δ > ). We then identify the value of MAPE and select the value of α when the MAPE value is the smallest. This value is the best background value coefficient. The MGM(1,4) model that optimizes the background value coefficients is denoted as Model 2. Its relative error is shown in Figure 4. In Figures 3-5, term B represents the relative error of the fit, and C represents the predicted relative error. Figure 3 shows that the maximum error is 41.33%. Further calculation shows that the average fitting error of the MGM(1,4) model is 8.62%, which does not meet the required accuracy requirements. Model improvements are needed to accommodate the current multivariate prediction of nanofiber diameter problems.
The selection of the background value coefficient α affects the sequence of system background values, thereby affecting the sequence of predicted values. The automatic optimization and weighting method is adopted to find the best background value coefficient from the internal modeling principle; this will improve the accuracy of the model. this will improve the accuracy of the model.
The automatic optimization-weight method is as follows: the initial value of α is 0, andα increases by a small amount of α Δ ( 0 α Δ > ). We then identify the value of MAPE and select the value of α when the MAPE value is the smallest. This value is the best background value coefficient. The MGM (1,4) model that optimizes the background value coefficients is denoted as Model 2. Its relative error is shown in Figure 4.  Figure 4 shows that the accuracy of Model 2 is much higher than that of Model 1, and the average relative error has reached the expected target value. It satisfies the requirements of general quality.
These models are all integer-order models. In fact, the sequence of fractional-order properties is more common, and the fractional-order grey prediction model can better reveal the essence of system The above algorithm was written and calculated in Matlab. The average relative error MAPE is the smallest when the background value coefficient α is 0.5 and the fractional order is 1.006. The relative error probability map is shown in Figure 5.  Figure 5 shows that the further improved Model 3 has a large improvement in the fitting error and the prediction error, and the error is further reduced. The optimized model is closer to the actual sequence and can meet the needs of higher precision.
We next compared the prediction results of the original model and the improved model ( Table  6)   The automatic optimization-weight method is as follows: the initial value of α is 0, and α increases by a small amount of ∆α (∆α > 0). We then identify the value of MAPE and select the value of α when the MAPE value is the smallest. This value is the best background value coefficient. The MGM (1,4) model that optimizes the background value coefficients is denoted as Model 2. Its relative error is shown in Figure 4. Figure 4 shows that the accuracy of Model 2 is much higher than that of Model 1, and the average relative error has reached the expected target value. It satisfies the requirements of general quality.
These models are all integer-order models. In fact, the sequence of fractional-order properties is more common, and the fractional-order grey prediction model can better reveal the essence of system sequences [38]. In fact, in order to obtain better prediction accuracy, it is more effective to optimize the grey model via the fractional order principle.
The fractional-order optimized multi-variable grey prediction model has the following principle: Here, X (0) is the original data sequence, X (r) is the r-order cumulative generation sequence of X (0) , and the expression and production formula of X (r) as follows: (2), · · · , x (r) (n)) (24) x (r) k = 1, 2, . . . , M; i = 1, 2, . . . , n; Γ is a gamma function. The differential equations of the fractional grey model FMGM(1,n) are expressed as: The intermediate transformation and calculation process is similar to the MGM (1, n) model and will not be described here. The r-order sequence of FMGM(1,n) is shown in Equation (27).
The simulation value of the original data can be generated by the r-order subtraction of Equation (27), and the accuracy of the model is evaluated by the average MAPE relative error.
It is known from Model 2 that selecting the appropriate background value coefficient α can improve the prediction accuracy. We next combined the fractional order and background value coefficient optimization as Model 3 and jointly explored a more accurate method. The steps are as follows: Step 1: Set the initial value r = 0 and accumulate the expression that generates X (r) .
Step 2: Set the initial value of the background value coefficient α = 0 and solve the model parameters by the least squares method.
Step 3: Solve the response function and further calculate the simulated value and the predicted value sequence.
Step 4: On the basis of 3, solve the average relative error MAPE.
Step 5: If r < 2 and α < 1, then add a small amount of ∆α and ∆r (∆r > 0) and repeat the above operation.
Step 6: The corresponding r and α values when the MAPE is minimal is the optimization target value.
The above algorithm was written and calculated in Matlab. The average relative error MAPE is the smallest when the background value coefficient α is 0.5 and the fractional order is 1.006. The relative error probability map is shown in Figure 5. Figure 5 shows that the further improved Model 3 has a large improvement in the fitting error and the prediction error, and the error is further reduced. The optimized model is closer to the actual sequence and can meet the needs of higher precision.
We next compared the prediction results of the original model and the improved model (Table 6).  The prediction error can be greatly reduced by selecting appropriate background-value coefficients during modeling ( Table 6). The fractional order provides an effective way to improve the prediction accuracy. This makes the improved model more in line with the actual sequence. The combination of background value coefficient optimization and fractional order optimization can reveal the essence of the system sequence. The method has strong applicability and will be useful for predicting fiber diameter in electrospinning.

Conclusions
In this paper, a new prediction method was proposed in the field of electrospinning fiber-diameter prediction; the method was verified with experimental data. The GM(1,1) model exhibits high prediction accuracy related to single-factor variation on fiber diameter (e.g., loading voltage, solution flow rate, and receiving distance). The average prediction error is within 2%. The error of GM(1,1) in predicting the mass fraction factor sequence is 2.41%. The GM(1,1) prediction model under this factor is improved. The results show that the improved DNGM(1,1) model has a prediction error of 1.62%, and the accuracy improvement is obvious.
The average fitting error of the MGM(1,n) model was 8.62% in terms of the multifactorial influence on fiber diameter. The model was optimized via the background-value coefficient, and the average fitting error was reduced to 1.01%; the average prediction error was reduced to 1.33%, which greatly improves the prediction accuracy. The best background-value coefficients and orders are found by combining fractional-order optimization with background-value-coefficient optimization. The results show that the average prediction error of the further optimized model is only 0.38%. The final improved model shows higher accuracy and provides an effective model for predicting the diameter of the electrospun fibers.
Both the single-variable prediction model and the multi-variable prediction model achieve higher accuracy predictions with less sample data. The multivariate model can achieve effective auxiliary control of the fiber diameter by adjusting multi-factor parameters. This theory has important guiding significance in actual production activities.
This work combined background-value coefficient optimization and fractional-order optimization to enhance the applicability of the grey model in a study of a small number of samples and some unknown factors. At the same time, it is also a successful application of the grey prediction theory in the prediction of the diameter of electrospun fibers. This work adds to the field of application of grey prediction theory.
Author Contributions: Q.Z., L.L. and G.C. participated in the design of this study. Q.Z., L.L. analyzed the data and wrote the paper. Q.Z., L.L., G.C. and Z.D. carried out the literature search, data acquisition, and manuscript preparation, and all authors read and approved the final manuscript.