Evolution of Ceftriaxone Resistance of Penicillin-Binding Proteins 2 Revealed by Molecular Modeling

Penicillin-binding proteins 2 (PBP2) are critically important enzymes in the formation of the bacterial cell wall. Inhibition of PBP2 is utilized in the treatment of various diseases, including gonorrhea. Ceftriaxone is the only drug used to treat gonorrhea currently, and recent growth in PBP2 resistance to this antibiotic is a serious threat to human health. Our study reveals mechanistic aspects of the inhibition reaction of PBP2 from the wild-type FA19 strain and mutant 35/02 and H041 strains of Neisseria Gonorrhoeae by ceftriaxone. QM(PBE0-D3/6-31G**)/MM MD simulations show that the reaction mechanism for the wild-type PBP2 consists of three elementary steps including nucleophilic attack, C–N bond cleavage in the β-lactam ring and elimination of the leaving group in ceftriaxone. In PBP2 from the mutant strains, the second and third steps occur simultaneously. For all considered systems, the acylation rate is determined by the energy barrier of the first step that increases in the order of PBP2 from FA19, 35/02 and H041 strains. Dynamic behavior of ES complexes is analyzed using geometry and electron density features including Fukui electrophilicity index and Laplacian of electron density maps. It reveals that more efficient activation of the carbonyl group of the antibiotic leads to the lower energy barrier of nucleophilic attack and larger stabilization of the first reaction intermediate. Dynamical network analysis of MD trajectories explains the differences in ceftriaxone binding affinity: in PBP2 from the wild-type strain, the β3-β4 loop conformation facilitates substrate binding, whereas in PBP2 from the mutant strains, it exists in the conformation that is unfavorable for complex formation. Thus, we clarify that the experimentally observed decrease in the second-order rate constant of acylation (k2/KS) in PBP2 from the mutant strains is due to both a decrease in the acylation rate constant k2 and an increase in the dissociation constant KS.


Introduction
Neisseria gonorrhoeae is a Gram-negative diplococci bacterium. It causes sexually transmitted infections, disseminated gonococcemia, septic arthritis, and gonococcal ophthalmia neonatorum. It took these bacteria 40 years to develop mechanisms of resistance to penicillin G. Further growth of resistance was more rapid [1,2]. In 2007, a list of recommended drugs for the treatment of infections caused by these bacteria was presented. However, after only three years, the only therapy that remained was ceftriaxone together with azithromycin. Emergence of resistance to azithromycin [3] resulted in ceftriaxone monotherapy as a major way of treatment in some countries. Still, failure of this therapy is already reported. Currently, there are several new compounds being clinically studied for treatment [4][5][6][7]. One of them, zoliflodacin, is in phase 3 trials and is considered as the most promising compound to replace ceftriaxone in case of further grown of Neisseria gonorrhoeae resistance.
There is a number of determinants associated with the mutation of genes responsible for emerging multiple resistance [8]. They occur in gene penA, which encodes an essential

Ceftriaxone Activation in the Active Site of Enzyme and Nucleophilic Attack
PBP2 is a serine hydrolase that shares common features with the enzymes of this type, EC 3. The active site carries a nucleophilic moiety that initiates acylation and an oxyanion hole that activates the substrate and stabilizes the tetrahedral intermediate [25][26][27]. In all considered PBP2 variants, the oxyanion hole is formed by NH groups of the main chains of Ser310 and Thr500 ( Figure 1B). These interactions of the substrate with the oxyanion hole result in polarization of the carbonyl group of the antibiotic and its activation as an electrophile. A Ser310 residue serves as a nucleophile, and the nearby Lys313 residue acts as a proton acceptor during the reaction ( Figure 1). Still, the active site of PBP2 FA19 differs from PBP2 from the mutant strains, which might be responsible for the k 2 /K s changes. The carboxyl group of ceftriaxone, which is conservative for β-lactam antibiotics, interacts with Thr498 and Ser362 residues in PBP2 FA19 ( Figure 1A). Both mutant PBP2 have G545S amino acid substitution that alters hydrogen bond partners of ceftriaxone carboxylate. In PBP2 from the mutant strains, Ser545 interacts with this group ( Figure 1A). This interaction shifts the position of the substrate closer to the exit from the binding pocket in PBP2 from the mutant strains. This influences the following reaction as shown below. formation. We perform electron density analysis of different states in ES complexes to explain how efficiency of substrate activation by enzyme affects the following nucleophilic attack. In addition, we perform dynamical network analysis to compare the antibiotic binding affinity by PBP2 from different strains.

Ceftriaxone Activation in the Active Site of Enzyme and Nucleophilic Attack
PBP2 is a serine hydrolase that shares common features with the enzymes of this type, EC 3. The active site carries a nucleophilic moiety that initiates acylation and an oxyanion hole that activates the substrate and stabilizes the tetrahedral intermediate [25][26][27]. In all considered PBP2 variants, the oxyanion hole is formed by NH groups of the main chains of Ser310 and Thr500 ( Figure 1B). These interactions of the substrate with the oxyanion hole result in polarization of the carbonyl group of the antibiotic and its activation as an electrophile. A Ser310 residue serves as a nucleophile, and the nearby Lys313 residue acts as a proton acceptor during the reaction ( Figure 1). Still, the active site of PBP2 FA19 differs from PBP2 from the mutant strains, which might be responsible for the k2/Ks changes. The carboxyl group of ceftriaxone, which is conservative for β-lactam antibiotics, interacts with Thr498 and Ser362 residues in PBP2 FA19 ( Figure 1A). Both mutant PBP2 have G545S amino acid substitution that alters hydrogen bond partners of ceftriaxone carboxylate. In PBP2 from the mutant strains, Ser545 interacts with this group ( Figure 1A). This interaction shifts the position of the substrate closer to the exit from the binding pocket in PBP2 from the mutant strains. This influences the following reaction as shown below. and PBP2 H041 (red) with bound ceftriaxone. Alignment is performed over nitrogen and carbon atoms of main chains of Ser310, Lys313 and Thr500. (B) The QM part of PBP2 H041 at the ES configuration; the leaving group of ceftriaxone (R 2 ) is highlighted by a yellow ellipse. The carbonyl group of the β-lactam ring and NH fragments of Thr500 and Ser310 residues forming the oxyanion hole are represented as balls and thick sticks. Ceftriaxone and other amino acid residues are depicted with thick and thin sticks, respectively. Here and in the following figures, carbon is green, oxygen is red, nitrogen is blue, hydrogen is white, and sulfur is yellow.
The efficiency of substrate activation in the enzyme-substrate complex can be quantified from molecular dynamics simulation with QM/MM potentials in different ways, including simple geometry criteria and more complicated electron-density-based descriptors. A detailed analysis of the geometry features in the active sites of the ES complexes of ceftriaxone with PBP2 from different strains demonstrates considerable variations of hydrogen bond distance between the NH fragments of Ser310 and the carbonyl The QM part of PBP2 H041 at the ES configuration; the leaving group of ceftriaxone (R 2 ) is highlighted by a yellow ellipse. The carbonyl group of the β-lactam ring and NH fragments of Thr500 and Ser310 residues forming the oxyanion hole are represented as balls and thick sticks. Ceftriaxone and other amino acid residues are depicted with thick and thin sticks, respectively. Here and in the following figures, carbon is green, oxygen is red, nitrogen is blue, hydrogen is white, and sulfur is yellow.
The efficiency of substrate activation in the enzyme-substrate complex can be quantified from molecular dynamics simulation with QM/MM potentials in different ways, including simple geometry criteria and more complicated electron-density-based descriptors. A detailed analysis of the geometry features in the active sites of the ES complexes of ceftriaxone with PBP2 from different strains demonstrates considerable variations of hydrogen bond distance between the NH fragments of Ser310 and the carbonyl oxygen atom of the substrate [28]. Practically, in PBP2 FA19 , the oxyanion hole is formed by two hydrogen bonds, and in PBP2 from the mutant strains, only one hydrogen bond with Thr500 is formed.
Here, we deepen the understanding of substrate activation by the enzyme, focusing on electron-density-based descriptors. We evaluate the Fukui electrophilicity index, f + , of the carbonyl carbon atom of ceftriaxone and Laplacian of electron density maps calculated in the plain of the carbonyl group of the substrate and a nucleophilic oxygen atom for sets of frames for each model system ( Figure 2B). oxygen atom of the substrate [28]. Practically, in PBP2 FA19 , the oxyanion hole is formed by two hydrogen bonds, and in PBP2 from the mutant strains, only one hydrogen bond with Thr500 is formed.
Here, we deepen the understanding of substrate activation by the enzyme, focusing on electron-density-based descriptors. We evaluate the Fukui electrophilicity index, f + , of the carbonyl carbon atom of ceftriaxone and Laplacian of electron density maps calculated in the plain of the carbonyl group of the substrate and a nucleophilic oxygen atom for sets of frames for each model system ( Figure 2B). First, we utilize a binary classifier, Laplacian of electron density, to evaluate the fraction of the reactive species that are characterized by the presence of the electron density depletion area on the carbonyl carbon atom in the direction of the nucleophilic attack [29][30][31]. The substrate activation is more efficient in PBP2 FA19 and equals 33%, whereas it decreases to 16% and 12% for PBP2 35/02 and PBP2 H041 , respectively. Figure 2B shows the violin plots of Fukui electrophilicity indices for ceftriaxone carbonyl carbon atom for PBP2 from three considered strains. For PBP2 FA19 -containing system, f + has a wider distribution with the larger average value compared with the other two PBP2. It is characterized by the 0.0117 ± 0.0053 a.u. mean value and standard deviation, whereas for PBP2 35/02 and PBP2 H041 , these values are 0.0047 ± 0.0018 and 0.0035 ± 0.0014 a.u., respectively. Larger Fukui atomic electrophilicity indices discriminate atoms with more pronounced electrophilic properties. Therefore, we can conclude that in PBP2 FA19 substrate, activation by the enzyme is considerably more pronounced. Thus, utilization of both spatially distributed Laplacian of electron density descriptor and atomic Fukui indices derives the consistent result that the substrate activation efficiency decreases in PBP2 from different strains in the following order: FA19 > 35/02 > H041.
We compare Gibbs energy profiles of the first step of the reaction to understand which particular features of the following reaction are affected by the efficiency of the substrate activation ( Figure 2A). The reaction coordinate at this step is a sum of distance between the proton transferring from the Ser310 residue side chain to the Lys313 residue side chain, d(NL…HS), and a distance of the nucleophilic attack, d(C…OS). The energy profile for PBP2 FA19 considerably differs from those for PBP2 from the mutant strains. The first intermediate is stabilized relative to the ES complex only in PBP2 FA19 . For PBP2 from both mutant strains, 2 kcal/mol destabilization is observed. The energy barrier is much lower for the system with PBP2 FA19 (4.6 kcal/mol), and it increases to PBP2 35/02 (8.4 kcal/mol) and PBP2 H041 (9.4 kcal/mol) in accordance with the decrease in the experimental k2/KS value. In addition, an important feature is the reaction coordinate at the ES minimum; this equals to 4.1 Å in PBP2 FA19 and around 4.5 Å in PBP2 from the mutant strains. This indicates that First, we utilize a binary classifier, Laplacian of electron density, to evaluate the fraction of the reactive species that are characterized by the presence of the electron density depletion area on the carbonyl carbon atom in the direction of the nucleophilic attack [29][30][31]. The substrate activation is more efficient in PBP2 FA19 and equals 33%, whereas it decreases to 16% and 12% for PBP2 35/02 and PBP2 H041 , respectively. Figure 2B shows the violin plots of Fukui electrophilicity indices for ceftriaxone carbonyl carbon atom for PBP2 from three considered strains. For PBP2 FA19 -containing system, f + has a wider distribution with the larger average value compared with the other two PBP2. It is characterized by the 0.0117 ± 0.0053 a.u. mean value and standard deviation, whereas for PBP2 35/02 and PBP2 H041 , these values are 0.0047 ± 0.0018 and 0.0035 ± 0.0014 a.u., respectively. Larger Fukui atomic electrophilicity indices discriminate atoms with more pronounced electrophilic properties. Therefore, we can conclude that in PBP2 FA19 substrate, activation by the enzyme is considerably more pronounced. Thus, utilization of both spatially distributed Laplacian of electron density descriptor and atomic Fukui indices derives the consistent result that the substrate activation efficiency decreases in PBP2 from different strains in the following order: FA19 > 35/02 > H041.
We compare Gibbs energy profiles of the first step of the reaction to understand which particular features of the following reaction are affected by the efficiency of the substrate activation ( Figure 2A). The reaction coordinate at this step is a sum of distance between the proton transferring from the Ser310 residue side chain to the Lys313 residue side chain, d(N L . . . H S ), and a distance of the nucleophilic attack, d(C . . . O S ). The energy profile for PBP2 FA19 considerably differs from those for PBP2 from the mutant strains. The first intermediate is stabilized relative to the ES complex only in PBP2 FA19 . For PBP2 from both mutant strains, 2 kcal/mol destabilization is observed. The energy barrier is much lower for the system with PBP2 FA19 (4.6 kcal/mol), and it increases to PBP2 35/02 (8.4 kcal/mol) and PBP2 H041 (9.4 kcal/mol) in accordance with the decrease in the experimental k 2 /K S value. In addition, an important feature is the reaction coordinate at the ES minimum; this equals to 4.1 Å in PBP2 FA19 and around 4.5 Å in PBP2 from the mutant strains. This indicates that in the wild-type PBP2, the ES complex is tighter, and this contributes to efficient substrate activation.
We analyze whether Fukui atomic electrophilicity indices depend on individual geometry parameters or are affected by their combination. We found no strong correlations between the individual interatomic distances and f + . Figure 2C depicts the scatter plot of a longer distance among two hydrogen bonds forming an oxyanion hole as a geometry parameter and Fukui index. The distributions of the distance of the longer hydrogen bond in the oxyanion hole practically do not overlap for PBP2 from the wild-type and mutant strains. Those are quantified by 1.97 ± 0.12, 2.92 ± 0.30 and 3.39 ± 0.22 Å for PBP2 FA19 , PBP2 35/02 and PBP2 H041 , respectively. We observe that the absence of a second hydrogen bond in the oxyanion hole almost absolutely eliminated states with high Fukui electrophilicity indices on the carbonyl carbon atom of the substrate as seen for PBP2 from the mutant strains ( Figure 2B,C). Still, the formation of two hydrogen bonds in the oxyanion hole does not guarantee efficient substrate activation, as there is no correlation between the f + and this interatomic distance for PBP2 FA19 .
To understand whether the geometry parameters together determine the extent activation of the substrate, we applied multiple linear regression analysis and random forest machine learning algorithms for Fukui indices and corresponding parameters for PBP2 FA19containing system. From multiple linear regression, we obtain the following relation This can estimate f + with the mean absolute value of 0.004 a.u. The distance of the nucleophilic attack comes with the positive coefficient. Distances of nucleophilic attack are relatively short in this system in different MD frames, and the influence is not obvious. For both hydrogen bond lengths, we find expected coefficients with negative signs; larger hydrogen bond distances decrease f + values.  [28], and its further elongation should have a more pronounced negative effect. The multiple linear regression demonstrates a relatively large mean absolute error of 0.004 a.u., indicating that the relation is far from linear. Therefore, we tried to obtain an implicit interrelation using a random forest approach. This allowed us to decrease the mean absolute error to 0.001 a.u. and estimate the importance of each geometry parameter. We found that the importance of all parameters is similar and ranges between 0.31 and 0.37. Utilization of only one geometry parameter in the random forest model leads to an increase in the mean absolute error to 0.002 a.u. for both distance of nucleophilic attack and hydrogen bond distances in the oxyanion hole. This is in line with the observations for the cysteine protease M pro from the SARS-CoV-2 virus such that both nucleophilic attack distance and oxyanion hole features should be considered [29].

Mechanism of Acyl-Enzyme Complex Formation in PBP2 from FA19 Strain
The pathway from the enzyme-substrate complex to the acyl-enzyme consists of three elementary steps in the case of PBP2 FA19 . The first step is a concerted process of proton transfer from the oxygen atom of the side chain of Ser310 residue to the side chain amino group of the Lys313 residue and formation of the covalent bond between the carbonyl carbon atom of ceftriaxone and an oxygen atom of the Ser310 side chain. The ES minimum on the Gibbs energy profile corresponds to the reaction coordinate 4.  Figure 3B). The third elementary step is cleavage of the covalent bond between the leaving group of ceftriaxone, R 2 , and the rest of ceftriaxone and proton rearrangements: the Ser362 proton returns to the amino group of Lys313, and the proton from the nitrogen of ceftriaxone returns to Ser362. Eliminating the leaving group provokes a redistribution of double bonds in the π-conjugated fragment of the substrate: a double bond between nitrogen and carbon of the thiazole ring is formed ( Figure 3C). The energy barrier value of the third step is 17 kcal/mol. d(HS362…OS362) is the proton transfer from an oxygen atom of Ser362. The Gibbs energy barrier of the second elementary step is 2.3 kcal/mol, and I2 is stabilized by 20 kcal/mol relative to the I1 ( Figure 3B). The third elementary step is cleavage of the covalent bond between the leaving group of ceftriaxone, R 2 , and the rest of ceftriaxone and proton rearrangements: the Ser362 proton returns to the amino group of Lys313, and the proton from the nitrogen of ceftriaxone returns to Ser362. Eliminating the leaving group provokes a redistribution of double bonds in the π-conjugated fragment of the substrate: a double bond between nitrogen and carbon of the thiazole ring is formed ( Figure 3C). The energy barrier value of the third step is 17 kcal/mol.

Mechanisms of Acyl-Enzyme Complex Formation in PBP2 from 35/02 and H041 Strains
As mentioned in Section 2.1 ( Figure 1A), in PBP2 from the mutant strains, the substrate bound to the active site is a different way that reflects the reaction mechanism as shown below. The overall acylation reaction occurs in two elementary steps in PBP2 from the mutant strains compared with three steps in PBP2 from the FA19 strain. The second step comprises simultaneous C-N bond cleavage of the β-lactam ring and elimination of the leaving group, R 2 , of the antibiotic. The reaction coordinate at this step is the sum of the distances d(C…N) and d(C…S). The energy barrier of the second elementary reaction step is 2.8 and 6.7 kcal/mol for PBP2 35/02 and PBP2 H041 , respectively ( Figure 4B). The overall acylation reaction leads to the stabilization of the I2 state by 12 and 4 kcal/mol for PBP2 35/02 and PBP2 H041 , respectively ( Figure 4B).

Mechanisms of Acyl-Enzyme Complex Formation in PBP2 from 35/02 and H041 Strains
As mentioned in Section 2.1 ( Figure 1A), in PBP2 from the mutant strains, the substrate bound to the active site is a different way that reflects the reaction mechanism as shown below. The overall acylation reaction occurs in two elementary steps in PBP2 from the mutant strains compared with three steps in PBP2 from the FA19 strain. The second step comprises simultaneous C-N bond cleavage of the β-lactam ring and elimination of the leaving group, R 2 , of the antibiotic. The reaction coordinate at this step is the sum of the distances d(C . . . N) and d(C . . . S). The energy barrier of the second elementary reaction step is 2.8 and 6.7 kcal/mol for PBP2 35/02 and PBP2 H041 , respectively ( Figure 4B). The overall acylation reaction leads to the stabilization of the I2 state by 12 and 4 kcal/mol for PBP2 35/02 and PBP2 H041 , respectively ( Figure 4B).

Dynamical Network Analysis of β3-β4 Loop
Molecular dynamic trajectories with a total length of 500 ns for each system were simulated to analyze conformational diversity introduced by amino acid substitutions. Dynamical network analysis was performed to identify communities, i.e., structural regions of the protein characterized by correlated motions. The method identified 12 communities in PBP2 FA19 and PBP2 35/02 and 11 communities in PBP2 H041 . Particular attention was paid on the β3-β4 loop that is supposed to be responsible for the stabilization of the noncovalent PBP2-ceftriaxone complex upon binding [24]. In PBP2 FA19  Step 1 Step 2

Dynamical Network Analysis of β 3 -β 4 Loop
Molecular dynamic trajectories with a total length of 500 ns for each system were simulated to analyze conformational diversity introduced by amino acid substitutions. Dynamical network analysis was performed to identify communities, i.e., structural regions of the protein characterized by correlated motions. The method identified 12 communities in PBP2 FA19 and PBP2 35/02 and 11 communities in PBP2 H041 . Particular attention was paid on the β 3 -β 4 loop that is supposed to be responsible for the stabilization of the noncovalent PBP2-ceftriaxone complex upon binding [24]. In PBP2 FA19 , the entire loop containing residues 501 to 515 is assigned to a separate community and shares common edges between the base of the loop (residues 501, 502, and 504) and the rest of the protein ( Figure 5). In PBP2 35/02 , residues 502 to 512 constitute a single community, and this community has six common edges with the protein core, including a common edge with a loop tip (residue 506). This shows that the movement of the β 3 -β 4 loop relative to the rest of the protein is limited. Moreover, at the tip of the β 3 -β 4 loop, there is a thick edge between residues 507 and 509, which indicates that the movement in the loop is highly correlated; the tip does not open and rotate as in PBP2 FA19 , but exhibits flat and stretched conformation. The β 3 -β 4 loop behavior is similar in PBP2 H041 and PBP2 35/02 . In PBP2 H041 , the difference from PBP2 FA19 is more pronounced, as a separate community is formed except in the 505 to 510 residues, and the edges are thicker, indicating a more pronounced correlation in this region ( Figure 5, PBP2 H041 , silver). In addition, the middle of the loop is assigned to another community, common with the protein core, and the tip of the loop has a common edge with the loop located under β 3 -β 4 ( Figure 5, PBP2 H041 , green). Thus, dynamical network analysis reveals that the β 3 -β 4 loop moves independently from the rest of the protein and demonstrates more twisted conformations in the in PBP2 FA19 in contrast to that from PBP2 from the mutant strains. Combining these data with experimental studies reposted in Ref [24], we can conclude that PBP2 FA19 should have higher affinity to ceftriaxone. network analysis reveals that the β3-β4 loop moves independently from the rest of the protein and demonstrates more twisted conformations in the in PBP2 FA19 in contrast to that from PBP2 from the mutant strains. Combining these data with experimental studies reposted in Ref [24], we can conclude that PBP2 FA19 should have higher affinity to ceftriaxone.

Discussion
According to the experimental studies, interactions of ceftriaxone and PBP2 are described by a three-step model: where E is the enzyme, S is an antibiotic, E·S is a noncovalent enzyme-substrate complex, E − S' is a covalent acyl-enzyme complex, and P is a hydrolyzed antibiotic. KS is a dissociation constant and characterizes the noncovalent complex formation, k2 and k3 are the rate constants of the acylation and deacylation steps, respectively [32]. According to the complex scheme of the experimental determination of k2/KS values for ceftriaxone reaction with PBP2 [32], k2 is referred to as the formation of the covalent complex between PBP2 and ceftriaxone. In the reaction mechanisms proposed here, the intermediates I1, I2 and I3 are in the system with PBP2 FA19 and I1, I2 are in the systems with PBP2 from the mutant strains, as a covalent adduct is already formed after the first elementary step. In other

Discussion
According to the experimental studies, interactions of ceftriaxone and PBP2 are described by a three-step model: where E is the enzyme, S is an antibiotic, E·S is a noncovalent enzyme-substrate complex, E − S' is a covalent acyl-enzyme complex, and P is a hydrolyzed antibiotic. K S is a dissociation constant and characterizes the noncovalent complex formation, k 2 and k 3 are the rate constants of the acylation and deacylation steps, respectively [32]. According to the complex scheme of the experimental determination of k 2 /K S values for ceftriaxone reaction with PBP2 [32], k 2 is referred to as the formation of the covalent complex between PBP2 and ceftriaxone. In the reaction mechanisms proposed here, the intermediates I1, I2 and I3 are in the system with PBP2 FA19 and I1, I2 are in the systems with PBP2 from the mutant strains, as a covalent adduct is already formed after the first elementary step. In other words, for the experimentally proposed mechanism, there is no particular composition of E-S with respect to the theoretically obtained intermediates. In PBP2 FA19 , all intermediates starting with I1 are stabilized relative to the ES; therefore, forward reactions occur with a higher rate than backward ones, and the first forward step should be considered when compared with the experimental k 2 value. For PBP2 from the mutant strains, I1 is destabilized relative to the ES; still, the following step resulting in the stabilized I2 state is characterized by a lower energy barrier. Therefore, for PBP2 35/02 and PBP2 H041 , the first forward step determines the overall process of acyl-enzyme complex formation. The energy barrier of the first forward step increases in the series of PBP2 FA19 , PBP2 35/02 and PBP2 H041 enzymes, indicating that the rate constant decreases in the same order.
Analysis of crystal structures further supports the proposed mechanisms and their differences between PBP2 from FA19 and mutant strains. For FA19, the crystal structure PDB ID 6P54 [33] presents a dimer, and in one monomer acyl-enzyme complex, there exists the leaving group R 2 (I2 in the mechanism on Figure 3), and in the other monomer, no R 2 fragment is observed (I3 in the mechanism on Figure 3). This is because the I2 state is stabilized and accumulates rapidly during the reaction. Then, it slowly proceeds to the I3 state without an R 2 fragment. I3 is less stable than I2; after the cleavage of the bond between the R 2 and the rest of ceftriaxone, leaving group dissociates to the solution, and thus, I2 to I3 transition may occur only in the forward direction. In contrast, in the PBP2 H041 crystal structure PDB ID 6VBD [22], the acyl-enzyme complex exists without a leaving group that is in line with the proposed mechanism ( Figure 4) with a direct formation of the stabilized state I2.
Explicit estimates of dissociation constants usually lead to large errors; therefore, we relied on experimental studies of the conformational diversity of PBP2 and performed molecular dynamic simulations. This allowed us to indirectly compare binding affinities of three considered PBP2. In long MD trajectories, we found that the behavior of the β 3 -β 4 loop differs in the considered systems. In wild-type PBP2, it mostly exists in conformation that facilitates the increase in binding affinity, whereas PBP2 from both mutant strains demonstrate conformations that are not favorable for efficient binding.
Thus, experimentally observed decrease of the k 2 /K S values in PBP2 from the mutant strains relative to the wild-type PBP2 is due to both decrease in the k 2 value and increase in the dissociation constant.

Methods
The complex of PBP2 FA19 with ceftriaxone was obtained from the crystal structure PDB ID 6P54 [33] with a resolution of 1.8 Å, and the complex of PBP2 H041 with ceftriaxone was obtained from the crystal structure PDB ID 6VBD [22] with a resolution of 1.8 Å. The crystal structure PDB ID 6VBL [22] of PBP2 35/02 apo-form with a resolution of 1.9 Å was used as the initial approximation for the coordinates. For all structures, amino acid residues 282-298 are not resolved in X-ray structures and therefore were reconstructed according to the primary amino acid sequences from UniProt. Ceftriaxone substrate was introduced in the active site of the apo-form of the crystal structure of PBP2 35/02 similarly to its complex with PBP2 H041 . Hydrogen atoms were added using the Reduce program [34] in such a way that the protonated forms of the amino acids with ionogenic groups corresponded to the neutral pH, except for the catalytic residue Lys313, which remained in the neutral form, being a proton acceptor from the catalytic residue Ser310 during the nucleophilic attack. Thus obtained full atomic models of the complexes were solvated by water molecules to rectangular boxes in which the distance from any atom of a protein to the cell boundaries was at least 12 Å. Then, chlorine anions were added to the system to neutralize the total charge of the cell. Complexes were equilibrated for 20 ns using the CHARMM36 [35,36] force field for the protein, TIP3P [37] for water molecules, and CGenFF [38] for ceftriaxone. Classical molecular dynamics (MD) simulations were performed at T = 300 K and p = 1 atm with the 1 fs integration time step. Calculations were performed using the NAMD program package [39]. The root mean-square deviation (RMSD) values along the obtained trajectories show that this length is enough for the complete relaxation of the system. Then, the representative frames from the last 2 ns of the MD runs were used to prepare the QM/MM models for calculation of reaction mechanisms.
The quantum subsystem includes the substrate molecule, catalytic residues Lys313 and Ser310, amino acid residues forming the oxyanion hole, interacting with substrate or catalytic residues, and water molecules ( Figure 1B). For PBP2 FA19 , they are Ser362, Asn364, Lys497, Thr498, and Thr500, and for PBP2 from mutated strains, they are Ser362, Asn364, Thr500 and Ser545. The quantum mechanical moiety was then described at density functional theory level in the Kohn-Sham variant with PBE0 hybrid functional [40] and D3 dispersion correction [41] in the 6-31G** basis set. The MD calculations were performed for 5 ps with the QM/MM potentials for each enzyme-substrate complex in the NAMD program package combined with the TeraChem [42] complex integrated for quantum chemical calculations [43]. The cutoff distance for point charges of the MM subsystem contributing to the QM Hamiltonian was 12 Å. PBE0 function is utilized, as it demonstrates consistent results for similar reactions initiated by the nucleophilic attack and comprising C-N bond cleavage [29][30][31]44].
We selected sets of 500 frames for each model system from 5 ps QM/MM MD simulations of the ES complex. We estimated reactivity of ceftriaxone using Laplacian of electron density maps and Fukui electrophilicity index, f + , of the carbonyl carbon atom. The spatial distributions of the Laplacian of the electron density were calculated in the plane of the carbonyl group of the substrate (C and O atoms) and the nucleophilic atom of Ser310, O S , to discriminate reactive and nonreactive moieties. Earlier works [29][30][31]44] demonstrate that this approach is a proper tool to visualize the substrate activation in nucleophilic addition reaction steps. More precisely, in nonreactive species, a carbonyl carbon atom is enveloped by the electron density concentration area. In reactive species, the carbonyl carbon atom has an electrophilic site that is an electron density deconcentration area located on the line of the nucleophilic attack. This deconcentration area is complementary to the nucleophilic site on the serine oxygen atom, being a lone electron pair oriented toward the electrophilic site. Apart from the spatially distributed Laplacian of electron density maps, an integral value of the Fukui atomic index can be calculated. The Fukui atomic index of electrophilicity, f + , is evaluated as a difference between Hirshfeld charges [45,46] calculated for the model system with N electrons (a system under consideration) and a system with the same geometry configuration, but with an additional electron. Larger Fukui atomic electrophilicity indices correspond to the more electrophilic species. The electron density analysis was performed in the Multiwfn program package [47].
The Gibbs energy profiles for each reaction step along the reaction pathway were calculated using the umbrella sampling approach [48,49]. The sets of 5-10 ps runs were performed with harmonic potentials centered at different values of reaction coordinates. The force constant of the harmonic potential 1 2 ·K·(ξ − ξ 0 ) 2 was usually set to 40 kcal/mol/Å 2 , and additional trajectories with the K = 80-120 kcal/mol/Å 2 in transition state regions were calculated in several runs. Harmonic potentials were centered every 0.2 Å along the reaction coordinates in regions close to minima and every 0.1 Å in transition state regions. Statistical analysis of reaction coordinate distributions was performed using both weighted histogram analysis method (WHAM) and umbrella integration (UI). The quality of distributions was monitored by their overlap and consistency of Gibbs energy profiles obtained in both WHAM and UI methods.
Classical MD simulations were carried out for 500 ns for the apo-form of proteins using the CHARMM36 force field for the protein and TIP3P for water molecules. Conformational diversity was analyzed using the dynamical network analysis, DNA [50]. According to this approach, a network is defined as a set of nodes connected by edges. In this work, every amino acid was represented by a single node. Any two nodes (except for the neighbors) were connected by an edge if the distance between any pair of atoms of the respective residues was less than 4 Å for more than 75% of the simulation time. Covariance and correlation matrices for dynamical network analysis were calculated with the CARMA program [51]. The MD calculations were carried out using the NAMD software package; trajectories analysis and visualization were performed in VMD program [52].

Conclusions
We propose molecular mechanisms of inhibition of PBP2 from wild-type FA19 and mutant 35/02 and H041 strains of Neisseria Gonorrhoeae by ceftriaxone due to the formation of the acyl-enzyme complex. In PBP2 FA19 , a reaction occurs via three elementary steps with rapid formation of the acyl-enzyme complex and slow elimination of the leaving group that is supported by X-ray data [33]. In PBP2 35/02 and PBP2 H041 from the mutant strains, a reaction happens in two elementary steps, as cleavage of the C-N bond in the β-lactam ring and elimination of the leaving group of the antibiotic occur simultaneously. Differences in the energy barrier of the first step of the reaction as well as subsequent stabilization of the tetrahedral intermediate are governed by the efficiency of the substrate activation in the ES complexes. It is quantified from the QM/MM MD trajectories of the ES complexes with subsequent electron-density-based analysis of features of the carbonyl carbon atom of ceftriaxone. Both Laplacian of electron density and atomic Fukui electrophilicity indices consistently demonstrate more efficient activation of ceftriaxone by PBP2 FA19 than by PBP2 from the mutant strains. Dynamical network analysis of classical MD trajectories reveals that ceftriaxone binding affinity is higher for PBP2 FA19 . Thus, we clarify that the increase in ceftriaxone resistance of PBP2 from the mutant strains evaluated from the decrease in k 2 /K S [22] can be explained by both decreases in the acylation rate constant, k 2 , and increases in the dissociation constant, K S .