Computational studies of AcrB, a Multidrug Efflux Pump
Dissertation
zur Erlangung des akademischen Grades eines Doktors der Naturwissenschaften (Dr. rer. nat.)
vorgelegt von
M. Sc. Biochem. Zhou Wenchang
an der
Mathematisch - Naturwissenschaftliche Sektion
Fachbereich Biologie
Tag der mündlichen Prüfung: 06.12.2012
1. Referent: Prof. Dr. Kay Diederichs 2. Referent: PD. Dr. Thomas E. Exner
,
Content
1. Zusammenfassung...1
2. Summary...2
3. Introduction...3
3.1 Antibiotic resistance... 3
3.2 Multi-drug efflux transporters... 3
3.3 RND/MFP/OMF efflux system... 4
3.4 An asymmetric structure of AcrB ... 6
3.5 Functional rotation mechanism... 7
3.6 Available experimental validation for the functional rotation mechanism ... 9
3.7 Why molecular simulations?...11
3.8 The aim of the project ... 12
4. Theories and Methods...14
4.1 Thermodynamics... 14
4.2 Statistical mechanics ... 15
4.2.1 Introduction... 15
4.2.2 Ensemble... 16
4.2.3 Monte Carlo simulation ... 17
4.2.4 Molecular Dynamics simulation ... 18
4.2.5 Free energy calculations... 22
4.3 Poisson-Boltzmann Calculations ... 25
4.4 GRIFFIN ... 27
4.4.1 Initial carving of the lipid bilayer ... 28
4.4.2 Calculation of the implicit protein force field... 28
4.4.3 Lipid bilayer equilibration with MD simulations ... 29
5. Using Poisson-Boltzmann Calculations to predict the pKa values...31
5.1 Motivation... 31
5.2 Methods... 32
5.2.1 Single-Site Titration Method... 32
5.2.2 Multi-Site Titration Method... 34
5.2.3 pKa calculations ... 35
5.3 Results and Discussion ... 36
5.3.1 Dielectric constant scanning for pKa prediction on X-ray structure ... 36
5.3.2 The probability of D407, D408 and K940 of being protonated... 38
6. Setting up a simulation system for AcrB in a lipid membrane...46
6.1 Motivation... 46
6.2 Methods... 46
6.3 Results and Discussion ... 48
7. Conventional MD simulations of the AcrB trimer...51
7.1 Motivation... 51
7.2 Methods... 51
7.2.1 Simulation system set up ... 51
7.2.2 Energy minimization and equilibration... 53
7.3 Results and Discussion ... 54
7.3.1 The equilibration of the trimer system... 54
7.3.2 Dynamical behavior of the trimer during production run... 56
7.3.3 Water distributions in the transmembrane domains ... 59
7.3.4 Water dynamics in the transmembrane domains... 62
7.3.5 Movement of the minocycline around the drug binding site in the T state ... 69
8. MD simulations on truncated protomers...72
8.1 Motivation... 72
8.2 Methods... 72
8.2.1 Simulation system set up ... 72
8.2.2 Energy minimization and equilibration... 73
8.3 Results and Discussion ... 74
8.3.1 Comparison between trimer and protomer simulations ... 74
8.3.2 Water distributions in the transmembrane domains ... 76
8.3.3 Water dynamics in the transmembrane domains... 78
9. Free energy perturbation calculations based on the equilibrated protomers85 9.1 Motivation... 85
9.2 Methods... 85
9.2.1 Dual-topology method ... 85
9.2.2 Simulation system set up ... 86
9.3 Results and Discussion ... 87
10. Structure analyses of the crystal structures...89
10.1 Motivation... 89
10.2 Methods... 90
10.3 Results and Discussion ... 91
10.3.1 Identifying unique structural elements in the transmembrane and periplasmic domain ... 91
10.3.2 Clash volumes between the periplasmic domains in adjacent protomers94 11. A mechanism of substrate/proton transport for AcrB...97
11.1 Mechanism of the substrate/proton antiporting developed within one protomer... 97
11.2 The proposed state transition paths ... 100
12. Conclusions and Perspectives...102
13. Acknowledgments...105
14. Bibliography...106
List of Figures & Tables
Fig.3.1 Proposed model of AcrA-AcrB-TolC tripartite efflux system . . . .5
Fig.3.2 Architecture of AcrB in symmetrical form . . . 5
Fig.3.3 Major structural differences of the AcrB protomers . . . .6
Fig.3.4 Side view superimposition of the AcrB L protomer and the PN1, PC2 subdomains as well as TM8 of the O protomer . . . .7
Fig.3.5 Schematic representation of the AcrB functional rotation mechanism of drug transport . . . .8
Fig.3.6 Positions of the engineered disulfide bridges introduced into AcrB . . . .10
Fig.4.1 The schematic protocol of GRIFFIN for optimization of a protein-membrane interface . . . 28
Fig.4.2 Schematic of force-redirection methodology . . . 30
Fig.5.1 The different types of geometrical objects mapped into the modeling system .35 Fig.5.2 The titration curves for residue D407, D408 and K940 in the L state . . . 36
Fig.5.3 The titration curves for residue D407, D408 and K940 in the T state . . . 37
Fig.5.4 The titration curves for residue D407, D408 and K940 in the O state . . . .38
Fig.5.5 The probability of D407 being protonated in the L, T and O state . . . .39
Fig.5.6 The probability of D407 being protonated in the L, T and O state . . . .39
Fig.5.7 The probability of D407 being protonated in the L, T and O state . . . .40
Fig.5.8 The titration curves for residue E346 and E417 in the L state . . . .41
Fig.5.9 The titration curves for residue E346 and E417 in the T state . . . .42
Fig.5.10 The titration curves for residue E346 and E417 in the O state . . . .43
Fig.5.11 The probability of E346 and E417 being protonated in the L, T and O state . . .44
Fig.6.1 The structure of the lipid bilayer before and after carving . . . .47
Fig.6.2 Number of atoms inside protein volume and the maximum depth of any atoms inside the protein volume . . . .49
Fig.6.3 The structure of the lipid bilayer after the equilibration . . . 50
Fig.7.1 Representation of the AcrB trimer simulation system . . . 52
Fig.7.2 The RMSD values of the L protomer during the equilibration . . . .55
Fig.7.3 The RMSD values of the T protomer during the equilibration . . . .55
Fig.7.4 The RMSD values of the O protomer during the equilibration . . . .56
Fig.7.5 The RMSD values of the L protomer during the 250 nanoseconds MD simulation . . . .57
Fig.7.6 The RMSD values of the T protomer during the 250 nanoseconds MD simulation . . . .57
Fig.7.7 The RMSD values of the O protomer during the 250 nanoseconds MD simulation . . . .58
Fig.7.8 The RMSF values of the L state for all the residues during the 250 nanoseconds
MD simulation . . . 58
Fig.7.9 The RMSF values of the T state for all the residues during the 250 nanoseconds MD simulation . . . 59
Fig.7.10 The RMSF values of the O state for all the residues during the 250 nanoseconds MD simulation . . . .59
Fig.7.11 Water distribution in the transmembrane domain of the L state for the AcrB trimer simulation . . . .60
Fig.7.12 Water distribution in the transmembrane domain of the O state for the AcrB trimer simulation . . . .61
Fig.7.13 Water distribution in the transmembrane domain of the T state for the AcrB trimer simulation . . . .61
Fig.7.14 Number of water molecules in L-I water wire during the MD simulation of the AcrB trimer . . . .63
Fig.7.15 Number of water molecules in L-II water wire during the MD simulation of the AcrB trimer . . . .63
Fig.7.16 Number of water molecules in L-III water wire during the MD simulation of the AcrB trimer . . . .64
Fig.7.17 Number of water molecules in L-IV water wire during the MD simulation of the AcrB trimer . . . .64
Fig.7.18 Number of water molecules in T-I water wire during the MD simulation of the AcrB trimer . . . .65
Fig.7.19 Number of water molecules in T-II water wire during the MD simulation of the AcrB trimer . . . .65
Fig.7.20 Number of water molecules in T-III water wire during the MD simulation of the AcrB trimer . . . .66
Fig.7.21 Number of water molecules in T-IV water wire (equivalent to water wire L-IV) during the MD simulation of the AcrB trimer . . . .66
Fig.7.22 Number of water molecules in T-IV water wire (equivalent to water wire O-IV) during the MD simulation of the AcrB trimer . . . .67
Fig.7.23 Number of water molecules in O-I water wire during the MD simulation of the AcrB trimer . . . .67
Fig.7.24 Number of water molecules in O-II water wire during the MD simulation of the AcrB trimer . . . .68
Fig.7.25 Number of water molecules in O-III water wire during the MD simulation of the AcrB trimer . . . .68
Fig.7.26 Number of water molecules in O-IV water wire during the MD simulation of the AcrB trimer . . . .69
Fig.7.27 Skeletal formulas of the minocycline and dodecyl-alpha-D-maltoside . . . .70
Fig.7.28 Top view and side view of the drug binding site in the T protomer . . . .70
Fig.7.29 Top view and zoom in view of the drug binding site in the T protomer . . . .71
Fig.8.1 The MD simulation setup of one protomer system . . . .73
Fig.8.2 RMSD comparison between trimer and protomer simulations . . . .74
Fig.8.3 RMSF comparison between trimer and protomer simulations . . . 75
Fig.8.4 Water distribution in the transmembrane domain of the L state for the AcrB
protomer . . . 76
Fig.8.5 Water distribution in the transmembrane domain of the O state for the AcrB protomer . . . 77
Fig.8.6 Water distribution in the transmembrane domain of the T state for the AcrB protomer . . . 77
Fig.8.7 Number of water molecules in L-I water wire during the MD simulation of the AcrB protomer . . . 78
Fig.8.8 Number of water molecules in L-II water wire during the MD simulation of the AcrB protomer . . . 78
Fig.8.9 Number of water molecules in L-III water wire during the MD simulation of the AcrB protomer . . . 79
Fig.8.10 Number of water molecules in L-IV water wire during the MD simulation of the AcrB protomer . . . .79
Fig.8.11 Number of water molecules in T-I water wire during the MD simulation of the AcrB protomer . . . .80
Fig.8.12 Number of water molecules in T-II water wire during the MD simulation of the AcrB protomer . . . .80
Fig.8.13 Number of water molecules in T-III water wire during the MD simulation of the AcrB protomer . . . .81
Fig.8.14 Number of water molecules in T-IV water wire (equivalent to water wire L-IV) during the MD simulation of the AcrB protomer . . . 81
Fig.8.15 Number of water molecules in T-IV water wire (equivalent to water wire O-IV) during the MD simulation of the AcrB protomer . . . 82
Fig.8.16 Number of water molecules in O-I water wire during the MD simulation of the AcrB protomer . . . .82
Fig.8.17 Number of water molecules in O-II water wire during the MD simulation of the AcrB protomer . . . .83
Fig.8.18 Number of water molecules in O-III water wire during the MD simulation of the AcrB protomer . . . .83
Fig.8.19 Number of water molecules in O-IV water wire during the MD simulation of the AcrB protomer . . . .84
Fig.9.1 Dual-topology description for an alchemical free energy perturbation . . . .86
Fig.9.2 Skeletal formulas of the aspartate acid, an analog of aspartate acid, acetate, lysine and its analog n-butyl-amine . . . .88
Fig.10.1 Top view of the two types of contact between the domain pairs . . . .90
Fig.10.2 Transmembrane domain of AcrB in the T state . . . .91
Fig.10.3 Drug binding site in the T state . . . .93
Fig.10.4 Top view of the two types of domain contacts between the two protomers A and B . . . .94
Fig.10.5 Clash volumes of the hypothetical trimers . . . .96
Fig.11.1 Schematic representation of the AcrB substrate/proton antiporting mechanism . . . .99
Fig.11.2 Proposed state transition paths for the AcrB trimer . . . .101
Table 3.1 Distances of residue pairs and levels of disulfide cross-linking . . . 11 Table 7.1 The protocol of applying restraints during equilibration . . . .53
Table.9.1 Summary of free energy calculations . . . .87
Table.10.1 RMSD values of repeat 1 excluding TM2 and repeat 2 excluding TM8 between the L, T and O state . . . .92 Table.10.2 RMSD values of repeat 2 excluding TM8 after superimposing on repeat 1 excluding TM2 . . . .92 Table.10.3 RMSD values of PN1, PN2, PC1 and PC2 domains between L, T and O state of the wild type . . . 93 Table.10.4 Two types of clash volume of different domain pairs in the A and B protomer . . . 95 Table.10.5 The clash volume comparison between hypothetical trimers and crystal structures of the LLL and TTT state . . . 96
Abbreviations
AcrB, Acriflavine resistance protein B MD, Molecular Dynamics
L, Loose T, Tight O, Open
PBE, Poisson-Boltzmann Equation GRIFFIN, GRId-based Force Field INput
CHARMM, Chemistry at HARvard Molecular Mechanics NAMD, NAnoscale Molecular Dynamics
POPC, 1-Palmitoyl-2-Oleoyl-sn-glycero-3-PhosphoCholine TM, TransMembrane
RMSD, Root Mean Square Deviation RMSF, Root Mean Square Fluctuation
1. Zusammenfassung
Der AcrA-AcrB-TolC Komplex spielt eine zentrale Rolle in der Resistenz von Escherichia coli gegen Antibiotika und andere Chemikalien. Dieser Proteinkomplex exportiert Gallensalze, Detergenzien, organische Lösungsmittel und einige Antibiotika.
Diese Chemikalien zeigen keine strukturelle Ähnlichkeit. AcrB ist dabei der Kern dieser Maschinerie. Es wurde gezeigt, dass AcrB ein asymetrisches Trimer ist und jedes Protomer während des Transportprozesses zyklisch eine von drei Konformationen einnimmt: offen (open, O), locker (loose, L) und eng (tight, T). In der Transmembranregion konnte ein Protonenfluss nachgewiesen werden. Mittels Titrationexperimenten wurden die Aminosäuren D407, D408 und K940 als mögliche Protonen-Akzeptoren bestimmt, was durch Mutationsstudien nachgewiesen werden konnte. In der vorliegenden Arbeit verwenden wir Molekulardynamische simulationen um diesen Protonenfluss näher zu untersuchen. Mit Hilfe von Poisson-Blotzmann-Berechnungen und “free-energy pertubation“, einer speziellen Simulationsmethode, haben wir zunächst den theoretischen pKa dieser Aminosäuren bestimmt. Die Resultate zeigen wie viele Protonen sich in der Transmembranregion aufhalten. Die konventionelle Molekulardynamische simulation zeigte ausserdem unterschiedliche Wasserwege durch die verschiedenen Protomere. Diese Wasserwege zeigen mögliche Pfade für Protonenbewegungen durch das Protein auf. Die Analyse bekannter Kristall-Strukturen ergab einzigartige strukturelle Elemente sowohl in der periplasmatischen als auch in der Transmembrandomäne. Durch Berechnung der strukturellen Konflikte in mehreren hypothetischen Trimeren schlagen wir mögliche Zustandsübergangspfade unter unterschiedlichen Wirkstoffkonzentrationen vor. Durch Kombination mit früheren strukturellen und Computerstudien und den Kristallstrukturanalysen wird in der vorliegenden Arbeit eine globale Arbeitshypothese für den Substrate/Protonen-Transport durch AcrB aufgestellt.
2. Summary
AcrA-AcrB-TolC plays a major role in drug resistance in Escherichia coli by extruding bile salts, detergents, organic solvents and many antibiotics that are not structurally related. AcrB, the core of this machinery, has been shown to be functional as an asymmetrical trimer, each protomer of which adopts a conformational state designated as loose (L), tight (T) and open (O) in the cycle of drug transport. Proton translocation has been shown to happen in the transmembrane domain where there are the three titratable amino acids D407, D408 and K940. Mutagenesis studies have shown that these three amino acids serve as the transient station for protons. In the present thesis we used molecular simulations to study the proton translocation across the transmembrane domain. At first, the pKa values of D407, D408 and K940 were predicted using Poisson-Boltzmann calculations and alchemical free energy perturbation simulations. The results show how many protons there are in the transmembrane domain of each state. From the conventional Molecular Dynamics simulations, different water wires were found in the transmembrane domains of the L, T and O state. These water wires revealed the possible path for proton across the transmembrane domains of AcrB. In the analyses of crystal structures, unique structural elements were identified both in the periplasmic and transmembrane domain. By calculating the clashes between protomers in several hypothetical trimers, a scheme of possible state transition pathways is proposed for AcrB working under different drug concentrations. Combining with the former structural studies, computational studies and crystal structure analyses in the current thesis, a global working hypothesis of substrate/proton transport is proposed for AcrB.
3. Introduction
3.1 Antibiotic resistance
Antibiotic resistance of pathogenic bacteria is an increasing threat for public health in the treatment of cancer or infectious diseases [1-4]. There have been many outbreaks of bacterial infections in the last few decades in Europe, Asia and America with many people affected and killed [5]. In dealing with drug resistance, on one hand there are many new improved antibiotics coming out and on the other hand more and more basic research move to this field in order to address problems related to antibiotic resistance [6-8]. As shown by recent studies [9-14], overexpression of multi-drug resistance efflux pumps plays an important role on the intrinsic drug resistance in Gram-negative bacteria; whereas in early studies, low permeability of outer membrane proteins contributed to the drug resistance entirely [15-18].
3.2 Multi-drug efflux transporters
Multi-drug efflux transporters have been described and classified in 5 different superfamilies [12], namely the ATP binding cassette (ABC) superfamily, the major facilitator superfamily (MFS), the resistance-nodulation-division (RND) superfamily, the drug/metabolite transporter (DMT) superfamily and the multi antimicrobial extrusion (MATE) superfamily. Compared with other superfamilies each one with a narrow range of substrate specificity, members of the RND superfamily [9] have been shown the ability to extrude a wide range of substrate including detergents, dyes, bile salts, and many antibiotics with the help of another two components from the family of membrane fusion protein (MFP) and outer membrane factor (OMF).
3.3 RND/MFP/OMF efflux system
RND/MFP/OMF efflux system has a unique architecture which spans both the inner and outer membrane so that substrates can be extruded directly from cellular membrane to the medium [9]. AcrA/AcrB/TolC is an important drug efflux system and well characterized in E. coli. It is composed of three parts (Fig.3.1): TolC, a member of the OMF family, which is located in the outer membrane; AcrB, a member of the RND family, which is the core component in the inner membrane which provides energy for substrate extrusion; and AcrA, a member of the MFP family, which serves as a periplasmic adaptor linking the two components in the outer membrane and inner membrane [19]. Although the crystal structure of this tripartite efflux system has not been observed yet, structural studies on single components have been carried out for the last decade taking the advantage of crystallography [20-21], especially the research on one component of this tripartite system, AcrB. This component is responsible for the broad substrate specificity and active transport of various substrates by using proton motive force, sustained by different proton concentrations on the two sides of the inner membrane [22-31]. Three protomers of AcrB are arranged symmetrically as a homotrimer in the first solved crystal structure.
Each protomer is composed of a TolC docking domain where TolC might directly dock into AcrB, a pore domain where the substrates bind and a transmembrane domain which is the place for the proton translocation. In the pore domain, there are four subdomains, PN1, PN2, PC1 and PC2, arranged in a clockwise order.
Fig.3.1 Proposed model of AcrA-AcrB-TolC tripartite efflux system [32]. Structures of AcrA and TolC are manually docked to AcrB according to engineered cysteine
cross-linking studies between AcrA-TolC [33] and AcrB-TolC [34].
Fig.3.2 Architecture of AcrB in symmetrical form [35]. Three protomers are colored separately in blue, green and red. Side view (left) and top view (right) of AcrB are
shown as ribbon representation.
3.4 An asymmetric structure of AcrB
An asymmetric structure of AcrB was first solved subsequently by two independent groups [28-29] by growing crystals in a new form which is different from the one for the symmetry structure solved previously [35]. Compared with the symmetric structure, the asymmetric structure has three protomers in different conformations referred to as loose (L), tight (T) and open (O), each conformation of which corresponds to one of the consecutive functional states participates in the cycle of the substrate transport. In the upper periplasmic domain, they found three unique tunnels structurally formed by different orientations of the PN1+PC2 domain and the PN2+PC1 domain.
Fig.3.3 Major structural differences of the AcrB protomers [28]. (A) Top view of three AcrB protomers shown as cylinder representation in blue (L), yellow (T) and red (O).
The three asymmetric protomers are superimposed onto the symmetric trimer shown in grey. A hydrophobic pocket in each protomer is defined by phenylalanines 136, 178, 610, 615, 617, and 628; valines 139 and 612; isoleucines 277 and 626; and tyrosine 327 at the PN2/PC1 interface. (B) Structural changes in the putative proton translocation site, viewed from the cytoplasm. In the L and T protomer, the same conformation is observed for conserved residues D407, D408 and K940; whereas in the O protomer, the side chain of K940 flips away from the two aspartates 407 and 408.
In the lower transmembrane domain, three titratable residues (K940, D407 and D408) form the proton binding site which is thought to be the transit stop for proton transfer.
In the third protomer, the configuration of proton binding site is different from the one in the other two protomers, with the side chain of K940 flipping away from D407 and D408. Another difference in the transmembrane domain is the different orientations of TM2 and TM8 which connect the transmembrane and the periplasmic domain.
Fig.3.4 (A) Side view superimposition of the AcrB L protomer (grey) and the PN1, PC2 subdomains as well as TM8 of the O protomer (red) [28]. (B) The close-up view of the boxed region of (A) the N-terminal part of TM8 (residues 859-880) and the PC2 subdomain (residues 679-721 and 822-858) are superimposed. The structures colored in blue, yellow and red represent the conformations of the TM8 in the L, T and O protomers, respectively.
3.5 Functional rotation mechanism
A functional rotation mechanism was proposed by Pos and co-workers based on their findings from the asymmetrical crystal structure [28]. In the upper periplasmic domain where the substrate binding sites is, they found three unique tunnels in three protomers, each protomer of which had a different tunnel penetrating the periplasmic domain. These three tunnels gradually connect the pathway from lateral clefts to deeper substrate binding sites and then to the TolC docking domain which docks with TolC, the upper port for this tripartite efflux system. It is proposed that each protomer has to go over the three states to accomplish one cycle of substrate transport. The
three protomers are called loose (L), tight (T) and open (O), according to the different situations for substrate binding.
Fig.3.5 Schematic representation of the AcrB functional rotation mechanism of drug transport [28]. The three unique conformational states loose (L), tight (T) and open (O) are colored in blue, yellow and red, respectively. (A) The side-view schematic representation of two of the three protomers of the AcrB trimer. AcrA and TolC are indicated in light green and light purple colors. The proposed proton translocation site (D407, D408 and K940) is indicated in the membrane part of each protomer. (B) The lateral grooves in the L and T protomer indicate the substrate binding sites. The different geometric shapes reflect low (triangle), high (rectangle), or no (circle) binding affinity for the transported substrates. The PN1 subdomains (including the pore helices) located in the middle of the model are highlighted and form the corners of an asymmetric triangle (white) to indicate the communication between the protomers.
This mechanism also described how proton transfer is coupled with substrate transport as they also found three close-contact titratable residues which had different orientations in state T and O. In the upper periplasmic domain, what happened in L state is one protomer binds a substrate; subsequently the substrate is transported to the deeper binding site in the T conformation and finally extruded towards the TolC docking main along a new tunnel in the O conformation. The state transition from T to O is suggested to be coupled with proton transfer happening in the transmembrane
domain where TM8 zips on and off reversibly induced by the side chain rearrangement of the three titratable residues.
The conformational cycling and functional rotation mechanism adopted by AcrB, is very similar to the sequence of conformational changes of ATP synthesis by the F1Fo ATPase [36]. This similarity implies that the conformational state of each protomer is dependent on the conformational states of the adjacent protomers, which presents a scheme of cooperation between the three AcrB protomers.
3.6 Available experimental validation for the functional rotation mechanism
Cross-linking and site-specific mutagenesis experiments have been done independently by several groups [24, 37-38] to verify the proposed functional rotation mechanism which suggests that each protomer has to cycle through the three different states L, T and O. These three states have differences in the transmembrane and the subdomains PN1, PN2, PC1 and PC2 which form the pore domain of the AcrB protomers (Fig.3.3). Therefore, it was expected that the distances between the residues of the moving subdomains were considerably variable. According to the structural information of the asymmetric trimer, residues were carefully selected and mutated to cysteines which can form disulfide bond between pair of cysteines when the Sγ atoms are within at a distance of 6.4 Å. These cysteines were introduced at interfaces of subdomains which were considered to have a relative movement during the substrate transport by AcrB. An increased susceptibility to noxious compounds was observed for several pairs of mutations as the result of cross-linking, which regain the ability to transport substrate after adding the reducing agent DTT.
Among those selected cysteine pairs, there were several pairs at the interface PN2-TM1 and PN1-PN2 (Table 3.1) where the observed amounts of cross-linking were considered to be higher if all the three unique protomers were present at the same time. This suggested that more than one protomer per trimer were cross-linked.
Comparing the distances of residue pairs that were mutated to cysteine between the L,
T and O state, it was expected that the cysteine pair S132C_A294C was cross-linked only in the L state, and the cysteine pair V32C_N298C was cross-linked only in the T state, assuming that cross-linking occurs only at distances within the theoretical value of 6.4 Å between the Sγ atoms. The higher amount of cross-linking products indicated that there must be more than one protomers that are in the L state or T state.
Fig.3.6 Positions of the engineered disulfide bridges introduced into AcrB [24]. The three AcrB protomers shown as cylinder representations in blue (L), yellow (T) and red (O) are superimposed onto the structure of the L protomer shown in transparent blue. a depicts a side view of the entire AcrB trimer, and b represents a top view onto the pore domain from the periplasm. Enlarged views of the boxed regions (dashed lines) are shown at the left and right borders of the figure. The locations of the disulfide bridges are indicated by encircled arrows (green). Below the close-up views, the subdomains harboring the cysteine residues and the respective substitutions are indicated.
Table 3.1 Distances of residue pairs and levels of disulfide cross-linking [24]. Relative amount of cysteines involved in disulfide cross-linking is determined by quantitative MALDI-TOF analysis. R558C_E839C cysteine pair substitution is chosen as the control since these two cysteine residues are not expected to be cross-linked because of the long-range distance between the two residues in all three protomers.
Distance between Sγ atoms (Å) Protomer
Linked cysteine residues L T O Disulfide cross-linking (%)
R558C_E839C (PC2-TM7) 16.2 14.0 10.3 -9.5±4.7
S132C_A294C (PN1-PN2) 6.3 17.5 11.1 41.6±0.6
V32C_N298C (PN2-TM1) 7.2 3.5 7.0 41.3±1.2
V32C_A299C (PN2-TM1) 9.5 4.8 11.9 18.1±1.0
Q229C_T583C (Loop-PC1) 5.8 7.4 6.4 69.4±0.5
Q229C_R586C (Loop-PC1) 5.5 7.8 6.7 46.4±0.5
3.7 Why molecular simulations?
Although a functional rotation mechanism has been understood by the evidence we got from structural biology and crystallography, many aspects of the details are still not understood. One of them is how protons are transferred through the transmembrane domain during drug efflux. While the three titratable residues are thought to be the media for proton binding and release, which is confirmed by mutagenesis experiments [23, 27-29], the information about protons is missing in the crystal structures since the resolution is not high enough to see protons.
Molecular simulation has various and well-established theoretical frameworks for modeling interactions and dynamical properties that are not easy to access by experiments [39-46]. Powered by techniques of high-performance computing and development of computer hardware, nowadays we can extend our simulations to a microsecond timescale which is a tremendous progress compared with the first all-atom MD simulation carried out more than three decades ago [47-48]. In the last
few years, simulations applied on systems with tens of thousands of atoms have become feasible and practical [47, 49]. In the meantime various systems, ranging from soluble protein, DNA, RNA, virus and membrane protein, are studied by molecular simulations [50-57].
MD simulations are now widely used to study membrane protein mechanisms at atomic level, which can provide novel interpretations and working hypotheses to understand the structural basis of the protein function [39, 58]. There are many examples to show how MD simulations can help to understand the relationship between structure and function of membrane proteins at an energetic and thermodynamic level. These examples include the simulations on family of aquaporin water channels and glycerol uptake facilitator (GlpF) [59-63]. Extensive MD simulations together with the umbrella sampling technique were done on aquaporin 1 (AQP1) and GlpF to compare the different selectivity mechanism of the two channels.
Simulations on aquaporin 0 (AQP0) revealed the structural basis for its much slower water transportation than other aquaporin channels. Other examples are simulations on the potassium channel [64-65], the ClC transporter [66-67], the c ring of ATPase [68-69], the sodium/proton antiporter [57], the ABC transporters [70-71] etc.
3.8 The aim of the project
The aim of the project is to understand how the proton transfer works in a way that powers the substrate efflux. Further, we would like to propose a working model for AcrB based on the structural analysis of the recently solved crystal structures.
Determining the protonation states for titratable residues is a way to assign protons in different states so that we can figure out how many protons are transferred during drug transportation. Therefore, we apply the Poisson-Boltzmann calculations to each protomer in chapter 5 to predict the pKa values for titratable residues in the transmembrane domain. In chapter 9 we revisit the protonation states for the triad residues by using alchemical free energy perturbations after having a long MD simulation in which we obtain an equilibrated system of AcrB.
After having a clear picture of the protonation states for the triad residues, we advance to explore the water wires in different protomers which are suggested to harbor different number of protons based on our pKa predictions. In chapter 7 we perform extensive conventional MD simulations to see different water wires which are thought to be the media for proton translocation.
In chapter 10 we explore the state transitions to which the trimer AcrB is likely to go or unlikely to go based on our structural analyses. We also find a unique structural element in the transmembrane domain where there are two repeats moving as two independent rigid bodies. Summarized with all the information we get from our simulations and also structural analysis, we have a global working hypothesis for AcrB in the end of the chapter 11.
4. Theories and Methods
This chapter includes the methodological framework that is used throughout this project. We start with the example of receptor-ligand binding to step into the field of thermodynamics. The basic concepts of thermodynamics, statistical mechanics, MD simulation, Poisson-Boltzmann equation, GRIFFIN and free energy perturbation methods are briefly described, while a description of more specific methods will be provided inside each chapter.
4.1 Thermodynamics
Thermodynamics is originally an experimentally based concept when it was developed in the 1800s [72-73]. Like other fundamental laws, thermodynamics is also embodied many great scientists [74] - Otto von Guericke, Robert Boyle, Robert Hooke, Sadi Carnot, Albert Einstein, James Clerk Maxwell, Ludwig Boltzmann, Max Planck, Rudolf Clausius and J. Willard Gibbs, to name just a few. According to their different basis or different types of systems, there are three branches in thermodynamics namely classical thermodynamics, statistical thermodynamics and chemical thermodynamics. Classical thermodynamics describes the states and processes of thermodynamic systems, using macroscopic properties like temperature or pressure which can be directly measured in the laboratory. A microscopic description of these can be provided by statistical thermodynamics which supplement classical thermodynamics with the development of atomic and molecular theories.
The third type of thermodynamics, chemical thermodynamics, is the description of energy relationships in chemical reactions.
Now let us consider the practical example of receptor-ligand binding, which is a very basic phenomenon related to enzymes, cell adhesion, immune response and many other biological functions [75-78]. The term binding affinity is a measurement to characterize the strength with which the ligand binds to the receptor. Higher affinity of a ligand normally corresponds to a greater tendency to bind to the receptor relative
to its dissociation from the receptor. In the laboratory, ligand binding is often characterized in terms of concentration of ligand at which half of the binding sites are occupied. In the microscopic world, there is an equilibrium between ligand binding and unbinding. The measurement of how many ligands bind to the receptor binding sites is actually a probability of the chance for the ligands to be bound to the binding sites. After we get this probability from experiments, then we can find out how many ligands bind to the binding sites which will be the concentration of the receptor-ligand complex. The concentration of the complex then could be used to calculate the dissociation constant which is normally reported as Kd in experiments. This dissociation constant, the macroscopic property measured in experiments, has a microscopic interpretation connecting to the free energy of ligand binding by means of statistical thermodynamics, which will be covered in the next section.
The measurement of pKa [79-81] on amino acids in proteins is another example which is very similar to the case of receptor-ligand binding. In this case, proton is the ligand binding to single amino acid. Similarly, the probability of how big is the chance for the proton to bind one single amino acid can be measured by experimental methods. As not all the measurements of pKa can be done in pure experiments, statistical mechanics opens ways to accomplish this task by employing the framework of thermodynamics.
4.2 Statistical mechanics
4.2.1 Introduction
Statistical mechanics provides a connection between the microscopic and macroscopic world [72]. In order to formulate this connection, we need to introduce the “ensemble” concept which describes a collection of all accessible microstates for some given macroscopic conditions (e.g. the volume, density, temperature and pressure). Averages performed over an ensemble yield many thermodynamics properties of the system and also many other equilibrium and dynamical behaviors.
As we want to calculate the averages on an ensemble, the next question is how we can do the sampling on an ensemble so that we can calculate those thermodynamical and dynamical properties. One way to accomplish this is the Monte Carlo simulation which relies on generation of random numbers for numerical experiments. Another method is the MD simulation, which is the main method we use in this thesis to explore the phase space, and which will be covered in one of the next few sections.
The free energy is a particularly important quantity in statistical mechanics as we can compare the results obtained from computations with the results from experiments. Normally we are less interested in absolute free energies than the free-energy differences between different thermodynamic states. The difference will tell us which state is more probable so that we know the possible protonation state of one specific titratable group, one major part of the thesis, or we know which ligand binds tighter namely a higher binding affinity in the case of receptor ligand binding.
In the next few sections, I will cover the basic and fundamental concepts and techniques, including ensemble, Monte Carlo simulations, MD simulations, force field and methods for Poisson-Boltzmann calculations and alchemical free energy perturbation.
4.2.2 Ensemble
As mentioned earlier, ensemble is a collection of systems which share the same set of microscopic interactions and macroscopic quantities. Each system represents a possible state that the real system might be in. So the ensemble of one specific system can be seen as a collection of points in the phase space. We can always connect a macroscopic observable directly to a microscopic function of coordinates and the momenta of the system according to the theory of classical mechanics [82-83]. For an ensemble with N members, the connection between a macroscopic equilibrium observable A and its microscopic function a(x) can be set up by an averaging procedure
〉
〈
≡
=
∑
=
a x N a
A
N
i i 1
) 1 (
. [1]
So A can be obtained by carrying out the average 〈a〉 over the ensemble.
There are several very important equilibrium ensembles according to their different macroscopic environmental constraints [72]. Microcanonical ensemble or NVE ensemble describes the thermodynamic properties of an isolated system with constant total energy, volume and number of particles. Canonical ensemble or NVT ensemble represents systems with constant volume, temperature and number of particles. Compared with NVE ensemble, NVT ensemble allows the system to exchange energy with its reservoir. Grand canonical ensemble also allows the energy exchange with its reservoir, but apart from energy exchange, the number of particles may also be changed. Another ensemble called Isothermal–isobaric ensemble or NPT ensemble plays a critical role in chemistry as chemical reactions are normally carried out under constant pressure. This ensemble is more practical as most experiments are performed under a standard, constant temperature and pressure.
4.2.3 Monte Carlo simulation
Originally, the Monte Carlo techniques are based on games of chance which would generate the solutions to particular problems when played many times. In sampling a specific ensemble, Monte Carlo sampling relies on trial moves. We take NVT ensemble as an example to see how Monte Carlo sampling works in general. The partition function of a system which has N particles with coordinates X(r1,…, rN) and potential energy U(r1,…, rN) is given by
∫
⋅⋅⋅ −= 3 1 (1,..., )
! ) 1 , ,
( N U r rN
N dr dr e
T N V N
Q β
λ (λ = βh2 2πm) [2]
.We need a criterion to judge whether one random trial move X' is accepted or not. This criterion is the acceptance probability. In this case it is simply
] ,
1 min[
)
|'
(X X e [U(X') U(X)]
A = −β − . [3]
So we can say the acceptance probability is determined uniquely by the change in the
potential energy resulting from the trial move. The trial move will be accepted with probability 1 if the potential energy decreases. If the potential energy increases, the move will be accepted with a probability which decreases exponentially with the change in potential energy.
A common example of trial move is particle displacement. A trial move starts a random selection of a particle, followed by a displacement operation on the Cartesian coordinates of the particle. The question here is, how do we select the appropriate value for the displacement ∆R? If it is too small, one may not explore the phase space at given a reasonable amount of computational time. If it is too big, a trial move will be rarely accepted according to the acceptance probability. Normally a rule of thumb is to set the displacement so that the acceptance probability of a new trial move is around 30%.
4.2.4 Molecular Dynamics simulation
Although the classical dynamics was originally applied to study the motion of planets, stars and other objects in a very large scale, it turns out to be remarkably successful at the molecular level, especially for the case where the computational methodology
“molecular dynamics” is applied to microscopic systems, such as the protein and nucleic acid in Biology for instance. In this section we will start Newton laws of motion which serve as the fundamental laws in molecular dynamics theories.
Newton’s laws of motion
Newton’s laws of motion describe the motion of an object in space at time t by giving the Cartesian position vector r(t), namely
)).
( ), ( ), ( ( )
(t x t y t z t
r = [4]
While v(t) is the first time derivative of the position and the acceleration a(t) is the first time derivative of the velocity, we can derive the acceleration which is actually the second time derivative of the position, here overdot notation is used for
differentiation with respect to time
= •
=dr dt r t
v( )
•
=•
=d r dt r t
a( ) 2 2 . [5]
Therefore, Newton’s second law can be expressed as
•
= •
= mr
dt r md
F 2
2
. [6]
Given a set of initial conditions, we could solve the equations of motion, normally numerically because of the complexity of the system, to describe the motion of the object in space against time. Nowadays, with the power of computer and high performance computing, we can design algorithms to solve the equation of motion numerically. Those algorithms could strongly influence the quality of calculation and sampling to obtain reliable averages. Here I will introduce one basic and simple example algorithm called Verlet algorithm [84]. In this algorithm, the position of one particle at time t+∆t is expressed as a function of the position, velocity and acceleration at time t
) 2
2 ( ) 1 ( ) ( )
(t t r t r t t r t t
r +∆ = + • ∆ + •• ∆ [7]
which is a Taylor expansion without all terms higher than second order in ∆t. Since )
( ) (t v t
r• = and r(••t)=F(t) m, the equation [7] could be rewritten as )
2 ( ) ( ) ( ) (
2
t m F t t t v t r t t
r ∆
+
∆ +
=
∆
+ . [8]
Following this scheme, we can obtain the position of the object at time t-∆t and derive the velocity of the particle at time t
) 2 ( ) ( ) ( ) (
2
t m F t t t v t r t t
r ∆
+
∆
−
=
∆
− [9]
) ( )
( ) ( 2 ) (
) ( )
( 2 ) ( ) (
2 2
t m F t t t r t r t t r
t m F t t r t t r t t r
+∆
∆
−
−
=
∆ +
+∆
=
∆
− +
∆ +
t t t r t t t r
v ∆
∆
−
−
∆
= +
2
) ( ) ) (
( . [10]
In this algorithm we could not obtain the velocity and position of the object at the same time. This issue is not a problem as it is already solved by many other algorithms developed in the last century [85-86]. So I will proceed to talk about several practical issues in MD simulations.
Assignment of initial conditions for MD simulations
As this thesis is devoted to investigate a membrane protein system, we focused this practical part on biological systems. For a biological system, the initial positions of all the atoms are obtained from the coordinates of experimentally determined structures deposited in Protein Data Bank. It might be necessary to add missing information like missing loops in proteins and in most cases the hydrogen atoms which normally can not be determined by experiments. It is also often necessary to solvate the macromolecules into a bath of water. Initial velocities are assigned generally by Maxwell–Boltzmann distribution at a specific temperature chosen for the MD simulations.
Once we specify the initial conditions, all the information is available to start a MD simulation.
Force Field
In molecular modeling, a force field refers to forms and parameters of mathematical functions which are used to describe the potential energy of a system of particles, typically molecules and atoms. A general form for total energy can be expressed as
vdw tic electrosta nonbonded
dihedral angle
bond bonded
nonbonded bonded
total
E E
E
E E
E E
E E
E
+
=
+ +
=
+
=
. [11]
Additionally, “improper torsion” term which maintains the planarity of aromatic rings for example, is also included for some force fields. The most expensive and intensive
computations are the force evaluations for the non-bonded term since it includes many more interactions per atom. A typical form of potential energy employed by CHARMM force field is
∑
∑
∑
∑
− + − + + − + −=
impropers imp dihedrals
angles bonds
b b b K K n K
K r
U( ) ( 0)2 θ(θ θ0)2 χ[1 cos( χ δ)] (ϕ ϕ0)2
ij l
j i
nobonds ij ij
UB
UB r
q q r
R r
R S
S
K ij ij
ε − + ε
+
−
+
∑
( 0)2∑
[( min )12 ( min )6] . [12]In this type of force field, 6 different terms are included in the potential energy.
Kb, Kθ, Kχ, Kimp and KUB are the force constants for bond, angle, dihedral, improper and Urey-Bradley (additional potential function for angle in order to restrain the motions of bonds involved in the angles). b, θ, χ, φ, S are the bond length, angle, dihedral angle, improper angle and Urey-Bradley 1-3 distance. The reference values are the values with subscript 0. Coulomb and Lennard-Jones 6-12 terms contribute to the non-bonded interactions with εfor the Lennard-Jones well depth and Rij for the distance between atom i and j. q is the partial charge and εl is the effective dielectric constant.
A typical parameter set, which is derived from experimental gas-phase geometries, vibrational spectra and high level quantum mechanical calculations [87-90], includes values for atomic mass, van der Waals radius and partial charge for single atom, equilibrium values for bond length, bond angle, dihedral angle and the corresponding spring constant for each potential energy.
Thermostat and barostat
Compared with NVE ensemble which maintains the total energy of the system constant, NVT ensemble is more practical as many experiments are performed under constant temperature instead of constant energy. In this case, the total energy of the system is not constant but fluctuates so that a Boltzmann distribution exp[-βH(r,p)]
can be generated. In order to maintain the temperature of the system, we need to mimic a thermal reservoir so that the system can communicate or exchange the energy
with this thermal reservoir. If the total energy goes higher, this thermal reservoir will take away the extra heat from the system. So actually the total energy of the system always fluctuates around a fixed temperature when the system reaches the equilibrium.
Of all the methods that are designed to control the temperature, by far the most popular ones are the “extended phase space” methods. These techniques supplement the physical phase space by additional variables that mimic the effect of a heat bath.
Most popular thermostats used in molecular dynamics are the Nosé–Hoover thermostat [72, 82, 91], the Berendsen thermostat [72, 92] and the Langevin dynamics [72].
Ensemble NPT is another practical and very important ensemble in which the pressure and temperature of the system are conserved since many reported physical properties including enthalpies, entropies and free energies of ligand binding were measured under constant temperature and pressure. In this scheme, the volume of the system is allowed to fluctuate since the number of particles, temperature and pressure are all conserved. In order to maintain a constant pressure for the system, we would mimic a similar effect of a piston that compresses or expands in response to the internal pressure fluctuation of the system so that the internal pressure would fluctuate around a specific external applied pressure. These popular barostats used in molecular dynamics include Berendsen barostat, Nosé–Hoover barostat and Langevin dynamics.
In this current thesis, the MD program we used, NAMD, is able to perform temperature coupling and pressure control, in which the Langevin equation is used to generate the Boltzmann distribution for canonical (NVT) ensemble simulations [93-94].
4.2.5 Free energy calculations
The free energy of thermodynamics describes the amount of work that a thermodynamic system can perform. It is a thermodynamic state function described by statistical mechanics. Compared with absolute free energy of the system, we are
more interested in free energy changes when the system goes from one thermodynamic state to another. In this section, we will introduce two popular methods which are more related to the current thesis in which we predicted the free energy difference between two states, namely different protonation states of the titratable residues.
The partition function is a quantity which is centered in the field of equilibrium statistical mechanics. It measures the number of microstates in phase space accessible for a given ensemble. Each ensemble has a particular partition function related to the macroscopic variables that describe the ensemble. As the free energy is the logarithm of the partition function, we will take the NVT ensemble as an example to illustrate the two methodologies.
The first method we want to introduce here is the free energy perturbation.
Considering the transformation between thermodynamic state A to thermodynamic state B, normally we characterize these two states in terms of their own potential energy function UA(r1,…,rN) and UB(r1,…,rN). In a case of acid dissociation, state A could be the charged state while state B would be the protonated state. The free energy difference between the two states is simply defined as ∆AAB=AB-AA. Each of the free energies is given by their own canonical partition function Q, free energy
Q kT
A=− ln while partition functions Q for the two states are
N A N
i
N A
i i N
N N
A N
T V N r Z
r m U
r p pd d C
Q 3
1
1 2
! ) , , ]} (
) ,..., 2 (
[
exp{−β + = λ
∫
=
∑
=
[13]
N B N
i
N B
i N i
N N
B N
T V N r Z
r m U
r p pd d C
Q 3
1
1 2
! ) , , ]} (
) ,..., 2 (
[
exp{−β + = λ
∫
=
∑
=
. [14]
Then the free energy difference is
A B A
B A
B
AB Z
kT Z Q
kT Q A
A
A = − =− ln =− ln
∆ . [15]
ZA and ZB are the configurational partition function
) ,...,
(1 N
A r r U N
A d re
Z =∫ −β , ZB =∫dNr e−βUB(r1,...,rN). [16]
In practice it will be easier to calculate the free energy difference if the above
equations could be expressed as terms related to the ensemble averages which we could directly obtain from the sampling either by Monte Carlo simulation or by Molecular Dynamics simulation. So we need to reformulate our configurational partition function as
) ,...,
(1 N
B r r N U
B d re
Z =∫ −β
) ,..., ( )) ,..., ( ) ,..., ( (
) ,..., ( ) ,..., ( ) ,..., (
1 1
1
1 1
1
N A N A N B
N A N A N B
r r U r r U r r N U
r r U r r U r r U N
e e
r d
e e
e r d
β β
β β
β
−
−
−
−
−
∫
=
∫
= . [17]
Now the ratio of ZB ZA becomes
) ,..., ( )) ,..., ( ) ,..., (
( 1 1 1
1 N UB r rN UA r rN UAr rN
A A
B d re e
Z Z
Z −β − −β
∫
= . [18]
A r r U r r
UB N A N
e 〉
〈
= −β( (1,..., )− (1,..., )
The angle bracket above denotes the average taken with respect to the canonical configurational distribution of state A. Now we can get the equation to calculate the free energy difference practically
A U U AB
A
e B
kT
A =− 〈 〉
∆ ln −β( − ) . [19]
The above equation, which is known as the free energy perturbation formula, describes a practical way to calculate the free energy difference between state A and state B. While this is easily implemented in practice, it will be challenging if these two states do not have enough overlap in configurational space. If there is little overlap, the free energy difference between the two states will be large; meaning that UB-UA is big and e-β(UB-UA) becomes very small so that most of the configurations have low weight in the ensemble average. This will introduce a problem of convergence which can be solved by introducing intermediate states between the two transforming states. In practice we divide the perturbation into a series of smaller windows so that the calculations could converge very well. We define these smaller windows as coupling parameter λ which is carefully selected ranging from 0 and 1.
State A corresponds to λ=0 and state B corresponds to λ=1. The free energy perturbation formula after introducing the coupling parameter is then
