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 are 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/molecular mechanics (MM) calculations: an ONIOM-based model as well as QM/MM/molecular dynamics (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 the 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.