Molecular Modeling and Simulation of Membrane Transport Proteins

Membranes fulfill the essential need of all living species to separate different compartments. On the other hand, in a cell the homeostatic environment can only be maintained by the cellular membrane acting as a selective ‘filter’, which allows the cell to continuously communicate with other cells. Mechanisms which facilitate the translocation of materials across the membrane regulate the entrance and disposal of ions, amino acids, nutrients, and signaling molecules.

This selective transport across cellular membranes is carried out by two broad classes of specialized proteins, which are associated with or embedded in those lipid bilayers: channels and transmembrane transporters. They work by different mechanisms: Whereas channels catalyze the passage of ions (or water and gas in the case of the aquaporin channel) (Agre, 2006) across the membrane through a watery pore spanning the membraneembedded protein, transporters are working via a cycle of conformational changes that expose substrate-binding sites alternately to the two sides of the membrane (Theobald & Miller, 2010).
If we regard the force that drives the transport process there is also a huge difference in the way ion channels and transporters act. Channels assist a downhill movement along a concentration gradient (passive diffusion), whereas in transporters it is usually directed against a concentration gradient of the substrate. Thus, in order to comply with their business, transporters are dependent on another source of the cellular energy. Secondary active transporters rely on ionic gradients. In the case of primary active transporters ATP is the driving force (Wang et al., 2010).
A comprehensive list of all annotated transport proteins is freely available online on the TCDB website (http://www.tcdb.org). This Transporter Classification Database uses an International Union of Biochemistry and Molecular Biology (IUBMB) approved system of nomenclature for transport protein classification. The TC system is analogous to the Enzyme Commission (EC) system for classification of enzymes, except that it incorporates both functional and phylogenetic information (Saier et al., 2006;Saier et al., 2009).
According to the TCDB system Membrane Transport proteins can be classified as follows (List of families and subfamilies of the TC system): www.intechopen.com 1. Pores and channels a. Helical channels b. Strand porins c. Pore-forming toxins d. Non-ribosomally synthesized channels e. Holins 2. Electrochemical-potential-driven transporters a. Transporters or carriers (uniporters, symporters and antiporters) b. Non-ribosomally synthesized transporters 3. Primary active transporters a. P-P-bond-hydrolysis-driven transporters b. Decarboxylation-driven transporters c. Methyl-transfer-driven transporters d. Oxidoreduction-driven transporters e. Light-driven transporters 4. Group translocators 5. Transmembrane electron carriers 6. Accessory factors involved in transport 7. Incompletely characterized transport systems Our special interest focuses on transmembrane transport proteins ('transporters'). Excellent manuscripts on membrane channels have been provided by other groups (Gumbart et al., 2005;Schmidt et al., 2006).
To date more than 400 membrane proteins have been annotated in the human genome. Two major superfamilies -which are also intensively investigated in our group -are the ATP-binding cassette transporters (ABC, e.g. P-glycoprotein), and the solute carrier superfamily (SLC). They will be discussed into more detail in section 3 of this chapter.

Transporters as pharmacological targets
Transport proteins are playing important roles in the whole drug discovery and development process. They regulate absorption, distribution, and excretion of drugs and therefore influence drug disposition, therapeutic efficacy and adverse drug reactions in the human body. This has to be taken into account in pharmacological studies (Giacomini et al., 2010).
It is estimated that transporters account for about 50% of drug targets. However, their modes of (selective) transport are only poorly understood. This is due to difficulties in membrane protein purification, expression, and crystallization (Caffrey, 2003), which is still in its childhood. As a consequence there exists a striking disproportion between the number of entries of resolved structures of soluble proteins and membrane proteins in the protein databank (PDB, http://www.pdb.org/). To date, only about 2% (1462 by Sept 2011) of the structures are from transmembrane proteins (75594 structures in total). Out of these there are only 302 unique structures (proteins of same type but from different species are included) (Irvine). Moreover, a significant number of the membrane protein structures determined are at relatively low resolutions (Lindahl & Sansom, 2008).
However, there are tremendous efforts as to ameliorate the methods in order to obtain atomic resolution structures of membrane protein molecules (Newby et al., 2009). Thus, for www.intechopen.com the past two decades, the number of available structures of membrane proteins has been climbing the exponential foot of the growth curve, since it is doubling every three years (Theobald & Miller, 2010).
For a medicinal chemist, the availability of a growing number of structures, paves the way for further in silico studies. Some are very promising in a way that they can be used as templates in order to build up a homology model of the membrane protein of interest. Those models, which may further be applied for Molecular Dynamics (MD) Simulation and Docking Studies, give us the opportunity to gather new insights into the molecular structure and function of the protein under investigation and the behaviour of certain ligands in the binding site.
As translocation always involves a dynamic process, which cannot easily be studied by mere experimental techniques, above all, the application of long-term MD simulations should be implemented into the whole process of drug discovery and development. Due to significant increase in computational power and improvements in parallelization techniques, nowadays simulations of membrane transport proteins may stretch up to microseconds -that is, to physiologically relevant time scales. In this review we are describing the theory and methodology related to computational techniques used in the modeling of transporters and we will outline the recent developments in the field of ABC transporters and neurotransmitter transporters.

Basic concepts
Despite the enormous increase of published structural data for proteins, the particular availability for a protein of interest can vary from the sheer presence of the amino acid sequence to a multitude of high-resolution X-ray structures. Fortunately, Mother Nature was not too generous in providing unique structural folds for functionally related proteins, as the structural arrangement within a family of homologous proteins is much higher conserved than the respective amino acid sequences (Lesk & Chothia, 1980). Thus, in many cases the combination of a primary sequence on the one hand and one or more reasonably well-resolved homologue structures on the other hand can result in homology models surprisingly well representing the molecular reality, paving the way to successful comparative modeling studies. The process of predicting the 3D structure of a protein can be achieved by four main steps: fold assignment, target-template sequence alignment, building and evaluation of the models (Cavasotto, 2011).

Assignment of the basic fold and sequence alignment
The first step towards a good model is the identification and careful selection of structurally related template proteins. Although generally a high percentage of global sequence identity is a good indicator for the model quality to be expected, it must be kept in mind that the identity in the area of interest, i.e. the binding site(s), can differ significantly from overall values. One should strive to put the main focus during template selection on facilitating a maximum of accuracy in modeling those vital regions.
Maybe the best source for structural templates is the Protein Data Bank (Berman et al., 2000). As mentioned earlier, it offers the coordinates of structurally resolved proteins, including large amounts of surplus information like primary sequence, experimental settings or cobound ligands and ions. Search tools like BLAST and FASTA (Altschul et al., 1997;Pearson, 1990) usually do reasonably well in identifying the correct fold of a protein. The second, even trickier step is the subsequent alignment step of target and template sequence, as it is possible with T-Coffee or CLUSTAL W (Notredame et al., 2000;Thompson et al., 1994). Conserved residues and regions of experimentally determined proximity need to be aligned as accurately as possible. Multiple alignments of sequences belonging to the same gene family can significantly enhance the performance of the search for conserved residues or even regions, but thorough literature search and manual adjustment of the alignment are inevitable in order to achieve best results.
A good example for the importance of taking a comprehensive look at the research subject is the meanwhile annual GPCR Dock competition (Kufareva et al., 2011), where prior to the release of newly resolved G-Protein Coupled Receptors (GPCR) 3D structures, modeling groups get the chance to submit their best efforts of predicting the correct receptor and ligand conformation. It turned out, that advanced modeling tools and human intervention contributed about equally to the success of the individual approaches.

Building and refinement of the models
Once it is assumed that the alignment meets all available experimental data, it can be started to calculate coordinates for the target residues. Although automated homology modeling methods exist, the yielded models tend to lack accuracy, especially in cases of low sequence identity (Dalton & Jackson, 2007).
Usually, the crude model is built by aligning the basic backbone framework, the so-called structurally conserved regions (SCRs). Conserved secondary structural elements like αhelices or β-sheets are inherited, being responsible for the general shape of the model. Several homology modeling approaches also try to include information about known ligands into the binding site construction in order to meet its particular geometry (Evers & Klebe, 2004;Sherman et al., 2006). Subsequently, assignment of the side chain conformations needs to be done according to steric and energetic constraints. Identical residues usually can be considered to be oriented similarly, likewise highly similar amino acids. For non related residues, rotamer libraries can provide initial geometrical guesses (Schrauber et al., 1993), although other effects like packing energies may lead to significant deviations. Up to 30% of side chain conformations in X-ray structures do not correspond to usual rotamers, yet adopting energetically allowed conformations. Naturally, selecting the most probable side chain orientation solely according to statistical criteria is problematic, so methods including structural features of the local environment have been developed (Deane & Blundell, 2001). Still, some cases require manual adjustment, for instance the incorporation of known disulfide bridges, specific internal hydrogen bonds or ion binding pockets.
The major challenge in comparative modeling is the treatment of structurally variable regions (SVRs). Especially flexible loop regions lacking a structural template are difficult to predict, since the calculation time increases nearly exponentially with the degrees of www.intechopen.com freedom added by every flexible residue. There are several strategies to meet this issue. Knowledge-based strategies try to find structural guesses by automated database search for related sequence sections in other proteins, possibly not even close to being genetically similar. From a computational point of view, conformational searches by ab initio calculation of the desired region are more costly, but recent approaches yielded reasonably good results for a loop length up to 17 residues (Mehler et al., 2002;Zhao et al., 2011). For significantly longer loops, in case that the problematic region is remote from the actual binding site(s) and not considered being directly linked to the binding process, it can be viable to leave it up to the standard modeling software how to build the respective flexible region, and hope for subsequent MD simulations to find a near-to-native conformation (Amaro & Li, 2010).
The model building process bears numerous sources of unfavorable steric strain energies, calling for an appropriate minimization procedure. As one can imagine, this step has to be carefully balanced in order to overcome steric clashes without compromising tediously elaborated side chain orientations or, even worse, entire conserved regions. Instead of global optimization attempts good minimization protocols start with local treatment of clashes with initially fixed backbone atoms. Thus, solvent molecules, ions and hydrogen atoms possibly responsible for large initial forces can adopt energetically more favorable positions. Gradually, initial tethering forces are reduced, avoiding artificial distortions (Höltje et al., 2008). This can be facilitated by molecular mechanics calculations using different force fields like CHARMM , OPLS (Jorgensen & Rives, 1988) or AMBER (Weiner et al., 1984). In contrast to force fields for small molecules, they have to handle huge systems, therefore being usually somehow simplified regarding the treatment of long-distance nonbonding interactions or non-polar hydrogen atoms, called united-atom models.
Energy minimized models may be further optimized by molecular dynamics (MD) simulations, as reported in-depth in section 2.3.

Model evaluation
Predictions including as many degrees of freedom as homology models desperately need reliable tools to estimate their quality, as the accuracy of a 3D model is responsible for the amount of information that can be gained by it (Marti-Renom et al., 2000). Several evaluation programs exist; many of them are available online on server-basis. The SWISS-MODEL (Arnold et al., 2006) (http://swissmodel.expasy.org) and the SAVES server (http://nihserver.mbi.ucla.edu/SAVES), for instance, offer a variety of local and global quality estimation tools (Benkert et al., 2011;Hutchinson & Thornton, 1996;Zhou & Zhou, 2002). It is important to look at both, as they are not necessarily mutually related.
A comprehensive stereochemistry check can be carried out using the Procheck suite (Laskowski et al., 1993), searching for geometrically unusual residues in a given protein structure by comparison with stereochemical parameters of high-quality benchmark structures. Likewise, the What_Check module of the WHATIF package and the VADAR web server do a similar job (Vriend, 1990;Willard et al., 2003). The features examined include the planarity of aromatic ring system and peptide bonds, bond lengths and basic checks like C α -chirality.
Once the yielded quality statistics are acceptable within the limitations of the possible, a model can be considered ready for the further use in, for instance, docking studies. www.intechopen.com

Molecular docking
Molecular Docking is a versatile tool in structure based drug design. This technique is able to predict possible orientations of one molecule to another. In this section we will focus on protein-ligand docking, describing the interaction of a small molecule in a binding pocket of the protein of interest.
In principle molecular docking is comprised of three consecutive steps: i) the definition of the binding site, ii) the placement of the ligand inside the defined site, and iii) the ensuing evaluation of this placement, called scoring.

Binding site identification
The right determination of the binding site of the ligand is essential for the subsequent docking process. If the active site is not known there are several algorithms that are able to detect potential binding pockets (extensively reviewed in (Henrich et al., 2010)). These programs scan the protein surface for cavities that fulfill certain geometrical constraints, which mark them as possible ligand binding sites. While the program LIGSITE uses a grid for the surface scan (Hendlich et al., 1997), the PASS algorithm utilizes layers of spheres that should describe buried cavities (Stouten & Brady, 2000).
In order to consider also physico-chemical criteria in binding site detection the surface of the protein can be scanned with fragments of ligands with subsequent calculation of their complementarity. Another approach to get an idea about potential active sites can be achieved by comparing the query protein with homologues, as proteins with related function share similar binding sites. For this purpose the program CAVBASE (Kuhn et al., 2007), which relies on the LIGSITE algorithm, has been developed.
As soon as the binding site is known, it has to be characterized in order to get information about specific binding possibilities through non-covalent interactions. To detect these hot spots atom probes, ligand fragments or whole small molecules are positioned inside the binding pocket. The program GRID is able to detect interactions and solvation effects by calculating the interaction energy between grid points of the binding pocket and certain atom probes (Reynolds et al., 1989). On the other hand, the multiple copy simultaneous search (MCSS) method places thousands of probe copies inside the pocket. After energy minimization those probes cluster at certain local minima defining the hot spots (Caflisch et al., 1993).

Search algorithms
The role of the search algorithm is the correct placement of the ligand in the binding pocket. Ideally it should therefore consider all possible degrees of freedom, which leads to higher accuracy. However, due to limitations regarding computer power, this is penalized in favor of higher speed by reducing the number of the degrees of freedom (Sousa et al., 2006).
Although in protein-protein docking the rigid-body approximation is still applied (Kuntz et al., 1982), in protein-ligand docking the small molecule is treated flexible.
Approaches that try to explore all degrees of freedom of the ligand systematically comprise conformational search methods, fragmentation methods or database methods. By applying www.intechopen.com conformational search methods, every rotatable bond of the ligand is rotated in fixed increments. As this can easily lead to a combinatorial explosion, this technique can only be applied for small or rigid ligands. More prevalently used are the so-called fragmentation methods that place fragments of the ligand in the binding pocket, which are subsequently fused. Depending on the fragmentation and placing of the ligand place and join and incremental approaches can be distinguished. Popular docking programs utilizing this type of search algorithms include LUDI (Bohm, 1992), FlexX (Rarey et al., 1996), DOCK (Ewing et al., 2001), ADAM (Mizutani et al., 1994) or Hammerhead (Welch et al., 1996).
A computationally efficient way to search for possible orientations forms the database method. For this protocol a conformational library of the ligand is prepared which is docked rigidly into the binding site.
Besides systematic search algorithms there are programs that prefer stochastic principles for binding mode prediction. At this the flexibility of the ligand is provided by random conformational changes that are either kept or rejected on basis of a direct evaluation of the conformation. Among others genetic algorithms present a convenient tool for this purpose. With this optimizing procedure a random population of possible ligand poses is generated, where the characteristics (degrees of freedom) of each are stored in its genetic code (chromosome). By applying genetic operations, like cross-over or mutation, new poses are generated and subsequently scored. Depending on this fitness score the pose is either rejected or it replaces the least fit member of the population. This procedure is c o n d u c t e d o v e r t h o u s a n d s o f c y c l e s w h i c h e n d s u p i n h i g h l y o p t i m i z e d l i g a n d orientations. This protocol is included in the popular docking programs GOLD (Jones et al., 1995;Verdonk et al., 2003) and AUTODOCK (Goodsell et al., 1996;Goodsell & Olson, 1990;Olson et al., 1998).
Another possibility to consider ligand flexibility is presented by molecular dynamics simulation of the ligand in the binding pocket. However, this is mainly used in combination with other search algorithms (Kitchen et al., 2004).
In the last years not only the flexibility of the ligand but also protein movements due to ligand binding gained more and more importance (B-Rao et al., 2009;Cozzini et al., 2008). Although it is known that some proteins undergo large structural changes, even domain rearrangements, upon ligand binding, by now it is not possible to cover that in reasonable time and effort. However, since docking a ligand into the right conformation of the binding site is extremely important for the quality of the resulting orientations, efficient workarounds have been developed. Soft docking is one possibility to account small movements of the protein side chains during docking (Jiang & Kim, 1991). For this technique soft potentials are applied on certain side chain atoms in the binding pocket, which therefore tolerate overlap with ligand atoms. The merit of this technique is the easy implementation, as only scoring parameters have to be adapted. On the other hand only small changes can be considered and there might be a bias towards the starting structure.
With the help of rotamer libraries, movements of side chains are included in the search algorithm (Leach, 1994). Depending on the size of the library this method calls on moderate computational power and is able to adapt to certain ligand conformations. Nevertheless, as the backbone is kept rigid large structural movements cannot be covered.
Docking into multiple protein structures (MPS) is therefore highly appreciated as they allow flexibility of the protein during the docking process. Different protein conformations (X-ray structures or taken from MD simulations) are selected and multiple docking runs are performed. As this approach is extremely costly, the more efficient method of ensemble docking should be used preferentially. Therefore an average receptor grid is generated and used for docking (Knegtel et al., 1997).
A hybrid technique that is commonly used to encounter protein flexibility is the induced fit docking protocol of the Schrödinger Suite (Sherman et al., 2006). This method turns major attention on the ligand-induced conformational changes of the protein residues surrounding the binding site. Therefore, the ligand is docked into the rigid binding pocket, amino acid residues that are within a certain radius of the resulting poses are removed and rebuilt using the Schrödinger homology modeling program Prime. After energy minimization of the complex the ligand is redocked into the modified binding pocket.

Scoring functions
The application of a scoring function is important to assess the quality of ligand orientations in the binding pocket that resulted from docking experiments. Basically there are three areas of use for scoring functions. In order to understand the interaction between a defined molecule and the target protein the scoring function needs to be able to identify the true pose among the plethora of orientations, generated by the search algorithm. For lead optimization in particular a scoring function should correctly determine the affinity between the ligand and the protein. However, for virtual screening of large compound databases scoring should provide correct ranking. As there are still limitations regarding computer power, the right balance between accuracy and speed has to be chosen, which is strongly dependent on the field of application (reviewed in ).
Force field based scoring functions use terms that describe the free energy of binding for evaluating binding poses. In that regard bond stretching, angle bending and dihedral angle forces for the ligand, but also non-bonded VDW and electrostatic interactions with the protein are calculated . Furthermore the accuracy of these methods depends on their treatment of the solvent. More accurate techniques, like thermodynamic integration or free energy pertubation, treat water molecules explicitly. As these are the most expensive affinity prediction methods, more simplified and computationally less expensive versions are linear interaction energy (LIE) models, where two additional empirical parameters can be used to reduce the number of simulations needed. On the other hand, MM/PBSA and MM/GBSA methods gain speed by using implicit solvent models.
However all of these methods are still not applicable for virtual screening as they are computationally too expensive.
Empirical scoring functions are therefore a fast alternative. They assess the quality of binding by a number of weighted terms that are derived by fitting data of complexes to known affinities (Bohm, 1994;Bohm, 1998). Numerous commonly used scoring functions belong to this group, including ChemScore (Eldridge et al., 1997) and X-Score (Wang et al., 2002). Nevertheless, a disadvantage of this method would be the dependence on the training set, as complexes with binding affinity are essential.
Thus, knowledge-based scoring functions may be preferred in this regard. These scoring functions make use of the statistical occurrence of protein-ligand interactions of complex databases. In contrast to empirical functions they do not aim at reproducing bindingaffinities, but experimentally determined structures, wherefore a much larger training set can be used (Tanaka & Scheraga, 1976). Representatives of this group of scoring functions are among others ITScore  and DrugScore (Gohlke et al., 2000). A further development of the ITScore by Zou et al. ITScore/SE managed to include solvation and entropic effects into the scoring function (Huang & Zou, 2008), which lead to a strong increase in scoring accuracy.
As the choice of the scoring function strongly depends on the research query, the combination of several functions, so-called consensus scoring, has been suggested (Charifson et al., 1999).

Molecular dynamics
It is obvious that the mechanism of action by which certain nutrients or drugs are translocated by a transporter implicates the protein to be flexible. In order to be able to allow for a sufficient comprehension of the dynamics of the transport protein, we can not only rely on experimental techniques. In addition, biomolecular simulations can provide a detailed description of particles in motion as a function of time. Thus, they are an important tool for understanding the physical basis of the structure and function of proteins, and biological macromolecules in general. However, experimental validation should always serve to test the accuracy of the calculated results and also to provide a basis for improving the methodology (Karplus & McCammon, 2002).
It is almost 35 years ago, since for the first time McCommon, Gelin and Karplus have studied the dynamics of the pancreatic trypsin inhibitor by solving the equations of motion for the atoms with an empirical potential energy function (McCammon et al., 1977). In this very beginning of Molecular Dynamics simulations, the calculations were still restricted to the picosecond timescale. However, according to Moore's Law computer power is doubling approximately every two years (Moore, 1965). Thus, MD simulations of biomolecules now are able to stretch up to microseconds. For the study of biological relevant phenomena like enzyme catalysis or even protein folding, MD has become a standard tool -always complementary to experimental techniques.

Theory, fields of application, strengths and limitations of MD simulations
By integrating the Newtonian Equations of Motion, Molecular Dynamics simulations are able to describe the behavior of particles in a certain system within the observed period of time. The interaction of the atoms is described by the potential energy function of the given force field [e.g. Amber (Cornell et al., 1995), CHARMM , GROMOS (Scott et al., 1999), OPLS (Jorgensen & Rives, 1988)]. Nowadays, there is an ongoing effort to ameliorate these parameters in a need for models being as less artificial as possible.
The field of application of biomolecular simulations is manifold. It reaches from validation and optimization of previously built homology models, refinement of crystal structures, to the prediction of protein-ligand, and protein-protein interactions, to the study of functional properties of biological systems at the atomic level (e.g. protein-folding, destabilization or structural change of a protein upon mutation), to even de novo design of proteins (Park et al., 2005).
Despite obvious drawbacks of classical MD simulations, which include limitations in time scales that can be studied but also certain inaccuracies of the force field, for instance with respect to polarization effects, the ability of bringing molecular structures alive also allows the researcher to sample the conformational space. This is especially interesting in liganddocking applications. On one hand, various extracted snapshots from a previously MDequilibrated protein-ligand complex may serve in order to perform an ensemble docking which is said to outperform docking into only one sample structure (Knegtel et al., 1997). Secondly, MD may also be used in order to refine certain poses and study the ligand-protein interactions on a molecular basis as a function of time.

Simulations in a membrane
The setup of a simulation system, which includes a protein embedded into a lipid bilayer requires additional efforts in comparison to a system with a soluble protein. There are different choices the researcher has to make regarding to the nature of the phospholipid bilayer used, the temperature at which the simulations should be performed (this also depends on the nature of the bilayer), the force field, the water model (e.g. SPC, SPC/E, TIP3P, TIP4P, TIP5P; this also depends on the choice of the force field), and many more.
One of the most challenging parts is the correct parameterization of the ligands. According to the force field that has been selected there are diverging approaches, ranging from a pure manual assignment of partial charges and force constants, to the use of scientific software like Gaussian (www.gaussian.com), to an automated procedure by taking advantage of platforms such as the Automated Topology Builder and Repository (Malde et al., 2011). However, it has to be stated clearly that a manual inspection and refinement of suchlike obtained topologies will always be needed.
Membrane proteins should be placed in a bilayer which is as similar as possible to its native environment. There is a diverse spectrum of phosphlipid bilayers available -differing mainly in the charges of their polar head groups, lengths and saturation of their acyl chains. If lipids play key roles in the proteins function, different combinations of lipids will probably better represent the in vivo conditions. It should always be kept in mind that in order to simulate the membrane in a liquid-crystalline state the temperature of the simulation needs to be above the melting temperature of the chosen lipids (phase-transition temperature).
The protocols for setting up MD simulations of membrane proteins are manifold. In any case, however, one needs a pre-equilibrated bilayer, which can be retrieved from different groups around the world (e.g. Peter Tieleman, Scott Feller, Helmut Heller, Mikko Karttunen) or an individual bilayer may be generated and equilibrated with regard to the respective size and nature of the protein to be studied.
When it comes to the insertion of the respective protein into the pre-equilibrated bilayer, again no standard protocol has been established up to now. In any case, it is of utmost importance to obtain a system with a tightly packed bilayer around the protein, so that the www.intechopen.com consecutive equilibration time for the membrane can be kept quite short. Protocols like inflategro (Kandt et al., July 2009;Kandt et al., 2007) or g_membed (Wolf et al., 2010) seem most suitable. Whereas, inflategro inflates the lipid bilayer, insert the protein and then deflate the lipid bilayer again, g_membed does it the other way around. It grows a protein into an already hydrated and equilibrated lipid bilayer during a short MD simulation. A special case of insertion procedure certainly is the use of coarse-grained simulations. Here the lipids are able to self-assemble around the protein. However, as this type of simulations use a very simplified description of interactions, for a lot of investigations the relevant information might not be captured.
An idea of a general protocol for the set up of a MD simulation can be found here: 1. Choose a force field for which you have parameters for the protein and lipids. 2. Insert the protein into the membrane. 3. Solvate the system and add ions to neutralize excess charges and adjust the final ion concentration 4. Energy minimize. 5. Let the membrane adjust to the protein. Typically run MD for ~5-10ns with restraints on all protein heavy atoms. 6. Equilibrate without restraints (gradually release the protein). 7. Run production MD. 8. Analysis.
As seen from this overview, after the insertion of the membrane protein it is inevitable to properly equilibrate the lipid bilayer again. This is done by restraining the protein (plus eventually existing ligands conserved water molecules, ions, and cofactors) during a MD run where the membrane is able to adjust to the protein. Subsequently, the whole system has to undergo an extensive equilibration procedure. The end point of the equilibration phase and simultaneous starting point for the MD production run can be determined mainly by evaluation of system parameters (e.g. total energy, temperature) and parameters concerning the protein (e.g. backbone root mean square deviation). A production run for membrane proteins typically resides somewhere in between 50 ns and hundreds of nanoseconds.

Enhanced sampling techniques
As already mentioned in chapter 2.3.1, classical MD simulations are confronted with their limitations in time scales. The limiting factor is the maximum timestep that can be used for the integration, determined by the fastest motion in the system (e.g. bond vibrations). Thus, it is not able to study 'slow' biological processes without taking advantage of enhanced sampling techniques. This includes of course always a method, which works at the expense of fidelity.
A s o u t l i n e d i n a n e x c e l l e n t r e v i e w o f Christen and van Gunsteren (Christen & van Gunsteren, 2008) we have to distinguish three different types of search and sampling enhancement techniques: deformation or smoothening of the potential energy surface (e.g. Coarse-graining the model by reducation of the number of interaction sites), scaling of system parameters (e.g. simulated temperature annealing), and multi-copy searching and sampling (e.g. replica-exchange algorithm).

www.intechopen.com
If we want to study membrane proteins and especially their interactions with ligands sampling along transition pathways will be needed. Such pathways are often characterized by high-energy barriers separating meta-stable states along the ligand/substrate transition. Here, methods like pulling or steered MD (SMD) and targeted MD (TMD) may be used in order to drive the sampling to a specific direction. In the SMD approach external forces are applied on certain atoms in order to accelerate processes that are otherwise too slow to model (Isralewitz et al., 2001). Subsequently, the potential of mean force required to induce the transition can be used to estimate free energy barriers. This method is well established and has been used in many applications (Isralewitz et al., 2001;Lu & Schulten, 1999).
In the TMD method, the reaction coordinate is defined by a single mass-weighted root mean-square 'target distance' between a known initial structure and a fixed final (target) structure. By gradually reducing the constrained target distance to zero, the system is driven from the reactant to product state without explicitly defining the reaction pathway (Schlitter et al., 1994).

ABC Transporters and multidrug resistance
ABC (ATP binding cassette) transporters are ubiquitous proteins that are expressed by prokaryotic and eukaryotic organisms. About 50 human ABC transporters are known, which are divided into seven different subfamilies, designated A-G.
Depending on ABC subfamily substrates include among others drugs, lipids, bile salts, peptides, ions and amino acids. Additionally some ABC proteins are known for transporting a broad variety of chemically diverse molecules, which are therefore referred as multidrug transporters. Besides their physiologically important protecting function of exporting xenotoxins, these efflux pumps affect pharmacokinetic profiles of many drugs. Furthermore the acquisition of multidrug resistance (MDR) can often be traced back to elevated expression of multidrug transporters in the affected cells.
The three ABC transporters mostly associated with MDR are P-glycoprotein (P-gp, ABCB1), multidrug resistance protein 1 (MRP1, ABCC1) and breast cancer resistance protein (BCRP, ABCG2). P-gp is encoded by the mdr1 gene and is expressed in epithelial cells of the blood brain barrier, liver, kidney and intestine, where it is located at the apical side of the membrane (Szakacs et al., 2008) (Fig. 1).
The cells of the blood brain barrier (BBB) are closely linked by tight junctions, which practically prevent hydrophilic molecules to diffuse between the cells into the central nervous system (CNS). However, as hydrophobic substances might diffuse through the membrane, it is the role of P-gp to keep those out as well (Neuhaus & Noe, 2009). The protecting function of P-gp at the BBB has been observed with mdr1 knock-out mice and the dog breed collie, which naturally lacks functional P-gp because of a mutated mdr1 gene. Collies are extremely susceptible to neurotoxic drugs and thus show dramatic adverse reactions after treatment with the antiparasitic drug ivermectin (Mealey et al., 2001). This detoxifying role of P-gp can be observed at other barriers as well (e.g. the fetal-maternal barrier), but also in faster clearance of administered drugs as it exports substrates form the hepatocytes into the bile, and from the intestinal epithelium into the intestinal lumen (Schinkel & Jonker, 2003). However, in drug research P-gp poses a large problem, since it highly influences pharmacokinetic properties of drugs. Because of the efflux behavior in the intestinal epithelium the oral bioavailability of drugs is hindered. There are a number of compounds that are able to modulate P-gp activity, which results in modified P-gp concentrations in the target tissue. As a consequence this can lead to adverse drug-drug interactions, when therapeutics are administered at the same time. Furthermore, elevated expression of P-gp, (as it is the case in cancer cells), is one major reason for the acquisition of MDR. One way to overcome these negative effects associated with P-gp activity would be the development of P-gp inhibitors that should restore sensitivity to therapeutics. Already 30 years ago, the reversal of resistance against the vinca alkaloids vincristine and vinblastine by the calciumchannel blocker verapamil was identified (Tsuruo et al., 1981). However, since then no inhibitor reached the market so far. This can be explained by its important physiological functions, rendering them rather antitargets than targets ).
Another ABC transporter that is highly associated with MDR belongs to the ABCC subfamily. MRP1/ABCC1 is located at the basolateral membrane of epithelial cells of the lung, kidney and the intestine (Fig. 1).
Although the substrate specificity of MRP1 shows some overlap with P-gp especially in terms of hydrophobic substances, MRP1 preferably binds to anionic substances in contrast to the positive charged substrates of P-gp (Borst & Elferink, 2002). Furthermore MRP1 is known for the export of hydrophilic substances as glutathione (GSH) conjugates. Therefore it is not only responsible for preventing xenotoxins entering the cell, but MRP1 also effluxes toxic metabolic compounds, which is highly important for faster clearance.
The role of MRP1 in the acquisition of MDR has particular impact on non-small cell lung carcinoma, a very aggressive cancer type, where high concentrations of MRP1 could be detected in the cancer cells.
The development of MRP1 specific inhibitors faces immense problems as MRP1 substrates and inhibitors show anionic properties, which lack good cell penetration properties (Schinkel & Jonker, 2003).
In 1998 Doyle et al. identified another ABC transporter that conferred resistance to the anthracenedione mitoxantrone, which is a poor substrate for P-gp and MRP1 (Doyle et al., 1998). BCRP belongs to the G or white subfamily of ABC transporters and received the name BCRP because of its isolation from a breast cancer cell line.
As P-gp, BCRP is located at the apical membrane of epithelial cells in the intestine, kidney and placenta (Schinkel & Jonker, 2003) (Fig. 1). Regarding substrate profiles, BCRP shows some overlap with P-gp and MRP1 but does not confer resistance to taxols, cis-platin and verapamil, or vinca alkaloids and anthracyclines. On the other hand, BCRP is known for transporting positively and negatively charged drugs (Sharom, 2008).
A specific inhibitor of BCRP is the tremorgenic mycotoxin fumitremorgin C (FTC). FTC blocked mitoxantrone transport by BCRP without affecting P-gp or MRP1-mediated drug resistance (Rabindran et al., 2000). However, due to neurotoxic effects in vivo application is still not possible.

Structure of ABC transporters
The minimal functional unit of an ABC transporter consists of two (pseudo)-symmetric halves that comprise a transmembrane (TM) and a nucleotide binding domain (NBD) (Fig.  2). In the case of P-gp and most other eukaryotic ABC transporters, these subdomains are fused to one polypeptide chain. On the contrary BCRP is a so-called half-transporter (Fig. 2). Half-transporters express each protein half separately and thus need to homo-dimerize to yield functional full transporters. The NBDs are responsible for the binding and hydrolysis of ATP, which is needed for the active transport of substrates. On each NBD sequence the characteristic ABC domain consisting of the Walker A and B region, as well as the "signature" or C motif, can be found (Fig. 3). One ATP molecule is supposed to be sandwiched between the Walker A and B of one NBD and the C motif of the other NBD (Fig. 3). As they are highly conserved, the NBDs show large sequence identity among ABC transporters. However, substrate binding and transport occurs at the TMDs. Each TMD consists of six TM helices, although this number varies between ABC transporters. The TMDs are much less conserved which leads to a large diversity in substrate profiles among ABC transporters.
During drug transport P-gp and its homologues undergo large conformational changes, converting an open-inward drug-binding state into an open-outward drug-releasing state (Rosenberg et al., 2001). This assumption was confirmed by cryo-electron microscopy and biochemical experiments, where P-gp was trapped in different states of the catalytic cycle (using the non-hydrolysable ATP analog AMP-PNP and ADP-Vi). The detailed mechanism of the energy driven drug transport, rendering the high-affinity into a low-affinity binding site, is currently hypothesized in two different ways and has been extensively reviewed in (Seeger & van Veen, 2009).

Homology modeling of P-gp
As already described in the introduction of this chapter, entries in the protein data bank (PDB) raise exponentially, but the structure determination of membrane proteins is still problematic and only relatively few structures have been resolved up to now. Thus, homology modeling is essential for performing docking or MD studies on most of the ABC transporters.
In 2001 the publication of X-ray structures of E. coli MsbA (PDB code: 1JSQ, resolution: 4.5Å), a lipid A transporter, raised a lot of interest in the ABC-transporter community (Chang & Roth, 2001).
The lipid flippase MsbA is an ABC protein that is responsible for the transport for lipid A and lipopolysaccharide (LPS). A non-functional MsbA leads to accumulation of lipopolysaccharide and phospholipids in the inner membrane of gram-negative bacteria.

www.intechopen.com
According to the X-ray structure published in 2001 the association of the two TMDs was interpreted as a chamber that provides alternating access for potential ligands during the catalytic cycle. The theory of MsbA switching between different conformations was confirmed by the subsequent publications of the X-ray structures of Vibrio cholerae MsbA in 2003 (PDB code: 1PF4, resolution: 3.80Å) (Chang, 2003) and of Salmonella typhimurium MsbA in complex with ADP·Vi (PDB code: 1Z2R, resolution: 4.20Å) (Reyes & Chang, 2005). The former presents the apo protein in a closed state, while the latter captures the Protein/ADP·Vi complex in the posthydrolytic state.
At this time these structures were the only source for structure-based design on MsbA and its homologues. Numerous ABCB1 homology models were generated relying on these MsbA templates Seigneuret & Garnier-Suillerot, 2003;Shilling et al., 2003;Stenham et al., 2003;Vandevuer et al., 2006). However, with the publication of the X-ray structure of the Staphylococcus aureus transporter Sav1866, an MsbA homologue (Dawson & Locher, 2006), the previous MsbA and two additional EmrE structures had to be retracted (Chang et al., 2006). According to Chang, an error in the in-house software that should process the crystallographic data resulted in a sign change and therefore to a momentous misinterpretation of the data (Matthews, 2007). This incident became the center of numerous discussions, often referred to as the "pentaretraction" (Davis et al., 2008;Penders et al., 2007). In contrast to the retracted MsbA models, the architecture of Sav1866 shows a helix arrangement that is analogous to domain swapping in other enzymes. Thus, TM helices of one TMD are in close contact with the opposite NBD via so-called coupling helices.
One year later, in 2007, Ward et al. published the corrected MsbA structures (Ward et al., 2007), which are in agreement with the SAV1866 architecture.
As SAV1866 (PDB code: 2HYD, resolution: 3.00Å) is one of the best resolved ABC exporters it has been often used as template for further modeling studies. In addition, this structure also fulfills most of the structural restraints that were obtained by cross-linking studies. However SAV1866 was crystallized in the nucleotide-bound conformation, which represents the ligand-releasing state of the protein. Thus the suitability of this template for docking studies can be questioned. In this respect Stockner et al. generated a data-driven homology model on the basis of SAV1866 that should represent the ligand-binding state of the protein by applying structural restraints in TM helices 6 and 12 obtained by crosslinking data on the model (Stockner et al., 2009).
The corrected MsbA coordinates cover different catalytic states, including a nucleotide-free ligand-binding conformation. Unfortunately these structures are resolved at resolutions far from being suitable for docking experiments, with some templates only represented by C atoms. Models on basis of the MsbA structures were therefore mainly used for exploring the conformational changes during the catalytic cycle.
With the publication of murine P-gp in March 2009 the first X-ray structure of a eukaryotic ABC exporter was available (Aller et al., 2009) (PDB code: 3G5U, resolution: 3.8Å). Two additionally published structures that include co-crystallized enantiomeric cyclic peptide inhibitors (CPPIs; QZ59-RRR and QZ59-SSS) highlight the binding-competence of these conformations and thus their great value for further docking studies. Furthermore the high sequence identity of 87% with human P-gp highly facilitates the modeling process. www.intechopen.com

Docking and MD studies
The definition of a binding site is an essential preparation step for docking studies. Regarding P-gp and other ABC transporters, we face the problem that hardly any binding sites for known P-gp ligands have been identified unambiguously. So far, it has been assumed that there is a large binding cavity in the transmembrane region (Loo & Clarke, 1999), which comprises distinct active sites. Furthermore, cysteine-scanning mutagenesis studies showed that the protein is able to bind at least two different molecules simultaneously (Loo et al., 2003). By using biochemical techniques a more detailed characterization of concrete binding sites for distinct substrates was possible (extensively reviewed in (Crowley et al., 2010;Loo & Clarke, 2008;Seeger & van Veen, 2009)). This led to the characterization of the interaction regions of Rhodamine 123 and Hoechst 33342, named R-and the H-site (Loo & Clarke, 2002;Qu & Sharom, 2002), together with a regulatory site, which binds prazosin/progesterone (Shapiro et al., 1999). Furthermore, the release of the Pgp/CPPI-complexes presented another step forward in elucidating drug/P-gp interactions.
Since the co-crystallized enantiomers showed distinct binding patterns, this information raised the assumption of stereoselectivity of P-gp in its ligand binding quality (Aller et al., 2009). Stereoselectivity has also been shown for flupentixol (Dey et al., 1999) and propafenone derivatives (Jabeen et al., 2011). On the other hand there are also ample reports on equal activity of enantiomeres. Thus, as for niguldipine and verapamil both enantiomers showed equivalent activities (Hollt et al., 1992;Luurtsema et al., 2003), the distomers with respect to cardiovascular activity were used for clinical studies.
As the resolution of the hitherto available templates used for constructing protein homology models is quite low, only very few docking studies have been conducted so far. Shortly after the publication of mouse P-gp, Pajeva et al. docked quinazolinones into a homology model of human P-gp based on the murine homologue, which is in complex with the cyclopeptide SSS-QZ59 (Pajeva et al., 2009). The binding site they used was defined by the co-crystallized ligands and was extended by 14Å. The results suggested interaction with TM helices 5, 6 and 11 and were further confirmed by a pharmacophore model. Becker et al. performed docking studies of the P-glycoprotein modulators colchicine, rhodamine B, verapamil and vinblastine into a homology model based on the nucleotidefree corrected MsbA structure (Becker et al., 2009). The resultant poses predicted that all ligands were able to interact with residues that were experimentally identified as important for ligand binding, strongly involving TM helices 5, 6, 7, 11 and 12. However, none of the drugs was able to contact every identified residue, which favors the hypothesis of distinct interactions sites forming one binding cavity.
Recently Dolghih et al. published a docking approach that was able to discriminate between known P-gp binders and non-binding metabolites (Dolghih et al., 2011). In this study there was a major interest in considering the high flexibility of P-gp. Therefore the induced fit protocol of the Schrödinger Suite was applied (Sherman et al., 2006). However, the discrimination between binders and non-binders can be more efficiently performed on basis of physicochemical properties than different binding mechanisms.
In our group, docking into a homology model based on mouse P-gp was used for explaining the stereoselective P-gp modulating activity of tricyclic benzopyranooxazines (Jabeen et al., www.intechopen.com 2011). Besides from activity differences, compounds with 4aS,10bR configuration showed a clear logP-activity correlation (r 2 =0.96), which was not the case for the 4aR,10bS series. This characteristic could be partly explained by the received binding hypotheses. The analysis of the docking poses by agglomerative hierarchical clustering resulted in distinct clusters for the different diastereomers. Therefore it has been hypothesized, that activity differences of the diastereomers is due to their different binding modes in the P-gp binding cavity. In addition, molecules with 4aR,10bS chirality were found close to the entry path of the protein, wherefore activity is primarily affected by the molecules' partition coefficient. On the other hand compounds of the 4aS,10bR series also showed docking poses at an active site in the binding pocket of P-gp, thus suggesting that the activity is dependent on multiple factors.
Furthermore, we were able to propose reliable binding hypotheses of propafenone analogs in P-gp by applying a knowledge driven docking protocol . Based on our extensive data from SAR studies on propafenones (Ecker et al., 2008;, we selected a small set of compounds for docking into a homology model based on mouse P-gp. As propafenone analogs show a clear SAR we assumed a similar binding mode of the docked propafenone derivatives. In that sense the resultant docking poses were clustered according the RMSD of their common scaffold. The clusters were prioritized according a combination of SAR data and protein-ligand interaction fingerprint information. With this protocol a high number of docking poses could be reduced to two reliable binding modes. Key interactions formed by these two clusters were formed with amino acids of TM helices 5, 6, 7 and 12 which were shown previously to be involved in ligand binding (Loo & Clarke, 2008;Seeger & van Veen, 2009).
In contrast to the compounds investigated above steroidal compounds are assumed to bind to the NBD rather than the TMDs. Several docking studies could show ATP-like binding of flavonoids, flavones and chalcones at the ATP-binding site, which is extensively reviewed in (Klepsch et al., 2010).
Regarding P-gp's high flexibility MD simulation represents a convenient technique to consider structural changes of the protein. Unfortunately, a number of MD studies were conducted relying on homology models based on the retracted MsbA X-ray structures Omote & Al-Shawi, 2006;Vandevuer et al., 2006) and are therefore partly no longer valid.
By now MD methods were mainly used for functional investigations of the protein. In order to determine the mechanisms of ATP hydrolysis numerous studies were conducted on isolated NBDs (Campbell & Sansom, 2005;Jones & George, 2007;Jones & George, 2009;Newstead et al., 2009;Wen & Tajkhorshid, 2008), as this comprises the sequence motives essential for ATP-binding. However, recently also studies considering the behaviour of the whole protein upon ATP hydrolysis were published (Becker et al., 2010;Gyimesi et al., 2011;Oliveira et al., 2011). All of those studies relied on the SAV1866 crystal structure, which represents the ligand-releasing and therefore open-outward state of the protein.
While Oliveira et al. (Oliveira et al., 2011) were able to show that replacing both ATP molecules in the NBDs by ADP structural changes in the protein occurred, Gyimesi et al. (Gyimesi et al., 2011) observed structural rearrangements already by exchanging one ATP molecule. This could be of great relevance for heterodimeric ABC proteins like P-gp, where an asymmetric ATP hydrolysis might be possible. In addition movements in TM helices 3 and 6 could be identified, which is in agreement with MD studies conducted by Becker et al. (Becker et al., 2010). Both groups observed closure of the TM domains after ATP hydrolysis.
The investigation of the drug-binding open-inside conformation of P-gp by MD simulation still faces numerous problems, due to instability of the mouse P-gp structure. In that sense the validity of this model is somewhat doubted (Gyimesi et al., 2011;Loo et al., 2010).

Biological background of the SLC-6 family
The concerted release and reuptake of transmitter substances is a basic principle of proper signal transduction in the nerve cells. In order to terminate a synaptic signal after neural firing, transporter proteins have to remove about 10 5 -fold of basal concentrations (Chen et al., 2004;Gouaux, 2009). The transporters practically have to act as selective molecular vacuum cleaners to deal with such huge loads of neurotransmitters in order to re-establish pre-signaling conditions within milliseconds. A major ion gradient serves as driving force and patron for the protein class: the Neurotransmitter:Sodium Symporters (NSS).
Synonymously called the solute carrier 6 family (SLC-6), NSS members include the sodiumand chloride-dependent transporters for GABA, dopamine, serotonin, norepinephrine and glycine, but also just sodium-dependent transporters of amino acids. Thus, the protein family is of particular medical importance, as many CNS diseases like depression, anxiety or epilepsy can be targeted by inhibiting transporters (Iversen, 1971).
They share a basic scaffold consisting of 12 transmembrane regions (TMs), segments 1-5 and 6-10 forming two pseudosymmetric domains housing the substrate and ion binding sites in partially unwound regions half-way across the membrane (Kanner & Zomot, 2008).
The crystal structure of LeuT, a bacterial orthologue of the eukaryotic NSS members, became available in an occluded state conformation in 2005 and in the open to out orientation in 2008 (Singh et al., 2008;Yamashita et al., 2005), thus revealing first detailed insights into the binding site topology. Furthermore, very recently a double mutant stabilized in an inward-open conformation was published (Krishnamurthy & Gouaux, 2012). These crystallographic snapshots fortify the so-called alternating access model for neurotransmitter membrane transport (Jardetzky, 1966). Various attempts have been made to clarify the exact molecular transport mechanism (Forrest et al., 2008;Shi et al., 2008), yet many questions remain unanswered. Concerning the quaternary structure, it is generally assumed that neurotransmitter:sodium symporters form constitutive oligomers (Forrest et al., 2008;Sitte et al., 2004). Despite a comparably poor average overall sequence identity between eu-and prokaryotic SLC-6 members of slightly above 20%, these structures paved the way to comparative modeling studies. Predominantly the monoamine transporters DAT, NET and SERT, but also GAT, have been modeled and studied extensively. For a comprehensive summary of the state of knowledge about the SLC-6 family, the reader is referred to the recent review by Kristensen et al. (Kristensen et al., 2011).

Examples of studies on the hSERT
As mentioned earlier, especially when dealing with low template-target sequence identity, a very careful sequence alignment including all possible experimental knowledge is crucial www.intechopen.com for the construction of reliable homology models. For the main members of the SLC-6 family a lot of effort has been put into this work, resulting in the comprehensive alignment of NSS sequences with the LeuT published by Beuming et al. in 2006(Beuming et al., 2006. Since then, some new structural insights into the protein class have been gained leading to slightly altered regions, but still the alignments can be considered a good starting point for experiments with NSS models. In the case of the hSERT, the recent work of Sarker et al. (Sarker et al., 2010) provides a good example for the cumulative value of combining molecular modeling methods with mutagenesis experiments in order to verify in silico elaborated hypotheses. For investigating the binding mode of tricyclic antidepressants (TCAs) in the serotonin transporter, comparative modeling marked the starting point for subsequent studies. Using the Beuming alignment, homology models of hSERT were built based on the previously mentioned high-resolution open-to-out structure of the LeuT published in 2008 (PDB code 3F3A). Subsequent docking studies of imipramine resulted in three pose clusters of potential binding modes, showing interactions to previously reported key residues (Andersen et al., 2009;Chen & Rudnick, 2000;White et al., 2006). A diagnostic Y95F mutation, a candidate residue for hydrogen bonding with the imipramine diaminopropyl moiety, significantly decreased imipramine affinity without affecting serotonin binding, ruling out one cluster. Further uptake and docking assays demonstrated that carbamazepine, structurally a truncated and slightly more rigid derivative of imipramine, was able to bind mutually non-exclusive with the substrate serotonin, whereas binding of its large-tailed relative is mutually exclusive. This led to the following conclusions: a) the tricyclic ring system of TCAs binds in an outer vestibule, and b) the basic side chain of imipramine points into the actual substrate binding site. Fig. 4. Molecular dynamics simulations of SERTThr-81 mutants reveal models favoring inward facing states. A, snapshot of wild type SERT after 16 ns of MD simulation. The Thr81 side chain forms a stable H-bond with the backbone carbonyl of Tyr350 in IL3. B, snapshot of SERTT81A after 6 ns of MD simulation; the H-bond is not formed between Ala81 and Tyr350 during the course of the simulation. C, snapshot of SERTT81D after 6 ns of MD simulation; no H-bond is formed between Asp81 and Tyr350 during the course of the simulation. (taken from ).
As an example for a more functional study on the SERT, the work of Sucic et al.  can be mentioned. As it was analogously reported for the DAT (Guptaroy et al., 2009), the important role of a highly conserved phosphorylation site at the N-terminus of the transporter in mediating the action of amphetamines was studied. Amphetamines are said to induce substrate efflux, but the way they do so is not well understood. Sucic et al. reported that mutating the highly conserved N-terminal residue T81 (a candidate site for phosphorylation by protein kinase C), to alanine or aspartate leads to subsequent fail of the transporter to support amphetamine-induced efflux. As it was also confirmed by molecular dynamics simulations of the wild type transporter, the in silico mutated SERT T81A and SERT T81D , the data suggested that by phosphorylation or in silico mutation of T81 the conformational equilibrium of the serotonin transport cycle alters towards the inward facing conformation. As seen in the MD studies, this happens due to a loss of a hydrogen bond network of T81 with Y350 in IL3 by these mutations. Furthermore, an increased distance between the C terminus (i.e. the most distal point of TM12) and the N terminus after in silico mutation was observed. This example nicely indicates how functional MD studies might aid in elucidating biological relevant phenomena.

Studies on hGAT models
The four Na + -and Cl --dependent GABA transporters, GAT-1-3 and BGT-1 (SLC6A1, A16, A11, A12), provide a similar percentage of sequence identity to the LeuT. The subtype showing the highest quantity in the CNS is GAT-1. It is also the best-investigated, and the only one currently targeted by a marketed drug, the second-line antiepileptic tiagabine (Gabitril®). Accordingly, systematic synthesis studies in order to discover even more selective compounds have been performed mainly on GAT-1. Nevertheless, other subtypes should not be ignored, as they may be the key to a less side-effect afflicted antiepileptic therapy, as tiagabine efficacy as anticonvulsant is limited, and its use was connected to several adverse effects like sedation, agitation, or even seizure induction. Neuronal GABA reuptake, mainly done by GAT-1, leads to subsequent recycling of the transmitter substance. On the contrary, astroglial uptake of GABA leads to degradation, suggesting subtypes predominantly present in glia cells being an interesting target for enhancing overall GABA levels. For example, the lipophilic GABA analog EF-1502, characterized by GAT1 and GAT2 (BGT-1) selectivity, showed synergistic anticonvulsant activity, when administered with tiagabine (Schousboe et al., 2004), although BGT1 levels in the CNS are about 1000-fold lower, and even a recent study with BGT-1 knockout mice did not show any change in seizure susceptibility (Lehre et al., 2011).
In the search for potent selective non-GAT-1 inhibitors, GABA mimetic moieties (like Rnipecotic acid in tiagabine, β-alanine or THPO [4,5,6,7-Tetrahydroisoxazolo(4,5-c)pyridin-3ol]) were systematically combined with large aromatic side chains, both in order to increase the affinity and to make the compounds blood-brain barrier permeable (Andersen et al., 1993;Andersen et al., 1999;Clausen et al., 2005;Knutsen et al., 1999;Kragler et al., 2008). Unfortunately, up to now no truly selective tools for the evaluation of non-GAT-1 inhibition are available, although the GAT-1/BGT-1 inhibitor EF1502 and SNAP-5114, showing a certain GAT-2/GAT-3 selectivity, mark a good starting point (Madsen et al., 2010). Thus, further insights into the molecular basis of ligand binding are sought by the aid of in silico methods.
GAT-1 has been subject of several comparative modeling studies. Initial studies predominantly aimed at clarifying the GABA binding mode in the occluded transporter state, which is quite well documented so far (Pallo et al., 2007;Wein & Wanner, 2009). Though, compounds with large aromatic tails cannot be accommodated in the occludedstate active site, as the entrance to the binding pocket is barred by the two extracellular gate residues R69 and D451, as well as the F294 side chain, forming the binding site "roof". In order to study tiagabine-like ligands, constructing open-to-out models seemed inevitable, as it was done by Skovstrup et al. (Skovstrup et al., 2010). Structures of both states were modeled and refined exhaustively, as described in section 2.1. The combined use of docking and molecular dynamics simulation was chosen to investigate binding of GABA, its analogue (R)-nipecotic acid and the high active (R)-enantiomer of tiagabine. The results for GABA binding were in line with the earlier mentioned experiments. In case of tiagabine, MD simulations helped to distinguish between the cis-and trans-conformer, both being possible states due to the protonated state of tiagabine at physiological pH. During the MD, the trans-conformer immediately stirred away to the extracellular space, whereas the other one remained stable in the binding site. Summing up, GABA and (R)-tiagabine turned out having two different binding modes, sharing the orientation of the carboxy group towards one of the co-transported sodium ions as a common feature.
For the other GAT subtypes, things are a bit more complicated. Looking at the residues corresponding to LeuT substrate binding site, just a few candidate residues differ significantly, being somehow unlikely to be fully responsible for subtype selective binding. So far, molecular modeling studies have been performed, but highly similar binding sites and the lack of selective ligand data limited their explanatory power (Pallo et al., 2009). Thus, a huge field of activity remains to be explored on the way to fully understand the differences between the GABA subtypes, in silico methods being a valuable tool for stepwise adding pieces of information to the big puzzle.

Concluding remarks
Membrane transport proteins are responsible for one of the most important processes in living cells: directed transport across barriers. They comprise about 30% of known proteomes and constitute about 50% of pharmacological targets. Although, due to difficulties in expression, purification and crystallization, only about 2% of the high resolution crystal structures in the Protein Data Bank (PDB) are transporters. Thus, computational methods have been utilized extensively to provide significant new insights into protein structure and function. Above all, molecular modeling and molecular dynamics (MD) simulations may deliver atomic level details to reveal the molecular basis of e.g. drugtransporter interactions. As shown on basis of recent research examples, in silico methods in many cases can provide additional information to biological experiments, either underpinning pharmacological results or they may even lead to new insight, not being biologically accessible.

Acknowledgments
The authors gratefully appreciate financial support provided by the Austrian Science Fund (FWF), grant SFB3502 and SFB3506.  Over the recent years, medicinal chemistry has become responsible for explaining interactions of chemical molecules processes such that many scientists in the life sciences from agronomy to medicine are engaged in medicinal research. This book contains an overview focusing on the research area of enzyme inhibitors, molecular aspects of drug metabolism, organic synthesis, prodrug synthesis, in silico studies and chemical compounds used in relevant approaches. The book deals with basic issues and some of the recent developments in medicinal chemistry and drug design. Particular emphasis is devoted to both theoretical and experimental aspect of modern drug design. The primary target audience for the book includes students, researchers, biologists, chemists, chemical engineers and professionals who are interested in associated areas. The textbook is written by international scientists with expertise in chemistry, protein biochemistry, enzymology, molecular biology and genetics many of which are active in biochemical and biomedical research. We hope that the textbook will enhance the knowledge of scientists in the complexities of some medicinal approaches; it will stimulate both professionals and students to dedicate part of their future research in understanding relevant mechanisms and applications of medicinal chemistry and drug design.