Detection of dynamic protein complexes through Markov Clustering based on Elephant Herd Optimization Approach

The accessibility of a huge amount of protein-protein interaction (PPI) data has allowed to do research on biological networks that reveal the structure of a protein complex, pathways and its cellular organization. A key demand in computational biology is to recognize the modular structure of such biological networks. The detection of protein complexes from the PPI network, is one of the most challenging and significant problems in the post-genomic era. In Bioinformatics, the frequently employed approach for clustering the networks is Markov Clustering (MCL). Many of the researches for protein complex detection were done on the static PPI network, which suffers from a few drawbacks. To resolve this problem, this paper proposes an approach to detect the dynamic protein complexes through Markov Clustering based on Elephant Herd Optimization Approach (DMCL-EHO). Initially, the proposed method divides the PPI network into a set of dynamic subnetworks under various time points by combining the gene expression data and secondly, it employs the clustering analysis on every subnetwork using the MCL along with Elephant Herd Optimization approach. The experimental analysis was employed on different PPI network datasets and the proposed method surpasses various existing approaches in terms of accuracy measures. This paper identifies the common protein complexes that are expressively enriched in gold-standard datasets and also the pathway annotations of the detected protein complexes using the KEGG database.

As PPI network continuously transforms with respect to the environment, time and various phases of the cell cycle, the clustering analysis on static PPI does not emulate these dynamic attributes and it is far from optimal solution. Thus, in recent times, various attempts on the clustering process of dynamic PPI network has been initiated along with the gene expression data to enhance the protein complex detection. Also, many evolutionary approaches were employed for analysing the clustering process of the PPI network such as ant colony optimization ACC-DPC 1 , ACO-MCL 10 , cuckoo search optimization (CSO) 11 , BiCAMWI using genetic algorithm 12 , Soft Regularized-MCL 13 , particle swarm optimization (PSO-MCL) 14 and artificial fish school algorithm (AFA-MCL) 15 . The firefly optimization was employed along with Markov Clustering (F-MCL) on the dynamic PPI network for predicting complexes. The execution time for F-MCL is higher as all the fireflies (proteins) in the population (network) tries to reach the optimal solution (cluster). There are few proteins that are not eligible to come under the cluster and take more iterations to reach the cluster, which may take a long time 16 .
The above-mentioned approaches were effective, but they do not promise a global outcome since they suffer from the effect of unwanted clusters which leads to time consuming. In order to discard the drawbacks of the above-mentioned approaches, a novel approach was proposed to detect the dynamic protein complexes through Markov Clustering based on Elephant Herd Optimization Approach. One of the most important advantages for EHO is that it is the most computationally efficient and has less time consuming compared to F-MCL and other approaches. This is because the unwanted noisy data (unclustered proteins) will be removed from the clan separating operation of EHO approach. The remaining sections of this paper is ordered as follows: Section 2 discusses briefly about the methodology of the proposed approach. Section 3 illustrates the experimental results with various performance measures, Section 4 deliberates about the implementation and discussion of the proposed method in detail and finally Section 5 concludes the paper and recommends for the future enrichments.

Methods
For detecting the protein complexes, initially, the proposed method divides a static PPI network into a sequence of subnetworks below diverse time points by combining gene expression data to form dynamic model. In order to build a dynamic model, the static PPI network is integrated with gene expression data, which declare the level of gene expression, as well as protein expression. As a protein does not always becomes active at a cell cycle, it is assumed that a protein was active at the time points with its highest expression level 17 . The expression level of a protein will be increased before its expression and will be decreased once the protein has completed its function, and the time points are identified with its expression level, which are higher than a threshold.
Given is a static PPI network P P = (P ver , P Edg ), where P ver , is a set of proteins and P Edg , is a set of interactions between these proteins. In gene expression data, there is a series of T time stamps coming with |P ver | × (T * TR) matrix M, where TR is the number of repetitions of the time series. Each element M(P ver , j) of this matrix represents the level of gene expression.
The three-sigma principle is employed to determine if a gene is expressed in a single stamp. For each gene P ver , the gene expression is defined as given in the following Eqs (1)(2)(3)(4)(5) = ∑ + × − = Ev(P ) M(P , i T (tr 1)) TR (1) i ver tr 1 TR ver where Ev i (P ver ) is the mean of the expression value of gene P ver at timestamp i, UE(P ver ) is the mean of its expression values over times ranging from 1 to T, σ(P ver ) is the standard deviation of its expression values, Fl(P ver ) is used to show fluctuation of the expression curve of gene P ver . Suppose that the gene expression data is governed by a normal distribution, then S 1 (P ver ) and S 2 (P ver ) are the associated mean and three-sigma value, that is S 1 (P ver ) = UE(P ver ) and S 2 (P ver ) = UE(P ver ) + 3σ(P ver ). In virtue of three-sigma principle, the probability that a value greater than S 2 (P ver ) is not an active point is less than 0.1%. AT(P ver ) is the active threshold of gene P ver . Consider the gene (P ver ) at timestamp i. If Ev i (P ver ) > AT(P ver ), then the gene P ver is expressed and the gene product exists 16 .
In the clustering procedure of every subnetwork, the proposed method starts with constructing the initial protein clusters depending on the protein complexes attained at the prior time point. The initial clusters constructed in the first generation have three steps The procedure for constructing initial clusters has three steps: seed node selection, attachment nodes addition and finally refining 1 . To clearly demonstrate the three steps, a 1. Selecting seed nodes: This step first computes the clustering coefficient of every node. Then it selects the nodes whose clustering coefficients are greater than a given threshold λ c as seed nodes, and puts them into the set of seed nodes at the current time point t, denoted by S t . The seed nodes are considered as the candidate clustering centers and represent different clusters of protein complexes. The clustering coefficient of any node i is defined in Eq. (6): where Neigh(i) = {j є P ver t |(i · j) є P Edg t } represents the neighbor nodes of node i, and |Neigh(i)| is the number of neighbor nodes of node i, n i t is the number of links between neighbour nodes of i at the time point t. 2. Attachment nodes addition: For any seed node i (i є S t ) of current time point t, if it is also the seed node of previous time point (t − 1), then the nodes which are in the cluster i at the previous time point (t − 1) and also exists in the subnetwork P p t at the current time point t are put into the cluster i of current time point t. In this way, initial clusters are built. However, some clusters may be too sparse since that not all proteins of previous time point (t − 1) exist at the current time point t. Thus, a refining step is needed to be carried out on the initial clusters. 3. Refining: For any initial cluster of protein complex c i t at the current time point t, if its density is smaller than a given threshold λ d all the nodes in c i t are sorted in a descending order according to their clustering coefficients, and the node with the smallest clustering coefficient is iteratively removed until the density of cluster c i i t is not smaller than the given threshold λ d . The density of a protein complex c i t is computed by Eq. (7): where n i and l i are number of nodes and edges in cluster c i t respectively 1 . Now, the clustering analysis of the remaining generations is employed by utilizing the Markov Clustering technique along with the EHO algorithm on every subnetwork. The matrix is constructed that depicts the probabilities of transition of a Markov Chain (random walk) based on the graph. The MCL procedure comprises of two activities such as expansion and inflation, which was applied to the matrix that was constructed. The construction of matrix M at for a graph description and the process of Markov clustering method is briefly described 18 .
Let P p = (P ver , P Edg ), where P ver , is a set of proteins and P Edg , is a set of interactions between these proteins. Denote a node in P ver by p vi and an edge between p vi and p vj in P Edg by (p vi , p vj ), in which i and j are the indexes of the corresponding nodes 16 . W(p vi , p vj ) is the weight of edge (p vi , p vj ), which represents the confidence level of the interaction in a weighted PPI networks. Adj is the adjacency matrix of a weighted graph given as Eq. (8) A canonical flow matrix M at is an n × n (n = |P ver |) matrix that shows the probabilities of transition of a random walk defined on the graph. M at (i, j) represents the probability of a transition from node p vi to p vj . The transition probability from p vi to p vj is referred to as the stochastic flow from p vi to p vj . All the elements in each column of M at will sum up to 1 and the matrix is expressed as given in Eq. (9) at k 1 n

= ∑ =
The three crucial parameters of MCL are inflation constant (ic), balance (b) and penalty proportion (P p ), where ic defines the size of each cluster, b defines the user-specific balance constant that is employed for penalizing higher-propensity neighbours and P p defines the penalty ratio of the protein nodes, which is also user-specified 16 . The clustering process using EHO algorithm is briefly explained here for clustering protein complexes. The overall flowchart of the proposed method is shown in Fig. 1. elephant herd optimization. One of the contemporary swarm intelligence technique is the elephant herd optimization which was projected in 2016 19 . This algorithm was stimulated by the herding characteristics of elephants. In general, elephants are social mammals with the composite social group comprising of numerous clans under the guidance of a matriarch. A clan comprises of one or more female elephant with their calves. Female desires to live in domestic clusters while male elephants prefer to live alone and they will exit from the clan when they grow with each generation 20 . The characteristics of the clans signifies exploitation and leaving elephants signifies the exploration of the population. (2019) 9:11106 | https://doi.org/10.1038/s41598-019-47468-y www.nature.com/scientificreports www.nature.com/scientificreports/ The characteristics of an elephant are measured using two main operators, namely clan updating and clan separating operators that are used for producing better clustering of proteins. Here, the elephant population is referred to as the static PPI network, each clan is referred to as the dynamic PPI subnetwork, and the elephants inside each clan is represented as proteins.
Clan updating operator. The static PPI is initially separated into k dynamic PPI. Each dynamic PPI is headed by the individual protein, which represents the best solution of the dynamic PPI. In each generation, protein e of dynamic PPI cl i moves towards the p best cl , i which has the best fitness in dynamic PPI cl i . The fitness of the dynamic PPI is computed by employing the accuracy values of the protein complex. For new protein e in dynamic PPI cl i , the position is updated by following Eq. (10).
new cl e cl e b est cl c l e , , , , , is the new position of protein e in dynamic PPI cl i and p cl e , i denotes the position in previous generation. p best cl , i signifies dynamic PPI cl i which has the best fitness, α is the scale factor that determines the influence of best fitness and rand is the random variable employed to enhance the diversity of the populations and defined in the range (0, 1).
The movement of a protein e for best fitness can be updated using Eq. (11).
p p (11) best cl e c enter cl , , , where β belongs to (0, 1) which is a scale to regulate the effect of p center cl is the centre of dynamic PPI cl i and for the di th dimension it can be computed using the Eq. (12).
The centre of the dynamic PPI cl i is computed through DI computations using Eq. (12). The pseudocode for the dynamic PPI updating operator is depicted in Algorithm 1.
Clan separating operator. To enhance the search capacity of the proposed method, the unclustered proteins and clusters with the lowest fitness will exit in every generation as given in Eq (13) 19 .
where p max and p min are the upper and lower bound of the single protein.
is the protein or complex with the lowest fitness. The rand is the random variable that has stochastic and uniform distribution in the range (0, 1). The pseudocode for the clan separating operator is given in Algorithm 2.
Depending on the clan updating and separating operator, the module of the proposed algorithm is framed as given in Algorithm 3.  www.nature.com/scientificreports www.nature.com/scientificreports/ The relationship between the DMCL-EHO and the protein complex is given in the Table 1.

experimental Results
Datasets. In this experiment, the datasets which consists of interactions for both Saccharomyces cerevisiae and Homo Sapiens are DIP 21 , BioGRID 22 and STRING 23 . The benchmark PPI datasets employed only for Saccharomyces cerevisiae are Gavin2 and Gavin6 24 , Krogan-core and Krogan-extended 25 , Collins 26 , and WI-PHI 27 . The Gavin + Krogan dataset was generated by merging Gavin and Krogan Core datasets. The PPI datasets employed only for Homo Sapiens are HPRD 28 , HPID 29 and PIPs 30 . Table 2 shows the list of datasets used in this experiment.
The gene expression data used in this study for Saccharomyces cerevisiae (GSE3431) 31 and Homo Sapiens (GSE3933) 32 are taken from the GEO database.
The predicted complexes are compared to gold standard benchmark databases such as CYC2008 33 , MIPS 34 , SGD 35 for Saccharomyces cerevisiae organism and PCDq 36 benchmark dataset for Homo sapiens organism. The percentage of overlapping interactions among the datasets in Gavin2 is 32%, Gavin6 is 53%, Krogan-core is 46%, Collins is 56%, HPRD is 23%, PIPs is 57%, DIP is 2%, BioGRID is 55% and STRING is 47% 37,38 . performance measures. To evaluate and compare the clustering results of predicted protein complexes, the generated complexes were compared and matched with the gold standard benchmark protein complexes. Assume P r (V Pr ,E Pr ) and B e (V Be ,E Be ) be the set of vertices (proteins) and edges (interaction) of a predicted protein complex and benchmark protein complexes.
Complex similarity score (CSS). CSS is defined as the closeness of two protein complexes namely predicted (P r ) and benchmark (B e ) protein complexes and they are computed based on Eq. 14.  where V pr and V Be denotes the set of proteins in predicted and benchmark protein complexes. If CSS(P r , B e ) is equal to 0, it denotes that the predicted and benchmark protein complexes do not have any common protein complexes. On the contradictory, if CSS(P r , B e ) is equal to 1, then the predicted complex P r (V Pr , E Pr ) has the same equal nodes as the benchmark complex B e (V Be , E Be ). Here, if CSS(P r , B e ) > 0.2, it is considered as the predicted and benchmark protein complexes are identical 39 .
Recall. Recall is defined as the accuracy of benchmark protein complexes that are identical to the predicted complexes. If the recall value is high, it indicates that the predicted complex has a good number of coverage of the proteins in the gold standard complexes. The recall of the protein complexes is computed based on Eq. (16).
where N Pc is denoted as the number of predicted complexes which match at least one recognized benchmark complex, N Bc is denoted as the number of recognised benchmark complexes which match at least one predicted complex, Predicted set is denoted as the set of complexes predicted by the proposed approach and Known set is denoted as the set of recognised gold standard benchmark protein complexes.
Coverage ratio (CR). CR is defined as the fraction of proteins in benchmark complex V Be found in predicted complex V pr and they are computed based on Eq. (17).

Elephant
The temporary proteins in dynamic subnetwork.
where V Be is denoted as the set of proteins in benchmark protein complexes. T i, j is denoted as the common number of proteins between V pr and V Be .
F-Measure. F-Measure is defined as the harmonic mean, i.e., a rational mixture of both precision and recall and it is computed based on Eq. (18).     Number of Clusters. The number of clusters is defined as the total quantity of clusters formed from the PPI network after the clustering process has been completed.
The performance measures such as coverage ratio, the number of clusters, precision, recall, f-measure and accuracy of the proposed method for Saccharomyces cerevisiae are compared with various datasets and existing algorithms against CYC2008 benchmark database and the graphical representation of the comparison is depicted in Figs 2-7. Also, the performance measures such as coverage ratio, the number of clusters, precision, recall,     www.nature.com/scientificreports www.nature.com/scientificreports/ From Figs 2 and 8, it is inferred that the number of clusters in the proposed method is less when compared to FOCA, AFA-MCL and ACO-MCL as they try to get solution from all the proteins in the network. These methods will not discard the undesirable proteins which may result in false positives. But in the proposed method, the clusters which has less than three proteins are discarded. Hence the precision, recall, F-Measure and accuracy are high for the proposed method.
From Figs 3 and 8, it is observed that the proposed method has more coverage ratio than the existing methods since it employs the iterated clustering approach. This enhances the coverage of proteins in the network as the proteins in the benchmark complexes are highly found in the predicted complexes. From Figs 4-7 and 9 it is observed that the precision, recall, F-Measure and accuracy shows fluctuations for PSO-MCL, ACO-MCL, AFA-MCL, F-MCl, FOCA and EHO-MCL. The mean of these measures for all the datasets shows that the proposed method performs better than the existing methods because it has employed the dynamic PPI along with EHO.

implementation and Discussion
The computational issue of attaining a solution with a high accuracy solution for protein complex detection from dynamic PPI is still a challenging task. In this paper, the elephant herd optimization algorithm along with Markov clustering technique is combined to solve the protein complex detection problem. The proposed method provides an enhancement of the results compared to all the other popular existing methods. This work was executed on 2.00 GHz Intel i3 with 8GB of memory running on Windows 10.
The number of clusters is small in an average when compared to other existing methods, due to the deletion of proteins without interactions. Here, the minimum number of proteins inside a cluster should be three or more and that are considered as a protein complex. The protein cluster with less than three proteins are removed. The proposed method was evaluated based on the removal of noise, insertion and deletion of random protein interactions, large PPI network, namely WI-PHI, various parameter analysis, statistical significance and finally with biological significance. evaluation by noise removal. The PPI networks are obtained from high-throughput experiments, the large coverage of the PPI network comprises of noise in the format of false positive interactions and redundant data. The main challenge of clustering these PPI networks is present in the PPI networks itself. In this method, after the clustering process is accomplished, the proteins that do not present in any of the clusters is also considered as a noise. These solitary proteins that do not interact with any other proteins will not provide any valuable information. The minimum number of proteins inside the cluster is set to be three in this work. Thus, the isolated proteins and clusters with below three proteins are considered as a noise and they are removed by the clan separating operator by the elephant herd optimization method. Many evolutionary approaches are inheriting the undesirable proteins from one generation to another which may lead to loss of accuracy, but EHO approach will discard the undesirable proteins from the population in the clan separating operator that leads to the optimal solution. The comparison of EHO with other existing methods is depicted in  evaluation by adding and removing random protein interactions. The testing of the proposed method is accomplished by inserting and deleting the random interactions of the PPI network to evaluate its performance. The noise can also be any missing information (false negatives) or added noise (false positives) in the PPI network. The DIP dataset is used for evaluation of adding and removing random interactions. The missing information of PPI network is processed by removing the proportion of edges randomly (0%, 20%, 40%, 60%,   ) and the false positive information of PPI network is processed by adding the proportion of edges randomly (0%, 20%, 40%, 60%, 80%, 100%). The performance of the proposed method by adding and removing the random interactions are depicted in Figs 10 and 11. From the Figs 10 and 11, it is observed that even though the random insertions and deletion of the protein interactions are employed on the dataset, the proposed method performs better than other existing approaches.    www.nature.com/scientificreports www.nature.com/scientificreports/ evaluation by parameter analysis. Generally, every metaheuristic approach is based on certain stochastic dissemination. Hence, diverse runs will produce various diverse results. This work implements 500 independent runs in order to score optimal solution. In general, 20 numbers of clans were employed as per literature. The execution process will be terminated, if the best result generated in each iteration remains interchangeable for 100 successive iterations or the maximum number of generations is attained. The assignment of parameter values was adjusted based on the experimental results. It was identified that the parameters of the proposed method that has values of α = 0.5 and β = 0.1 produced better solution among different values and hence were allocated. It was observed that the optimal solution was identified after 315 th generation. For all the performance measures, there were fluctuations during the first 10 runs of the experiment and in the future runs reliability was observed. Figures 2-9 shows the average outcome of performance measures for the above parameter values of the proposed method. Table 3 shows the various parameter values for the proposed approach and the other existing approaches of protein complex detection. evaluation by statistical significance. The proposed method was also assessed by utilizing non-parametric test such as, Wilcoxon Matched-Pair Signed-Rank Test among each pair of approaches that produces the statistical consequence. The discrepancy between the F-Measure and Accuracy for every entry in Figs 6 and 7 was tested based on the confidence level of 1% (p-value < 0.01). The p-value less than 0.01 are assumed as highly significant and the values greater than 0.01 are assumed as insignificant values. The scores of F-Measure and Accuracy is alone considered as they are computed based on precision and recall. The Statistical Significance  of the proposed and existing approaches based on F-Measure and Accuracy is depicted in Table 4. The scores of upper right positions of the table are attained from F-Measure of proposed and various existing algorithms based on DIP dataset against CYC2008 benchmark database. The scores of lower left positions of the table are attained from Accuracy of proposed and various existing algorithms based on DIP dataset against CYC2008 benchmark database. From Table 4, it is shown that the proposed method is statistically significant in nature compared to all the existing methods.
Evaluation by biological significance. Many of the existing methods solve the protein complex detection problem based on the topological similarity. But to obtain some useful biological information, the computational methods should be biologically significant in nature. This proposed method is evaluated in the biological significance test. The gold standard benchmark databases are manually annotated based on the information from biologically experimental analysis. Thus, the detected protein complexes obtained from the proposed method is compared and matched with the benchmark databases. Few benchmark databases such as CYC2008, MIPS, SGD databases for Saccharomyces cerevisiae and the PCDq database for Homo sapiens are employed for assessing the proposed method. Table 5 displays the common protein complexes between the CYC2008 benchmark database and the proposed method for DIP and Krogan-extended. Also, the common biological process, molecular function and the cellular component of the obtained protein complexes are displayed. Correspondingly, the common pathway annotations of the predicted protein complexes are obtained from the KEGG database are displayed.
The predicted complex gene ontology and KEGG pathway enrichment analysis were predicted by using the DAVID gene function classification online tool. The overall predicted complex enrichment score and the respective gene ontology elements and KEGG pathway enrichment scores are displayed in Table 5. The pictorial representation of the common RNA Polymerase KEGG Pathway of the predicted protein Complex on Krogan-extended dataset and common Oxidative Phosphorylation KEGG Pathway of the predicted protein complex on DIP dataset is exhibited in the supplementary information. The RNA polymerase is essential for nucleolar assembly and for high polymerase loading rate. Oxidative phosphorylation is the metabolic pathway in which cells use enzymes to oxidize nutrients, thereby releasing energy which is used to produce adenosine triphosphate (ATP) [40][41][42] . The pictorial representation of the Top 5 common protein complexes, gene ontology functions and KEGG pathways of the predicted complexes of proposed method is given as Venn diagram in Figs 12 and 13. execution time. Besides the accuracy, the time required to detect the dynamic protein complexes is also an important factor. Processing the various benchmark datasets with various numbers of proteins and different interactions requires more time complexity due to stochastic optimization methods. Subsequently, not all methods were available under the same platform, the execution of many of the approaches were done on virtual machines, which prohibited us from accomplishing an exact comparison of their relative execution times. Thus, here the average execution time of SR-MCL, ACO-MCL, PSO-MCL, AFA-MCL, F-MCL, FOCA AND EHO-MCL is displayed in Fig. 14.
From Fig. 14, it is observed that in this research, the proposed algorithm has less execution time when compared to other algorithms, due to the clan separating operator of EHO approach. It is inferred that the proposed EHO-MCL is efficient for detecting dynamic protein complexes. In future, the EHO-MCL can be further optimized in multicore CPU.

conclusion
The Protein Complex detection is an exposed problem for scientists. The solution for the complex problem should be recurrently improved as they are important in the analysis of the biological process. The volume of PPI networks has also been increased due to high-throughput experiments, the lack of accurate computational model for protein complex detection exists. Many of the existing researches were employed on the static PPI data that do not provide accurate biological results. Thus, in this proposed method initially, the static PPI data is converted into dynamic PPI data by integrating the gene expression data. Later, every dynamic subnetwork was clustered based on the popular clustering technique MCL along with the elephant herd optimization method for exploring and exploiting the better solution. The proposed method was employed on various 11 widespread datasets and the predicted complexes were compared with 4 different benchmark databases. Also, the proposed method was evaluated based on noise removal, insertion and deletion of random protein interactions, using the large PPI dataset, various parameter analyses, statistical significance and biological significance. On every evaluation phase, the proposed method was outperforming all other existing approaches and identified the common protein complexes, Gene Ontology functions and KEGG pathways of predicted protein complexes. As a future work, additional information on the unknown protein complexes predicted by the proposed method is to be addressed with the help of biological experts. The proposed method can also be applied and analyzed on weighted PPI networks. Also, various other diseased databases can be used to experiment.