•Overview of the development of research in natural and artificial photosynthesis.•Relation between structure and protonation and redox states of Mn-complexes.•Results on pKA computations of the Mn-cluster in photosystem II are reported.•The theoretical background of ab initio pKA computation is summarized.•pKA computations on the hexa-aqua manganese complex from different labs are compared.•The dependencies on atomic radii and charges are discussed in detail.•New results on pKA computations of di-Mn complexes are compared with published data.•Problems of ab initio pKA computations and ways out of it are discussed.We report state-of-the-art pKA computations of manganese model systems. The employed methods can be used to perform pKA computations on the oxygen evolving Mn4CaO5 cluster (OEC) in photosystem II (PSII). These computations are very useful for understanding the function of molecular machines and nano-structures considered for artificial photosynthesis. We first summarize the activities going on in natural and artificial photosynthesis. Next, it is outlined how OEC and manganese model structures can vary with different protonation and oxidation states. We also discuss the role of the dielectric constant of proteins in electrostatic energy computations. Then the results obtained so far on pKA computations of the OEC are summarized. Subsequently, we outline the theoretical background behind performing ab initio pKA value computations. In the application part, we first consider pKA computations for the simple hexa-aqua manganese complex. There it is discussed in detail how computed pKA values depend on atomic radii, charges and the electron leakage effect. A presentation of results follows including a critical discussion of the computed pKA values of di-manganese model complexes. Possible reasons for deviations from measured pKA values and techniques for improvement are discussed and summarized in a final section.Download high-res image (97KB)Download full-size image

We introduce a method that requires only moderate computational effort to compute pKa values of small molecules in different solvents with an average accuracy of better than 0.7 pH units. With a known pKa value in one solvent, the electrostatic transform method computes the pKa value in any other solvent if the proton solvation energy is known in both considered solvents. To apply the electrostatic transform method to a molecule, the electrostatic solvation energies of the protonated and deprotonated molecular species are computed in the two considered solvents using a dielectric continuum to describe the solvent. This is demonstrated for 30 molecules belonging to 10 different molecular families by considering 77 measured pKa values in 4 different solvents: water, acetonitrile, dimethyl sulfoxide, and methanol. The electrostatic transform method can be applied to any other solvent if the proton solvation energy is known. It is exclusively based on physicochemical principles, not using any empirical fetch factors or explicit solvent molecules, to obtain agreement with measured pKa values and is therefore ready to be generalized to other solute molecules and solvents. From the computed pKa values, we obtained relative proton solvation energies, which agree very well with the proton solvation energies computed recently by ab initio methods, and used these energies in the present study.

To speed up the generation of an ensemble of poly(ethylene oxide) (PEO) polymer chains in solution, a tetrahedral lattice model possessing the appropriate bond angles is used. The distance between noncovalently bonded atoms is maintained at realistic values by generating chains with an enhanced degree of self-avoidance by a very efficient Monte Carlo (MC) algorithm. Potential energy parameters characterizing this lattice model are adjusted so as to mimic realistic PEO polymer chains in water simulated by molecular dynamics (MD), which serves as a benchmark. The MD data show that PEO chains have a fractal dimension of about two, in contrast to self-avoiding walk lattice models, which exhibit the fractal dimension of 1.7. The potential energy accounts for a mild hydrophobic effect (HYEF) of PEO and for a proper setting of the distribution between trans and gauche conformers. The potential energy parameters are determined by matching the Flory radius, the radius of gyration, and the fraction of trans torsion angles in the chain. A gratifying result is the excellent agreement of the pair distribution function and the angular correlation for the lattice model with the benchmark distribution. The lattice model allows for the precise computation of the torsional entropy of the chain. The generation of polymer conformations of the adjusted lattice model is at least 2 orders of magnitude more efficient than MD simulations of the PEO chain in explicit water. This method of generating chain conformations on a tetrahedral lattice can also be applied to other types of polymers with appropriate adjustment of the potential energy function. The efficient MC algorithm for generating chain conformations on a tetrahedral lattice is available for download at https://github.com/Roulattice/Roulattice.

Structural data of the oxygen-evolving complex (OEC) in photosystem II (PSII) determined by X-ray crystallography, quantum chemistry (QC), and extended X-ray absorption fine structure (EXAFS) analyses are presently inconsistent. Therefore, a detailed study of what information can be gained about the OEC through a comparison of QC and crystallographic structure information combined with the information from range-extended EXAFS spectra was undertaken. An analysis for determining the precision of the atomic coordinates of the OEC by QC is carried out. OEC model structures based on crystallographic data that are obtained by QC from different research groups are compared with one another and with structures obtained by high-resolution crystallography. The theory of EXAFS spectra is summarized, and the application of EXAFS spectra to the experimental determination of the structure of the OEC is detailed. We discriminate three types of parameters entering the formula for the EXAFS spectrum: (1) model-independent, predefined, and fixed; (2) model-dependent that can be computed or adjusted; and (3) model-dependent that must be adjusted. The information content of EXAFS spectra is estimated and is related to the precision of atomic coordinates and resolution power to discriminate different atom-pair distances of the OEC. It is demonstrated how a precise adjustment of atomic coordinates can yield a nearly perfect representation of the experimental OEC EXAFS spectrum, but at the expense of overfitting and losing the knowledge of the initial OEC model structure. Introducing a novel type of penalty function, it is shown that moderate adjustment of atomic coordinates to the EXAFS spectrum limited by constraints avoids overfitting and can be used to validate different OEC model structures. This technique is used to identify the OEC model structures whose computed OEC EXAFS spectra agree best with the measured spectrum. In this way, the most likely S-state and protonation pattern of the OEC for the most recent high-resolution crystal structure of PSII are determined. We find that the X-ray free-electron laser (XFEL) structure is indeed not significantly affected by exposure to XFEL pulses and thus results in a radiation-damage-free model of the OEC.

For a benchmark set of 194 measured pKa values in 13 proteins, electrostatic energy computations are performed in which pKa values are computed by solving the Poisson–Boltzmann equation. In contrast to the previous approach of Karlsberg+ (KB+) that essentially used protein crystal structures with variations in their side chain conformations, the present approach (KB2+MD) uses protein conformations from four molecular dynamics (MD) simulations of 10 ns each. These MD simulations are performed with different specific but fixed protonation patterns, selected to sample the conformational space for the different protonation patterns faithfully. The root-mean-square deviation between computed and measured pKa values (pKa RMSD) is shown to be reduced from 1.17 pH units using KB+ to 0.96 pH units using KB2+MD. The pKa RMSD can be further reduced to 0.79 pH units, if each conformation is energy-minimized with a dielectric constant of εmin = 4 prior to calculating the electrostatic energy. The electrostatic energy expressions upon which the computations are based have been reformulated such that they do not involve terms that mix protein and solvent environment contributions and no thermodynamic cycle is needed. As a consequence, conformations of the titratable residues can be treated independently in the protein and solvent environments. In addition, the energy terms used here avoid the so-called intrinsic pKa and can therefore be interpreted without reference to arbitrary protonation states and conformations.

ProPairs is a data set of crystal structures of protein complexes defined as biological assemblies in the protein data bank (PDB), which are classified as legitimate protein–protein docking complexes by also identifying the corresponding unbound protein structures in the PDB. The underlying program selecting suitable protein complexes, also called ProPairs, is an automated method to extract structures of legitimate protein docking complexes and their unbound partner proteins from the PDB which fulfill specific criteria. In this way a total of 5,642 protein complexes have been identified with 11,600 different decompositions in unbound protein pairs yielding legitimate protein docking partners. After removing sequence redundancy (requiring a sequence identity of the residues in the interface of less than 40%), 2,070 different legitimate protein docking complexes remain. For 810 of these protein docking complexes, both docking partners possess corresponding unbound structures in the PDB. From the 2,070 nonredundant protein docking complexes there are 417 which possess a cofactor at the interface. From the 176 protein docking complexes of the Protein–Protein Docking Benchmark 4.0 (DB4.0) data set, 13 differ from the ProPairs data set. Twelve of them differ with respect to the composition of the unbound structures but are contained in the large redundant ProPairs data set. One protein docking complex of the DB4.0 data set is not contained in ProPairs since the biological assembly specified in the PDB is wrong (PDB id 1d6r). For one protein complex (PDB id 1bgx) the DB4.0 data set uses a fabricated unbound structure. For public use interactive online access is provided to the ProPairs data set of nonredundant protein docking complexes along with the source code of the underlying method [http://propairs.github.io].

Biomolecular systems are commonly exposed to a manifold of forces, often acting between multivalent ligands. To understand these forces, we studied mono- and bivalent model systems of pyridine coordination complexes with Cu2+ and Zn2+ in aqueous environment by means of scanning force microscopy based single-molecule force spectroscopy in combination with ab initio DFT calculations. The monovalent interactions show remarkably long rupture lengths of approximately 3 Å that we attribute to a dissociation mechanism involving a hydrogen-bound intermediate state. The bivalent interaction with copper dissociates also via hydrogen-bound intermediates, leading to an even longer rupture length between 5 and 6 Å. Although the bivalent system is thermally more stable, the most probable rupture forces of both systems are similar over the range of measured loading rates. Our results prove that already in small model systems the dissociation mechanism strongly affects the mechanical stability. The presented approach offers the opportunity to study the force-reducing effects also as a function of different backbone properties.

We present a robust method for the calculation of electrostatic potentials of large molecular systems using tetrahedral finite elements (FE). Compared to the finite difference (FD) method using a regular simple cubic grid to solve the Poisson equation, the FE method can reach high accuracy and efficiency using an adaptive grid. Here, the grid points can be adjusted and are placed directly on the molecular surfaces to faithfully model surfaces and volumes. The grid point density decreases rapidly toward the asymptotic boundary to reach very large distances with just a few more grid points. A broad set of tools are applied to make the grid more regular and thus provide a more stable linear equation system, while reducing the number of grid points without compromising accuracy. The latter reduces the number of unknowns significantly and yields shorter solver execution times. The accuracy is further enhanced by using second order polynomials as shape functions. Generating the adaptive grid for a molecular system is expensive, but it pays off, if the same molecular geometry is used several times as is the case for pKA and redox potential computations of many charge variable groups in proteins. Application of the mFES method is also advantageous, if the molecular system is too large to reach sufficient accuracy when computing the electrostatic potential with conventional FD methods. The program mFES is free of charge and available at http://agknapp.chemie.fu-berlin.de/mfes.

A first step toward three-dimensional protein structure description is the characterization of secondary structure. The most widely used program for secondary structure assignment remains DSSP, introduced in 1983, with currently more than 400 citations per year. DSSP output is in a one-letter representation, where much of the information on DSSP’s internal description is lost. Recently it became evident that DSSP overlooks most π-helical structures, which are more prevalent and important than anticipated before. We introduce an alternative concept, representing the internal structure characterization of DSSP as an eight-character string that is human-interpretable and easy to parse by software. We demonstrate how our protein secondary structure characterization (PSSC) code allows for inspection of complicated structural features. It recognizes ten times more π-helical residues than does the standard DSSP. The plausibility of introduced changes in interpreting DSSP information is demonstrated by better clustering of secondary structures in (φ, ψ) dihedral angle space. With a sliding sequence window (SSW), helical assignments with PSSC remain invariant compared with an assignment based on the complete structure. In contrast, assignment with DSSP can be changed by residues in the neighborhood that are in fact not interacting with the residue under consideration. We demonstrate how one can easily define new secondary structure classification schemes with PSSC and perform the classifications. Our approach works without changing the DSSP source code and allows for more detailed protein characterization.

Cytochrome c oxidase (CcO) is a central enzyme in aerobic life catalyzing the conversion of molecular oxygen to water and utilizing the chemical energy to pump protons and establish an electrochemical gradient. Despite intense research, it is not understood how CcO achieves unidirectional proton transport and avoids short circuiting the proton pump. Within this work, we analyzed the potential role of Glu286 as a proton valve. We performed unconstrained MD simulations of CcO with an explicit membrane for up to 80 ns. Those MD simulations revealed that deprotonated Glu286 (Glu286-) is repelled by the negatively charged propionic acid PRD of heme a3. Thus, it destabilizes a potential linear chain of waters in the hydrophobic cavity connecting Glu286 with PRD and the binuclear center (BNC). Conversely, protonated Glu286 (Glu286H) may remain in an upward position (oriented toward PRD) and can stabilize the connecting linear water chain in the hydrophobic cavity. We calculated the pKa of Glu286 under physiological conditions to be above 12, but this value decreases to about 9 under increased water accessibility of Glu286. The latter value is in accordance with experimental measurements. In the time course of MD simulation, we also observed conformations where Glu286 bridges between water molecules located on both sides (the D channel being connected to the N side and the hydrophobic cavity), which might lead to proton backflow.

Higher organisms defend themselves against invading micro-organisms and harmful substances with their immune system. Key players of the immune system are the white blood cells (WBC), which in case of infection move in an extravasation process from blood vessels toward infected tissue promoting inflammation. This process starts with the attachment of the WBC to the blood vessel wall, mediated by protein pair interactions of selectins and counter-receptors (C-R). Individual selectin C-R binding is weak and varies only moderately between the three selectin types. Multivalency enhances such small differences, rendering selectin-binding type specific. In this work, we study selectin C-R binding, the initial step of extravasation. We performed electrostatic energy computations based on the crystal structure of one selectin type co-crystallized with the ligating part of the C-R. The agreement with measured free energies of binding is satisfactory. Additionally, we modeled selectin mutant structures in order to explain differences in binding of the different selectin types. To verify our modeling procedures, surface plasmon resonance data were measured for several mutants and compared with computed binding affinities. Binding affinities computed with soaked rather than co-crystallized selectin C-R structures do not agree with measured data. Hence, these structures are inappropriate to describe the binding mode. The analysis of selectin/C-R binding unravels the role played by individual molecular components in the binding event. This opens new avenues to prevent immune system malfunction, designing drugs that can control inflammatory processes by moderating selectin C-R binding.

Extensive quantum chemical DFT calculations were performed on the high-resolution (1.9 Å) crystal structure of photosystem II in order to determine the protonation pattern and the oxidation states of the oxygen-evolving Mn cluster. First, our data suggest that the experimental structure is not in the S1-state. Second, a rather complete set of possible protonation patterns is studied, resulting in very few alternative protonation patterns whose relevance is discussed. Finally, we show that the experimental structure is a mixture of states containing highly reduced forms, with the largest contribution (almost 60%) from the S–3-state, Mn(II,II,III,III).(1)

The mutual information (MI) expansion is applied to two molecular systems to probe algorithms that serve to estimate conformational entropy differences more precisely. The individual terms of the MI expansion are evaluated with a histogram method. Internal coordinates are used to avoid spurious correlations, which would require higher order terms in the MI expansion. Two approaches are applied that compensate for systematic errors that occur with a histogram method: (1) Simulation data are balanced by using the same number of coordinate sets (frames) for both conformer domains considered for the entropy difference computation. Balancing puts fluctuations of the histogram bin contents on the same level for both conformer domains, allowing efficient error cancellation. (2) Bias correction compensates for systematic deviations due to a finite number of frames per bin. Applying both corrections improves the precision of entropy differences drastically. Estimates of entropy differences are compared to thermodynamic benchmarks of a simple polymer model and trialanine, where excellent agreement was found. For trialanine, the average error for the estimated conformational entropy difference is only 0.3 J/(mol K), which is 100 times smaller than without applying the two corrections. Guidelines are provided for efficiently estimating conformational entropies. The program ENTROPICAL, used for the computations, is made available, which can be used for molecular dynamics or Monte Carlo simulation data on macromolecules like oligopeptides, polymers, proteins, and ligands.

A first step toward predicting the structure of a protein is to determine its secondary structure. The secondary structure information is generally used as starting point to solve protein crystal structures. In the present study, a machine learning approach based on a complete set of two-class scoring functions was used. Such functions discriminate between two specific structural classes or between a single specific class and the rest. The approach uses a hierarchical scheme of scoring functions and a neural network. The parameters are determined by optimizing the recall of learning data. Quality control is performed by predicting separate independent test data. A first set of scoring functions is trained to correlate the secondary structures of residues with profiles of sequence windows of width 15, centered at these residues. The sequence profiles are obtained by multiple sequence alignment with PSI-BLAST. A second set of scoring functions is trained to correlate the secondary structures of the center residues with the secondary structures of all other residues in the sequence windows used in the first step. Finally, a neural network is trained using the results from the second set of scoring functions as input to make a decision on the secondary structure class of the residue in the center of the sequence window. Here, we consider the three-class problem of helix, strand, and other secondary structures. The corresponding prediction scheme “SPARROW” was trained with the ASTRAL40 database, which contains protein domain structures with less than 40% sequence identity. The secondary structures were determined with DSSP. In a loose assignment, the helix class contains all DSSP helix types (α, 3–10, π), the strand class contains β-strand and β-bridge, and the third class contains the other structures. In a tight assignment, the helix and strand classes contain only α-helix and β-strand classes, respectively. A 10-fold cross validation showed less than 0.8% deviation in the fraction of correct structure assignments between true prediction and recall of data used for training. Using sequences of 140,000 residues as a test data set, 80.46% ± 0.35% of secondary structures are predicted correctly in the loose assignment, a prediction performance, which is very close to the best results in the field. Most applications are done with the loose assignment. However, the tight assignment yields 2.25% better prediction performance. With each individual prediction, we also provide a confidence measure providing the probability that the prediction is correct. The SPARROW software can be used and downloaded on the Web page http://agknapp.chemie.fu-berlin.de/sparrow/.

Experiments show that a ligand–receptor complex formed by binding a bivalent ligand (D) in which the two ligating units are joined covalently by a flexible polymeric spacer (S) can be orders of magnitude more stable than the corresponding complex formed with monomeric ligands. Although molecular models rationalizing this “enhancement effect” have been proffered, they ignore spacer–receptor (S–R) interactions, which can substantially influence the relative stability of complexes. Here, the results of a computational study designed to assess the impact of S–R interactions in the prototypic bivalent complex are presented and compared to results of experiments. The S–R interactions mimicking general features of biological systems are modeled by contoured R surfaces with hills (or depressions) at the binding sites. In the fictitious limit of vanishing S–R interactions, the enhancement is pronounced, as observed in experiments. For strictly repulsive S–R interactions (hard R surface), the enhancement vanishes, or even reverses. This is particularly the case if the R surface is convex (i.e., rising between the binding sites), while the enhancement is only moderately reduced if the R surface is concave. Alternatively, a weak S–R attraction close to the R surface can increase the enhancement. It is concluded that large enhancement should be observed only if both features are present: a concave R surface plus a weak S–R attraction. The latter occurs for spacer material such as polyethylene glycol (PEG), which is weakly hydrophobic and thus attracted by protein surfaces. It is shown that the enhancement of bivalent binding can be characterized by a single key parameter, which may also provide guidelines for the design of multivalent complexes with large enhancement effect.

Continuum electrostatic theory was applied to compute redox potentials of rubredoxin (Rd) proteins. We used multiple side chain conformers of Rd crystal structures, optimized geometries of salt bridges, mutated residues, and residues in the neighborhood of the iron−sulfur complex (FeS complex) self-consistently for given solvent pH and redox potential. The following contributions to Rd redox potentials are discussed: side chain conformations, H-bond geometries of the FeS complex, dielectric environment, charged residues, and salt bridges. We considered 15 different Rd's (of different species/strains and mutants) with available crystal structures whose redox potentials vary between −86 mV and +31 mV. The computed redox potentials deviated by less than 16 mV, root-mean-square deviation (RMSD), from measured values. The amide H-bond geometry is considered to be crucial for the variation of Rd redox potentials. To test this assumption, we considered 14 mutant Rd's for which we modeled the structures based on Rd from WT Clostridium pasterianum (Cp) leaving the amide H-bond geometry of the FeS complex invariant. Here, we obtained an RMSD of only 14 mV with measured values demonstrating that the amide H bond geometries cannot be a major factor determining Rd redox potentials. We analyzed the factors determining the Rd redox potentials of a mesophilic and a thermophilic Rd differing by nearly 90 mV. We found that half of the difference is due to sequence and half is due to backbone variations. Albeit salt-bridge networks vary considerably between these two Rd's and are considered to be responsible for differences in thermostability, their overall influence on Rd redox potentials is small.

We constructed an accurate polyether force field for implicit solvent (IS) molecular dynamics (MD) simulations that matches local and global conformations of 1,2-dimethoxy-ethane (DME) and polyethylene glycol (PEG), respectively. To make appropriate force field adjustments for IS models of PEG, we used long-term MD simulation data of 1 μs in explicit solvent (ES) based on the most recent CHARMM35 ether force field that includes adjustments for PEG in explicit water. In IS models, competition of attractive van der Waals (vdW) interactions between solute−solute and solute−solvent atom pairs is often not considered explicitly. As a consequence, the attractive vdW interactions between solute atom pairs that remain in IS models explicitly can yield equilibrium structures that are too compact. This behavior was observed in the present study comparing MD simulation data of the DME and PEG ES model with corresponding IS models that use generalized Born (GB) electrostatics combined with positive surface energy terms favoring compact structures. To regain balance of attractive vdW interactions for IS models, we considered the IS generalized Born with simple switching (GBSW) model in detail, where we turned off surface energy terms and reduced attractive vdW interactions to 90%, or we used alternatively even slightly negative surface energies. However, to obtain quantitatively the same local and global distributions of PEG conformers as in ES, we needed additional force field adjustments involving torsion potentials and 1−4 and 1−5 atom pair Coulomb interactions. This CHARMM ether force field, specifically optimized for IS simulation conditions, is equally valid for dimeric and polymeric ethylene glycol. To explore the conformational space of PEG with MD simulations, an IS GBSW model requires 2 orders of magnitude less CPU time than the corresponding ES model. About a factor of 5 of this gain in efficiency is due to the lack of solvent viscosity in IS models.

A multivalent ligand can bind to a multivalent receptor to form a multivalent complex that is orders of magnitude more stable than its “monovalent” counterpart formed by the binding of monovalent ligands, each of which is equivalent to the corresponding ligating unit of the multivalent ligand. This “enhancement effect” is frequently rationalized qualitatively in terms of the concept of “effective concentration” borrowed from polymer physics. The purpose of this article is to present the details of a fundamental theory of the enhancement effect that was only summarized recently. Classical statistical mechanics is used to compute closed analytic expressions for the binding constants of the various ligand−receptor complexes associated with the prototypal system, divalent ligand binding to a divalent receptor. A simple analytic expression for a defined “enhancement factor” reveals insight into the mechanism of enhancement. The enhancement factor is proportional to a generalized effective concentration of one ligating unit with respect to the other, which depends not only on the separation between the binding sites of the receptor but also on the thermodynamic state of the system. The theory is related to typical experiments and is used to analyze data on the activation of ion channels (tetravalent receptors) in cell membranes by polymer-linked divalent ligands.

The energetics of redox states in different models of rubredoxin-like iron−sulfur complexes (ISC) was computed using a combination of density functional and electrostatic continuum theories. In agreement with experiment, the calculated redox potential for the small ISC model [Fe(SCH2CH3)4]1−,2− in acetonitrile was −813 mV [Galstyan, A. S.; Knapp, E. W. J. Comput. Chem. 2009, 30, 203−211] as compared to the measured value of −838 mV. Surprisingly the experimental values for rubredoxin (Rd) are much higher ranging between −87 and +39 mV. These large variations in redox potentials of ISC models and ISC in Rd are due to specific conformational symmetries adopted by the ligands due to both the protein environment and type and the number of H-bonds, and the dielectric environment. In a dielectric environment corresponding to proteins (ε = 20), the computed ISC redox potentials shift positive by about 64 mV for Fe−S···H−N and 95 mV for Fe−S···H−O H-bonds, correlating well with data estimated from experiments on ISC proteins. In aqueous solutions (ε = 80), a positive shift of 58 mV was computed for Fe−S···H−O H-bonds (using a model with the same ISC conformation as in Rd) in agreement with a measured value for Rd with partially solvent exposed ISC. The latter demonstrates the dependence of the ISC redox potentials on the environment (solvent or protein). For a model whose chemical composition is analog to the relevant part of ISC in a specific Rd, the computed redox potential of the model agrees with the measured value in Rd. This study allows to understand redox potential shifts for small ISC models and ISC in proteins.

In the absence of structural knowledge on the target protein, the bound ligand conformer (BLC) can be constructed approximately by an indirect drug-design approach that uses a set of ligands binding to the same target. Once the bound ligand conformer (BLC) is known, different strategies of drug design can be pursued. The indirect drug-design approach of the present study is based on the assumption that a set of ligands with chemically different architecture binding to the same target protein carry hidden information of their corresponding true BLCs. It is shown how this information can be extracted by pairwise flexible structure alignment (FSA) using molecular dynamics (MD) simulations with attractive intermolecular interactions that derive from the molecular similarity of the ligands and allow the ligands to adopt the same space. The FSA approach is performed with a newly designed module overlap in the experimental CHARMM-29a1, which soon will become publicly available. Combining the conformations obtained from FSA of different ligand pairs yields consensus ligand conformers (CLCs) that should be similar to the BLCs. This procedure was validated on HIV-1 protease (HIV-P), where at present 44 crystal structures with bound ligands of sufficiently diverse chemical composition are available. The FSA approach identifies four different clusters of HIV-P BLCs. These clusters are consistent with the H-bond patterns of the ligands bound to HIV-P in the crystal structures exhibiting four different binding modes. The cluster-specific CLCs are indeed very similar (rmsd ≈ 2 Å) to the corresponding BLCs from the crystal structure, demonstrating the feasibility of the present approach.

Insights in structural biology can be gained by analyzing protein architectures and characterizing their structural similarities. Current computational approaches enable a comparison of a variety of structural and physicochemical properties in protein space. Here we describe the automated detection of rotational symmetries within a representative set of nearly 10 000 nonhomologous protein structures. To find structural symmetries in proteins initially, equivalent pairs of secondary structure elements (SSE), i.e., α-helices and β-strands, are assigned. Thereby, we also allow SSE pairs to be assigned in reverse sequential order. The results highlight that the generation of symmetric, i.e., repetitive, protein structures is one of nature’s major strategies to explore the universe of possible protein folds. This way structurally separated ‘islands’ of protein folds with a significant amount of symmetry were identified. The complete results of the present study are available at http://agknapp.chemie.fu-berlin.de/gplus, where symmetry analysis of new protein structures can also be performed.

Cd2+ binding at the bacterial photosynthetic reaction center (bRC) from Rhodobacter sphaeroides is known to inhibit proton transfer (PT) from bulk solvent to the secondary quinone QB. To elucidate this mechanism, we calculated the pKa for residues along the water channels connecting QB with the stromal side based on the crystal structures of WT-bRC and Cd2+-bound bRC. Upon Cd2+ binding, we observed the release of two protons from His-H126/128 at the Cd2+ binding site and significant pKa shifts for residues along the PT pathways. Remarkably, Asp-L213 near QB, which is proposed to play a significant role in PT, resulted in a decrease in pKa upon Cd2+ binding. The direct electrostatic influence of the Cd2+-positive charge on these pKa shifts was small. Instead, conformational changes of amino acid side chains induced electrostatically by Cd2+ binding were the main mechanism for these pKa shifts. The long-range electrostatic influence over ≈12 Å between Cd2+ and Asp-L213 is likely to originate from a set of Cd2+-induced successive reorientations of side chains (Asp-H124, His-H126, His-H128, Asp-H170, Glu-H173, Asp-M17, and Asp-L210), which propagate along the PT pathways as a “domino” effect.

Numerous selective estrogen receptor modulators (SERMs) have been synthesized and assayed in recent years. The focus of this study is to apply coarse-grain molecular docking procedures coupled with fine-grain all-atom force field optimization strategies to shed light on the binding mechanisms of currently available estrogen receptor-active compounds. Although the mechanics of ligand binding in estrogen receptors is generally well understood, there is room for surprises. In this paper computational evidence corroborating the experimentally observed type I agonistic binding mode for estradiol (E2) and diethylstilbesterol (DES) and the type II antagonistic binding mode for 4-hydroxytamoxifen and raloxifen is presented. Included in this type I agonistic mode are the DES derivatives, transstilbene and 1,2-diaryldiaminoethane. In addition, a novel ‘type II agonistic’ binding mode for 2,3-diarylimidazolines, 4,5-diarylimidazoles, 2,3-diarylpiperazines is introduced. This mode is stabilized by suggesting alternative hydrogen bond anchor points in the ligand binding domain as potential leads for future drug design.Molecular structures of the alkene (1), diamine (2), imidazole (3), imidazoline (4), and piperazine (5) test compounds. Ar is an abbreviation for a substituted benzene ring, for which X can be either chlorine or fluorine, and R is either a hydroxyl or methoxy group.

We present a method for deriving energy functions for protein folding by maximizing the thermodynamic average of the overlap with the native state. The method has been tested by using the pairwise contact approximation of the energy function and generating alternative structures by threading sequences over a database of 1,169 structures. With the derived energy function, most native structures: (i) have minimal energy and (ii) are thermodynamically rather stable, and (iii) the corresponding energy landscapes are smooth. Precisely, 92% of the 1,013 x-ray structures are stabilized. Most failures can be attributed to the neglect of interactions between chains forming polychain proteins and of interactions with cofactors. When these are considered, only nine cases remain unexplained. In contrast, 38% of NMR structures are not assigned properly.

To understand the rapid diffusion mechanism of water molecules in the crystal lattice of hydrated β-cyclodextrin (β-CD), molecular dynamics (MD) simulations of crystalline β-CD were performed at five different relative humidities corresponding to hydration states ranging from β-CD-9.4H2O to β-CD-12.3H2O, and in aqueous solution. The trajectories for the crystalline β-CD hydrates had lengths of 4 ns each, whereas the simulation in aqueous solution extended to 2 ns. Transport of water molecules in the crystal was characterized in terms of a spatially varying diffusion constant and the main direction of diffusion, which were evaluated using newly developed algorithms. The main diffusion pathway winds through the cavities of adjacent doughnut shaped β-CD molecules and is slightly slanted with respect to the crystallographic b-axis. Water molecules outside the β-CD cavities have access to the main diffusion pathway. The diffusion constant for transport of water molecules along the main pathway calculated from the MD simulation data adopts 1/30 of the value in bulk water at room temperature. This is in agreement with estimates that can be made from experimental data on the adjustment of a β-CD crystal to changes in relative humidity.

A key enzyme in aerobic metabolism is cytochrome c oxidase (CcO), which catalyzes the reduction of molecular oxygen to water in the mitochondrial and bacterial membranes. Substrate electrons and protons are taken up from different sides of the membrane and protons are pumped across the membrane, thereby generating an electrochemical gradient. The well-studied A-type CcO uses two different entry channels for protons: the D-channel for all pumped and two consumed protons, and the K-channel for the other two consumed protons. In contrast, the B-type CcO uses only a single proton input channel for all consumed and pumped protons. It has the same location as the A-type K-channel (and thus is named the K-channel analog) without sharing any significant sequence homology. In this study, we performed molecular-dynamics simulations and electrostatic calculations to characterize the K-channel analog in terms of its energetic requirements and functionalities. The function of Glu-15B as a proton sink at the channel entrance is demonstrated by its rotational movement out of the channel when it is deprotonated and by its high pKA value when it points inside the channel. Tyr-244 in the middle of the channel is identified as the valve that ensures unidirectional proton transfer, as it moves inside the hydrogen-bond gap of the K-channel analog only while being deprotonated. The electrostatic energy landscape was calculated for all proton-transfer steps in the K-channel analog, which functions via proton-hole transfer. Overall, the K-channel analog has a very stable geometry without large energy barriers.