IDENTIFICATION OF PHARMACOLOGICAL TARGETS COMBINING DOCKING AND MOLECULAR DYNAMICS SIMULATIONS

Studies that include both experimental data and computational simulations ( in silico ) have increased in number because the techniques are complementary. In silico methodologies are currently an essential component of drug design; moreover, identification and optimization of the best ligand based on the structures of biomolecules are common scientific challenges. Geometric structural properties of biomolecules explain their behavior and interactions and when this information is used by a combination of algorithms, a dynamic model based on atomic details can be produced. Docking studies enable researchers to determine the best position for a ligand to bind on a macromolecule, whereas Molecular Dynamics (MD) simulations describe the relevant interactions that maintain this binding. MD simulations have the advantage of illustrating the macromolecule movements in more detail. In the case of a protein, the side chain, backbone and domain movements can explain how ligands are trapped during different conformational states. Additionally, MD simulations can depict several binding sites of ligands that can be explored by docking studies, sampling many protein conformations. Following the previously mentioned strategy, it is possible to identify each binding site that might be able to accommodate different ligands through atomic motion. Another important advantage of MD is to explore the movement of side chains of key catalytic residues, which could provide information about the formation of transition states of a protein. All this information can be used to propose ligands and their most probable site of interaction, which are daily tasks of drug design. In this review, the most frequent criteria that are considered when determining pharmacological targets are gathered, particularly when docking and MD are combined.


INTRODUCTION
Experimental techniques of molecular biology can be used to explore the intrinsic mechanisms of storage and transmission of information within the cell. In particular, the cloning and purification of a protein can permit its study by Nuclear Magnetic Resonance (NMR) or X-ray studies. A useful result of these techniques is the structural chemistry of the protein (tridimensional models, 3D), which helps to elucidate its biological Science Publications AJABS properties. Using the 3D structure of a protein, molecular recognition studies between protein-protein or proteinligand can be achieved and used to explain biological events. Signal transduction and enzyme activation are now explained by complex structures that are determined by the aforementioned techniques and the design and production of new drugs depend largely on understanding how the functions of the proteins are inhibited or activated (Klaic et al., 2012;Baud et al., 2012). Furthermore, in addition to binding studies, there are other drug properties considered when designing drugs, such as absorption, distribution, metabolism, excretion and toxicity, among others. These can be studied using computational approaches (Alonso et al., 2006).
There are theoretical methods that include algorithms to solve the problem of molecular recognition. Initially, docking simulations considered the protein molecule and the ligand to be rigid elements with a specific geometry, for which energetic functions were evaluated to find an optimal coupling (Yue, 1990;Aqvist and Tapia, 1992;DesJarlais et al., 1986). However, various recent methods for the docking of flexible proteins and ligands have been developed Andrusier et al., 2008). Docking simulations depend on initial positions of the atoms of participating molecules and the sampling of the conformational space is deficient. Therefore, improvements in these docking methods and even their conjugation with other algorithms are strongly recommended. The incorporation of flexibility properties in the molecular recognition process increases the likelihood of finding suitable complexes, i.e., the refinement of the atomic positions can be obtained as a result of the mutual interdependence of participating atoms. This description is offered by Molecular Dynamics (MD) simulations (Feher and Williams, 2012;Coupez and Lewis, 2006;Alonso et al., 2006). MD simulations of biomolecules permit the construction of a hypothesis of molecular mechanisms that are involved in a biological phenomenon, explaining them by the behavior of their constituent atoms (Rahman et al., 2012;Lin et al. 2012;Tsai et al., 2012;Rosas-Trigueros et al., 2011). MD uses a numerical method for solving the Newtonian equations for the atoms of biomolecules. The solution gives the consecutive positions and velocities of particles subjected to a potential function derived from an empirical model that approximates atomic interactions in terms of classical mechanics. MD simulations also offer an ensemble of conformations that provide an interpretation of sampling in the scope of mechanical statistics, where the mean values of observables describe a system (Cuendet, 2006;Lee et al., 2009). As with docking techniques, computational processing costs limit the MD simulations; moreover, this could be critical if the size and the number of molecules grow considerably. Consequently, the algorithms, codes and methods available are constantly improved (Gotz et al., 2012;Bauer et al., 2011;Pool et al., 2012;Zhao et al., 2012). This study reviews the potential of strategies for computational molecular recognition that combine docking and MD simulations to gain insight into molecular behaviors and conditions to identify pharmacological targets.

Methods Used for Performing MD Simulations
As stated above, several docking methods have yielded promising results by considering binding models that take into account the atomistic motions of receptors and ligands. Without doubt, it is particularly beneficial to visualize how various proteins move and modify their shape, atom-by-atom, as a function of time, as they perform their functions (e.g., catalytic activity) (Roux, 2010). MD simulation methods simulate atomic motions, yielding molecular trajectories. This information can be used to calculate thermodynamic quantities and estimate binding affinities and kinetic rates. Unlike other computational methods such as Monte Carlo methods, MD methods allow us to directly examine dynamical processes, driven by the finite temperature of the simulation (Tilocca, 2012). MD simulations have also been used to examine the cooperativity in DNA-drug recognition and Virtual Screening (VS) of snapshots from MD simulations has been successfully used for position prediction and ranking of compound libraries. Molecular flexibility and binding properties can help to identify specific protein-ligand complexes at varying frequencies along typical MD simulations and the sampling of molecular movements can be performed on different time scales (Nichols et al., 2011;Harris et al., 2001). Additionally, the reduction of cost associated with these methods has made their intensive use for drug discovery attractive in recent times (Harvey and Fabritiis, 2012;Ou-Yang et al., 2012).
Conventional MD simulation methods evaluate the time evolution of a system by numerically integrating Newton's equations of motion. A molecule is considered to be a collection of spheres corresponding to atoms with a fixed electronic distribution. The molecular mechanical model considers interactions between bonded atoms where bonds are modeled by springs with a Hookean pairwise potential. Interactions between atoms connected by consecutive bonds usually include triplets (angles) and quadruplets (dihedral and improper angles). Together with nonbonded interactions (Coulomb and van der Waals), this model yields the potential energy for a Science Publications AJABS fixed conformation, where the negative gradient provides the force needed for the numerical calculation of the system trajectories. To approximate the physical behavior of real molecules in motion, the energy terms mentioned above are parameterized to fit quantummechanical calculations and experimental data. The parameters thus obtained, together with the equations, are collectively called a "force-field" and although the parameterization is performed with a limited training set of molecules, a force-field is expected to be transferable to similar molecules. However, the reliability of forcefields is often questioned (Paton and Goodman, 2009). Commonly used force fields for MD simulations of biochemical systems include OPLS-AA (Jorgensen and Tirado, 1988), CHARMM (MacKerell et al., 1998) and AMBER (Cornell et al., 1995), which were developed or optimized for the simulation of proteins (Mackerell, 2004) and have yielded similar results in simulations of peptidic chains with 165 amino acids or less (Price and Brooks, 2002). MD methods include numerical algorithms that manipulate the pressure, temperature and volume of the system. These algorithms allow researchers to simulate molecules under different conditions (Schlick, 2010).
Molecular simulation protocols (docking and MD simulations) start with an initial 3D model of the system, which can be obtained from NMR or crystallographic data, collected in the Protein Data Bank (PDB) (Berman et al., 2000). When such information is not available for a peptidic chain, a protein structure prediction method can be used to obtain a 3D model that can later be used as input for the MD simulation. Among protein structure prediction methods, homology modeling has become very popular for obtaining the initial coordinates for a protein, providing coordinates based on homologous sequences with an experimentally determined structure; the obtained models can be further refined using MD simulations (Nurisso et al., 2012). When homology modeling is not feasible, several other protein structure prediction methods are available, as has been shown in the biannual CASP experiment (Kinch et al., 2011). With the protein coordinates as a starting point, a set of boundaries are defined for the system and the space is filled with solvent molecules. Once the initial model is ready, the trajectory is calculated in a stepwise manner. In each step, the forces acting on each atom are calculated and each atom is moved according to those forces in a small time step (1 or 2 fs). The positions and velocities of the atoms are stored for subsequent analyses, together with other relevant data such as energy and pressure (Durrant and McCammon, 2011). The initial model often presents interactions with artificially high potential energy (e.g., overlapping atoms). To alleviate this problem, the system is typically subjected to one or more rounds of structural optimization. The desired temperature of the system is then slowly reached, typically within the NVT ensemble (constant number of particles, volume and temperature). An equilibration stage follows, often in the NPT ensemble (constant number of particles, pressure and temperature), to allow the system density to converge and for the structure to relax. The desired relaxation is evaluated by checking the convergence of timedependent system properties such as energy, density, temperature, pressure and root mean square deviation (RMSD) to the initial structure. The equilibrated system can now enter a production phase, in which the goal is to generate enough representative conformations in a trajectory to satisfy the ergodic hypothesis, which states that the average values over time of physical quantities that characterize a system are equal to the statistical average values of these quantities. If enough representative conformations are sampled, the relevant biophysical properties can then be calculated (Nurisso et al., 2012).
Although MD simulations of systems with up to 1 million atoms for over 50 ns have been reported (Freddolino et al., 2006), considerable efforts continue to be devoted to the development of theories, algorithms, software and hardware for the purposes of reducing the cost of performing MD simulations. These contributions facilitate application to systems of larger sizes and performance of simulations of a larger number of systems of interest. Also, these improvements permit to extend the simulation time to the millisecond range, which has already been reached for a 58-residue protein (Shaw et al., 2010). Among these efforts, Coarse-Graining (CG) reduces some of the accuracy in speeding up the calculations. In CG, some sort of vast approximation is made to greatly increase the simulation speed. For instance, certain united atom methods of CG represent groups of atoms with one large pseudo-atom for representing the overall properties of the represented atoms (Balabin et al., 2009). A large number of calculations are employed on solvent molecules. Implicit solvent models attempt to reduce explicit solute-solvent interactions to their mean field characteristics, which are expressed as a function of the solute configuration alone (Feig and Brooks, 2004). While implicit solvent modeling has been shown to study well with proteins (Chopra et al., 2008), these methods may not be adequate to study systems where solvent molecules have been shown to be crucial for binding affinity, such as the case of Small-Molecule (SM) inhibitors of blood coagulation reported by Abel et al. (2011). The

AJABS
Accelerated Molecular Dynamics (AMD) paradigm is another method available, where a non-negative bias potential is added to the potential energy when the latter lies below a certain threshold (Hamelberg et al., 2004). This increase in the potential energy accelerates the exchange between low-energy conformational states while still maintaining the essential details of the potential energy surface. This paradigm has been shown to yield accurate free energy statistics and its potential application in studying functional dynamics in biomolecules has been evaluated by Mangoni and McCammon (2011). In the context of preventing MD simulations from getting trapped in local minimumenergy states, the method of Replica-Exchange Molecular Dynamics (REMD) can be of help. In a REMD study, several MD simulations of the same system (replicas) are started at different temperatures. Pairs of replicas corresponding to neighboring temperatures are exchanged periodically, which results in a random walk in the "temperature space" that allows the simulation to escape from local minima-energy states (Sugita and Okamoto, 1999). The combination of REMD generalized ensemble sampling with ensemble docking and free energy pathway analysis has been recently proposed as a novel research protocol for the simulation of protein-ligand induced-fit recognition (Park and Li, 2010). A powerful technique for reconstructing the Free-Energy Surface (FES) as a function of few selected degrees of freedom is known as metadynamics, first introduced by Laio and Parrinello (2002). In metadynamics, the exploration of the FES is guided by forces along the selected degrees of freedom and these forces are scaled by the estimated size of the FES basin, which would allow the CG dynamics of the system to escape the local minima. After their calculation, the forces are replaced by a history-dependent term that discourages the system from revisiting points, thus encouraging an efficient exploration of the FES. The difficulty of choosing an appropriate set of degrees of freedom has been the focus of recent developments in metadynamics (Barducci et al., 2011) and the applications of this technique have expanded to binding profile determination, as in the study by Fidelak et al. (2010). Another MD based protocol worth mentioning is Steered Molecular Dynamics (SMD), developed with the aim of providing atomic level descriptions of the underlying events in single-molecule measurement techniques. SMD applies forces that are external to the force field model to investigate the mechanical properties of biopolymers and accelerate processes that are otherwise too slow to simulate (Isralewitz et al., 2001). SMD can be used to pull the ligand from a ligand-receptor complex to obtain the irreversible mechanical work necessary for the undocking and to discern active from inactive enzyme inhibitors, as reported by Colizzi et al. (2010).
Other commonly used methods with lower computational cost are continuum electrostatic calculations and Brownian Dynamics (BD), a variation of MD in which the use of approximations makes long timescale calculations possible (Dodson et al., 2008). Another time-saving approach used to study the dynamics of a biomolecule is Normal Mode Analysis (NMA), in which the simple harmonic motions of the molecule about a local energy minimum are calculated by means of a Hessian matrix built from a Hookean potential model of the system. This method has been suggested as a preliminary step in drug design (Floquet et al., 2006). Although NMA does not produce time-dependent trajectories, it can provide insight into the large-scale and long-time conformational motions of proteins (Bahar and Rader, 2005). Another problem with MD simulations is that, despite the continuous enhancement of force-field parameters over time (Lindorff-Larsen et al., 2012), the native conformation may not be the lowest free energy state of the system (Nurisso et al., 2012). In particular, some force fields tend to favor commonly seen secondary structures such as alpha helices and beta sheets (Best et al., 2008;Patapati and Glykos, 2011;MacKerell et al., 1998).
Different approaches have aimed to expand the scope of physical and chemical phenomena that can be represented by conventional MD models and thus allow them to simulate processes such as charge transfer between atoms and breaking/making bonds. These proposals include polarizable models, which simulate induced point dipoles Wang et al., 2011) and reactive force-fields, which simulate bond breaking and formation (Farah et al., 2012), as well as methods that include quantum mechanics information (Braga, 2012). Despite the advantages of these enhanced MD simulations, their high computational cost tends to make them unattractive for the simulation of biomolecular systems.

Methods for Achieving Docking Studies
Every bimolecular reaction begins with a recognition event for which shape complementarity, as well as the chemical surfaces involved, are crucial for reaching the Science Publications AJABS most stable ligand-protein complexes. The intermolecular forces of recognition between molecules are often weak and of short range and include mechanical anchoring, hydrogen bonds, metal bonds, salt bridges and aromatic stacking (Rebek, 2009). Molecular recognition is useful in biology because biological processes such as the formation of the double helical structure of DNA (Leblanc, 2006;D'Abramo et al., 2012), the ligandreceptor interaction (Okazaki and Takada, 2008;Razzaghi-Asl et al., 2012) and the enzyme-substrate interaction (Vitorovic-Todorovic et al., 2010;Mazur et al., 2010) are driven by this phenomenon.
Docking methods can depict the ligand-receptor interaction, which can be divided into three general strategies: (1) rigid ligand and receptor; (2) flexible ligand and rigid receptor; and (3) flexible ligand and flexible receptor. Although limited, the interaction between rigid bodies in a three-dimensional coordinate space is computationally affordable. A certain amount of flexibility is usually initially calculated for the ligand, generating a set of ligand conformations before intensive calculations are initiated (Meng et al., 2011). In the case of a flexible ligand and rigid receptor, the structures involved are considered complementary following the Koshland's Induced-Fit (IF) model, where the fit occurs only after the structural changes, induced by the ligand itself, take place (Koshland, 1963;Hammes, 2002) (Fig.  1). This suitable coupling, however, is improved by including receptor flexibility, although a high computational cost is incurred. It is worth mentioning that MD simulations are an alternative to acquire the conformations of the ligand-receptor system, although the technique requires high computational resources. Perhaps the screening of a large chemical database could make MD simulations unaffordable.
Because there are several different positions of atoms into an interface ligand-protein, the possible binding modes are sampled for. Therefore, computational costbenefit methods are designed to find these molecular interactions: geometry-based algorithms (Fischer et al., 1993;Norel et al., 1994), fragment-based and docking incrementally (Rarey et al., 1996;DesJarlais et al., 1986), fragment-based methods for the de novo design (Miranker and Karplus, 1991;Eisen et al., 1994;Bohm, 1992), stochastic searches in the form of Monte Carlo or genetic algorithms (Goodsell et al., 1993;Hart and Read, 1992;Jones et al., 1997;Oshiro et al., 1995) and MD simulations (Cornell et al., 1995;Weiner et al., 1984;Brooks et al., 1983) are some of the strategies reported. MD simulations are used for further protein refinement before and after docking simulations (Cornell et al., 1995;Weiner et al., 1984;Brooks et al., 1983). Once a general strategy to simulate molecular recognition has been selected, a function capable of evaluating the mutual interaction between the ligand and the receptor will guide the evaluation of the proposed complexes. Scoring functions can be based on forcefields or on the knowledge of a type of interactions (Kitchen et al., 2004).
Some mathematical expressions of scoring functions can be obtained by taking into account the binding free energy values, which are calculated by the sum of the non-bonded interactions such as electrostatics and van der Waals forces and other restrictions such as bond angles. Along with the rest of the parameters for these equations, classical force-field-based scoring functions have been defined (Aqvist et al., 2002;Carlson and Jorgensen, 1995). Additionally, some software programs include extensions of force-fields that consider the hydrogen bonds, solvations and entropic contributions (Kuntz et al., 1982;Verdonk et al., 2003;Morris et al., 1998). Moreover, when the binding energy function is explicitly dependent on the hydrogen bonds, ionic interactions, the hydrophobic effect and the binding Science Publications AJABS entropy, the scoring function is catalogued as empirical, whereby a ponderation of each component in the function is employed (Bohm, 1998;Verkhivker et al., 2000;Gehlhaar et al., 1995). In contrast, knowledgebased scoring functions consider that the most favorable interaction presents the greatest frequency of occurrence; thus, a statistical analysis of the complex's crystal structure is performed, focusing on inter-atomic contacts and distances. For the screening of large databases and the modeling of particular interactions, as sulfuraromatic or cation-π, scoring functions are used (Muegge and Martin, 1999;Ishchenko and Shakhnovich, 2002;Feher et al., 2003). In addition, some physics-based scoring functions are used to assess the solvation effect (Kollman et al., 2000;Srinivasan et al., 1998;Still et al., 1990;Guimaraes and Mathiowetz, 2010). Another strategy termed consensus scoring, which unites the results of various schemes of docking, also offers virtual screening to improve the predictions of bound conformations (Charifson et al., 1999;Feher, 2006). Finally, because the ligand-binding process is cooperatively driven by enthalpic and entropic effects, the scoring functions should address the limited resolution of crystallographic models, the inherent flexibility permitted in the simulation, the binding dependency of conformational changes and the influence of water molecules. Thus, docking methods are an attractive area of study because of the complexity of the theory (Gorse and Gready, 1997;Kitchen et al., 2004;Mangoni et al., 1999;Hildebrandt et al., 2007;Roy and Mandal, 2008).

Combining Docking and MD Simulations
The chemical structure of a molecule is a critical factor that determines its physical properties and mechanisms of action. Although recent biological investigations have focused on discovering the functionality of biomacromolecules, SMs have demonstrated their vital role in regulating metabolism, biosynthesis and signaling in cellular networks. Moreover, SMs are widely used as therapeutic products because these drugs inhibit or activate enzymatic reactions (Chepelev and Dumontier, 2009).
The conformational variations of a molecule offer different options for how they interact with a biomacromolecule. The combinatorial chemistry of SMs has benefited from a collection of algorithms that produce several dispositions of the atoms in a molecule. For short polypeptide chains, the thermodynamics and kinetics are distinct from those in proteins because of the different dimensionality of their free energy hypersurfaces (Daura et al., 2002).
The role of a SM, as a counterpart of a larger molecule, is evaluated by algorithms that simulate this interaction. For instance, a protein-ligand docking program has two essential components, sampling and scoring. First, sampling is the process used to generate the putative ligand binding orientations/conformations close to the binding site of a protein; the ligand sampling and flexibility of the protein are taken into account. Second, the prediction of the binding tightness for each ligand pursuant to its orientations/conformations is termed scoring. In this calculation, the physical or empirical energy function is used, as mentioned above. Finally, after sampling and the scoring evaluation, the lowest energy score of the orientation/conformation is used to predict a binding mode (Huang and Zou, 2010).
A classification of strategies to obtain putative conformers has proposed four methods: a Key-Lock model (KL) using rigid-backbone docking, a Conformer Selection model (CS) using a novel ensemble docking algorithm, an Induced Fit model (IF) using energy gradient-based backbone minimization and a Combined Conformer Fit model (CS/IF) (Koshland, 1958). In a comparative study of these methods, a set of 21 complexes of unbound crystal structures were analyzed. The steps required to achieve this goal considered the backbone flexibility only for the smaller partner of the complex, an algorithm to generate the structural ensembles and a docking procedure by local perturbations around the complexed conformation. As a result, the lowest energy in the complexes showed more than 30% of native contacts for KL, CS, IF and CS/IF docking. Even when 15 targets using NMR ensembles of the smaller protein were studied, a similar result was obtained (Chaudhury and Gray, 2008).
Given that changes in backbone conformation influence the intra-and intermolecular energies of putative complexes, more information to distinguish near-native structures is needed. Thus, both backbone conformational sampling and discrimination of contacts should be considered together in a flexible docking. Koshland's Induced-Fit (IF) model affirms that a protein recognizes a ligand to form a complex because of the structural modifications in the binding site required for the interaction with the ligand. In other words, the conformation of a protein is the result of the presence of the partner in the complex and the backbone conformation should be obtained during molecular recognition by evaluating the local energetics of the interface (Koshland, 1958). Therefore, docking algorithms currently incorporate flexibility in the ligand and to a lesser extent, in the protein; the internal energy Science Publications AJABS of the protein continues to be evaluated as a parameter of interest and the analysis of the interaction is also restricted to selected residues (Alonso et al., 2006).
Variations in the interacting rules, for instance, decreasing the van der Waals repulsion to increase the atom-atom proximity between the receptor and ligand, result in flexibility in the receptor (Jiang and Kim, 1991;Gschwend et al., 1996). Furthermore, the use of rotamer libraries provides a set of side chain conformations with experimental origins that improve the sampling and prevent the minimization barriers, consequently conferring flexibility to the system (Leach, 1994). An ensemble of proteins is also considered as a conformational preselected set that provides more options for adjustability in the molecular recognition, although it contains rigid targets (Knegtel et al., 1997;Cavasotto and Abagyan, 2004). Changes in the protein backbone, giving rise to alternative conformations, can also be used to generate an average structure that maintains its most conserved features, which can be considered when the docking is subsequently performed. In this pretreatment of the set of coordinates, some loop movements that are involved reveal their participation (Alonso et al., 2006). Some methods focusing around the binding site have been developed that form an ensemblebased grid or employ precalculated two-body potentials to determine the interaction energy of the ligand. Particularly, ensemble-based grids reduce the effect of steric clashes in the interaction and therefore, lead to the selection of reliable conformations (Alonso et al., 2006).
The modeling of the ligand-receptor interaction should consider the mobility of both backbone and side chain flexibility. A refinement stage primarily includes the choice of atoms and the method by which these will result in an increase in mobility, with specific spatial restrictions (Andrusier et al., 2007). FireDock, a program devoted to refine and re-score rigid-body protein-protein docking, restricts the side chain flexibility exclusively in the clashing interface residues, smoothing the atomic radii of the partners. In this method, the scoring of the candidates is based on softened van der Waals interactions, Atomic Contact Energy (ACE), electrostatic and binding free energy calculations. In summary, algorithms dedicated to refining docking include at least three procedures: side chain prediction, rigid-body optimization and ranking of the candidates (Andrusier et al., 2007).
Other computational techniques that involve MD, energy minimization and gradient-based methods in Monte Carlo Minimization (MCM) have also been used to model backbone flexibility (Dominguez et al., 2003;Smith et al., 2005;Krol et al., 2007;Vries et al., 2007). Moreover, a high refinement of protein models is required because slight variations (1-2 Å) of atomic positions can regulate the formation of hydrogen bonds or steric clashes. Previous optimization of the participating molecules will prevent inappropriate docking (Katritch et al., 2012).
MD simulations have been shown to be helpful to study the dynamic behavior of the bound conformation of proteins complexed to synthetic ligands. For instance, Novak et al. (2009) used MD simulations to highlight the importance of protein flexibility in the specificity of Bcl-xL to bind different inhibitors. MD, as a very popular simulation approach, is often unable to overcome highenergy barriers within reasonable simulation time periods. Therefore, when both docking and MD are coupled to simulate the protein-ligand interaction, the spatial disposition of the molecules can only be a representation of an energetic local minimum. Different temperatures and alternative starting positions of the ligand are evaluated to address this inconvenience and even Monte Carlo-type algorithms are proposed to enhance the strategy (Kitchen et al., 2004;Joannis et al., 2011).
The use of MD for simulating the interactions of proteins with synthetic ligands as organic and peptide molecules is a common tool. MD simulations have enabled researchers to rationalize experimentally measured properties, to analyze the ligand-receptor interactions and to refine models of biomolecules determined by X-ray or NMR methods (Alonso et al., 2006). A properly constructed MD simulation, in which a solvated protein and its unbound ligand are subjected to physical laws, can yield a stable protein-ligand complex (Brooks et al., 2009). Using MD, it is possible to theoretically characterize how protein structure and stability are affected by an explicit solvent. Furthermore, time-averaged properties such as density, conductivity, dipolar moment, thermodynamic parameters, energies and entropies can be assessed by a sampling of the conformational space (Alonso et al., 2006). The solvent also impacts molecular recognition because water molecules can shield the interactions. Fortunately, at least a parametrization of the solvent is also included in the above mentioned force-fields (Brooks et al., 2009). Furthermore, when the effect of mutations in the receptor on protein-ligand interactions is evaluated, the inclusion of flexibility is essential because subtle variations in the content of atoms should be determined (Gorse and Gready, 1997;Kitchen et al., 2004;Mangoni et al., 1999). An analysis of the atomic behavior, induced fit Science Publications AJABS effects, the role of explicit solvent and stability of the complex over time can be achieved by MD simulations. Additionally, the final optimized structures can be utilized for calculating the binding-free energies (Alonso et al., 2006). A total failure, even if the ligand has been positioned into the putative binding pocket, occurs when the docking calculation results in unfavorable steric overlaps between ligand and receptor. Some homology models can result in imprecision that effect an appropriate interaction between molecules; therefore, a previous refinement of the model is always suggested. In some cases, the formation of the complex induces local conformational adjustments involving changes in the secondary structure; in this context, full flexibility of the receptor protein might lead to unrealistic complexes because of force field limitations and the omission of solvent molecules. These problems can be addressed by MD simulations, taking into account the limitations of the force-field and the requested level of resolution (Zacharias, 2004).
MD simulations methods can be combined with docking protocols to predict reliable protein-ligand complexes. The particularities of both techniques are complementary: the rigidity and driven strategy of some docking methods and the force-field-dependent flexibility of MD simulations can be combined for a common goal. The positioning of a ligand within a binding site is predicted by a docking calculation, thereby yielding the energy-dependent location and conformation. Once the ligand is in the most probable site, the MD simulation models the movements of the atoms involved in the interaction (Alonso et al., 2006).
The MM-PB/SA method combines Molecular Mechanics (MM) and the Poisson-Boltzmann/solventaccessible Surface Area (PB/SA) continuum solvent approaches to estimate binding energies; this MM-PB/SA function is also proposed as a post-docking filter during the virtual screening of compounds. It must be considered that the ligand affects the structure of the binding site and the dynamic equilibrium between distinct conformational states of the protein, providing information to identify the most likely complex ligand-protein (Alonso et al., 2006). A similar method that uses the Generalized Born equation (GB) to estimate the electrostatic contribution to the solvation free energy is known as MM-GB/SA and has been recently applied to SRC Kinase inhibitor potency prediction (Kohlmann et al., 2012).

CONCLUSION
The use of MD before and after docking is an appropriate way to study the conformational space of the protein-receptor complex (Fig. 1). Thus, MD simulations of the final docked structures in an aqueous environment can help in rationalizing the dynamics of molecular recognition. MD simulations are an attractive option for structural refinements of docked complexes, incorporating the freedom of both ligand and receptor, improving interactions and enhancing complementarity and thus accounting for the induced fit. Additionally, time-dependent evolution provides a dynamic picture of the complex and helps to distinguish the correctly docked conformations from the unstable ones (Alonso et al., 2006). Table 1 shows some examples where a molecular recognition technique has been combined with MD simulations.

Glycolitic enzymes as Drug Targets Analyzed by Combining MD Simulations and Docking Strategies
Glycolytic enzymes are targets for drug design, particularly against those organisms that depend mainly on glycolysis for ATP production. We present a study of Triosephosphate Isomerase (TIM).
This enzyme is an established target for drugs against various organisms, particularly protozoan parasites from the genera Giardia, Entamoeba, Plasmodium and Trypanosoma (Rodriguez and Rodrik, 2001;Joubert et al., 2001;Enriquez-Flores et al., 2008;Olivares-Illana et al., 2006;Ogungbe and Setzer, 2009). The enzyme is a 23 kDa protein that catalyzes the isomerization between glyceraldehyde 3-phosphate and dihydroxyacetone phosphate in glycolysis and gluconeogenesis, with a turnover number approaching the diffusion limit (Blacklow et al., 1998). The crystal structure of TIMs from several species have been solved by X-ray crystallography. This enzyme is a periodical arrangement of alternate α-helices and β-strands. The β-strands of this arrangement form an inner cylinder, whereas the αhelices form an outer cylinder. This topology is known as (β/α)8 barrel or TIM barrel (Banner et al., 1975) and approximately 10% of the proteins whose structure is known share this scaffold. TIM is only active as a dimer, even though each monomer possesses its own complete catalytic site. Exceptions are found in archaea and some thermophilic bacteria, in which TIM forms a tetramer (Kohlhoff et al., 1996). The TIM dimer is held by an extended interface in which loop 3 of one subunit aids in the arrangement of the positions of the catalytic amino acids of the other subunit, forming a favorable placement for catalysis. This implies that alterations in the intersubunit contacts of the dimer should bring about the abolition of catalysis (Tellez-Valencia et al., 2004).  (Cerqueira et al., 2009) Mohanty, 2012) The TM4/TM5 dimerization HADDOCK (Dominguez et al., 2003) CHARMM (Gorinski et al., 2012) interface of the serotonin 5-HT1A I-TASSER (Roy et al., 2010) Binding modes of flavonoid derivatives with VEGA ZZ (Pedretti et al., 2004) AMBER (Lu and Chong, 2012) the neuraminidase of the 2009 NAMD H1N1 influenza A Virus (Phillips et al., 2005) Autodock vina (Trott and Olson, 2010) Prediction of the human EP1 GROMACS (Lindahl et al., 2001) Gromos87 (Zare et al., 2011) receptor binding site HyperChem (Froimowitz, 1993 Derivatives of peptide epoxyketone and Surflex (Jain, 2003) AMBER (Liu et al., 2011) tyropeptin-boronic acid as inhibitors against the β5 subunit of human 20S proteasome Regarding drug design, the overall structure and the catalytic site are highly conserved among TIMs. Therefore, the efforts for drug design must focus on inhibitors that interact with non-catalytic residues. The interface represents a particularly attractive region to focus on because of the loss of activity upon subunit dissociation.

AJABS
The first studies mediated by docking and/or MD simulations were inspired by the interest in structurally describing the interaction of potential inhibitors found by in vitro analysis. In 2004, a low-molecular-weight compound, 3-(2-benzothiazolylthio)-1-propanesulfonic acid (compound 8), was found to bind to the dimer interface of the triosephosphate isomerase from Trypanosoma cruzi (TcTIM) and to abolish its function with a high level of selectivity (Tellez-Valencia et al., 2004). In this study, it was hypothesized that compound 8 would likely bind Cys 15, a conserved residue within several parasites (Trypanosoma brucei (Garza-Ramos et al., 1996), Leishmania Mexicana (Garza-Ramos et al., 1998, Plasmodium falciparum (Maithal et al., 2002), Entamoeba histolytica (Rodriguez-Romero et al., 2002) but not in Homo sapiens (Tellez-Valencia et al., 2004). To further describe the possible binding sites of benzothiazoles at the interface of tripanosomal TIM, fully flexible benzothiazoles were docked onto the dimer interface.
It was found that dimer disruption did not occur via Cys 15 but instead through the unstabilization of π-π interactions of two aromatic clusters present at the interface (Espinoza-Fonseca and Trujillo-Ferrara, 2004). Later, the same research group presented the docking of seven benzothiazoles into the interface of both human and trypanosomal triosephosphate isomerases using the program AutoDock. Structural and energetic analysis of the complexes showed that large benzothiazoles could form more stable complexes with the trypanosomal triosephosphate isomerase than with the human triosephosphate isomerase (Espinoza-Fonseca and Trujillo-Ferrara, 2005).
From this study, it was concluded that the distribution of the residues forming the aromatic clusters at the enzyme's interface as well as the size of the inhibitors ligands may play crucial roles in the selective inhibition of TcTIM. This information was largely improved when the same research group performed a series of combined docking/molecular dynamics simulations to determine the factors that play a role in the selectivity of certain benzothiazoles over parasite TIMs. The interaction of the compound 6,6'-bisbenzothiazole-2,2' diamine depicted V7 with TIMs from Trypanosoma cruzi, Trypanosoma brucei, Entamoeba histolytica, Plasmodium falciparum, yeast and humans was analyzed. It was found that different accessibilities of the protein's interface of TIMs are a key determinant of the inhibitory activity of benzothiazoles on the enzyme. It was found that V7 directly interacted with both aromatic clusters located at the interface of the TIM from T. cruzi. These aromatic clusters are formed by Phe75 from one monomer and Tyr102 and Tyr103 from the adjacent monomer.
Similarly, V7 had direct contact with Tyr101 and Tyr102, which, together with Phe74 from the adjacent monomer, constitute the aromatic clusters of TIM from T. brucei. In contrast, it was found that V7 does not interact very tightly with TIMs from E. histolytica, P. falciparum, yeast and human TIMs due to the reduced accessible surfaces of interfaces and to the packing of the aromatic clusters that did not allow for the formation of a Science Publications AJABS well-defined binding site for V7 (Espinoza-Fonseca and Trujillo-Ferrara, 2006).
Undoubtedly, combining docking and molecular dynamics simulation strategies has provided a complete view on the dynamics of the complexes between antitrypanosomatid agents and TIMs from different species. This improves our understanding of how parasite TIMs could be effectively and specifically inhibited, leading to better rational drug design.

Identification of Neuropharmacological Targets and their Importance
The activities of single neurons, neural networks and neural centers have dynamic behavior that can be addressed using dynamical systems theory (Nowacki et al., 2012;Ghorbani et al., 2012;Serletis et al., 2011). The analysis of the spatiotemporal characteristics of brain cells offers a perspective in which even biomolecules as proteins can be incorporated to obtain a multivariable mathematical model. Computational neuroscience and computational methods in neuropharmacology address how to distinguish normal information processing from pathological information processing, with the purpose of finding therapeutic alternatives for neurological and psychiatric disorders. The combined use of the above disciplines is called computational neuropharmacology (Aradi and Erdi, 2006). In an attempt to expand the scope of this interdisciplinary approach to neuropharmacology, diverse computational methods such as cheminformatics and bioinformatics can be included, with the aim of improving drug design. Thus, information on the three-dimensional structure of the target macromolecule and its binding molecules, to model the receptor-ligand interactions, can be enhanced with computational simulations of brain signals produced by the information of integrative physiology of neurons. These latter signals are compared with electrophysiological measures (Veselovsky and Ivanov, 2003;Schneider and Fechner, 2005).
Approximately 140 types of voltage-gated channels and even more ligand-gated channels have been identified (Yu et al., 2005;Novere and Changeux, 2001). These channels are tissue-specific and they are associated with different phases of development. Furthermore, diverse diseases, including epilepsy, cystic fibrosis and some forms of diabetes, are related to dysfunctional channels (Kass, 2005).
When a compound binds a regulatory site of a neural receptor to regulate the efficiency of the binding, a socalled allosteric modulator enhances or suppresses the electrical signal and consequently alters the synaptic transmission. The description of the mechanism of this modulation is still insufficient; however, the computational methods that integrate conductance-based techniques promise a conceptually new perspective for the computational design of drugs. Searching for novel target-specific drugs can be aided by the simulations and experimental results of spatiotemporal neural activity patterns, whereby pathological and normal dynamical states can be identified by a previous calibration of the integrated system (Aradi and Erdi, 2006). Meanwhile, a wider integration system of information is still improved and molecular recognition and molecular dynamics simulations are the most widely used techniques to design drugs to regulate the activity of the nervous system.
An efficient presynaptic transport is necessary to upload the neurotransmitters into small vesicles at the axon terminals. In particular, the amino acid glutamate, the main excitatory neurotransmitter, exhibits an uploading system driven by Vesicular Glutamate Transporters (VGLUTs) (Takamori, 2006). To study the stability of the human VGLUT1 protein, a structural model was built based on a bacterial homologue, the glycerol-3-phosphate transporter GlpT from Escherichia coli. This model was analyzed by docking and molecular dynamics techniques. Furthermore, the latter was simulated into a lipid bilayer (Almqvist et al., 2007). This simulation confirmed that the VGLUT1 model stably maintains all transmembrane helices and that these structures display the lowest RMSD fluctuations over the simulated 10 ns. Furthermore, to draw a conclusion about the orientation of the amino acids embedded into the membrane, docking studies with the VGLUT substrates such as L-glutamate and inorganic phosphate were performed.

Concluding Remarks
Drug design must consider the intrinsic flexibility of the proteins when any strategy devoted to model the ligand-receptor interaction is implemented. Moreover, understanding the flexibility of biomolecules in diverse contexts is an area of interest. Pharmaceutical targets can be favored by backbone and side chain flexibility increasing or decreasing the ligand binding. A future challenge is developing conjugate algorithms such as docking and MD as well as hybrid schemes, where a combination of efficient modeling and computational cost could guide an exhaustive search of molecular interactions with pharmacological applications.