Probing molecular docking problem by an improved quantum-behaved particle swarm optimization algorithm

Abstract

The main objective of molecular docking is to find a model of interaction between a protein and ligand with a minimum binding energy. This process is driven by intricate algorithms and scoring functions. This paper mainly concentrates on the search algorithm used for solving the docking problem. Here, a new approach is proposed for the molecular docking problem that utilizes a hybrid algorithm that combines an improved quantum-behaved particle swarm optimization algorithm (QPSO) and the Solis and Wets algorithm. The improved QPSO algorithm that is based on individual particle evolutionary processes is known as individual particle evolutionary particle swarm optimization (IEQPSO). The IEQPSO algorithm was tested and compared with particle swarm optimization, QPSO, and its variants with a suite of benchmark functions. The results indicated the superiority of the proposed approach according to benchmark test functions. Then, the hybrid algorithm based on the IEQPSO algorithm was used for optimizing the energy function of the molecular docking problem and was compared with the classical Lamarckian genetic algorithm used by molecular docking software. Molecular docking and molecular dynamics simulation experiments revealed the effectiveness and feasibility of the proposed algorithm in solving the molecular docking problem.

Keywords

Quantum-behaved particle swarm optimization individual particle evolutionary process swarm intelligence molecular docking molecular dynamics simulation

Introduction

Among many evolutionary algorithms, the quantum-behaved particle swarm optimization (QPSO) algorithm¹ is more efficient than other conventional algorithms. It is based on quantum mechanics and trajectory analysis of particle swarm optimization (PSO).² The main difference between QPSO and PSO is that the QPSO algorithm requires no velocity vectors for particles, and it has fewer parameters to adjust. The QPSO algorithm has been shown to successfully solve a wide range of optimization problems. Moreover, many efficient strategies have also been proposed to improve algorithm performance.^3–4 However, in the QPSO algorithm, the mean best position is simply the average of the personal best position of all particles, which means that each particle is considered equal and exerts the same influence on the mean best position. In fact, it is not appropriate to consider each member equal. In this article, to improve the performance of QPSO in highly dimensional complex problems, we proposed an improved QPSO based on an individual particle evolutionary process (individual particle evolutionary quantum-behaved particle swarm optimization (IEQPSO)) algorithm to improve the global search capability and the convergence performance of the optimization process. The benchmark functions were used to test the performance of the IEQPSO algorithm. The ultimate aim of the optimization algorithm is to solve real-world problems. Here, we use the IEQPSO algorithm to optimize the energy function used for the molecular docking problem.

Molecular docking is an important tool in drug discovery, and many docking methods have been developed and widely used.⁵ It can be used for predicting the correct conformation between a protein and ligand, for identifying novel lead compounds, for preliminary lead optimization,⁶ and for investigating the mechanisms of action of biologically active compounds.⁷ The quality of molecular docking depends on two factors: the scoring function and the search algorithm.⁸ The scoring function aims to evaluate the conformations between a protein and a ligand to detect those conformations with the minimum binding energy. The search algorithm is used to identify the most favorable ligand conformation when the ligand binds to the receptor protein.

Various optimization algorithms and docking software have been proposed for use in solving different aspects of molecular docking problems. For the docking software, a popular tool used for molecular docking is AutoDock,⁹ which uses the C++ software package. In addition, there are many other excellent software programs, such as DOCK,¹⁰ GOLD,¹¹ FlexX,¹² Glide,¹³ and PythDock,¹⁴ and so on. With the development of docking software, a number of different metaheuristics have been proposed as search methods to solve the docking problem. Although all the above mentioned works are helping in solving molecular docking problems, due to the complexity of docking free energy landscapes with rugged funnel shapes, molecular docking is still a very challenging problem.

Based on above description, the paper highlights the following: (1) We proposed a new variant version of the QPSO algorithm, called the IEQPSO algorithm. The modification involves the use of the individual particle evolutionary process to calculate the mean best position to strike a balance between the needs of the global search process and the local search (LS) process. (2) To solve the molecular docking problem, a hybrid strategy based on combining the IEQPSO algorithm with the Solis and Wets algorithm¹⁵ was proposed for dealing with the energy optimization problem. (3) The protein 1B58¹⁶ was used as an experimental receptor model, and the detailed mechanisms involved in the protein–ligand binding interactions of the compounds were elucidated from the results of the molecular docking and molecular dynamics (MD) simulations.

The remaining exposition is as follows. The upcoming section introduces the related state-of-the-art stochastic optimization methods for molecular docking problems. Then, a subsequent section describes the proposed method. The results are discussed in penultimate section, and the final section concludes the paper.

Related state-of-the-art stochastic optimization methods

Many researchers have proposed various methods for predicting molecular docking interactions. The search algorithms mainly include simulated annealing (SA), Monte Carlo (MC), genetic algorithm (GA), PSO, differential evolution (DE), and so on, and the many variants of these algorithms. For example, the GA and SA were first used in AutoDock software, and Morris et al.⁹ proposed the Lamarckian genetic algorithm (LGA), which is a hybrid of GA and LS, to minimize the binding energy. The LGA search method is more efficient in escaping from a local minimum. In addition, there are many other remarkable studies that use GA to solve protein–ligand problems. Kang et al.¹⁷ presented an improved adaptive GA to solve the problem of protein–ligand docking. Li et al.¹⁸ used an entropy-based GA with multi-population evolution to solve the docking problem and created a coefficient adaptive scoring function that greatly improves the docking accuracy. Guan et al.¹⁹ presented a novel and robust optimization algorithm based on the LGA. This algorithm applied a population evolution direction-guided model of genetics and could find solutions with lower energy for flexible protein–ligand docking problems. Leonhart et al.²⁰ developed the Biased Random Key Genetic Algorithm (BRKGA-DOCK) as a sampling strategy for the protein–ligand docking problem. The proposed method has been compared to other existing tools and obtains the best average root-mean-square deviation (RMSD). Guan et al.²¹ proposed a novel GA based on a crossover elitist preservation mechanism (crossover elitist preservation genetic algorithm mechanism (CEPGA)) for solving protein–ligand docking problems. The algorithm employed crossover elitist preservation to keep the elite individuals of the last generation and make the crossover more efficient and robust. The algorithm was compared with GA, LGA, and SODOCK, and the results indicated that CEPGA could obtain the lowest energy and the highest accuracy.

The MC algorithm has also been studied for the molecular docking problem. For example, in the study by Liu and Wang,²² an MC simulation approach was used for instances of flexible ligands. Lee et al.²³ applied conformational space annealing (CSA) and the MC algorithm with the minimization (MCM) method to generate global optimization methods to perform molecular docking simulations. The simulation results showed that the CSA method is more efficient and accurate than the MCM method.

The DE algorithm is an alternative for the molecular docking problem. Koohi-Moghadam and Rahmani²⁴ used an enhanced DE algorithm with an LS algorithm and a pseudo-elitism operator to address the docking problem. Guan et al.²⁵ introduced an efficient ABC_DE-based hybrid algorithm for protein–ligand docking, integrating artificial bee colony algorithm and DE algorithm.

More recently, Swarm Intelligence techniques were used to tackle molecular docking problems. For example, Liu et al.²⁶ introduced a fully informed swarm optimization algorithm (FIPSDock) for solving flexible protein–ligand docking problems. The García-Nieto et al.²⁷ utilized multiobjective particle swarm optimization (MPSO) variants based on different archiving and leader selection strategies to minimize the intermolecular energy and the RMSD in the predicted ligand conformations. In the study by Dai et al.,²⁸ a new MPSO algorithm based on decomposition (MPSO/D) was proposed for the molecular docking problem. The results showed that the proposed algorithm had good performance in terms of convergence and diversity. García-Nieto et al.²⁹ proposed a new variant of speed-constrained multi-objective PSO (SMPSO) algorithm, which is characterized by the use of archive strategies. In the study, the authors analyzed several SMPSO variants based on different archiving strategies for molecular docking.

Based on reports in the literature, we have observed that GA-based methods have been prominent when solving the molecular docking problem. At the same time, PSO methods have been applied to these problems more frequently in recent years.

Methods

Quantum-behaved particle swarm optimization

We assume that during QPSO with M individuals, at the nth iteration, particle i flies in the N-dimensional quantum space. $X_{i, n} = (X_{i, n}^{1}, X_{i, n}^{2}, \dots, X_{i, n}^{N})$ represents the position vector of particle i at the n^th iteration. The jth dimension of the position of particle i at the (n + 1)th iteration can be expressed by the following

X_{i, n + 1}^{j} = p_{i, n}^{j} \pm β | X_{i, n}^{j} - C_{n}^{j} | l n (1 / μ_{i, n + 1}^{j})

(1)

where

p_{i, n}^{j}

is the local attractor in the canonical PSO,² which reflects the LS for the particle. It can be defined by the following

p_{i, n}^{j} = \frac{c_{1} r_{i, n}^{j} P_{i, n}^{j} + c_{2} R_{i, n}^{j} G_{i, n}^{j}}{c_{1} r_{i, n}^{j} + c_{2} R_{i, n}^{j}}

(2)

The parameter β is the contraction-expansion (CE) coefficient. $μ_{i, n + 1}^{j}$ is a random number uniformly distributed in (0, 1), $μ_{i, n + 1}^{j} \sim U (0, 1)$ . In equation (1), $C_{n}^{j}$ is called the mean best (mbest) position, which is defined as the average of the personal best position ( $P_{i, n}^{j}$ ) of all particles.

\begin{matrix} C_{n}^{j} = (1 / M) \sum_{i = 1}^{M} P_{i, n}^{j} & 1 \leq j \leq N \end{matrix}

(3)

Individual particle evolutionary QPSO algorithm

In QPSO, the mean best position is introduced to make the algorithm more efficient than the original algorithm. By using equation (3), the mean best position is calculated using the average of the personal best position of all particles; however, the effect and influence of each particle are not exactly the same. Considering the characteristics of each particle in the evolutionary process, we designed a new control method to replace the mean best position in the QPSO in this section. A new parameter was introduced, which is based on the ratio between the fitness value of the global best solution $f (G_{n})$ and the value of the personal best position of the current particle $f (P_{i, n})$ and is given by equation (4)

α_{i, n} = | \frac{f (G_{n})}{f (P_{i, n})} |

(4)

Then, by introducing the new parameter control, the mean best position equation is rewritten as follows

\begin{array}{l} C_{n}^{j} = \frac{1}{M} \sum_{i = 1}^{M} α_{i, n}^{j} \cdot P_{i, n}^{j} \\ = (\frac{1}{M} \sum_{i = 1}^{M} α_{i, n}^{1} \cdot P_{i, n}^{1}, \dots, \frac{1}{M} \sum_{i = 1}^{M} α_{i, n}^{N} \cdot P_{i, n}^{N}) \end{array}

(5)

In equation (5), the parameter $α_{i, n}$ represents the importance of particles in the population when they are evolving and reflects the individual particle evolutionary rate. When the personal best fitness $f (P_{i, n})$ of two particles differ noticeably from each other, namely, when the swarm is far from the optimum solution, the value of parameter $α_{i, n}$ is decreased. Conversely, when the values are close to each other and differ little from $f (G_{n})$ , the search ability of the swarm is enhanced by increasing the value of the parameter. The exploration aptitude of the swarm is encouraged by incorporating the novel parameter control. The new improved QPSO algorithm has a stronger search capability than the QPSO algorithm and this allows the convergence accuracy to be higher.

Pseudocode of IEQPSO:

1	Initialize the population size, the positions and the dimensions of the particles;
2	For t = 1 to Maximum Iteration T
3	For i = 1 to population size M
4	If $f (X_{i}) < f (P_{i})$ then $P_{i} = X_{i}$ ; End if
5	$G = argmin (f (P_{i}))$ ;
6	Compute the mean best position $C_{n}$ by Eq. (5);
7	For j = 1 to N
8	$φ_{i, n}^{j} = rand (0, 1)$ ; $μ_{i, n}^{j} = rand (0, 1)$ ;
9	$p_{i, n}^{j} = φ P_{i, n}^{j} + (1 - φ) G_{n}^{j}$ ;
10	If ( $rand (0, 1) < 0.5$ )
11	$X_{i, n + 1}^{j} = p_{i, n}^{j} + β \| X_{i, n}^{j} - C_{n}^{j} \| l n (1 / μ_{i, n + 1}^{j})$ ;
12	Else
13	$X_{i, n + 1}^{j} = p_{i, n}^{j} - β \| X_{i, n}^{j} - C_{n}^{j} \| l n (1 / μ_{i, n + 1}^{j})$ ;
14	End if
15	End for
16	Evaluate the objective function value $f (X_{i, n + 1})$ ;
17	Update $P_{i, n}$ and $G_{n}$ ;
18	End for
19	End for

Results and analysis

Performance evaluation and comparison of benchmark functions

To evaluate the overall performance of the proposed IEQPSO during function optimization, the algorithm was tested with eight well-known benchmark functions (Table 1). For the performance comparison, PSO with inertia weight (PSO-In), the original QPSO algorithm and quantum-behaved particle swarm optimization with weighted mean best position (WQPSO)³⁰ were also tested with the benchmark functions. Parameter settings for these algorithms are listed in Table 2. The mean best fitness value and standard deviation based on 50 runs of each algorithm on each problem are presented in Table 3. The standard deviation values are shown in brackets.

Table 1.

Benchmark functions.

Function	Formulation	Feasible solution space
Sphere function	$f_{1} (x) = \sum_{i = 1}^{n} x_{i}^{2}$	(−100, 100)
Schwefel’s 222 function	$f_{2} (x) = \sum_{i = 1}^{n} \| x_{i} \| + \prod_{i = 1}^{n} \| x_{i} \|$	(−10, 10)
DeJong function	$f_{3} (x) = \sum_{i = 1}^{n} i x_{i}^{4}$	(−100,100)
Rosenbrock function	$f_{4} (x) = \sum_{i = 1}^{n - 1} (100 {(x_{i + 1} - x_{i}^{2})}^{2} + {(x_{i} - 1)}^{2})$	(−30, 30)
Ackley function	$f_{5} (x) = 20 + e - 20 e^{- \frac{1}{5} \sqrt{\frac{1}{n} \sum_{i = 1}^{n} x_{i}^{2}}} - e^{\frac{1}{n} \sum_{i = 1}^{n} \cos (2 π x_{i})}$	(−32, 32)
Rastrigin function	$f_{6} (x) = \sum_{i = 1}^{n} (x_{i}^{2} - 10 \cos (2 π x_{i}) + 10)$	(−5.12, 5.12)
Schaffer function	$f_{7} (x) = 0.5 + \frac{{(\sin (\sqrt{x^{2} + y^{2}}))}^{2} - 0.5}{{(1.0 + 0.001 (x^{2} + y^{2}))}^{2}}$	(−2.048, 2.048)
Griewank function	$f_{8} (x) = \frac{1}{4000} \sum_{i = 1}^{n} x_{i}^{2} - \prod_{i = 1}^{n} \cos (\frac{x_{i}}{\sqrt{i}}) + 1$	(−600, 600)

Table 2.

Parameter settings for the PSO-In, QPSO, WQPSO and IEQPSO algorithms.

Algorithm	Parameter settings
PSO-In	Population size 40	Maximum iterations 2000	Inertia weight w0.9 ∼ 0.4	c₁, c₂2
QPSO	Population size 40	Maximum iterations 2000	CE coefficient $β$ 1.0 ∼ 0.5
WQPSO	Population size 40	Maximum iterations 2000	CE coefficient $β$ 1.0 ∼ 0.5
IEQPSO	Population size 40	Maximum iterations 2000	CE coefficient $β$ 1.0 ∼ 0.5

PSO: particle swarm optimization; QPSO: quantum-behaved particle swarm optimization; WQPSO: QPSO with weighted mean best position; IEQPSO: individual particle evolutionary QPSO; CE: contraction-expansion.

Table 3.

Results of the benchmark functions.

Function	PSO-In	QPSO	WQPSO	IEQPSO
f1	1.192430e-010 (2.572647e-010)	2.606978e-029 (1.123051e-028)	1.550875e-079 (4.496400e-079)	9.525978e-086 (4.203196e-085)
f2	3.206959e-008 (3.499209e-008)	4.624757e-018 (1.591423e-017)	2.408737e-049 (3.648797e-049)	1.669145e-052 (3.084726e-052)
f3	8.048631e-010 (3.613984e-009)	6.201183e-037 (3.341504e-036)	2.146062e-109 (1.037800e-108)	4.094578e-092 (2.636037e-091)
f4	7.814176e+001 (4.647823e+001)	4.378334e+001 (3.173342e+001)	6.730166e+001 (4.943229e+001)	2.361887e+001 (2.212292e-001)
f5	2.038785e-005 (1.017194e-004)	1.541008e-009 (3.804749e-014)	1.540960e-009 (8.437228e-016)	1.540978e-009 (8.220215e-015)
f6	3.285426e+001 (8.221160e+000)	2.152986e+001 (5.855429e+000)	9.381218e+001 (2.941288e+001)	1.783571e-021 (1.354936e-020)
f7	3.886364e-004 (1.903922e-003)	5.830047e-004 (2.307388e-003)	3.902530e-004 (1.903624e-003)	5.5511e-018 (3.4528e-016)
f8	1.312913e-002 (1.435475e-002)	6.597326e-003 (8.014635e-003)	6.254397e-003 (1.140907e-002)	8.881784e-018 (6.217249e-017)

PSO: particle swarm optimization with inertia weight; QPSO: quantum-behaved particle swarm optimization; WQPSO: QPSO with weighted mean best position; IEQPSO: individual particle evolutionary QPSO. The optimal solution of each problem is presented in boldface in Table 3 for four algorithms.

For the shifted Sphere Function (f1), all the algorithms could obtain a solution near the global optimal solution, but the IEQPSO algorithm obtained better results than the other methods. The results for the shifted Schwefel’s 222 function (f2) show that IEQPSO could find a global optimum successfully and exhibited a better search ability than PSO-In and QPSO. For the DeJong function (f3), the WQPSO algorithm outperformed the other methods. In regard to benchmark f4, the Rosenbrock function is a multimodal function, which is too difficult to solve for many optimization problems. The IEQPSO algorithm showed a slightly better performance among all the algorithms, and its obvious advantage over other algorithms is shown by its standard deviation values. The fifth benchmark function, f5, is the Ackley function, and the results indicate that all the QPSO-based methods showed better performance for this problem than PSO. For function f6, the Rastrigrin function, the IEQPSO algorithm yielded the best result. For this function, the difference between QPSO and PSO-In was not significant. Based on the results obtained for f7, the Schaffer problem, IEQPSO had a significant advantage over the other methods. The results for the shifted, rotated Griewank’s function (f8) suggest that the IEQPSO algorithm was able to find a better solution than the others.

Figure 1 presents the convergence process of the algorithms for each benchmark function. It shows that the IEQPSO algorithm had a better convergence property than its competitors during the later stage of iteration, except for functions f3 and f5. According to the results shown above in the tables and figures, it can be seen that the IEQPSO algorithm was able to find solutions of higher quality for the functions compared to those obtained using the other methods, in most cases.

Figure 1.

Convergence process of the algorithms for each benchmark function.

Results for molecular docking

Because the energy search during molecular docking is complex, we used a hybrid strategy based on combining the IEQPSO algorithm with the Solis and Wets algorithm to address energy optimization that was referred to as the hybrid algorithm of IEQPSO and local search (LIEQPSO). The IEQPSO algorithm is used as a global search algorithm, and the Solis and Wets algorithm is used as an LS algorithm. The addition of an LS algorithm effectively maintains the diversity of the particles and prevents premature convergence of the algorithm. In addition, the LGA algorithm was used as the other search algorithm. It is the classical algorithm used in AutoDock 4, in which the GA algorithm incorporates the Solis and Wets algorithm during the LS. The main motivation for this is to analyze the docking performance of both algorithms to determine which of them can obtain better scoring values. To achieve this goal, we attempted to incorporate the new algorithm into the source code of AutoDock 4.

Problem description

AutoDock 4 software is used to evaluate a docked conformation. The scoring function of AutoDock 4 is the semiempirical free energy function. The total docked energy is expressed as the sum of the intramolecular energy and intermolecular energy of the protein and ligand. The intramolecular energy was estimated for the transition from the unbound state to the bound state of the protein and ligand. The intermolecular energy was estimated for the combination of the protein and the ligand in the bound conformation.

The force field includes six pairwise estimates (V) and the evaluation of the conformational entropy lost upon binding ( $Δ S_{conf}$ ):

\begin{array}{l} Δ G = (V_{bound}^{L - L} - V_{unbound}^{L - L}) + (V_{bound}^{P - P} - V_{unbound}^{P - P}) \\ (V_{bound}^{P - L} - V_{unbound}^{P - L} + Δ S_{conf}) \end{array}

(6)

where P refers to the “protein” and L refers to the “ligand” in the protein–ligand docking calculation.

Each of the pairwise energetic items consists of the following sources of energy: the first source is the Lennard-Jones 12–6 Van der Waals interaction energy, the second source is the 12–10 hydrogen bond potential energy, the third source is the electrostatic potential energy, and the final is the desolvation potential energy.

\begin{array}{l} V = W_{vdw} \sum_{i, j} (\frac{A_{i j}}{r_{i j}^{12}} - \frac{B_{i j}}{r_{i j}^{6}}) + W_{hbond} \sum_{i, j} E (t) (\frac{C_{i j}}{r_{i j}^{12}} - \frac{D_{i j}}{r_{i j}^{10}}) \\ + W_{elec} \sum_{i, j} \frac{q_{i} q_{j}}{ε (r_{i j}) r_{i j}} + W_{sol} \sum_{i, j} (S_{i} V_{j} + S_{j} V_{i}) e^{(- r_{i j}^{2} / 2 σ^{2})} \end{array}

(7)

Experimental results for protein–ligand docking

The protein and ligand docking tests were run on the Red Hat Enterprise Linux OS. The prepared protein and ligand structures were saved in the PDBQT file format. The AutoDock package includes two programs, the AutoGrid program, and the AutoDock program. The AutoGrid program performs the calculation of the energy grid maps. The grid size was set to 60 × 60 × 60 points with a spacing of 0.375 Å. The AutoDock program is responsible for the conformational search and energy evaluation. Table 4 shows the parameter settings of the search algorithms used for docking experiments. The other parameters used the default settings in AutoDock 4.

Table 4.

Parameter settings for protein–ligand docking.

Number of particles	50
Number of algorithm runs	30
Maximal number of functional evaluations	2.5 × 10⁵
Maximum number of generations	2.7 × 10⁴
Number of local search iterations	300
The probability of performing a local search	0.06

For the test case, PDB 1B58 was used as the experimental receptor protein and was obtained from the RSCB protein data bank (http://www.rcsb.org/). The sequence Lys-Tyr-Lys was used as the experimental ligand and contained 43 atoms.

There are two criteria for evaluating molecular docking results: the binding energy and RMSD of docking. The RMSD represents the comparison between atomic positions in a docked ligand and the corresponding positions in the ligand crystal structure. In general, we assumed that the lower the binding energy, the closer the docked complex was to the native state of the complex. The lower RMSD value is assumed to have higher docking accuracy. An RMSD value of 2 Å was usually used as the cutoff value. For the LIEQPSO docking complex, the best binding result had an energy = −19.68 kcal/mol and an RMSD = 1.63 Å. The best result for the LGA docking complex had an energy = −12.74 kcal/mol and an RMSD = 2.00 Å. According to the lowest energy and RMSD values, the LIEQPSO algorithm obtained a better docking result than the LGA algorithm.

Among the sources of energy involved in docking binding, Van der Waals energy and electrostatic interaction energy are the principal binding forces in protein–ligand interactions. For Van der Waals energy, both attractive and repulsive interactions play a vital role in controlling the binding phenomenon. The electrostatic interaction plays an important role in changing the enthalpy resulting from binding. The Van der Waals interaction energy and the electrostatic binding energy of the ligand residue atoms are plotted in Figure 2. A more negative value for binding energy indicates favorable binding of ligand with the protein receptor. Figure 3 shows examples of complex docking obtained by the two methods. Based on the above results, the LIEQPSO algorithm has a higher docking accuracy than the LGA algorithm for the test case.

Figure 2.

Van der Waals interaction and electrostatic interaction energy calculated for the ligand residue atoms of the complexes docked by the LIEQPSO and LGA algorithms.

Figure 3.

Examples of complex docking obtained by the (a) LIEQPSO algorithm and (b) LGA algorithm. The receptor protein is represented by a novel cartoon secondary structure style. The color of the crystal structure is pink, and the predicted structures are purple (a) and cyan (b).

MD simulation results for protein–ligand docking

The objective of docking is to obtain a model of the interaction between a ligand and a protein that is characterized by a minimum binding energy. For further detailed analysis of protein–ligand interactions, we performed MD simulations with the docking results to establish a more reliable method to illustrate protein–ligand interactions. The interaction of proteins with ligands is very important in biological processes in the context of drug design and pharmacology, since proteins usually accomplish their functions through binding to various molecules. In this section, we used NAMD software³¹ and the CHARMM27 force field³² to perform the MD simulations. The simulated temperature was set to 300 K, and the simulation time was set to 10 ns.

Complex stability analysis

To evaluate the overall structural stability and flexibility of the complex structure, we calculated the RMSD of the complex structure. The RMSD is defined as follows.

RMSD = \sqrt{\frac{\sum_{i = 1}^{N_{atoms}} {(r_{i} (t_{1}) - r_{i} (t_{2}))}^{2}}{N_{atoms}}}

(8)

where

N_{atoms}

is the number of atoms whose positions are being compared, and

r_{i} (t)

is the position of atom i at time t.

The RMSD of the backbone atoms in the complex versus the simulation trajectory for 10 ns was calculated. As shown in Figure 4(a), we noticed that the RMSD of the LIEQPSO docking complex structure fluctuated around a value of 2.5 Å after 6 ns. The average RMSD value is 2.05 Å for the backbone atoms. In the MD simulation of the LGA docked complex, the backbone RMSD fluctuated between 0.49 Å and 2.78 Å during the simulation. The average RMSD value is 1.99 Å for the LGA docked complex. The plots showing the RMSD of the backbone atoms could indicate the stability of the complexes in the MD simulation.

Figure 4.

(a) RMSD of the backbone atoms as a function of time in the docked complex and (b) RMSD as a function of the residue number in the docked complex at a temperature of 300 K.

To investigate whether there are any specific dynamic residues involved in the protein–ligand interaction, the RMSD of the residues was calculated. This can indicate the flexibility of the protein residues. It can also clearly represent the differences in some residues in both complexes. As shown in Figure 4(b), a comparison of the RMSD of the structures of the LIEQPSO complex and the LGA complex revealed that the RMSD curves exhibited great similarity in their fluctuation trends, with only a few residues showing a significant difference. For instance, the RMSD values of residue PHE44 according to the LGA docking results were 7.68 Å, whereas the RMSD was 2.88 Å in the LIEQPSO docking complex. Similar residues, including PRO56, VAL221, and LEU512, among others, also showed higher RMSD values in the LGA docking complex. In general, the flexibility of the residues in the LIEQPSO docking complex was relatively lower than those in the LGA docking complex.

Salt bridge between interfacial residues

In the interactions between proteins and ligands, electrostatic interactions play an important role due to their high specificity. We used a salt bridge to examine the contribution of electrostatic interactions. Here, we have only considered the salt bridge present at the interface of the complex. A salt bridge was identified according to a cutoff distance of 4.0 Å between any of the nitrogen atoms in the basic residues and the oxygen atoms in acidic residues. It was found that there were two salt bridges in the interfacial residues in the complex formed by docking by the LIEQPSO algorithm at 300 K (Figure 5(a)). LYS1 of the ligand forms a salt bridge with the ASP419 residue of the protein, and LYS3 of the ligand forms a salt bridge with the GLU229 residue of the protein. In the complex that was docked by the LGA algorithm, one salt bridge was found, and the LYS1 residues in the ligand formed a salt bridge with ASP419 in the protein (Figure 5(b)). As shown in Figure 5, the salt bridge LYS1-ASP419 in the complex docked by the LIEQPSO algorithm was not as stable as that in the complex docked by the LGA algorithm. The salt bridge ASP419-LYS1 in the LGA docked complex was mostly maintained at a distance of 4.0 Å. In comparison, the salt bridge ASP419-LYS1 in the LIEQPSO docked complex was ruptured and restored throughout the simulation process. Fortunately, the two salt bridges that existed in the complex docked by the LIEQPSO algorithm were believed to promote the stability of the complex.

Figure 5.

Plots showing changes in the salt bridge in the binding domain as a function of time. (a) Changes in the salt bridges ASP419-LYS1 and GLU229-LYS3 in the LIEQPSO docking complex and (b) changes in the salt bridge ASP419-LYS1 in the LGA docking complex.

Hydrogen bond in between interfacial residues

Hydrogen bonds are also important factors that influence the binding stability of protein ligands. The presence of hydrogen bonds is an important indication of the significance of protein–ligand interactions. In the complex, interfacial residues forming hydrogen bonds are shown in Figure 6. For the complexes docked by both algorithms, hydrogen bond residues within the protein are mainly concentrated near ARG413, CYS417, and VAL34. The obvious distinction is that the hydrogen bonds of both complexes have different occupancy times.

Figure 6.

Hydrogen bond map between protein and ligand interfacial residues in the complex.

For the complex after molecular simulation, a distance cutoff of 3.5 Å and an angle cutoff of 30° were used for the hydrogen bond calculation. According to the above mentioned hydrogen bond calculation criteria, we have listed the hydrogen bonds in interfacial residues that exist >50% of the time throughout the entire simulated trajectory at 300 K. We observed 10 hydrogen bonds in the LIEQPSO docked complex and 10 hydrogen bonds in the LGA docked complex. Among these hydrogen bonds, seven common hydrogen bonds, ARG413::LYS3, CYS417::LYS1, ASP419::LYS1, GLU32::TYR2, VAL34::TYR2, ARG404::TYR2, and LYS1::HSD161, were found in both complexes. Although ALA418::LYS1, ASN247::LYS3, and ASN246::LYS3 were also found in the LGA docked complex, they exist for a relatively short time during the MD trajectory. The results show that three hydrogen bonds, ARG413::LYS3, CYS417::LYS1, and ASP419::LYS1, were the most stable in both complexes. In general, the hydrogen bonds in the LIEQPSO docked complex have higher occupancy time in the simulated trajectory (Table 5).

Table 5.

The major hydrogen bonds in the binding domains of both complexes.

Docking algorithm	Hydrogen bonds
LIEQPSO	ARG413::LYS3, CYS417::LYS1, ASP419::LYS1, ALA418::LYS1, GLU32::TYR2, VAL34::TYR2,ARG404::TYR2, HSD161::LYS1, ASN246::LYS3, ASN247::LYS3
LGA	ARG413::LYS3, ARG404::TYR2, TRP416::LYS1, ASP419::LYS1, CYS417::LYS1, GLY415::LYS3, VAL34::TYR2, TYR485::LYS3, GLU32::TYR2, HSD161::LYS1

LGA: Lamarckian genetic algorithm; LIEQPSO: hybrid algorithm of IEQPSO and local search.

Conclusion

This paper describes an improved QPSO algorithm based on individual particle evolutionary processes. A hybrid algorithm of IEQPSO and the Solis and Wets method was used to study molecular docking problems. Finally, we presented an MD simulation investigation on the binding interactions between proteins and ligands.

To assess the performance of the IEQPSO algorithm, eight test cases with unimodal functions and multimodal functions were used for comparison. The results indicated that the IEQPSO algorithm significantly improved the search ability and convergence precision compared to other algorithms.

To solve the molecular docking problem, the ligand can be bound to the protein by the LIEQPSO and LGA algorithms. The MD simulation results indicated that there were complex binding interactions between the protein and the ligand. The results indicate that the proposed LIEQPSO algorithm outperforms the LGA algorithm in terms of optimal solutions. Based on the above analysis, the proposed LIEQPSO algorithm may be used as an alternative method for solving molecular docking problems.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the National Natural Science Foundation of China under Grant Nos. 61502203, the Natural Science Foundation of Jiangsu Province under Grant Nos. BK20150122, the Natural Science Foundation of the Jiangsu Higher Education Institutions of China under Grant Nos. 17KJB520039 and Nos. 18KJB520046, the 333 High-level Talents Training Project of Jiangsu Province under Grant Nos. BRA2018147, “Qing Lan Project” of Wuxi City College of Vocational Technology.

References

Sun

Feng

Particle swarm optimization with particles having quantum behavior. In: Proceedings of the 2004 congress on evolutionary computation (IEEE Cat. No.04TH8753), Portland, OH, 19–23 June 2004.

Clerc

Kennedy

The particle swarm-explosion, stability, and convergence in a multidimensional complex space. IEEE Trans Evol Computat 2002; 6: 58–73.

Tharwat

Hassanien

AE.

Quantum-behaved particle swarm optimization for parameter optimization of support vector machine. J Classif. Epub ahead of print 21 January 2019. DOI: 10.1007/s00357-018-9299-1.

Dereli

Köker

A meta-heuristic proposal for inverse kinematics solution of 7-DOF serial robotic manipulator: quantum behaved particle swarm algorithm. Artif Intell Rev. Epub ahead of print 30 January 2019. DOI: 10.1007/s10462-019-09683-x.

García-Godoy

López-Camacho

García-Nieto

, et al. Bio-inspired optimization for the molecular docking problem: State of the art, recent results and perspectives. Appl Soft Comput 2019; 79: 30–45.

Kawada

Ebiike

Tsukazaki

, et al. Lead optimization of a dihydropyrrolopyrimidine inhibitor against phosphoinositide 3-kinase (PI3K) to improve the phenol glucuronic acid conjugation. Bioorg Med Chem Lett 2013; 23: 673–678.

Sengupta

Obiorah

Maximov

, et al. Molecular mechanism of action of bisphenol and bisphenol A mediated by oestrogenreceptor alpha in growth and apoptosis of breast cancer cells. Br J Pharmacol 2013; 169: 167–178.

Aleksandrov

Myllykallio

Advances and challenges in drug design against tuberculosis: application of in silico approaches. Expert Opin Drug Discov 2019; 14: 35–46.

Morris

Goodsell

Halliday

, et al. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J Comput Chem 1998; 19: 1639–1662.

10.

Lang

Brozell

Mukherjee

, et al. DOCK 6: combining techniques to model RNA-small molecule complexes. Rna. 2009; 15: 1219–1230.

11.

Jones

Willett

Glen

, et al. Development and validation of a genetic algorithm for flexible docking. J Mol Biol 1997; 267: 727–748.

12.

Kramer

Rarey

Lengauer

Evaluation of the FlexX incremental construction algorithm for protein-ligand docking. Proteins 1999; 37: 228–241.

13.

Friesner

Banks

Murphy

, et al. Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J Med Chem 2004; 47: 1739–1749.

14.

Chung

Cho

Hah

JM.

A python-based docking program utilizing a receptor bound ligand shape: PythDock. Arch Pharm Res 2011; 34: 1451–1458.

15.

Solis

Wets

Minimization by random search techniques. Math Oper Res 1981; 6: 19–30.

16.

Sleigh

Seavers

Wilkinson

, et al. Crystallographic and calorimetric analysis of peptide binding to OppA protein. J Mol Biol 1999; 291: 393–415.

17.

Kang

Jiang

, et al. An improved adaptive genetic algorithm for protein-ligand docking. J Comput Aided Mol Des 2008; 23: 1–12.

18.

Zhuang

, et al. Adaptive molecular docking method based on information entropy genetic algorithm. Appl Soft Comput 2015; 26: 299–302.

19.

Guan

Zhang

Ning

EDGA: a population evolution direction-guided genetic algorithm for protein-ligand docking. J Comput Biol 2016; 23: 585–596.

20.

Leonhart

Spieler

Ligabue-Braun

, et al. A biased random key genetic algorithm for the protein-ligand docking problem. Soft Comput 2018; 23: 4155–4176.

21.

Guan

Zhang

Ning

Genetic algorithm with a crossover elitist preservation mechanism for protein-ligand docking. AMB Expr 2017; 7: 174.

22.

Liu

Wang

MCDOCK: A Monte Carlo simulation approach to the molecular docking problem. J Comput Aided Mol Des 1999; 13: 435–451.

23.

Lee

Czaplewski

Kim

, et al. An efficient molecular docking using conformational space annealing. J Comput Chem. 2004; 26: 78–87.

24.

Koohi-Moghadam

Rahmani

AT.

Molecular docking with opposition-based differential evolution. In: Proceedings of the 27th annual ACM symposium on applied computing, Trento, Italy, 26–30 March 2012, pp. 1387–1392. New York, NY: ACM.

25.

Guan

Zhang

Zhao

An efficient ABC_DE_based hybrid algorithm for protein-ligand docking. Int J Mol Sci 2018; 19: 1181.

26.

Liu

Zhao

, et al. FIPSDock: A New molecular docking technique driven by fully informed swarm optimization algorithm. J Comput Chem 2012; 34: 67–75.

27.

García-Nieto

López-Camacho

García-Godoy

, et al. Multi-objective ligand-protein docking with particle swarm optimizers. Swarm Evol Comput 2018; 44: 439–452.

28.

Dai

Wang

A new multi-objective particle swarm optimization algorithm based on decomposition. Inform Sci 2015; 325: 541–557.

29.

García-Nieto J, López-Camacho E, Godoy MJG, et al. A Study of Archiving Strategies in Multi-objective PSO for Molecular Docking. In: Dorigo M. et al. (eds) Swarm Intelligence. ANTS 2016. Lecture Notes in Computer Science, vol 9882. Springer, Cham, 2016.

30.

Sun

An improved quantum-behaved particle swarm optimization algorithm with weighted mean best position. Appl Math Comput 2008; 205: 751–759.

31.

Phillips

Braun

Wang

, et al. Scalable molecular dynamics with NAMD. J Comput Chem 2005; 26: 1781–1802.

32.

MacKerell

Bashford

Bellott

, et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins and nucleic acids. J Phys Chem B 1998; 102: 3586–3616.