First Step of the Transglutaminase Reaction Catalyzed by Activated

Subscriber access provided by UNIV OF LOUISIANA

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

The First Step of the Transglutaminase Reaction Catalyzed by Activated Factor XIII Subunit A, Hybrid Quantum Chemistry-Molecular Mechanics Calculations Gábor Balogh, Laszlo Muszbek, and Istvan Komaromi J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b00542 • Publication Date (Web): 19 Apr 2019 Downloaded from http://pubs.acs.org on April 23, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The First Step of the Transglutaminase Reaction Catalyzed by Activated Factor XIII subunit A, Hybrid Quantum Chemistry-Molecular Mechanics Calculations

Gábor Balogh1, László Muszbek*1, and István Komáromi1 1Division

of Clinical Laboratory Sciences, Department of Laboratory Medicine, Faculty of Medicine, University of Debrecen, Debrecen, Hungary

*Corresponding author: E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 36

ABSTRACT The A subunit of blood coagulation factor XIII belongs to the family of transglutaminase enzymes. Its active form (FXIIIa) catalyzes isopeptide bond formation between Glu and Lys residues of specific substrates. Little data is available on the mechanism of this reaction. In this work, the first step of the proposed two-step process was investigated using two different protocols of hybrid QM/MM calculations: an ONIOM-based model as well as QM/MM/MD metadynamics simulations in explicit TIP3P solvent with Gromacs, Plumed and a DFTB3 package. Based on calculations involving a truncated system derived from docking of a peptide substrate, our study confirms the higher stability of a zwitterionic form of the catalytic Cys and His residues in the Michaelis complex as well as the “resting” state of the enzyme. Potential energy surfaces, obtained by geometry optimizations with Gaussian, show a two-step reaction mechanism with a zwitterionic tetrahedral intermediate formation in the first and NH3 dissociation in the second step in case of our ONIOM system. In contrast, in QM/MM MD metadynamics simulations all three steps occurred in a concerted manner. As a conclusion, our model is able to provide insights into the reaction mechanism of this enzyme.

2


Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction

Blood coagulation factor XIII (FXIII) belongs to the transglutaminase enzyme family. The human transglutaminase family comprises eight members showing enzymatic activity: FXIII-A and transglutaminases 1-7.1-4 The erythrocyte membrane protein band 4.2 is also considered a member of this family, but it is catalytically inactive.1-4 The transglutaminase enzymes can catalyze the formation of an isopeptide bond between a glutamine and a lysine residue of different proteins, a type of posttranslational modification, or the hydrolysis of the Gln amide group in the absence of an amine-donor substrate.1-4 Similarly to serine and cysteine proteases, a two-step mechanism has been proposed for the transglutaminase family. First, a thioester intermediate is formed (referred to as acylation), which can either react with an amine donor or the intermediate can be hydrolyzed if a water molecule enters the active site (i.e. deacylation).2-4 FXIIIa is a multifunctional transglutaminase exhibiting both extracellular and intracellular functions. It plays an important role in the late stage of blood coagulation by formation of crosslinks between fibrin α and γ chains, and covalently binding α2-antiplasmin to the fibrin clot.4 The crosslinking reaction mechanically strengthens the fibrin clot, protects it from the shear stress of circulation blood and increases its resistance against fibrinolysis.4, 5. Unlike other members of the family, FXIII is a zymogen homodimer of potentially active catalytic A subunits (FXIII-A) in its intracellular form. In the plasma it is of tetrameric structure (A2B2) two carrier B subunits (FXIII-B) are bound to an A2 homodimer.3, 4 FXIII-A consists of four domains: a beta-sandwich, the catalytic core domain and two beta-barrel domains. At the N-terminal end an activation peptide (AP-FXIIIA) of 37 amino acids is connected to the beta sandwich. During activation of FXIII-A2B2 AP-FXIIIA is cleaved off by thrombin at the Arg37-Gly38 peptide bond.3, 4 Then, in the presence of Ca2+ the FXIII-B subunits dissociate from the heterotetramer and FXIII-A assumes a catalytically active conformation (FXIIIa). In the absence of the FXIII-B subunits, Ca2+ binding is sufficient to bring FXIII-A2 into the active conformation without thrombin cleavage, although this is a much slower 3



Page 4 of 36

process. FXIIIa is a transglutaminase, in which three residues of the core domain, Cys314, His373 and Asp396 form the catalytic triad, while the Trp279 side chain is involved in the stabilization of the oxyanion intermediate of the reaction3, 4. Multiple X-ray diffraction structures of inactive FXIIIA2 have been published. In these structures, the access of substrates to the active center is blocked by the Tyr560 residue in the beta-barrel domain of the same subunit and by the activation peptide (AP) of the other A subunit, the latter interaction is stabilized by a salt bridge between the Arg11 residue of the AP3, 4 and the Asp343 residue near the catalytic center. An X-ray diffraction structure of non-proteolytically activated (non-cleaved) FXIII bound to the irreversible peptide-like inhibitor ZED-1301 was published by Stieler et al. in 2013, showing FXIII-A in a monomeric state.6 The monomeric nature of both the thrombin-cleaved and non-proteolytically activated FXIII has been confirmed by an analytical ultracentrifugation experiment by Anokhin et al.7 Mechanistic studies on both the acylation and the deacylation step catalyzed by transglutaminase 2, another enzyme in this family, have been reviewed by Keillor et al.8 The data available for factor XIII are limited. Kinetic measurements with various acyl donor substrates, in the presence or absence of an amine substrate have been reported,9-12 without investigating the reaction mechanism in depth. Furthermore, we are not aware of theoretical studies in which the catalytic mechanism itself was investigated in the transglutaminase family. A previous important in silico study on transglutaminase 2 deals with the binding of covalent inhibitors to the enzyme, but not the conversion of a model substrate.13 For cysteine proteases, a related class of enzymes, both static energy-minimization based QM/MM calculations and QM/MM dynamics simulations have been published.14-19 In the QM/MM dynamics based studies a semi-empirical method, such as SCCDFTB, was used for treatment of the QM part. Hybrid methods allow description of chemical reactions in complex systems. According to this scheme, the catalytically relevant residues are treated by a high-level quantum chemical (QM) approach, while the rest of the system is described by molecular mechanics (MM). The QM is required to describe the bond breaking and formation events, while lower resolution treatment of 4


Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the protein enables extensive sampling of the protein environment.20-23 Hybrid methods have been shown to be powerful to provide quantitative models for enzymatic catalysis22, 23 as well as some complex physical processes. The ONIOM scheme is an example of such a hybrid method.24, 25 For the QM subsystem, density functional theory (DFT) based methods are used frequently as they can give results comparable to post-Hartree-Fock methods, in many cases at a significantly lower computational cost. “Newer” DFT-based methods containing either terms for computing London dispersion energy (for example ωB97XD26) or methods with an excessive parametrization, such as Minnesota functionals, eg. M062X27, are able to take in account the dispersion interactions.19 The aim of present work was to investigate the acylation step of the reaction mechanism of FXIIIa, common in both transamidation and amide hydrolysis catalyzed by transglutaminases, on a model system, using both ONIOM-type calculations and QM/MM MD simulations. It was intended to obtain deeper general insight in the mechanism of the first step of the transglutaminase reaction and to obtain special information on the reaction catalyzed by FXIIIa. It is believed that the results provide information that could be utilized in constructing transglutaminase inhibitors. Our ONIOM-type QM/MM calculations we determined the geometries with two DFT methods (ωB97XD, M062X) and with the high-level MP2 method. We also optimized the geometries of energy minima and saddle points and computed potential energy surfaces. Furthermore, we performed QM/MM MD simulations to calculate free energies using the collective variable (CV)based well-tempered metadynamics method28, where the catalytically most important residues are treated with DFTB3, a relatively recently developed variant of the SCC-DFTB method.29-31

5



Scheme 1. Scheme of the proposed reaction mechanism of transglutaminases, based on the publication of Iismaa et al.32 Here the numbering of residues corresponds to FXIII-A. The modified scheme was created using MarvinSketch, (www.chemaxon.com)

2. Methods

2a. Model system preparation

Multiple conformations of the N-terminal dodecapeptide of α2-antiplasmin (the N-terminally cleaved form, see review33), a glutamine-containing substrate of Factor XIII were obtained from an equilibrium MD simulation. The starting conformation of the peptide was generated using UCSF Chimera34 by setting all backbone ϕ torsional angles to -139° and ψ torsional angles to 135°. The NPT simulation was performed at T = 310 K with GROMACS35 using the AMBER03 force field36 and TIP3P explicit water model.37 6


Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The structure of non-proteolytically activated factor XIII subunit A published by Stieler et al (PDB: 4KTY) was chosen for simulations.6 The conformations at t = 8 ns, 9 ns and 10 ns of the peptide were docked using HADDOCK38-40. Residues 214, 215, 223, 279-283, 289, 290, 313-315, 317, 339, 342, 350-352, 360, 365-374, 398-400, 439, 441, 456, 459 and 460 of factor XIII (subunit A) and all residues of the dodecapeptide were defined as “active residues”. “Unambigous” restraint was applied between the Cys314 residue of factor XIII and the amide carbon atom of the peptide Gln2 side-chain. The molecular docking resulted multiple conformations in several clusters. The conformation with the optimal relative position of residues (including the hydrogen bond of Trp279) was chosen. The core domain from the chosen structure, with the peptide bound, was refined using a simulated annealing protocol of 10 cycles (0-2 ns: heating from 5 K to 400 K, 2-5 ns: simulation at 400 K, 5-10 ns: cooling to 5 K in each cycle, under NPT conditions), using the GROMACS software package35 and the AMBER 14SB force field. The peptide in the refined structure is shown in Figure 1. In this simulation, the v-rescale thermostat41 and the Parinello-Rahman barostat42, 43 was used. The simulation was performed in TIP3P explicit solvent.37 The Coulomb and van der Waals cutoffs were set to 1.0 nm and the PME method was used for calculating long-range electrostatics calculations44.

2b. Model systems used in the QM/MM calculations and QM/MM/MD simulations

We investigated the acylation step of the reaction mechanism of FXIIIa, common in both transamidation and amide hydrolysis catalyzed by transglutaminases, on a model system, using both ONIOM-type calculations and QM/MM MD simulations. The ONIOM scheme is an example of a hybrid method. In its simplest form, the system is separated into a “high-level” and a “low-level” part. The total energy of the system is calculated as the energy of the full system using the lower level method plus the energy of the smaller, model system, minus the energy of the model system 7



Page 8 of 36

obtained by the lower lever method (a “subtractive” scheme.) Three or more subsystems, calculated using different quantum mechanics or molecular mechanics methods, are also possible.24, 25 In this paper, the term “ONIOM-type calculation” refers to a scheme with two subsystems, QM and MM. Our model system, used in the ONIOM-type QM/MM calculations, was obtained by truncation from the “refined” structure (See figure 2A). The high-level part consisted of catalytic cysteine (Cys314), histidine (His373) as well a propionamide substrate (the truncated glutamine residue). The model system included further amino acids (residues 211-216, 219-227, 266-268, 274-279, 283-292, 311-313, 315-322, 334-342, 368-372, 374-376, 395-400, 430-434, 464-466), corresponding roughly to a 12 Å sphere around the active center cysteine, the energy calculation of these residues was done using the AMBER molecular mechanics force field, as implemented in the Gaussian software. An N-acetyl group was connected to all “new” N-terminals and methylamide groups to the “new” C-terminals where peptide bonds were cut. To investigate the protonation state of the catalytic Cys and His residues also in the resting state of the enzyme, a second model system was created which included the same residues of the protein but not the substrate. A negatively charged cysteine and a protonated histidine were used in the two structures. The Gaussian input files were prepared using the TAO package.45 To obtain the starting structure used our QM/MM MD simulations, the geometry of the Michaelis complex optimized at the M062X/6-31+G(d,p) level using our ONIOM-based protocol was chosen as starting structure. The system was solvated using TIP3P water molecules37 and it was neutralized by adding Na+ and Cl- ions so that the salt concentration was close to 0.15 mol/dm³. The QM system included the same residues as in ONIOM-type calculations, the semi-empirical DFTB3 method was applied to the QM atoms29-31. The MM part was calculated using the AMBER 14SB force field 46. The 3ob-3-1 parameter set was used in our DFTB3 calculations.47, 48 The link atoms were located on the line connecting the Cα and the Cβ atoms approximately 1 Å from the beta carbon atoms of amino acids Cys314 and His373, defined as virtual sites. In order 8


Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to avoid the over-polarization of the link atoms, the charge of the two carbons were set to zero. Therefore the total charge of the system differed slightly from zero.

2c. ONIOM-type hybrid QM-MM calculations

The ONIOM-type hybrid QM/MM calculations were performed using the Gaussian 03 and 09 software packages49-51. An unconstrained geometry optimization was performed to obtain a starting geometry of the Michaelis complex for further relaxed scans. To confirm the greater stability of the histidine-protonated state in both the Michaelis complex and the enzyme without the substrate, two relaxed scans were performed using the MP2/6-31G(d) method. The reaction coordinate was the distance of the proton from the sulfur atom of cysteine. Energy minimization was performed in grid points (0.2 Å, 0.1 Å near the transition state) to obtain potential energy surfaces as a function of two different pairs of reaction coordinates (Cys314 sulfur-substrate amide C (d1) and substrate amide C-amide N (d2) distances vs S-C (d1) and substrate amide N-H distances (d3)) (see figure 2B). (Energies of some points were not calculated in the d1: 2.8-3.2 Å, d2: 2.0-2.6 Å region due to convergence problems. These are expected to be high energy conformations.) Two different methods were applied for the high level, QM system for both scans (MP2/6-31G(d) vs M062X/6-31+G(d,p)). To study the breakage of the C-N bond starting from the tetrahedral intermediate, two further relaxed scans were performed using the MP2/631G(d) and M062X/6-31+G(d,p) methods. The Michaelis complexes, the transition states and the intermediate geometries were optimized using the MP2 method and different DFT methods (M062X, ωB97XD) and basis sets. Zero-point vibration corrected energies were also computed for the energies calculated with the two DFT methods. (This correction was not calculated at the MP2 level of theory due to the huge computational costs of MP2 vibrational frequency calculations on such systems.)

9



Page 10 of 36

Visualization of molecular structures was performed using UCSF Chimera.34 The twodimensional plots of potential energy and free energy surfaces were created using DPlot.52

2d. QM/MM MD simulations

The QM/MM MD simulations were carried out using the GROMACS package with the DFTB3 implementation 29-31, 35, 53 by Kubař et al. and PLUMED for the metadynamics 28, 54, 55 simulations. After steepest descent minimization of our solvated system (as described in the “Model systems” section), a 1 ns “conventional” molecular dynamics simulation under NVT conditions was performed, using a v-rescale thermostat41 and a time step of 0.5 fs. All residues except 312-319, 370, 372-375, 396, 398 and 400 were position restrained. A weaker restraint was applied to Trp279 side-chain and the propionamide substrate during the equilibration to prevent the breakage of the hydrogen bond between them and fix them in the correct, hydrogen-bonded position. The free energy surface as function of the CVs was obtained by a 4.4 ns long well-tempered metadynamics simulation (WT-MetaD) with bias factor set to 4528. During the simulation, adaptively sized (ADAPTIVE=DIFF) Gaussians were added every 500 MD steps56. The two collective variables were the sulfur-amide carbonyl carbon distance (d1) and the amide carbonyl carbon-amide nitrogen distance (d2) (same as the first pair of collective variables used in ONIOMtype calculations). The same position restraints as in the NVT equilibration were applied with the exception of the propionamide substrate. The conformational space accessible to the system was limited using lower and upper walls (CV1: 0.17-0.37 nm, CV2: 0.13-0.33 nm). A further restraint was applied on the histidine proton-amide nitrogen distance to prevent “flipping” of the His373 residue. To determine the protonation state of the active center Cys314 and His373 amino acids, another WT-MetaD simulation with biasfactor = 30 was performed using the same reaction 10


Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


coordinate as in the MP2 calculation as a CV (lower wall: 0.13 nm, upper wall: 0.25 nm). The MD simulation parameters were the same as above, with additional restraints to prevent His373 flipping and displacement of substrate from the binding tunnel.

3. Results and discussion

3a. Binding of substrate peptide to factor XIII

Figure 1: The α2-antiplasmin dodecapeptide bound to the core domain of factor XIII, refined structure. Labels of amino acids in the catalytic domain of FXIII-A are shown in black, key residues of α2-antiplasmin peptide involved in the interactions with FXIII-A* are depicted in purple.

Starting from an “initial” structure obtained by protein-protein docking, we performed a simulated annealing molecular dynamics simulation in order to obtain a refined FXIII-A-peptide complex structure (figure 1). Further details of these simulations are found in the Methods section. In the “final” conformation (t = 100 ns) extracted from the trajectory, the positively charged N-terminal of 11



Page 12 of 36

Asn1 forms an intramolecular salt bridge with Glu3. There is a hydrogen bond between the hydroxyl group of Tyr372 and the amide carbonyl oxygen atom of Asn1 side-chain. Pro7 occupies a small pocket between Val369, Phe339, Leu439 and Tyr441 while Leu8 and Leu10 side-chains form hydrophobic contacts with residues near the beta-sheet homologous to TG2 hinge regions. Lys12 is bound to Asp456 with a salt bridge. In a previous work from our group, the binding of truncated and mutated peptides to FXIIIA was investigated by a kinetic assay and NMR measurements12. Cleary et al. and Marinescu et al. have also published similar studies11, 57. These experiments have demonstrated the importance of Asn1 and residues 7-12 in binding of the dodecapeptide to activated FXIII-A and suggested a secondary binding site on FXIII-A. Our theoretical model of the FXIII-A core domain-α2-antiplasmin dodecapeptide complex is in qualitative agreement with the experimental NMR data regarding the importance of Asn1 and the C-terminal part of the peptide in the interaction. As our simulated annealing-based method yields a structure corresponding to the most probable conformation, while the NMR measurements provide information on an ensemble of structures, complete agreement between the results obtained by the two methods cannot be expected.

3b. ONIOM-type QM/MM calculations – protonation state of Cys314 and His373 in the Michaelis complex and “resting” state

12


Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: A. The model system used in the ONIOM-type calculations. B. The “characteristic” distances (d1-d4). d1:Cys314 sulfur-substrate amide C distance, d2: substrate amide C-amide N distance, d3: substrate amide N-proton distance, d4: His373 imidazole epsilon N-proton distance.

Figure 3: Energy values obtained in the relaxed scans of the Cys-His proton transfer (Michaelis complex, structure with ligand removed), as a function of the reaction coordinate (cysteine S-proton distance)

Proton transfer plays an important role in the catalytic mechanism of transglutaminases, as well as serine and cysteine proteases: the catalytic Cys residue must be deprotonated (the His residue serves as a general base to deprotonate the cysteine, making it a stronger nucleophile) for the nucleophilic attack, while protonation of the amide NH2 requires proton transfer from the protonated His373 13



Page 14 of 36

residue to the amide nitrogen. The protonation states of these residues in the resting state of cysteine proteases have been the topic of several studies using various methods, confirming the existence of a cysteine thiolate-protonated histidine ion pair58-60. The experimental data available for transglutaminases are even more limited, but Case and Stein have proposed a neutral Cys/His diad for guinea pig and human transglutaminase 2 in the resting state of the enzyme based on kinetic isotope effects measurements. They suggested a nucleophilic attack by cysteine thiolate in the ratelimiting step. We carried out ONIOM-type hybrid QM/MM calculations to determine the protonation states of the two residues in both the Michaelis complex and the “resting” state of the active enzyme. Two relaxed potential energy scans were performed on the two truncated model systems of the enzyme as described in the Methods section, one lacking the substrate. The potential energies obtained by energy minimization as the function of the reaction coordinate (the distance of the proton from the sulfur atom of Cys314), are shown in Figure 3. The energy minimum is located at 2.4 Å and 2.3 Å in the two cases, which corresponds to a zwitterionic form with protonated histidine and deprotonated cysteine. The cysteine-protonated form (~1.3-1.4 Å) has significantly higher energy (36.1 kJ/mol and 54.7 kJ/mol respectively) than the histidine-protonated form. In both cases, the histidine-protonated ion-pair form is predicted to be more stable than the cysteine protonated form. We also calculated the zero-point vibration corrected energies for both end states of the proton transfer reaction, using multiple basis sets and two DFT methods (ωB97XD, M062X). Results can be found in Table S1. The energy of the optimized His-protonated form is always lower than the Cys-protonated state, further supporting its higher stability. Assuming that a Cys-protonated state of the protein exists, the proton transfer would occur before the substrate binding, due to the large energy differences between the zwitterionic and the hypothetical Cys-protonated forms. Therefore the proton transfer and the nucleophilic attack could be considered as two fully separate reaction steps. This is in a good agreement with the mechanism proposed by Case and Stein, e.g. the nucleophilic attack of a cysteine thiolate. On the other hand, 14


Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


our calculations demonstrated that the ion pair form is significantly more stable also in the resting state of active FXIII, as it was verified previously for cysteine proteases.

3c. ONIOM-type QM/MM calculations – Mechanism of the first (acylation) step

Figure 4: A. The two potential energy surfaces, obtained using the d1 (C-S distance) and d2 (Camide N distance) reaction coordinates. B. The two potential energy surfaces, obtained using the d1 (C-S distance) and d3 (amide N-Histidine H distance) reaction coordinates. 15



Page 16 of 36

Regarding the acylation step of the cysteine protease reaction, different reaction mechanisms have been proposed, the various reaction paths were summarized in the recent paper of Fekete et al.19. Lowe suggested a scheme similar to serine proteases61, while theoretical studies have resulted in different reaction mechanisms. In these studies, it was proposed that either the histidine imidazole nitrogen – amide nitrogen proton transfer and S nucleophilic attack may occur in the same elementary step, or the proton transfer can precede the attack. In the same studies, the existence of an anionic tetrahedral intermediate (hemi-thioacetal anion) during one or both steps of serine or cysteine protease-catalyzed reactions was also investigated. Some of them have demonstrated the existence of such an intermediate, while others did not find a local minimum corresponding to it. According to Iismaa et al.,32 the existence of such an intermediate may be one of the possible explanations for the decreases in reaction rates when Trp241 in transglutaminase 2 (corresponding to Trp279 in FXIII-A) is substituted with other amino acids. In order to study the acylation step of the transglutaminase reaction using ONIOM-type QM/MM calculations, we performed two potential energy scans, applying two different pairs of reaction coordinates (d1 and d2 vs. d1 and d3). Using the first pair of reaction coordinates (d1 for S-C and d2 for C-N distances, see Methods section), we obtained two similar potential energy surfaces when applying two different methods to calculate the energies for the QM system (MP2/6-31G(d) and M062X/6-31+G(d,p)). Results can be seen in figure 4A. According to the results, the Michaelis complex (located at ~3.4, 1.35 Å), has been predicted to have the lowest energy using both methods. The plotted PES suggests a possible saddle point near d1 = 2.0 Å, d2 = 1.5 Å. A local minimum can be observed at d1 = 1.8 Å, d2 = 1.7 Å. This corresponds to a tetrahedral, zwitterionic intermediate (TI) where the amide NH2 group of the peptide is protonated. The potential energy surface in the region d1 = 1.8 Å, d2 = 2.4-2.6 Å is very flat and the energy barrier for breaking the C-N bond in the TI is predicted to be low. This region may contain a second saddle point. The main difference between the results obtained by the two methods is that M062X/6-31G(d) predicts a 16


Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


higher energy barrier than the MP2/6-31G(d) method and the potential energy of the intermediate is also higher. This corresponds to a mechanism, where the proton transfer from the histidine imidazole nitrogen to the amide nitrogen occurs in the same elementary step together with the nucleophilic attack of the thiolate on the amide carbonyl carbon, but the breakage of the C-N bond occurs later, in a following step. We conducted two similar potential energy scans along a second pair of reaction coordinates (d1 for S-C and d3 for N-H distance), see figure 4B. The results indicate that, consistently with the previous results, a single-step reaction mechanism can be proposed for the attack and the NH2 proton transfer, in which a similar tetrahedral intermediate is formed through a single transition state. The shape of the two potential energy surfaces are remarkably similar, but the M062X/631+G(d,p) energies are consistently higher, than the MP2/6-31G(d)-obtained values. The tendency is similar to the results got with the d1-d2 scans. No C-N bond breakage was observed during the geometry optimizations when applying the second pair of reaction coordinates.

Figure 5: Energy values obtained in the relaxed scans of the C-N bond breakage using two different QM methods, starting from the optimized structure of the tetrahedral intermediate, as a function of the reaction coordinate (1: TI, 2: TS 2, 3: thioester with NH3)

17



Page 18 of 36

The breakage of the C-N bond was further investigated by a relaxed scan using only one reaction coordinate, the carbonyl carbon-amide nitrogen distance (d3). Figure 5 shows the computed energies as a function of the distance. A local maximum can be observed at 2.2-2.3 nm in both cases. The barrier is large enough to prevent the spontaneous dissociation of ammonia during an unconstrained energy minimization of the TI, but it is likely a rapid reaction step. At the end point of the scan (d3 = 2.9-3.0 nm), a thioester intermediate is formed, but the ammonia molecule is still located in the active site of the enzyme, close to the active center. The energy of this “end” state is higher than that of the tetrahedral intermediate, this may be a consequence of the fact that the NH3 is still bound. Our method is capable for the prediction of the intermediates and the calculation of energy barriers, but it is probably not suitable for determining the energy of the ammonia release or the energy of the thioester form compared to the Michaelis complex. The correct description of the reaction path of NH3 release may require more complex reaction coordinates and taking the solvation of NH3 into account would be computationally very demanding for an ONIOM system with this size. On the basis of our calculations we can conclude that the nucleophilic attack leading to the tetrahedral intermediate and the proton transfer are expected to occur in the same elementary step, while the release of the ammonia takes place separately. However, the latter reaction step is expected to occur very fast due to the low energy barrier.

3d. ONIOM-type QM/MM calculations – Energies and geometries of the stationary points

18


Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: A. Geometries of the Michaelis complex (MC), transition state 1 (TS 1) and the tetrahedral intermediate, optimized at M062X/6-311+G(d,p) level for the QM subsystem. B. Geometries of the transition state of the ammonia release step (TS2) and the thioester intermediate with the NH3 molecule still bound, optimized at the M062X/6-311+G(d,p) level for the QM subsystem.

19



Method/Basis set MP2/6-31G(d) MP2/6-31G(d,p) MP2/6-31+G(d,p) MP2/6-311+G(d,p) ωB97XD/6-31G(d) ωB97XD/6-31G(d,p) ωB97XD/6-31+G(d,p) ωB97XD/6-311+G(d,p) M062X/6-31G(d) M062X/6-31G(d,p) M062X/6-31+G(d,p) M062X/6-311+G(d,p)

E(Tetrahedral i.) – E(TS) – E(Michaelis c) E(Michaelis c.) (kJ/mol) (kJ/mol) 83.88 62.42 79.16 65.18 89.09 72.16 85.35 69.12 107.15 92.37 105.81 96.31 114.53 100.87 118.18 103.80 93.58 82.44 93.25 86.52 102.38 92.12 105.67 94.57

Page 20 of 36

E(TS) – E(Michaelis c), ZP vibration corrected (kJ/mol) 93.39 91.54 101.47 105.12 80.84 80.49 89.55 92.76

E(Tetrahedral i.) – E(Michaelis c.), ZP vibration corrected (kJ/mol) 91.26 93.94 100.22 103.51 80.81 83.54 90.34 92.74

Table 1: Energies of the optimized transition state and the tetrahedral intermediates calculated using multiple DFT methods and basis sets, compared to the energy of the Michaelis complex calculated using the same method. The energies corrected for zero-point vibration are also shown. (Michaelis c.: Michaelis complex, Tetrahedral i.: tetrahedral intermediate)


E(TS 2) – E(Tetrahedral i.) (kJ/mol) 22.88 20.72 27.87 25.77 22.51 22.15 26.87 24.33 22.96 22.89 26.95 24.07

E(Thioester i.) – E(Tetrahedral i.) (kJ/mol) 12.82 8.50 17.34 13.20 16.71 14.68 18.27 13.35 16.17 14.49 17.73 12.75

E(TS 2) – E(Tetrahedral i.), ZP vibration corr. (kJ/mol) 16.19 17.20 20.36 17.67 17.82 18.39 21.61 18.00

E(Thioester i.) – E(Tetrahedral i.), ZP vibration corr. (kJ/mol) 8.36 6.45 8.29 3.00 9.29 8.47 10.58 5.43

Table 2: Energies of the optimized second transition state and the thioester intermediate (with NH3 still bound) calculated using multiple DFT methods and basis sets, compared to the energy of the tetrahedral intermediate calculated using the same method. The energies corrected for zero-point vibration are also shown. 20


Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The zero-point vibration corrected energies of the five stationary points we calculated (Michaelis complex, first transition state (TS), tetrahedral, zwitterionic intermediate, the second transition state (TS2) and the thioester form with the ammonia bound) using different DFT methods and basis sets are shown in Tables 1 and 2, along with the uncorrected energies derived from unconstrained optimizations. (The geometries of these states are presented in Figure 6A and 6B.) Regarding the uncorrected energies of the stationary points, the activation energies for the nucleophilic attack step calculated by the MP2 method are 23-33 kJ/mol lower than the energy calculated by the ωB97XD method, and 10-20 kJ/mol lower than that calculated by the M062X method, depending on the basis set. The energies of the tetrahedral intermediate calculated by the above two methods are 29-35 and 20-25 kJ/mol higher than the values obtained by MP2 method. The zero-point vibration corrected energies of the TS are 80.8-91.5 kJ/mol higher than that of the Michaelis complex. The energy difference is 80.8-93.9 kJ/mol in the case of the TI. The energies of the transition state and the intermediate are very close in all cases. We conclude that the intermediate either has very short half-life or it does not exist at all at room temperature. As for the dissociation of the tetrahedral intermediate, all three methods predict a 21-27 kJ/mol higher energy for the second transition state compared to the TI. The zero-point corrected energies for the transition state of the ammonia release step (TS2) are 4-6 kJ/mol below the uncorrected energies. This result is consistent with the low energy barrier of this reaction step as it was discussed before, in section “Mechanism of the first (acylation) step”. According to the transition state theory, the activation free energy values can be compared to kinetic constants using the Eyring-Polanyi equation. However, to our knowledge all available kinetic data for FXIII-A was obtained by enzyme kinetics studies of cross-linking or amide hydrolysis. In this case, the reaction rates depend on both steps of the reaction, not just the acylation 21



Page 22 of 36

step. It is not known whether the first or the second step is rate-limiting in case of factor XIII. Assuming the slowest, Cys “attachment” step of our proposed reaction mechanism as rate-limiting, the energy barrier was around 85 kJ/mol in our MP2-based calculations, the highest level method in our study. Unfortunately we were not able to calculate the zero point vibration correction at this level, but it may be in the range of -15 – -10 kJ/mol based on the DFT data. Assuming a corrected energy barrier of 70 – 75 kJ/mol, this would correspond to reaction rates 0.5 – 4 s-1 for the model compound. In contrast, the kcat values determined by kinetic measurements of FXIII-A and various substrates are in the range of 16 – 490 s-1. 11, 12 We can conclude that substrates with more than one residue are required for catalysis with higher efficiency.

Michaelis complex


d1 (Å) 3.436 3.406 3.630 3.435 3.525 3.507 3.511 3.515 3.370 3.364 3.368 3.370

d2 (Å) 1.352 1.351 1.347 1.350 1.349 1.349 1.346 1.344 1.350 1.349 1.347 1.346

d3 (Å) 2.874 2.856 2.962 2.857 2.824 2.827 2.840 2.852 2.684 2.688 2.715 2.730

Transition state 1 d4 (Å) 1.040 1.035 1.033 1.036 1.032 1.032 1.028 1.027 1.034 1.033 1.029 1.028

d1 (Å) 2.179 2.161 2.100 2.129 2.330 2.261 2.144 2.180 2.284 2.242 2.162 2.202

d2 (Å) 1.516 1.519 1.518 1.517 1.508 1.514 1.516 1.513 1.515 1.521 1.519 1.516

d3 (Å) 1.303 1.287 1.286 1.290 1.229 1.226 1.250 1.246 1.228 1.209 1.238 1.232

Tetrahedral intermediate d4 (Å) 1.293 1.292 1.294 1.292 1.363 1.363 1.335 1.340 1.373 1.395 1.360 1.366

d1 (Å) 2.016 2.017 1.985 1.997 2.040 2.044 2.004 2.012 2.033 2.041 2.002 2.009

d2 (Å) 1.595 1.594 1.583 1.585 1.592 1.587 1.580 1.583 1.595 1.588 1.583 1.586

d3 (Å) 1.060 1.061 1.057 1.061 1.060 1.065 1.060 1.058 1.063 1.072 1.065 1.060

Table 3: d1-d4 distances from optimized structures of the Michaelis complex, Transition state 1, Tetrahedral intermediate. (d1-d4: See figure 2B)

22


d4 (Å) 1.800 1.750 1.781 1.760 1.773 1.741 1.776 1.780 1.765 1.717 1.756 1.772

Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Transition state 2


d1 (Å) 1.854 1.859 1.845 1.854 1.852 1.855 1.848 1.857 1.862 1.867 1.856 1.866

d2 (Å) 2.170 2.139 2.180 2.154 2.233 2.170 2.215 2.179 2.205 2.174 2.179 2.144

d3 (Å) 1.030 1.027 1.026 1.029 1.028 1.028 1.026 1.026 1.029 1.029 1.028 1.027

Thioester + NH3 d4 (Å) 1.990 1.960 1.992 1.971 1.990 1.982 2.015 2.000 1.971 1.952 1.987 1.975

d1(Å) 1.801 1.800 1.796 1.799 1.810 1.810 1.804 1.807 1.810 1.809 1.802 1.805

d2 (Å) 2.972 2.979 2.991 2.987 2.824 2.850 2.878 2.900 2.912 2.984 3.005 3.052

d3 (Å) 1.023 1.019 1.019 1.020 1.024 1.023 1.022 1.020 1.023 1.021 1.020 1.019

d4 (Å) 2.270 2.266 2.266 2.269 2.177 2.195 2.226 2.241 2.221 2.253 2.287 2.317

Table 4: d1-d4 distances from optimized structures of the second transition state, and the thioester intermediate (with ammonia still bound to FXIII)

From our optimized geometries of the stationary points discussed previously, we analyzed the “characteristic” distances d1-d4 (Figure 2B). The calculated distances are shown in Tables 3 and 4. Regarding the Michaelis complex, a similar structure is predicted, independently from the basis set or DFT method. The transition state and the tetrahedral intermediate have a similar geometry, the difference in the d1 and d2 distances is less than 0.1 Å in all cases. In the Michaelis complex, two important hydrogen bond interactions are present between the protein and the substrate. The Trp279 residue (corresponding to Trp241 in TG2), involved in stabilization of the transition state of the reaction, forms a hydrogen bond with the carbonyl oxygen of the amide group in the substrate. The amide NH2 group of the molecule forms a hydrogen bond with the backbone carbonyl oxygen of Tyr342. This interaction is preserved in the tetrahedral intermediate and in both transition states. The unreleased NH3 molecule also interacts with the same oxygen atom. The (histidine) proton-amide nitrogen distance was significantly longer in the transition state (1.21-1.25 Å) than in the tetrahedral intermediate (1.06 – 1.07 Å). As for the tetrahedral intermediate, d1 is still longer than a typical C-S bond (2.00-2.04 Å versus ~1.82 Å). A slightly

23



Page 24 of 36

shorter C-S bond length is predicted when basis sets with diffuse functions were used. In the second transition state, the carbon-nitrogen distance is between 2.14-2.23 Å depending on the method and basis set used, the carbon-sulfur distance is only slightly longer than its length in the thioester intermediate. The d1 distance in the Michaelis complex and the d2 and d4 distances have a slightly higher variability but this is expected as the substrate and the NH3 molecule are noncovalently bound to the enzyme in the two cases. In general, these differences in distances (relatively small in most cases) and their dependence on the basis sets are difficult to explain due to the extensive parametrizations of the used DFT methods.

3e. QM/MM MD simulation of the histidine-cysteine proton transfer

Figure 7: Free energy as a function of the CV (cysteine S-proton distance), calculated from the output of the metadynamics simulation of this reaction step.

We conducted hybrid QM/MM MD simulations, besides our ONIOM-type QM/MM calculations. We used the DFTB3 semi-empirical method and the 3ob-3-1 parameter set for the QM system to study both steps of the reaction investigated in our ONIOM-type QM/MM MD calculations (proton transfer in the Michaelis complex as well as the nucleophilic attack-thioester formation step). 24


Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In our His-Cys proton transfer metadynamics simulation we could observe multiple proton transfers from His373 to the sulfur atom of Cys314, and back to the histidine nitrogen with 2 ns of simulation time. (This does not reflect the kinetics of the actual chemical reaction as this is a “biased” MD simulation for free energy calculation.) The “reverse” protonation reaction is faster, the system remains in the Cys-protonated state for a relatively short period of time in all cases according to our non-equilibrium MD simulation. (The collective variable, as a function of time, is shown on Figure S1.) In a good agreement with our ONIOM-type calculations the zwitterionic, Hisprotonated form is predicted to have lower free energy. The obtained energy difference between the two protonation states of the active center residues is in a good agreement with the MP2 calculations (~38 kJ/mol vs. 36.1 kJ/mol) (see Figure 7). The positions of unrestrained amino acid sidechains remain relatively stable during the simulation.

3f. QM/MM MD simulation of the nucleophilic attack step and the ammonia release

Figure 8: Free energy surface calculated from the output of the two-CV metadynamics simulation (CV1: C-S distance (d1); CV2: C-amide N distance (d2)) 25



Page 26 of 36

For the nucleophilic attack-ammonia release simulation, a QM/MM MD free energy surface as the function of (d1, d2) distances was obtained by summing the Gaussians from PLUMED output (see Figure 8). According to the FES, the reaction occurs in a single step, and the TS is located close to d1 = 0.22 nm, d2 = 0.18 nm. In contrast to the results obtained in our ONIOM-type calculations, there is no local minimum corresponding to an intermediate comparable to the tetrahedral, zwitterionic form observed using DFT or MP2 methods. On the basis of this simulation, the tetrahedral intermediate is unstable, the ammonia molecule is almost immediately released after the nucleophilic attack. The nucleophilic attack, the ammonia release and the proton transfer from histidine to amide/ammonia nitrogen all occur in a highly concerted manner. The energy barrier of the reaction is predicted to be ~100 (97-98) kJ/mol, which falls in the range of the energies calculated using the ONIOM method with various functionals and basis sets. The thioester intermediate (before the ammonia release) has ~40 kJ/mol higher energy than the Michaelis complex. During the simulation, the first nucleophilic attack and thioester formation occurs at about ~1.54 ns. The reverse reaction can be observed at ~3.05 ns and a second nucleophilic attack at ~4.10 ns. The system is mostly in a Michaelis complex-like state between approximately 0 and 1.54 ns and 3.05 and 4.10 ns, while a thioester intermediate-like form is present between these intervals and after 4.10 ns. The values of the two collective variables, as well as the amide nitrogen-proton distance, are shown in Figure S2. Figure S3 shows the “characteristic” conformations of the Michaelis complex and the thioester intermediate. The positions of unrestrained amino acid residues in the active center of the enzyme remain relatively stable during the simulation until around ~4.3 ns. A hydrogen bond, between the backbone N-H of the Cys314 amino acid and the amide carbonyl oxygen is relatively weak as judged by the O-H distance, but it was present during most of the simulation time. 26


Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


No significant displacement of the unrestrained residues can be observed in the first ~4.2 ns of the simulation. The Trp279-peptide substrate hydrogen bond breaks at ~4.3 ns due to the movement of the substrate in the binding pocket. The part of simulation after 4.3 ns was not used in the free energy calculations.

4. Conclusion

In this study we described the first step of the FXIIIa transglutaminase catalytic mechanism based on multiscale QM/MM in silico calculations. New aspects of the first step in the FXIIIa-catalyzed transglutaminase reaction were revealed. The findings might provide information important for future studies aiming to create transglutaminase inhibitors. We have found a significantly more stable Cys314 thiolate-protonated His373 ion pair in both ONIOM-type and QM/MM MD simulations with explicit solvent. This suggests a low occurrence of a form containing a neutral dyad. A mechanism where the deprotonation of the cysteine residue would occur in the same step as the nucleophilic attack is therefore unlikely. The energy profile of the nucleophilic attack step as obtained by our ONIOM-type calculations suggests that the Cys nucleophilic attack and the protonation of the amide NH2 in the substrate (formation of tetrahedral intermediate) probably occur in the same elementary step and NH3 dissociation in a second step. In contrast in QM/MM MD simulations, all three reaction steps occurred in a concerted manner. The (in part) different mechanism obtained by the two different QM/MM protocols may be explained by the higher level QM methods used in our ONIOM-type calculations compared to the semi-empirical DFTB3 in QM/MM MD and the presence of an explicit solvent on QM/MM MD simulations. Iismaa and coworkers carried out interesting complementary studies with transglutaminase 2. The mechanism 27



Page 28 of 36

we found for FXIII-A* differs from the scheme proposed for transglutaminase 2: in our study the intermediate is zwitterionic rather than oxyanionic with a neutral NH2 group (in the ONIOM-type calculations) or no such intermediate exists (in QM/MM MD simulations). In conclusion, our study provides a model of the acyl-enzyme formation during the transglutaminase reaction of FXIII-A, supported by hybrid QM/MM calculations.

Conflict of Interest The authors declare no conflict of interest with the content of this article.

Supporting Information

Table S1: Energies of the Cys-protonated (neutral Cys and His) states of the Michaelis complex and “resting state” compared to the zwitterionic (His-protonated) forms. Figure S1: The collective variable as a function of time, from the QM/MM MD metadynamics simulation of the Cys-His proton transfer. Figure S2: (A) The two collective variables and the N-H distances (d3) as a function of time, during the QM/MM MD metadynamics simulation of the nucleophilic attack and NH3 release steps. Figure S3: “Characteristic” conformations for the Michaelis complex and the thioester intermediate from the QM/MM MD simulation of the nucleophilic attack/ammonia release step.

Acknowledgements

28


Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The work was supported by the Hungarian Scientific Research Fund (OTKA-106294, 113097) and by the GINOP-2.3.2-15-2016-00050 project. The latter project is co-financed by the European Union and the European Regional Development Fund. Gábor Balogh is recipient of a scholarship from “Talentum Foundation” of Gedeon Richter PLC. We acknowledge [PRACE/NIIF] for awarding us access to supercomputer resource based in Debrecen, Hungary.

References

1.

Facchiano, A.; Facchiano, F., Transglutaminases and their substrates in biology and human

diseases: 50 years of growing. Amino Acids 2009, 36 (4), 599-614. 2.

Iismaa, S. E.; Mearns, B. M.; Lorand, L.; Graham, R. M., Transglutaminases and disease:

lessons from genetically engineered mouse models and inherited disorders. Physiol Rev 2009, 89 (3), 991-1023. 3.

Komáromi, I.; Bagoly, Z.; Muszbek, L., Factor XIII: novel structural and functional aspects.

J Thromb Haemost 2011, 9 (1), 9-20. 4.

Muszbek, L.; Bereczky, Z.; Bagoly, Z.; Komáromi, I.; Katona, É., Factor XIII: a

coagulation factor with multiple plasmatic and cellular functions. Physiol Rev 2011, 91 (3), 931972. 5.

Schroeder, V.; Kohler, H. P., Factor XIII: Structure and function. Semin Thromb Hemost

2016, 42 (4), 422-428. 6.

Stieler, M.; Weber, J.; Hils, M.; Kolb, P.; Heine, A.; Büchold, C.; Pasternack, R.; Klebe,

G., Structure of active coagulation factor XIII triggered by calcium binding: basis for the design of next-generation anticoagulants. Angew Chem Int Ed Engl 2013, 52 (45), 11930-11934. 7.

Anokhin, B. A.; Stribinskis, V.; Dean, W. L.; Maurer, M. C., Activation of factor XIII is

accompanied by a change in oligomerization state. FEBS J 2017, 284 (22), 3849-3861.

29



8.

Page 30 of 36

Keillor, J. W.; Clouthier, C. M.; Apperley, K. Y.; Akbar, A.; Mulani, A., Acyl transfer

mechanisms of tissue transglutaminase. Bioorg Chem 2014, 57, 186-197. 9.

Muszbek, L.; Polgár, J.; Fésüs, L., Kinetic determination of blood coagulation Factor XIII

in plasma. Clin Chem 1985, 31 (1), 35-40. 10.

Kárpáti, L.; Penke, B.; Katona, E.; Balogh, I.; Vámosi, G.; Muszbek, L., A modified,

optimized kinetic photometric assay for the determination of blood coagulation factor XIII activity in plasma. Clin Chem 2000, 46 (12), 1946-1955. 11.

Cleary, D. B.; Maurer, M. C., Characterizing the specificity of activated Factor XIII for

glutamine-containing substrate peptides. Biochim Biophys Acta 2006, 1764 (7), 1207-1217. 12.

Pénzes, K.; Kövér, K. E.; Fazakas, F.; Haramura, G.; Muszbek, L., Molecular mechanism

of the interaction between activated factor XIII and its glutamine donor peptide substrate. J Thromb Haemost 2009, 7 (4), 627-633. 13.

Jasim, M. H.; Rathbone, D. L., Reaction profiling of a set of acrylamide-based human tissue

transglutaminase inhibitors. J Mol Graph Model 2018, 79, 157-165. 14.

Harrison, M. J.; Burton, N. A.; Hillier, I. H., Catalytic mechanism of the enzyme papain:

Predictions with a hybrid quantum mechanical/molecular mechanical potential. J Am Chem Soc 1997, 119 (50), 12285-12291. 15.

Han, W. G.; Tajkhorshid, E.; Suhai, S., QM/MM study of the active site of free papain and

of the NMA-papain complex. J Biomol Struct Dyn 1999, 16 (5), 1019-1032. 16.

Ma, S.; Devi-Kesavan, L. S.; Gao, J., Molecular dynamics simulations of the catalytic

pathway of a cysteine protease: a combined QM/MM study of human cathepsin K. J Am Chem Soc 2007, 129 (44), 13633-13645. 17.

Mladenovic, M.; Fink, R. F.; Thiel, W.; Schirmeister, T.; Engels, B., On the origin of the

stabilization of the zwitterionic resting state of cysteine proteases: a theoretical study. J Am Chem Soc 2008, 130 (27), 8696-8705. 18.

Shokhen, M.; Khazanov, N.; Albeck, A., Challenging a paradigm: theoretical calculations

30


Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the protonation state of the Cys25-His159 catalytic diad in free papain. Proteins 2009, 77 (4), 916-926. 19.

Fekete, A.; Komáromi, I., Modeling the archetype cysteine protease reaction using

dispersion corrected density functional methods in ONIOM-type hybrid QM/MM calculations; the proteolytic reaction of papain. Phys Chem Chem Phys 2016, 18 (48), 32847-32861. 20.

Friesner, R. A.; Guallar, V., Ab initio quantum chemical and mixed quantum

mechanics/molecular mechanics (QM/MM) methods for studying enzymatic catalysis. Annu Rev Phys Chem 2005, 56, 389-427. 21.

Senn, H. M.; Thiel, W., QM/MM methods for biomolecular systems. Angew Chem Int Ed

Engl 2009, 48 (7), 1198-1229. 22.

Warshel, A., Computer simulations of enzyme catalysis: methods, progress, and insights.

Ann Rev Biophys Biomol Struct 2003, 32, 425-443. 23.

Kamerlin, S. C.; Haranczyk, M.; Warshel, A., Progress in ab initio QM/MM free-energy

simulations of electrostatic energies in proteins: accelerated QM/MM studies of pKa, redox reactions and solvation free energies. J Phys Chem B 2009, 113 (5), 1253-1272. 24.

Vreven, T.; Byun, K. S.; Komáromi, I.; Dapprich, S.; Montgomery, J. A.; Morokuma, K.;

Frisch, M. J., Combining quantum mechanics methods with molecular mechanics methods in ONIOM. J Chem Theory Comput 2006, 2 (3), 815-826. 25.

Chung, L. W.; Sameera, W. M.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.;

Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; et al. Morokuma, K., The ONIOM method and its applications. Chem Rev 2015, 115 (12), 5678-5796. 26.

Chai, J. D.; Head-Gordon, M., Long-range corrected hybrid density functionals with damped

atom-atom dispersion corrections. Phys Chem Chem Phys 2008, 10 (44), 6615-6620. 27.

Zhao, Y.; Truhlar, D. G., The M06 suite of density functionals for main group

thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other

31



Page 32 of 36

functionals. Theoretical Chemistry Accounts 2008, 120 (1), 215-241. 28.

Barducci, A.; Bussi, G.; Parrinello, M., Well-tempered metadynamics: a smoothly

converging and tunable free-energy method. Phys Rev Lett 2008, 100 (2), 020603. 29.

Slater, J. C.; Koster, G. F., Simplified LCAO Method for the Periodic Potential Problem.

Phys Rev 1954, 94 (6), 1498-1524. 30.

Elstner, M.; Seifert, G., Density functional tight binding. Philos Trans A Math Phys Eng Sci

2014, 372 (2011), 20120483. 31.

Gaus, M.; Cui, Q.; Elstner, M., DFTB3: Extension of the self-consistent-charge density-

functional tight-binding method (SCC-DFTB). J Chem Theory Comput 2012, 7 (4), 931-948. 32.

Iismaa, S. E.; Holman, S.; Wouters, M. A.; Lorand, L.; Graham, R. M.; Husain, A.,

Evolutionary specialization of a tryptophan indole group for transition-state stabilization by eukaryotic transglutaminases. Proc Natl Acad Sci U S A 2003, 100 (22), 12636-12641. 33.

Abdul, S.; Leebeek, F. W.; Rijken, D. C.; Uitte de Willige, S., Natural heterogeneity of α2-

antiplasmin: functional and clinical consequences. Blood 2016, 127 (5), 538-545. 34.

Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng,

E. C.; Ferrin, T. E., UCSF Chimera--a visualization system for exploratory research and analysis. J Comput Chem 2004, 25 (13), 1605-1612. 35.

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E.,

GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19-25. 36.

Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.;

Cieplak, P.; Luo, R.; Lee, T.; et al. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J Comput Chem 2003, 24 (16), 1999-2012. 37.

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.,

Comparison of simple potential functions for simulating liquid water. J Chem Phys 1983, 79 (2),

32


Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


926-935. 38.

Dominguez, C.; Boelens, R.; Bonvin, A. M., HADDOCK: a protein-protein docking

approach based on biochemical or biophysical information. J Am Chem Soc 2003, 125 (7), 17311737. 39.

van Zundert, G. C.; Bonvin, A. M., Modeling protein-protein complexes using the

HADDOCK webserver "modeling protein complexes with HADDOCK". Methods Mol Biol 2014, 1137, 163-179. 40.

van Zundert, G. C. P.; Rodrigues, J. P. G. L.; Trellet, M.; Schmitz, C.; Kastritis, P. L.;

Karaca, E.; Melquiond, A. S. J.; van Dijk, M.; de Vries, S. J.; Bonvin, A. M. J. J., The HADDOCK2.2 web server: User-friendly integrative modeling of biomolecular complexes. J Mol Biol 2016, 428 (4), 720-725. 41.

Bussi, G.; Donadio, D.; Parrinello, M., Canonical sampling through velocity rescaling. J

Chem Phys 2007, 126 (1), 014101. 42.

Parrinello, M.; Rahman, A., Polymorphic transitions in single crystals: A new molecular

dynamics method. J Appl Phys 1981, 52 (12), 7182-7190. 43.

Nosé, S.; Klein, M. L., Constant pressure molecular dynamics for molecular systems.

Molecular Physics 1983, 50 (5), 1055-1076. 44.

Darden, T.; York, D.; Pedersen, L., Particle mesh Ewald: An N⋅log(N) method for Ewald

sums in large systems. J Chem Phys 1993, 98 (12), 10089-10092. 45.

Tao, P.; Schlegel, H. B., A toolkit to assist ONIOM calculations. J Comput Chem 2010, 31

(12), 2363-2369. 46.

Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling,

C., ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J Chem Theory Comput 2015, 11 (8), 3696-3713. 47.

Gaus, M.; Goez, A.; Elstner, M., Parametrization and benchmark of DFTB3 for organic

molecules. J Chem Theory Comput 2013, 9 (1), 338-354.

33



48.

Page 34 of 36

Gaus, M.; Lu, X.; Elstner, M.; Cui, Q., Parameterization of DFTB3/3OB for sulfur and

phosphorus for chemical and biological applications. J Chem Theory Comput 2014, 10 (4), 15181537. 49.

Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J., A new ONIOM

implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives. Dedicated to Professor Keiji Morokuma in celebration of his 65th birthday.1. Journal of Molecular Structure: THEOCHEM 1999, 461-462, 1-21. 50.

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman,

J. R.; Montgomery, Jr., J. A..; Vreven, T.; Kudin, K. N.; Burant, J. C.; et al. Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT, 2004. 51.

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman,

J. R.; Scalmani, G.; Barone, V.; Petersson, G. A., et al. Gaussian 09, Revision B.01. Gaussian, Inc., Wallingford, CT, 2009. 52.

DPlot, https://www.dplot.com/index.htm. (accessed April 15, 2019)

53.

Kubař, T.; Welke, K.; Groenhof, G., New QM/MM implementation of the DFTB3 method

in the gromacs package. J Comput Chem 2015, 36 (26), 1978-1989. 54.

Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G., PLUMED 2: New

feathers for an old bird. Computer Physics Communications 2014, 185 (2), 604-613. 55.

Laio, A.; Parrinello, M., Escaping free-energy minima. Proc Natl Acad Sci U S A 2002, 99

(20), 12562-12566. 56.

Branduardi, D.; Bussi, G.; Parrinello, M., Metadynamics with adaptive gaussians. J Chem

Theory Comput 2012, 8 (7), 2247-2254. 57.

Marinescu, A.; Cleary, D. B.; Littlefield, T. R.; Maurer, M. C., Structural features

associated with the binding of glutamine-containing peptides to Factor XIII. Arch Biochem Biophys 2002, 406 (1), 9-20. 58.

Polgár, L., Mercaptide-imidazolium ion-pair: the reactive nucleophile in papain catalysis.

34


Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FEBS Lett 1974, 47 (1), 15-18. 59.

Sluyterman, L. A.; de Graaf, M. J., The fluorescence of papain. Biochim Biophys Acta 1970,

200 (3), 595-597. 60.

Johnson, F. A.; Lewis, S. D.; Shafer, J. A., Determination of a low pK for histidine-159 in

the S-methylthio derivative of papain by proton nuclear magnetic resonance spectroscopy. Biochemistry 1981, 20 (1), 44-48. 61.

Lowe, G., The structure and mechanism of action of papain. Philos Trans R Soc Lond B Biol

Sci 1970, 257 (813), 237-248.

35



TOC graphic

36


Page 36 of 36

First Step of the Transglutaminase Reaction Catalyzed by Activated

Recommend Documents