Prediction of active ingredients in Salvia miltiorrhiza Bunge. based on soil elements and artificial neural network

Sample collection and processing

Roots and rhizosphere soil of S. miltiorrhiza (cultivation or wild) were collected from 25 producing areas in eight provinces of China from mid-November to early December 2007 (Fig. 1, Table S1). The samples were collected by the “S” parallel sampling and multi-point mixing method, each sample was collected at 20 to 25 points, and the rhizosphere soil was collected by the root shaking method, and finally, the samples for analysis were obtained by the quadratic method, of which 2 kg/sample of rhizosphere soil and 5 kg/sample of medicinal parts were retained. After the samples were collected and quickly transported back to the laboratory, the herbs were processed routinely, and the soil samples were air-dried and prepared.

Determination of inorganic elements

The available potassium (K), copper (Cu), zinc (Zn), and manganese (Mn) were determined by atomic absorption spectrophotometry after ammonium bicarbonate-diethylenetriaminepentaacetic acid (AB-DTPA) extraction (Lu, 2000); the available nitrogen (N) was determined by alkaline diffusion method (Lin, 2004); available phosphorus (P) was determined by sodium bicarbonate extraction-molybdenum antimony anti-colorimetric method (Zhao et al., 2020).

Determination of the content of active ingredients

A high-performance liquid chromatographic method was used to determine the contents of water-soluble components and lipid-soluble components (Yang et al., 2010; Yang et al., 2011a). The working parameters for the determination of water-soluble components was as follows: the column was a phenomenex Gemini C¹⁸ column (250  × 4.6 mm, 5 µm, Guangzhou Philomen Scientific Instruments Co., Ltd.); the mobile phases were water-acetonitrile-formic acid (90:10:0.4) (phase A) and acetonitrile (phase B), with gradient elution, phase A: 0∼40 min, 100%∼70%; phase B: 0∼40 min, 0%∼30%; the detection wavelength was 280 nm; the flow rate was 1 mL/min; the column temperature was 25 °C; the injection volume was 20 µL. The working parameters for the determination of liposoluble components were as follows: the column was a Welchrom™ C¹⁸ column (Analytial 4.6  ×  250 mm, 5 µm, Welch Corporation, USA); the mobile phases were methanol (A phase) and water (B phase), with gradient elution, A phase: 0∼25 min, 67%∼67%; 25∼45 min, 67%∼90%; B phase: 0–25 min, 33%–33%; 25–45 min, 33%–10%. The detection wavelength was 270 nm; the flow rate was 1 mL/min; the column temperature was 25 °C; the injection volume was 20 µL.

BP (back propagation) neural network

Bootstrap is a statistical inference method based on resampling and data simulation. Due to the nonlinear nature of small sample sizes and the difficulty of characterizing the overall distribution, the Bootstrap method can be applied to improve the estimation accuracy of the model (Wang et al., 2018). Therefore, the Bootstrap method has the potential to be widely applied to modeling estimations of small sample sizes. The steps are described as follows: (i) perform resampling to select a certain number (given) of samples and to allow repeat sampling; (ii) calculate the given statistics T based on the given samples; and (iii) repeat the above steps N times to gain N number of statistics T (Wang et al., 2018). The back propagation neural network proposed by Rumelhart in 1986 is a multi-layer feedforward network trained according to the error back propagation algorithm and is one of the most widely used neural network models (Rumelhart, Hinton & Williams, 1986). The BP neural network model algorithm consists of two aspects: the forward propagation of the signal and the back propagation of the error. In other words, the actual output is calculated according to the direction from input to output. However, when the actual output contradicts the expected output, the back propagation of error is performed according to the direction from output to input, and the output error of each layer neuron is calculated layer by layer. Then, the weights and thresholds of each layer are adjusted according to the error gradient descent method to achieve the final output value that can be close to the expected value. The model-building process was described in detail as follows.

In this study, a BP neural network model was constructed using the neural fitting tool (nftool) in MATLAB 2019b mathematical software and trained with soil element parameters as input and the content of plant active ingredients as output. The input variables were as follows: Mn, Cu, Zn, N, P, and K; output variables: water-soluble components (danshensu, protocatechuic aldehyde, caffeic acid, rosmarinic acid, salvianolic acid B, salvianolic acid A), and lipid-soluble components (dihydrotanshinone I, cryptotanshinone, tanshinone I, and tanshinone IIA). In the solution process, 70% of the samples (n = 18) were used to obtain samples using the bootstrap method and build a BP neural network model. The remaining 30% of the samples (n = 8) were used to verify the model. First, to prevent the negative impact of different ranges of input variables on the model’s efficiency, the input variables in the model were therefore normalized for the specified range (0–1) (Eq. (1)): (1)${\displaystyle F_j=\frac{X_j-X_{min}}{X_{max}-X_{min}}.}$

In the equation, F_j, X_j, X_min and X_max are the standardized value, original value, minimum value and maximum value of the input variable respectively.

Second, the number of hidden layers and the number of neurons in each hidden layer in a BP neural network impact the overall neural network structure. Currently, many empirical formulas are applied to determine the number of neurons in the hidden layer, and one of these formulas is as follows (Eq. (2)): (2)${\displaystyle h=\sqrt{m+n}+a.}$

In the equation, m (m =6) is the number of nodes in the input layer, n (n =2) is the number of nodes in the output layer, and a (1 ≤a ≤10) is a constant.

According to Eq. (2), the number of nodes in the hidden layer was set as an integer between 4 and 12. Eventually, the number of neurons in the hidden layer was determined to be 8 by iterative trials. A neural network model with 6-8-2 structure was finally established (Fig. 2), in which the input layer consisted of 6 neurons corresponding to the 6 input variables, and the output layer had 2 neurons representing the content of active components in the model.

Third, another problem in establishing neural network models is the choice of network learning or training algorithms. Since the Levenberg–Marquardt algorithm minimizes the sum of the error function of the form (Eq. (3)), thus the Levenberg–Marquardt algorithm has the best performance (R-value) compared to other training algorithms (i.e., Bayesian regularization (BR) and scaled conjugate gradient (SCG)) (Mahmoudi & Mahmoudi, 2014). Therefore, the Levenberg–Marquardt algorithm will be used to train the network. The training epochs, learning rate, and minimum performance gradient were set as 1000, 0.01, and 1e⁻⁷. (3)${\displaystyle E=\frac{1}{2}\sum k{\left(e_k\right)}^2=\frac{1}{2}{∥e∥}^2}$

where e_k is the error in the k^th exemplar or pattern and e is a vector with element e_k.

During the BP training process, a sigmoid function (Eq. (4)) is used to describe the nonlinear relationship between the input and output of each neuron in the hidden layer as follows (Yi et al. 2007). (4)${\displaystyle f\left(x\right)=\frac{1}{1+e^x}.}$

The output y_j of the hidden layer neuron j is calculated by Eq. (5): (5)${\displaystyle y_j=\varnothing \left(\sum _{i=1}^n{W}_{ij}{X}_i+{\theta }_j\right).}$

In the equation, ø(x) represents the activation function of the hidden layer, n is the number of neurons in the input layer, X_i is the input of the input layer neuron i, W_ij is the weight from the input layer neuron i to the hidden layer neuron j, and θ_j is the threshold value of the hidden layer neuron j.

The output z_k of the output layer neuron k are calculated by Eq. (6): (6)${\displaystyle z_k=f\left(\sum _{j=1}^m{V}_{jk}{y}_j+{\alpha }_k\right).}$

In the equation, f (x) represents the activation function of the output layer, and m is the number of neurons in the hidden layer, V_jk is the weight of the hidden layer neuron j to the output layer neuron k, and α _k is the threshold value of the output layer neuron k.

The network weights adjustment is defined as follows (Eq. (7)) (Huang et al., 1996): (7)${\displaystyle \Delta \omega \left(t\right)=-\alpha \frac{\partial E}{\partial \omega \left(t\right)}+\beta \Delta \omega \left(t-1\right)}$

where α and β are assumed constants, called the learning rate and momentum factor, respectively, E is the error function (in which MSE was used), ω is the weight vector, and t is the iteration number (epoch in MATLAB).

Finally, the steps of training can be summarized as follows: (i) applying the input vectors, (ii) calculating the output of the network and comparing it with the corresponding target vectors, (iii) feeding the difference (error) through the network, and (iv) changing the weights according to the algorithm, which tends to minimize the error. The vectors of the training set are applied sequentially. This process is repeated several times using the entire training set until the error is within acceptable criteria or until the output does not change significantly.

Predictive performance evaluation

To evaluate the predictive power of the model, the MSE (mean square error) and the R² (coefficient of determination) were used to evaluate the overall performance of the model (Eqs. (8)–(9)): (8)${\displaystyle R^2=1-\frac{\sum _{i=1}^n{\left(Y_i-Y_j\right)}^2}{\sum _{i=1}^n{\left(Y_i-\overline{Y_i}\right)}^2}.}$ (9)${\displaystyle MSE=\frac{1}{n}\sum _{i=1}^n{\left(Y_i-{\hat{Y}}_i\right)}^2.}$

In the equation, Y_i is the experimental value of the evaluation model, ${\hat{Y}}_i$ is the corresponding predicted data, ${\bar{Y}}_i$ is the average of the experimental data and n is the number of experimental data.

Data analysis

In this study, correlation plots of the active ingredients of S. miltiorrhiza and rhizosphere soil elements were created using the “corrplot” package (version 0.9; https://cran.r-project.org/web/packages/corrplot/index.html) in R. The Pearson’s correlation test was used to demonstrate the correlation between the active ingredients of S. miltiorrhiza and rhizosphere soil elements. P < 0.05 was considered to indicate a statistically significant difference. BP neural network model was created using the neural fitting tool (nftool) in MATLAB 2019b mathematical software (MathWorks Corporation, America). Scatter and regression plots were created using the “ggplot2” package (version 3.3.5; https://cran.r-project.org/web/packages/ggplot2/index.html) in R. The study used the kappa value to evaluate the multicollinearity of characteristic indicators by the “Kappa” package in R. The map plot was created using the “plyr” package (version 1.8.6; https://cran.r-project.org/web/packages/plyr/index.html), “maptools” package (version 1.1-2; https://cran.r-project.org/web/packages/maptools/index.html), and “ggplot2” package (version 3.3.5; https://cran.r-project.org/web/packages/ggplot2/index.html) in R.

Correlation analysis of active ingredients and rhizosphere soil elements

In this study, the active ingredients and rhizosphere soil elements of S. miltiorrhiza from 25 producing areas in eight provinces of China were determined, and the contents of 6 elements (Mn, Cu, Zn, N, P, K) (Fig. S1) and 10 the kinds of active ingredients (danshensu (D1), protocatechuic aldehyde (D2), caffeic acid (D3), rosmarinic acid (D4), salvianolic acid B (D5), salvianolic acid A (D6), dihydrotanshinone I (D7)), cryptotanshinone (D8), tanshinone I (D9), tanshinone IIA (D10)) (Fig. S2). From the perspective of linear relationship, Fig. 3 shows the correlation between the content of different active ingredients and soil elements. Most active ingredients are related to major and trace elements (r>0). In addition, the relationships between soil elements (Fig. 3) showed that Mn had a significant positive correlation with Zn (p = 0.004<0.01); Cu showed a weak correlation with Zn, K (r < 0.4); N showed a weak correlation with P (r < 0.4); Zn showed a weak correlation with K (r < 0.4). It is suggested that there were synergistic properties in the utilization of elements by S. miltiorrhiza. Then, the study used the kappa value to evaluate the multicollinearity of typical indicators. A k value below 100 was interpreted as low multicollinearity, and a k value exceeding 1000 indicates high multicollinearity (Ma et al., 2021). The kappa value (k = 480.153) was exceeded 100. These results indicated multicollinearity among fundamental soil indicators, and linear regression cannot be performed directly. From the perspective of a nonlinear relationship, Figs. 4 and 5 showed the scatter and regression plots between the active ingredients of S. miltiorrhiza and the soil elements. The results showed that there was not a simple increase or decrease between the active ingredients of S. miltiorrhiza and soil elements, but an inevitable fluctuation would follow, which means that as the content of soil elements increases, the content of active ingredients of S. miltiorrhiza will increase and decrease. These results indicated that the active ingredients of S. miltiorrhiza had a nonlinear relationship with multiple soil elements, and a particular element index cannot be used to express the effective ingredients of S. miltiorrhiza. The artificial neural network methods can learn and create nonlinear and complex relationships and can flexibly solve the complex problems of multiple interacting variables (Bayat et al., 2019). It is one of the best techniques for extracting information from inaccurate and nonlinear data (Caselli et al., 2009). Therefore, in this study, a BP neural network model with soil elements as input values was established to explore the relationship between soil elements and the active ingredients for predicting the active ingredients.

Predictive performance of BP neural network model

In this study, a BP neural network model was constructed using the neural fitting tool (nftool) in MATLAB 2019b mathematical software and trained with soil element parameters as input and the content of plant active ingredients as output. The input variables were as follows: Mn, Cu, Zn, N, P, and K; output variables: water-soluble components (danshensu, protocatechuic aldehyde, caffeic acid, rosmarinic acid, salvianolic acid B, salvianolic acid A), and lipid-soluble components (dihydrotanshinone I, cryptotanshinone, tanshinone I, and tanshinone IIA). Through repeated training, a neural network model with 6-8-2 structure was finally established (Fig. 2). In order to evaluate the predictive ability of the established BP neural network model, the MSE (mean square error) and the coefficient of determination (R²) were used to evaluate the overall performance of the model. Ideally. The closer the MSE value is to zero, the closer the R² value is to 1, which indicates that the average training and testing performance is appropriate. The BP model showed fast training and high simulation accuracy in model testing (Fig. 6A). The relationship between the predicted and measured active ingredients content was favorable, thereby indicating that this model has a basic consistency and a high degree of simulation. A scatter diagram was constructed between the predicted and real values of the inversion model. The coefficient of determination R2 between the predicted and real values of the inversion model was 0.93, the linear regression line between the measured and predicted values were close to 1 (i.e., linear), and the MSE (MSE =0.0203) was low. Therefore, the predictive ability of the model was relatively high and showed strong nonlinear fitting ability, indicating that the soil elemental content could be used for accurate inversion of the active ingredient content.

From the results obtained, the empirical equation based on the Levenberg–Marquardt algorithm for predicting active ingredient content in normalized form is Eq. (10): (10)${\displaystyle \left(D1,D2\right)=\sum _{j=1}^3\left[purelin\left\{\sum _{i=1}^8\sum _{j=1}^{6}logsig\left[\left(Mn{j}_1+Cu{j}_2+Zn{j}_3+N{j}_4+P{j}_5+K{J}_6\right)+b_i\right]\right\}\right.\left.\times L{w}_{i,j}+b_{k_i}\right].}$

The transfer function “purelin” correlated the linear relationship between the input and output variables, while “logsig” calculated the layer’s output from the network input. The variables j ₁, j ₂, j ₃, j ₄, j ₅, and j ₆ were input weights from the input layer to the hidden layer. Also, Lw ₁, Lw ₂ were hidden layer weights from the hidden layer to the output layer. The variables b_i and b_k were biases connected to the hidden layer neurons and output layer neurons, respectively. The values for these weights and biases for the model were in Table 1.

To further evaluate the model’s generalizability, the developed neural network model was tested using 41 new datasets. The new datasets were derived from the results of Zhang’s study (Zhang et al., 2018b). The correlation coefficient, mean square error, was used to evaluate the generalization ability of the model. A scatter diagram (Fig. 6B) was constructed between the predicted and real values of the inversion model. The coefficient of determination R² between the predicted and real values of the inversion model was 0.94, the linear regression line between the measured and predicted values were close to 1:1 (i.e., linear), and the MSE (MSE =0.0164) was low. Therefore, the predictive ability of the model was relatively high. Our results indicate that the BP neural network model based on the content of soil elemental can be a powerful tool for predicting the content of active ingredients of S. miltiorrhiza.