Research

We develop modern methods of statistical thermodynamics, implements and apply them in computational biophysics. These methods are based on novel combinations of several different approaches, in particular classical molecular dynamics simulations, efficient quantum chemical calculations, extended-sampling schemes, as well as constructs of machine learning. Such approaches make it possible to elucidate the mechanisms of chemical reactions taking place in complex biomolecular systems which would otherwise be difficult to access due to their low rates or their coupling to conformational transitions.

Development & implementation

We aim to uncover the underlying mechanisms, thermodynamics and kinetics of chemical reactions occurring in biomolecular complexes, possibly in coupling to conformational transitions, with no or little a priori assumptions on the mechanism. The general strategy is to derive the properties of the biomolecular complex or of the process from a free energy surface, generated as a potential of the mean force (PMF) from a molecular simulation. The real challenge is the slow rates of many of the involved processes, so that a naive approach based on free (unbiased) MD simulation would have to cover extremely long time periods. As that would be impossible with any computational resources realistically available, the only practicable procedure seems to be some combination of hybrid QM/MM MD simulations with efficient quantum chemical methods and extended sampling approaches. Each of those additional components is designed to increase the computational efficiency while, at the same time, not compromising the accuracy beyond whatever is deemed acceptable in the relevant applications.

In our specific strategy, these components are

• the hybrid QM/MM approach itself, with the QM region constituting of the active site of the protein and its closest surroundings,
• the semi-empirical density-functional method DFTB3, benchmarked and/or parametrized for the reaction under investigation,
• optionally, a machine-learning calculation to improve the accuracy of the QM calculation or even to replace it altogether,
• an extended-sampling procedure such as metadynamics to make the rare events of interest occur more frequently in simulations while keeping access to unbiased thermodynamics and kinetics, and
• optionally, a trivially parallel variant of that procedure, such as multiple-walker metadynamics, to make better use of computational resources and approach microsecond MD sampling.

We have been implementing our methodical developments in a combination of established open-source software packages for bio/molecular modeling such as Gromacs, Plumed and DFTB+. They were described in two publications [8, 11], and our code is freely available for download and use from a public repository.

Proton transfer & proton-coupled electron transfer

One of the reactions that we focus on is the proton transfer (PT) occurring over long distances in protein complexes, such as those participating in biomolecular energy transduction.

Bacteriorhodopsin

The transmembrane protein bacteriorhodopsin (BR) is an archetypal light-activated proton pump, which uses the energy of an absorbed photon to translocate a proton through the cellular membrane outwards. We aimed to characterize a feature that still remained mechanistically unclear – the final proton transfer (PT) step, in which the side chain of Asp85 passes the proton to the proton-release group (PRG) consisting of a couple of Glu side-chains, bringing the protein from the so-called O-state back to the resting state. This process is exoergic but slow, passing through a detectable intermediate (protonated Asp212) and several water molecules, and is affected by several floppy side-chains on the way. [7]

We refined a structural model for the previously poorly characterized initial state, and identified main structural distinctions between the different states (initial O, intermediate O*, and final, ground state BR). Then, we optimized the simulation setup to make it possible to generate a 3D PMF of the entire long-range PT reaction, involving one collective variable (CV) describing the PT itself, and two CVs describing other motions in the protein. The rate-limiting barrier observed in the PMF corresponds to the experimentally estimated kinetics, and the computed reaction free energy agrees with the experimentally estimated differences of acidity constants, supporting the validity of our results.

We observed that PT in BR proceeds via a proton hole mechanism, as opposed to the textbook-like Grotthuss’ mechanism involving an excess proton. This result is in line (i) with the chloride pumping ability of a single-amino-acid mutant of BR, and (ii) with the similarities of mechanisms of BR and halorhodopsin, pointed out previously. Considering BR as an example of a more general phenomenon, we concluded that an analysis of PT reactions in heterogeneous environments shall in general consider both ionic forms of water. Also, we characterized the structural changes accompanying the long-range PT: an Arg side chain swinging between two distinct orientations and mechanically breaking the path of potential backward PT, the pair of Glu residues forming the PRG approaching closely, and additional water molecules flowing out of the active site.

Proton-coupled electron transfer

We have developing a new concept of free energy calculations for chemical reactions by means of extended-sampling MD simulations that involve biasing potentials applied to electron density [3, 9]. The novel feature is that the biasing potential is applied on a CV constructed from Mulliken atomic charges, which appears to be a convenient CV to describe and bias an electron transfer reaction (ET) in many cases. The simulation may also involve one or more other CVs consisting of atomic coordinates, which may be used to bias a PT.

The calculation of biasing force in such simulations additionally requires the gradients of QM charges to be calculated (which is not the case in unbiased QM/MM MD simulations). These are obtained by solving coupled-perturbed equations within DFTB3, which were originally developed by others to serve different purposes (to calculate polarizabilities or vibrational frequencies). Our own development consists of an additional set of equations to calculate the derivatives of QM charges with respect to the coordinates of MM atoms. We demonstrated that these are necessary to obtain converged free energies in QM/MM situations.

The method was tested on PCET reactions taking place in a model system mimicking two Tyr side-chains, in several different spatial orientations, immersed in a water nano-droplet. Two CVs were introduced naturally, one for the PT the other for the ET. Stable metadynamics simulations were obtained, and yielded PMF that showed the expected topography and energy barriers, in agreement with extended unbiased simulations performed for reference. Hence, this method has the potential to distinguish between concerted and sequential mechanisms of PCET, and can be applied to larger molecular complexes.

To provide reference data for PCET processes taking place in a range of varied biomolecular complexes, we have analyzed PCET occurring between aromatic amino-acid side-chains in several different biomimetic peptides [0]. The PMFs for the reactions were obtained by tracking two CVs along the simulation: the difference of proton donor–proton distances for the PT, and the difference of summed Mulliken charges for the ET. This was performed in two different ways, (i) from unbiased MD simulations as well as (ii) from 1D metadynamics simulations with biasing the PT CV, followed by performing a reweighting procedure. Our simulations show that the geometry, in which the residues are arranged relative to each other, strongly influences the electron transfer, and that an increased solvent exposure affects the PCET mechanism, altering its free energy landscape.

Multi-scale simulation models of reactions in proteins

We participate in two different projects within the framework of the DFG research training group (RTG) 2450 Tailored Scale-Bridging Approaches to Computational Nanoscience. In the first funding period, we have investigated the thiol-disulfide exchange and phosphoryl transfer reactions in selected proteins, which occur in coupling to conformational transitions.

Thiol–disulfide exchange

We started off by investigating the thiol–disulfide exchange reaction in a minimal model system as well as in small peptides [10]. The PMFs generated by way of multiple-walker QM/MM metadynamics provided the heights of free energy barriers opposing the spontaneous exchange. The structure of the transition state as well as the barrier heights were very similar in all cases. Also, the comparison of the different small peptides showed sterically unfavorable disulfide bonds in CXC fragments, while CXXC disulfide bonds were favored over longer-range disulfide bonds along the peptide backbone, corroborating with the high abundance of CXXC motifs in redox proteins.

Then, we analyzed the effect of mechanical stretching stresses on the mutated immunoglobulin domain I27* [6], which was investigated by our collaboration partners in the Gräter lab (HITS Heidelberg) using more simplified methods previously. Specifically, the aim was to uncover the effects of structural factors and electrostatic interactions with the environment on the thiol–disulfide exchange reactions in I27*. We relied on generating an extensive ensemble of unbiased QM/MM MD trajectories for a total sampling of multiple microseconds. This trivially parallel approach was enough to observe a significant number of thiol–disulfide exchanges, and the free thiolate preferred to attack one of the Cys side-chains involved in the disulfide bond (Cys55) over the other, in agreement with previous findings by others.

We found that the reaction is additionally directed by the electrostatic interactions of the involved sulfur atoms with the molecular environment. The relationships of atomic charges, modulated by those electrostatic interactions, leads to the kinetic preference of the attack on Cys55. To quantify these electrostatic effects in more detail, we set up a small model system with varied artificial external electric potentials, and employed our combined QM/MM metadynamics scheme to study the thiol–disulfide exchange in it. The observed changes in reaction kinetics were of the same magnitude as in I27*. This finding confirmed the role that the electrostatic interactions play in the regioselectivity of the thiol–disulfide exchange reactions in the protein.

While the accuracy of the standard parametrization of DFTB3 for the description of thiol–disulfide exchange proved sufficient to make qualitative comparisons such as those above, it turned out to be insufficient for drawing more quantitative conclusions [5]. We worked to alleviate this problem by developing a specific reaction parametrization (SRP) of DFTB3 for the thiol–disulfide exchange, in two different ways: (i) by generating a new S–S repulsive function, and (ii) by creating an artificial neural network (ANN) implementation in the DFTB+ software and training on an extended set of molecular structures relevant for the reaction. We applied both approaches to investigate the thiol–disulfide exchanges in a solvated model system and then in a blood protein. Highly accurate PMFs were obtained efficiently, as the ANN component only brings about a small additional computational cost.

A particularly interesting class of enzymes supporting thiol–disulfide exchange reactions are glutaredoxins, which catalyze the oxidation and reduction of protein disulfide bonds. Their active site has either one or two Cys residues, and that leads to different catalytic reaction cycles, which have been a subject of considerable research interest. We investigated a proposed (but not fully uncovered) mechanism for the reduction of the disulfide bond in the protein HMA4n by a mutated monothiol human glutaredoxin and the co-substrate glutathione, which involves three successive thiol–disulfide exchanges between the three participating molecules [4]. First, we estimated the regioselectivity of the different attacks by analyzing classical MD trajectories. Then, we generated PMFs for each reaction with multiple-walker QM/MM metadynamics, which involved the ANN correction to improve the accuracy. We showed the same regiospecificity as observed in the experiment, and the obtained barrier heights for the different reaction steps confirmed the proposed pathway.

Phosphoryl transfer

We investigated the auto-phosphorylation of histidine kinases (HK) – a complex process that consists of a large transition from an inactive to an active conformation, followed by an actual phosphoryl transfer reaction. This project has been undertaken in collaboration with the Schug lab (FZ Jülich), who focused on the conformational transition, and we initially concentrated on the chemical step. The significance of such studies is given as HKs are some of the main prokaryotic signaling systems. The HK structure consists of two structurally conserved catalytic domains (the dimerization histidine phosphotransfer domain DHp, and the catalytic ATP-binding domain CA), which enable autokinase, phosphotransfer,
and phosphatase activities.

The phosphorylation of a His side-chain consists of a P–N bond being created, therefore, this process needs to be described accurately with the QM method used. Since this interaction pattern was not really considered in the original parametrization of DFTB3 (our method of choice in the simulations), this parameter set performs poorly. Therefore, we developed an SRP consisting of a modified P–N repulsive function for DFTB3, to improve its performance for reactions in which a P–N bond is replaced by a P–O bond or vice versa [2]. We showed that QM/MM metadynamics simulations using the new parameterization describe biochemical phosphorylation and hydrolysis reactions involving P–N bonds accurately. The parameter set was benchmarked on a model reaction mimicking the auto-phosphorylation of His, and also was applied to study the auto-phosphorylation of a genuine HK protein.

Then, we focused on a detailed mechanistic understanding of the functional cycle of WalK HK, i.e., the activation followed by the phosphoryl transfer [1]. We employed a multi-scale simulation approach, consisting of classical as well as hybrid QM/MM MD simulation. We found that the low rates of the process result from the high rate-limiting barrier of the conformational transition between inactive and active states, which has to be overcome by the energy transduction from the sensory domain in the naturally occurring process. Concentrating on the chemical step, we found that the auto-phosphorylation inside DHp proceeds via a penta-coordinated transition state to a protonated phosphohistidine intermediate, which is consequently deprotonated by a suitable nearby base. The reaction energetics are controlled by the final proton acceptor and the presence of an Mg²⁺ cation. The phosphorylation step exhibits a lower barrier and down-the-hill energetics, which represents the first portion of energy released during the long phosphorylation cascade in the naturally occurring process. Thus, our work offers a detailed mechanistic model for the auto-phosphorylation of HK.

Selected publications

[0] K. Spies; L. Pfeffer; M. Elstner; T. Kubař; N. Gillet: Multiscale simulations of proton-coupled electron transfer in biomimetic peptides, 2024, manuscript in preparation.
[1] F. Idiris; M. Kansari; H. Szurmant; T. Kubař; A. Schug: Mechanism of activation and autophosphorylation of a histidine kinase, Commun. Chem. 2024, 7(1), 196.
[2] M. Kansari; L. Eichinger; T. Kubař: Extended-sampling QM/MM simulation of biochemical reactions involving P–N bonds, Phys. Chem. Chem. Phys. 2023, 25, 9824–9836.
[3] D. Maag; J. Böser; H. A. Witek; B. Hourahine; M. Elstner; T. Kubař: Mechanism of proton-coupled electron transfer described with QM/MM implementation of coupled-perturbed density-functional tight-binding, J. Chem. Phys. 2023, 158(12), 124107.
[4] J. Böser; T. Kubař; M. Elstner; D. Maag: Reduction pathway of glutaredoxin 1 investigated with QM/MM molecular dynamics using a neural network correction, J. Chem. Phys. 2022, 157(15), 154104.
[5] C. L. Gómez-Flores; D. Maag; M. Kansari; V.-Q. Vuong; S. Irle; F. Gräter; T. Kubař; M. Elstner: Accurate free energies for complex condensed-phase reactions using an artificial neural network corrected DFTB/MM methodology, J. Chem. Theory Comput. 2022, 18(2), 1213–1226.
[6] D. Maag; M. Putzu; C. L. Gómez-Flores; F. Gräter; M. Elstner; T. Kubař: Electrostatic interactions contribute to the control of intramolecular thiol–disulfide isomerization in a protein, Phys. Chem. Chem. Phys. 2021, 23(46), 26366–26375.
[7] D. Maag; T. Mast; M. Elstner; Q. Cui; T. Kubař: O to bR transition in bacteriorhodopsin occurs through a proton hole mechanism, Proc. Natl. Acad. Sci. USA 2021, 118(39), e2024803118.
[8] B. Hourahine; B. Aradi; ... T. Kubař; ...T. Frauenheim: DFTB+, a software package for efficient approximate density functional theory based atomistic simulations, J. Chem. Phys. 2020, 152(12), 124101.
[9] N. Gillet; M. Elstner; T. Kubař: Coupled-perturbed DFTB-QM/MM metadynamics: Application to proton-coupled electron transfer, J. Chem. Phys. 2018, 149(7), 072328.
[10] M. Putzu; F. Gräter; M. Elstner; T. Kubař: On the mechanism of spontaneous thiol–disulfide exchange in proteins, Phys. Chem. Chem. Phys. 2018, 20(23), 16222–16230.
[11] T. Kubař; K. Welke; G. Groenhof: New QM/MM implementation of the DFTB3 method in the Gromacs package, J. Comput. Chem. 2015, 36(26), 1978–1989.

ORCID: 0000-0002-2419-6912