A matrix M is formed by n-dimensional feature vectors x₁, …, x_m in space, where each row x_i (x_i1, x_i2, …, x_in) is regarded as a data point, and we assume that we can gather the attribute parts of other data points to get the current one. Adjacent data points are usually viewed as possessing similar properties. Hence, the neighbors can be selected as the contributing force for the reconstruction of the point, whereas other irrelevant data points participate in the calculation but are assigned a weight of 0. For x_i, the common calculation method for selecting neighbors is Euclidean distance. Considering that various features, such as GC content and MFE, characterize the components of the data point in these dimensions and that similar data points are assumed to eventually form multidimensional vectors with similar directions, the Cosine distance is then chosen to select neighbors N₁(x_i) (a total of n₁), and then the Euclidean distance is chosen to select nearer neighbors N₂(x_i) (a total of n₂) on the basis of the Cosine distance to select more exact neighbors. The latter is a subset of the former, ensuring that neighbors are nearer in both direction and position. The order of the two can be swapped since the final selected neighborhood is the same. The percentages of neighbors are expressed by K₁ and K₂ where K₁ = n₁/m and K₂ = n₂/n₁.

To minimize the reconstruction error for m data points, we propose the objective function:

where M is an m × n matrix in the feature space, and both C₁ and C₂ are indicator matrices that separately indicate whether they are neighbors on the basis of the Cosine distance and nearer neighbors on the basis of the Euclidean distance. G and W are both weight matrices. Here, μ is a weight parameter, and e is an m × 1 column vector with all elements being 1. ‖•‖_F is to obtain the Frobenius norm of a matrix. ‖•‖₂ is to obtain the 2-norm of a vector. The first term of the objective function is to get the optimal weight matrix to minimize the reconstruction error of all data points, whereas the second term is to reduce overfitting during the reconstruction. For the weight matrices G and W, the elements are non-negative.

Given that the first neighbors are required to find the second ones, the objective function must be decomposed:

The Lagrange multiplier method is used to solve Equation (2), which then has the following form:

Then

where λ₁ and λ₂ are Lagrange factors.

According to the Karush–Kuhn–Tucker condition (Kjeldsen, 2000), the following conditions must be satisfied to determine the optimal value:

The partial derivative of L is determined with respect to G:

Then

If λ₁^T = 0, then there will be

In that case,

If λ₁^T ≠ 0, then there will be G = 0. Thus, G_ij = 0.

Given the relevance of data point-based reconfiguration that G_ij ≠ 0 when x_j∈N₁(x_i), the solution is

Thus far, the iterative form with unknown parameters has been obtained. We inscribe its equivalent form for Equation (2) to obtain the following parameter:

where Grⁱ is the gram matrix. If x_j∈N₁(x_i) and x_k∈N₁(x_i), then Gr_j, _k = (x_i – x_j)^T(x_i – x_k). Otherwise, Gr_{j, k} = 0. Solving Equation (9) using the Lagrange multiplier method yields

By taking the partial derivatives of g and λ, λ_i can be obtained:

where the reconstruction error is close to none, i.e., 1/2ϑiTGriϑi≈0. According to the Lagrange multiplier method, e^Tϑ_i−1 = 0. Thus, λ_i = μ. If λ = μ × e, G can be represented as

The similarity matrix based on the Cosine distance is then acquired through iteration until convergence. Similarly, W is obtained with respect to G:

The final similarity matrix can be acquired by iterating until convergence or the maximum rounds.

Given a set of l lncRNAs, l₁,…, l_i,…, l_l, and a set of m miRNAs, m₁,…, m_j,…, m_m, whose interaction is represented by a matrix Y of l × m, if there exists an interaction between lncRNA l_i and miRNA m_j, then Y_ij = 1; otherwise, Y_ij = 0.

Four features of lncRNAs and miRNAs are extracted (110 features in total). The frequency of l lncRNAs with 4^k long contiguous subsequences is calculated. Herein, we assumed k = 1, 2, 3 to obtain the lncRNA-related feature vector. The sequence similarity between pairs of l lncRNAs is calculated to yield the similarity matrix of l × l, which is denoted as S_lncS. In the same manner, the similarity matrix of m miRNAs is denoted as S_miS of m × m. The interaction profiles of lncRNAs and miRNAs are derived from the interaction matrix. For lncRNA l_i, the interaction profile indicates whether it interacts with each miRNA, matching the i-th row of Y, i.e., Y(i,:). Similarly, for miRNA m_j, it matches the j-th row of Y, i.e., Y(:, j). The similarity between two interaction profiles of l lncRNAs is calculated as the matrix S_lncP of l × l, whereas the similarity between two interaction profiles of m miRNAs is computed as the matrix S_miP of m × m.

Label propagation (Kato et al., 2009) assigns labels to previously unlabeled data points. During label assignment, the labels of labeled data points are propagated to unlabeled data points. The core idea of the label propagation algorithm is that similar nodes should have similar labels. It involves two stages, namely, calculating the similarity matrix and propagating the labels.

The edge from node i to node j represents the similarity of these nodes. All edge weights constitute a weight matrix, where the higher the similarity the larger the weight. Herein, ILNS is adopted to construct the similarity matrix and calculate the Cosine-distance neighbors and the Euclidean-distance neighbors of each node until convergence. Nearer neighbors of each data point are fixed to a certain proportion, and the weights of others are 0. The weight matrix is actually a sparse matrix. The labels are propagated through the edges between the nodes. The larger the weight of the edge, the more similar the nodes are to each other and the easier to propagate the labels (Zang and Zhang, 2012; Zhang et al., 2016). For m data points x₁,…, x_m, an m × m probability transfer matrix P is defined as

where P_ij represents the transferring probability and w_ij is the weight. The propagation involves three steps. First, a unique label is allocated to each node, i.e., label one for node one and label i for node i, where the labels are different from each other. Second, for node j, all nodes are traversed to discover their neighboring nodes and obtain their labels to obtain the label with the most occurrences. If more than one label satisfies the largest number of occurrences, then one is randomly selected to replace the current label. Finally, if the label of node j no longer changes after this round of relabeling or the pre-set number of rounds is reached, then the iteration is stopped. Otherwise; step 2 is repeated.

The sequence and interaction profiles of lncRNAs and miRNAs are captured to develop our model. The workflow of MINLP is shown in Figure 2. The specific steps are as follows:

Step 1: the sequence feature similarity S_lncS and S_miS are calculated using the ILNS algorithm.

Step 2: the interaction profiles of all lncRNAs and all miRNAs are exported according to the interaction matrix of lncRNAs and miRNAs.

Step 3: the interaction profile similarity S_lncP and S_miP are calculated using the ILNS algorithm.

Step 4: S_lncS and S_lncP are combined to obtain lncRNA integrated similarity, and S_miS and S_miP are combined to obtain miRNA integrated similarity.

Step 5: for the two integration similarities, the linear neighborhood propagation method is used to generate the prediction of lncRNA and the prediction of miRNA.

Step 6: the weighted sum of the two prediction matrices is calculated, and the final interaction prediction matrix is determined.

The effects of the different parameters are visualized in Figure 3. First, α and β are fixed. Theoretically, the neighbors are set twice to find a more accurate batch faster. However, serious analysis reveals that a change in K₂ is logical within a certain range of K₁; otherwise, even when K₂ is 100%, it will not be very helpful. Let K₂ = 1.0, as K₁ changes from 0.1 to 0.9, we find that K₂ = 0.9 is the optimal value, as shown in Figure 3A. AUC initially decreases and then increases with K₁ and reaches its lowest point at 0.4. The AUC values are always above 0.975. K₁ = 0.9 is then fixed, and K₂ is changed. As K₂ becomes larger, AUC tends to increase globally and reaches the maximum at 0.9 (Figure 3B). The most pronounced increase is evident from 0.7 to 0.8, clearly demonstrating that 0.8 is a cut-off point. Finally, the impact of α and β is investigated. The contour and concentration plots in Figures 3C,D, respectively, demonstrate the variations in AUC with α and β. The gradient from blue to yellow is set to indicate the increase in AUC values. Herein, the scenario with β = 0 is dropped to detect subtle variations in the other cases owing to our prior observation that the AUC value is as low as 0.4 in the case of β = 0, thereby forming a cliff-like change from others. Both plots show that the best results are achieved at α of 0.2–0.8 and β of 0.15–0.9, which are extremely close to 0.98. The yellowest areas appear locally as a result of drawing tool error, but it does not affect the fact that the results are roughly consistent. Figures 3E,F present the grid plots of AUC with respect to variations in α and β, respectively. Figure 3E is the case with β = 0 removed from the corresponding plots in 3(C) and 3(D). Figure 3F includes all cases and validates the previous interpretation that the inclusion makes distinguishing the changes from the others difficult. The models of all ranges show that the best parameters are K₁ = 0.9, K₂ = 1.0, α = 0.7, and β = 0.35 when AUC at this point is 0.9797. This tells that the performance of our model is attributed to Cosine distance, determining a more accurate neighborhood which is preferable to applying Euclidean distance only.

To demonstrate the superiority of our model, we calculate a similarity network with a single information source and build linear neighborhood propagation models to compare with MILNP. The results are presented in Table 3. The optimal values obtained are set in bold typeface. We construct these models based on the optimal combination of parameters. Except for the difference in information sources, all the other processes are guaranteed to be the same. The model using only sequence similarity with ILNS as the core algorithm is named MILNP-I. The model using only interaction profile similarity with ILNS as the core algorithm is named MILNP-II. Overall, both of them are inferior to MILNP with integrated information. MILNP-II is generally very close to our model, indicating that the interaction information is a key contributor throughout the prediction. Thus, with specific optimization, MILNP-II could be credible for predicting isolated lncRNAs or miRNAs. The performance of MILNP-I is not as good as that of MILNP-II, but its AUC is barely satisfactory. We also tried another method for calculating similarity, LNS (Zhang et al., 2018b), to obtain new models, named the SLNPM-series, as shown in Table 3. With the information source controlled, a comparison of SLNPM-II and MILNP-II reveals that ILNS is more accurate. All these results clearly validate the superiority of MILNP in terms of information integration and similarity calculation.

As far as we know, very few studies have investigated lncRNA–miRNA interactions in plants. We select Pmlipred (Kang et al., 2020) and CIRNN (Zhang et al., 2020) as reference methods. Both methods predict plant lncRNA–miRNA interactions. PmliPred (Kang et al., 2020) builds a prediction model by using a machine learning approach combined with a deep learning approach, and the final prediction results are made from fuzzy decisions of the two components. We use the publicly available source code on GitHub to implement Pmlipred (Kang et al., 2020). CIRNN (Zhang et al., 2020) builds integrated deep learning models with both a CNN and an IndRNN, where the former is used to automatically extract gene sequence functional features, whereas the latter is utilized to obtain sequence feature representations and dependencies. We replicate it in the detailed description of CIRNN (Zhang et al., 2020). The results of the comparison are summarized in Table 4. The optimal values obtained are set in bold typeface. The performance of the two methods is clearly good, but not as good as that of our model. In particular, the AUC and ACC values have a relatively large gap. We notice that both methods stretch to deep learning. Thus, we implement similar models on their basis to measure their effectiveness. We construct a CNN–Gate Recurrent Unit (GRU) combinatorial model. Sequential features are extracted from the original data by CNN and compressed into a one-dimensional vector in the flattened layer to input into GRU. This process is well-suited for processing sequential information. We consider GRU instead of others because GRU has fewer parameters and reduces overfitting. We first use a three-layer CNN and a single-layer GRU mixed with RF to obtain CNNRF1. We then add another layer of CNN on top of that to obtain CNNRF2. Gladly, although our MILNP is simple and built by deriving mathematical formulas via a top–down approach and layer by layer, the results are satisfactory. For the baseline SLNSM (Zhang et al., 2018b) that is originally created for animal prediction, we run our dataset and observe that our method is slightly better. That is because we focus on improving the computational procedure for linear neighborhood similarity by adding the spatial direction restriction. The optimal parameter combination shows that such performance is attributed to Cosine distance, determining a more accurate neighborhood, which is preferable to the approach of the baseline.

Top-rank prediction is an important way of visualizing the performance of the models. We examine the top-rank predictions from 200 to 2,000 and identify the percentage of interactions that are truly correct. As shown in Figure 4, an average of 186 positive interactions per prediction is reached in the top 200 predictions, whereas 1,380 real interactions are determined in the top 2,000 predictions. The results demonstrate the good performance of our model.

Furthermore, we predict the interactions of isolated lncRNAs and miRNAs with MILNP. For the isolated lncRNA or miRNA, only sequence-dependent information can be used. In separate cases gma-miR395a and lcl| Gmax_Glyma.18G279100.1 are taken as examples. We validate the prediction of the selected miRNA and lncRNA with respect to RNAhybrid2.1.2 (Rehmsmeier, 2004). All predictions are sorted in descending order of probability. For miRNA gma-miR395a, 4 of the top 10 are correctly predicted, as shown in Table 5, thereby confirming the predictive power of MILNP. However, the list reveals that the fifth and eighth detected lncRNAs in the prediction of miRNA “gma-miR395a” belong to Medicage truncatula, as evidenced by their nomenclature. The situation is worth contemplating. On the one hand, this indicates that the selected samples may affect the performance of MILNP. On the other hand, it inspires us to further explore cross-species linkages and assume that the remaining uncertified interactions are possible. Likewise, we make predictions for lncRNA lcl| Gmax_Glyma.18G279100.1. To our surprise, five of the results happen to be identified by the tool as having interactions, a result that is very encouraging. We conjecture that the remaining ones predicted by our model can be possible. As demonstrated by the results of the comparison of the two sets of predictions, the fact that sample selection has a great influence on the prediction results should not be ignored. The association between the selected sample and other samples also affects the results. For those that have a similarity with many samples, the prediction results may be more accurate. If a sample has little similarity with other samples, then predicting its potential interactions will be difficult.

The different results also prompt us to reconsider the results of previous experiments. We find that, although MILNP achieves good AUC and ACC, its PRE is relatively low, which may be attributed to the model itself and the distribution of the dataset. Some very similar samples may have confused the model and that induces it to arrive at a wrong judgment. Nevertheless, it could also be a new revelation that suggests these possible associations. This assumption warrants further biological laboratory validation.