Evaluating the covalent binding of carbapenems on BlaC using noncovalent interactions

Carbapenems, as irreversible covalent binders and slow substrates to the class A β-lactamase (BlaC) of Mycobacterium tuberculosis, can inhibit BlaC to hydrolyze the β-lactam drugs which are used to control tuberculosis. Their binding on BlaC involves covalent bonding and noncovalent interaction. We introduce a hypothesis that the noncovalent interactions dominate the difference of binding free energies for covalent ligands based on the assumption that their covalent bonding energies are the same. MM/GBSA binding free energies calculated from the noncovalent interactions provided a threshold with respect to the experimental kinetic data, to select slow carbapenem substrates which were either constructed using the structural units of experimentally identified carbapenems or obtained from the similarity search over the ZINC15 database. Combining molecular docking with consensus scoring and molecular dynamics simulation with MM/GBSA binding free energy calculations, a computational protocol was developed from which several new tight-binding carbapenems were theoretically identified.


Introduction
Mycobacterium tuberculosis, a pathogenic bacterial species in the family of Mycobacteriaceae, causes tuberculosis which is one of the top 10 causes of death and makes millions of people fall sick each year [1]. This is complicated by the emergence of multidrug-resistant and extensive drug-resistant strains of Mycobacterium tuberculosis, which reveals the drawback of the current treatment strategies for tuberculosis [1][2][3]. The βlactam class of antibiotics such as penicillin was effective in the treatment of tuberculosis due to their inhibition of bacterial cell wall synthesis [4]. However, the efficacy is reduced due to the resistance caused by class A β-lactamase (BlaC) which inactivates the β-lactam antibiotics via serine acylation followed by hydrolytic diacylation [5]. This involves the nucleophilic attack of the serine residue against the β-lactam carbonyl to form a covalent acyl-enzyme complex while the β-lactam ring is opened [6] followed by the hydrolysis of the ester bond to release the inactive ring-opened product [7]. In order to prevent this mechanism of resistance, it is necessary to find novel compounds to inhibit BlaC hereby being of utility in β-lactam/ β-lactamase inhibitor combinations such as meropenem/ vaborbactam [8] and meropenem/clavulanate [9].
Carbapenems also contain the β-lactam ring and can be used as antibiotics targeting penicillin-binding proteins [10]. They are poor substrates with slow hydrolysis on BlaC and could act as inhibitor of BlaC or an active partner with inhibitor clavulanate [9] such that they may act as a template for the development of novel inhibitor combinations. The crystal structures of BlaC-carbapenem complexes have been reported for some carbapenems including meropenem (MERO) [9], doripenem (DORI), and ertapenem (ERTA) (Fig. 1a) [11], which provide structural basis for the design of new carbapenems using computational approaches. Kinetic parameters k cat and K m have been measured for some carbapenems including the above three carbapenems, which are poor substrates with small k cat and K m [12]. The structures of these carbapenems contain two parts, as shown in Fig. 1b, the core structure in blue and the R group in black. All carbapenems have the same core structure but different R groups.
The core structure of the carbapenems contains two rings, the β-lactam ring with the substituent 6-α-1R-hydroxyethyl and the pyrroline ring to which the R group is attached via the thioether sulfur atom on C3 atom. There are four chiral centers in the core structure, i.e., the carbon atom C4 (the stereoisomer R) in the pyrroline ring, C5 (S) at the joint of the pyrroline ring and the β-lactam ring, C6 (S) in the β-lactam ring and C8 (R) in the substituent 6-α-1R-hydroxyethyl group (Fig. 1b). The core part completely lies in the binding site while, the R group spreads out of the binding site and is flexible (Fig. 1a). The β-lactam ring sterically blocks the binding site of BlaC and prevents its attack on the β-lactam class of antibiotics [9,11]. The binding contains both covalent and noncovalent interactions. Here, it should be pointed out that Ser84 is numbered according to the crystal structures [9,11], while Ser70 was used to indicate the same residue in the literature according to the Ambler position [6,7,13]. The noncovalent interactions involve the residues Ile117, Ser142, Glu182, Thr232, and Thr253 of BlaC.
The R groups cause the differences between these carbapenem molecules since they have the same core structure. Some carbapenems with different R groups were experimentally evaluated [12], which may be used to build models to develop more efficient inhibitors based on variations of the R group.
There can be many ways to construct the R group. In general, arbitrary permutations and combinations of any chemical structural units can generate a large number of structures. However, it might be restricted by chemical synthesis. In the present paper, we show how to build new R groups using the structural units of these experimentally identified carbapenems. Meanwhile, similarity search is used to find carbapenem-like structures. All these structures are then compared via analyzing the BlaC-carbapenem intermolecular interactions using computational tools with our filtering and scoring treatments.
Several computational methods already address covalent binding. CovalentDock [14] forms the covalent bond via a dummy atom, then computes the covalent bonding using a Morse potential which was fitted to the enthalpy change on bond length based on quantum chemical calculations. CovalentDock uses the AutoDock scoring function for the noncovalent interaction. Another covalent docking tool CovDock [15] forms the covalent bond using the OPLS3 force field for initial ligand binding poses generated by Glide docking sampling, then optimizes the linked structures using Prime, and ranks the binding poses using VSGB2.0 scoring function which is based on OPLS force field and Generalized Born solvation model. The covalent binding affinity can be computed in two individual states separately, i.e., the covalent binding and noncovalent binding states, using the free energy perturbation coupled with λ-exchange molecular dynamics method, while a QM/MM calculation of warhead was also performed to address the noncovalent interactions. The calculations of a series of αketoamide analogues relative to a common warhead scaffold revealed that covalent and noncovalent binding states of an inhibitor do not necessarily exhibit the same selectivity [16].
In the present work, we show how to use the noncovalent interactions to evaluate the covalent ligands and help design new covalent carbapenems. First, MD simulations were applied to the BlaC-carbapenem covalent systems to evaluate the role of R groups by analyzing the structural difference and energy contributions. Then, some new compounds were built using the structural units in the R groups. Meanwhile, a similarity search with respect to the core structure of carbapenem, performed over ZINC15 database [17], provided more compounds which were screened using molecular dockings. Furthermore, all selected compounds were consequently evaluated using MD simulations and the MM/GBSA binding free energy calculations. Finally, some compounds will be suggested for further experimental validation.

Structure preparation
The crystal structures of BlaC-carbapenem complexes were obtained from the Protein Data Bank (http://www.rcsb.org/) Fig. 1 a Superimposition of the three BlaC-carbapenem crystal structures. BlaC is shown in cartoon while the carbapenems are shown as sticks (DORI carbons in sky blue, ERTA carbons in plum and MERO carbons in tan). The binding site is expressed in surface using the structure of BlaC-MERO (PDB ID: 3DWZ). b Structural change during the covalent binding, i.e., the opening of β-lactam ring and the isomerization of pyrroline ring. c Structure used in molecular dockings [18], including the complexes with DORI (PDB ID: 3IQA) [11], ERTA (PDB ID: 3M6H) [11], and MERO (PDB ID: 3DWZ) [9], respectively. The root-mean-square deviations (RMSDs) between any pair of the three structures are very small (0.10~0.15 Å) with respect to either the Cα atoms of the whole backbone or the atoms of the binding site residues within 6 Å around the ligand (Table S1). Therefore, any of them, here 3DWZ, can be used as a representative to compute the binding of carbapenem ligands. The binding of carbapenem to BlaC causes a structural change which involves the formation of a covalent bond with the hydroxyl oxygen of Ser84, the opening of the β-lactam ring, and the isomerization of the pyrroline ring (Fig. 1b). The resulting structure is the conformation of carbapenem in the crystal and was used in our MD simulations.
The protonation of the protein was prepared using Tleap module in AMBER 14 [19]. The histidine tautomers were determined using our in-house program which use criteria similar to those used in literature such as that by Baker and Hubbard [20]. We use two criteria to judge if a histidine N atom will be protonated. The first criterion is the possibility to form a hydrogen bond (H-bond) with an H-bond acceptor at a proper position. The second criterion is if there is a spacious cavity to avoid repulsion with other possible hydrogen atom on an H-bond donor such as the amide nitrogen. A complementary convention is that the delta position is usually protonated if both N positions have no chance to form H-bond and the cavities are spacious. It should be very strict to assign the protonation to both N atoms unless they can form very strong H-bond with the surrounding H-bond acceptors or if the residue is adjacent to a negative, acidic group.

Parametrization
Parametrization for the covalently bonded residue-ligand structure was performed using tools implemented in the AMBER package. There is a covalent bond between the residue Ser84 and ligand carbapenem. Therefore, this residue should be a new type, named SEP herein. The new residue SEP does not have the hydroxyl hydrogen atom on its sidechain but leaves a linking site for forming a bond with ligand carbapenem. The partial atomic charges for SEP and ligand carbapenem were derived from the RESP (restrained electrostatic potential) charge model using the antechamber module implemented in AMBER, on the basis of the molecular electrostatic potential (ESP) calculated at HF/6-31G* level using Gaussian 09 [21] for a molecular model which contains the ligand carbapenem and the covalently bonded serine.
Here, this serine is capped with CH 3 CO-at its N terminal and with -NHCH 3 at its C terminal. The above ESP computation and RESP charge fitting were repeated for each individual SEP-ligand system. More details of the parameterization for these carbapenems and the covalently bonded serine SEP are given in Table S2.

MD simulation and binding free energy calculations
The protein-ligand systems were submitted to MD simulations using AMBER 14 package [19] with the AMBER ff14SB force field for protein [22] and the General Amber Force Field for ligand [23]. The process of MD simulation is described in Table S3. The MD trajectories were analyzed using CPPTRAJ module [24] implemented in AMBER 14 to derive the information about RMSD, root-mean-square fluctuation (RMSF), and H-bonds. The conformations in the trajectories were extracted for clustering using NMRCLUST [25] to obtain representative structures for further analysis of the binding modes in the active site. The details of clustering are given in Table S4.
Conventional binding free energy calculations via MD simulations cannot be directly applied to the covalently bonded system. Fortunately, all carbapenems have the same covalent bonding. Thus, the comparison between different carbapenems can be estimated by considering the noncovalent interactions. Therefore, we introduced a hypothesis that the difference of binding free energy for different carbapenems is determined by the noncovalent interactions if the covalent bonding energies are the same to all carbapenems. The binding free energy ΔG bind can be expressed as the following equation.
In the AMBER MM/GBSA binding free energy calculation, the receptor and ligand can be normally specified, no matter of covalent bonding or noncovalent interaction. By this way, the binding free energy for noncovalent interactions is a relative value and can be calculated using AMBER as follows: where ΔE internal is the change of the internal energy including the contributions from bonds, angles, and dihedrals. ΔE vdw is the non-bonded van der Waals interaction. ΔE ele is the electrostatic interaction. ΔG GB and ΔG SA represent the polar and nonpolar solvation energy, where GB stands for the generalized born model and SA for the solvent accessible surface area according to the MM/GBSA model [26] implemented in AMBER 14 [19],, using the model option GB2 [27,28].

Similarity search and molecular docking
For finding more similar compounds, a similarity search, with respect to the ring part of the core structure of carbapenem (SMILES: O=C1CC2CC=CN12, InChl Key in the ZINC scaffold set: YZBQHRLRFGPBSL-UHFFFAOYSA-N), was performed online (http://zinc15.docking.org/) using the default similarity criteria over the ZINC15 database which contains more than 750 million purchasable compounds [17]. A total of 243 similar molecules were obtained and further screened via a preliminary visual inspection to check if a similar compound contains the same core part and the same chiral atoms of carbapenem (Fig. 1b).
In order to estimate the noncovalent interactions, it is better to use noncovalent docking tools. The compounds obtained from the above design and similarity search were computed using two docking tools AutoDock Vina [29] and LeDock (http://www.lephar.com/) [30] with their default docking and scoring options. In order to avoid spatial blocking during molecular docking, the -C7 = O part of ligand was removed (Fig. 1c). The details of the noncovalent dockings are described in Table S5.
For balancing the scores obtained from multiple docking and scoring systems, normalization and consensus scorings were usually employed [31]. We constructed two consensus scoring functions as follows: where E Vina is the normalized AutoDock Vina score, E LeDock the normalized LeDock score, N the number of the ligand heavy atoms, and N the average of all N numbers. Score1 is an average of the two scores. Score2 is a scaled value with respect to the ligand size. Since the original scores from the docking outputs are negative, the normalization can be expressed as follows: where E is the score from the docking output for the considered ligand and E min and E max the minimum and the maximum of all scores for the calculated ligands. The factor 10 is used to tune the value to be close to the original scores. An ideal docking pose which should receive better scores simultaneously from the two consensus scoring functions will be selected for further MD simulations and binding free energy calculations.

Results and discussion
The flowchart in the present work is illustrated in Fig. 2. The basic process is that MD simulations and binding free energy calculations performed on the BlaC-carbapenem systems (Fig. 3) can provide a threshold for selecting hits. Then, the similarity search and combinations of structural units can be used to build new ligands which will be estimated using a series of computational treatments including molecular docking with consensus scoring and MD simulations with binding free energy calculations. The MD simulations were performed on the protein-ligand complexes with covalent bond which guarantees the stability of ligand in the binding site. The molecular dockings were performed on the noncovalent protein-ligand complexes, in which the covalent bond between protein and ligand was omitted and the atoms involved in the covalent bond were removed to avoid spatial repulsion during the noncovalent docking. The noncovalent interactions are mainly evaluated using MD simulation and binding free energy, while the docking scores can be used as preliminary criteria for filtering a larger pool of molecular structures.

Carbapenems
Based on the three known carbapenems, DORI, ERTA, and MERO (Fig. 3a), we built four compounds, named SL1, SL2, SL3, and SL4 (Fig. 3b), by combining the common parts of the R groups. SL1 contains only the core part of the carbapenem, without the R group. It will be used as a basis for comparison with those with different R groups. SL2 contains the pyrrolidin-3-ylthio moiety that is a common structural unit in the three R groups. SL3 is the common part of MERO and ERTA. SL4 contains the main part of DORI but removing the sulfamoxil-aminomethyl tail. Two more ligands Laura_7a and Laura_7h (Fig. 3c) were taken into consideration since they are also slow substrates with small k cat and K m , even smaller than those of the above three carbapenems [12]. The common part of Laura_7a and Laura_7h forms a new compound SL7 which contains the core part linked via a two-methylene chain to a triazole ring (Fig. 3d). SL8 and SL9 are formed from SL7 by deleting one-or two-methylene group, respectively (Fig.  3d). Another six compounds (SL10-SL15) (Fig. 3e) were built from the sub-structures of above twelve compounds.

Binding free energy and slow substrates
Three replicas of 50 ns MD simulations and the subsequent MM/GBSA binding free energy calculations were performed on each BlaC-ligand system. The averaged binding free energy derived from the three MD simulations will be used to evaluate the binding strength ( Table 1). The standard deviations for the binding free energy are small, indicating the convergence of the binding free energy calculations (Table S6). More details of the binding free energies are shown in Table S7.
In Table 1, four tight-binding slow substrates have more negative binding free energy, i.e., DORI (−44.14 kcal/mol), MERO (−36.44 kcal/mol), Laura_7a (−36.70 kcal/mol) and Laura_7h (−43.05 kcal/mol). The weaker binding substrate ERTA has a smaller binding free energy −17.95 kcal/mol. This indicates that a carbapenem analog can be judged using the binding free energy. A small change in the structure of the R group might cause a large change in the binding free energy. The R groups of DORI, ERTA, and MERO are similar. However, ERTA has a benzene ring with a carboxyl group and results at the less favorable binding free energy.
It should be pointed out that the MD simulations of DORI and ERTA were performed on their own crystal structures, while all other carbapenems were placed on the crystal structure from the BlaC-MERO complex. Therefore, DORI and ERTA were submitted again for three replicas of 50 ns MD simulations using the protein structure from the BlaC-MERO complex. The results are very close, i.e., − 43.83 vs. − 44.14 kcal/mol (Table 1)

SL1-SL4 and the binding free energy threshold
The structure of SL1 is just the core part of carbapenems (Fig.  3). The binding free energy of SL1 on BlaC is −36.26 kcal/ mol (Table 1), which can be used as a threshold to judge new carbapenems which contain different R groups (Fig. 3). An improved R group should receive a binding free energy more negative than this value. This value is already close to these values of the four good carbapenems. This means that the core part contributes the majority of binding free energy. SL2, SL3, and SL4 (Fig. 3b) are constructed from the structural units of DORI and MERO. However, their binding free energy values have not been improved, −31.47 kcal/mol for SL2, −34.55 kcal/mol for SL3 and − 33.81 kcal/mol for SL4, respectively (Table 1). It is interesting that these values are even less than the energy value of SL1. Comparing SL2-SL4 with DORI, it can be seen that the sulfamoyl group -NHSO 2 NH 2 or the sulfamoxil-aminomethyl group -CH 2 NHSO 2 NH 2 on DORI should be the important structural unit for contributing to the binding free energy.

SL7-SL9
We constructed SL7-SL9 (Fig. 3d) by removing the benzene ring and indole ring from the R groups of Laura_7a and Laura_7h (Fig. 3c). The triazole group in SL7-SL9 are kept. There is an improvement in SL7 which has a binding free energy of −37.83 kcal/mol (Table 1), more negative than that of SL1 (−36. 26 kcal/mol). Further, we reduced the size of R groups for SL8 and SL9 since a smaller group retains more space for further optimization which will usually increase the ligand size and present them here although the binding free energy values are temporarily not good for SL8 (− 34.00 kcal/ mol) and SL9 (−33.91 kcal/mol). Comparing the three Fig. 2 Flowchart of the carbapenem design and evaluation approach structures of similar size, SL8, SL9 and SL2, the energy values of SL8 and SL9 are better than that of SL2. It can be seen that the triazole group -C 2 H 2 N 3 on SL8 and SL9 should be an important structural unit, while the pyrrolidine group -C 4 H 9 N + on SL2 is not so favorable. This is definitely proven by Laura_7a and Laura_7h which both have the triazole group and favorable binding free energy.

SL10-SL15
On the basis of above structural analyses, we picked up two important structural units as building blocks for constructing new carbapenems. One structural unit is the sulfamoyl group or its extended unit sulfamoxil-aminomethyl group on DORI. Another is the triazole group on Laura_7a and Laura_7h. Considering three different lengths of the linker from the core Fig. 3 Structures of the carbapenems part, we can construct six new structures SL10-SL15 (Fig. 3e) via the combination of the triazole group and the sulfamoyl group or the sulfamoxil-aminomethyl group. The six carbapenem molecules will be evaluated later together with other compounds obtained from similarity search.

Molecular docking
The molecular dockings were performed on the above six new carbapenems SL10-SL15 (Fig. 3e) and 10 similar compounds Z1-Z10 (Fig. 4) which were obtained from the similarity search and visual inspection using the process described above. The names Z1-Z10 are the abbreviations of their ZINC IDs (Table S8). As references, the five experimental carbapenems were also docked, including DORI, ERTA, MERO, Laura_7a and Laura_7h. For each ligand, multiple top binding poses from each docking tool were extracted for further selection. One filter is the RMSD of the core part with respect to the X-ray structure of MERO (PDB code: 3DWZ), which should be smaller than 2 Å. Other filters are the consensus scores Score1 and Score2 which should be close to or even more negative than those of the identified MERO, i.e., Score1 < −4.93 kcal/mol and Score2 < −4.78 kcal/mol. Seven compounds meet the requirements, including SL10, SL11, SL12, SL13, SL14, SL15, and Z1 (Table 2), which will be further evaluated using MD simulations and MM/GBSA binding free energy calculations.

Stability, intermolecular interaction, and binding pose
Small RMSD values (Table S9) guarantee the stability of MD trajectory for binding free energy calculations. The RMSF values of ligands (Fig. S1~S3) indicate that the core structures are conserved in the binding site while the R groups are more flexible, associated with them being more exposed to the solvent.
During MD simulation, atoms are constantly moving within a small region and some inter-atomic distances keep changing. Thus, the number of H-bonds will constantly change. Therefore, we use the statistical number of H-bonds, which is the count of H-bonds divided by the number of snapshots (i.e., 5000). Table 4 provides the averaged numbers of Hbonds over the three MD simulations, which details for the three MD replicas are given in Table S10. The H-bonds can be categorized into three types, core-protein, R-protein, and ligand-water which mostly comes from the R-water part. A tight-binding carbapenem should have stable core-protein Hbonds, such as SL7, SL10, SL12, and SL15 which have about five core-protein H-bonds. A flexible carbapenem might have more H-bonds with water which tends to pull ligand out of the binding site, thereby reducing the binding free energy, such as ERTA which has about ten H-bonds with water. Figure S4 and Fig. S5 show the H-bonds and hydrophobic interactions for two examples SL15 and ERTA, based on their representative structures which are selected via clustering the conformations extracted from the MD trajectories. The detailed description of hydrogen bond and hydrophobic interaction for the two compounds SL15 and ERTA are shown in Tables S11 and S12, which are correlated to the binding poses in Fig. S4 and Fig. S5.
The designed compounds SL10-SL15 contain a sulfonyl group which S atom can act as a chalcogen bond donor. The chalcogen bond acceptor in protein can be a lone pair possessing atom-like O and N atoms, or a π system which can contain C atom. A chalcogen bond can be formed if the interatomic distance is less than the sum of the VDW radii and more than the sum of the covalent radii. We counted the  The binding poses of the eight selected carbapenems, SL7, SL10-SL15, and Z1, are expressed in Fig. 5 using the representative conformations which are superimposed with respect to the binding site residues. The locations of the ligands are marked using three residues, Ser84 in red, Thr231 in dodger blue and Glu292 in spring green, respectively. All core parts are stably lying at the central binding site near Ser84. The R groups, however, spread toward two different directions, i.e., SL11, SL13, and SL14 toward Thr231, but other five toward Glu292.

Alternate approach for atomic charge and molecular geometry
For the purpose of consistency, we used the method HF/6-31G* to compute the atomic charges (Table S2) since it was employed to generate the atomic charges of biomolecules in the traditional force fields. Molecular electrostatic potential and atomic charges are typically derived from the singledeterminant wavefunction using the HF method. It should be reasonable as long as the basis sets are sufficiently large. It is also interesting to apply hybrid density functional methods such as CAM-B3LYP [33] to the computation of molecular electrostatic potential for further deriving out atomic charges. Therefore, we calculated the RESP atomic charges of one sample model SEP-SL15 using the two methods HF/6-31G* and CAM-B3LYP/6-31G*. The two sets of atomic charges seem to be consistent (Tables S14 and S15). Based on the two sets of atomic charges, the binding free energy values of compound SL15 were calculated. These values are close to each other in the averages over the three replicas of 50 ns MD simulations (Table S16). In addition, for molecular docking, we also compared the binding poses of the designed compounds SL10-SL15 for two kinds of initial molecular geometries optimized using MMFF94x and CAM-B3LYP/6-31G* (Table S5). These docking poses possess similar orientations ( Fig. S6-Fig. S11).  Fig. 5 Binding poses of the eight compounds, SL7, SL10-SL15 and Z1. The figure was created using Chimera [32] Conclusions This research provides a feasible approach to design covalent drug, based on two issues: (1) a hypothesis that the noncovalent interactions can be used to compare the binding strength if the covalent bonding is the same to all ligands and (2) a computational protocol combining MD simulations, binding free energy calculations, similarity search, molecular docking, and consensus scoring. MD simulations and MM/ GBSA binding free energy calculations were first applied to the experimentally identified carbapenems to obtain a threshold for selecting novel carbapenem ligands. The structural units extracted from the experimental carbapenems were used to construct new carbapenem compounds. These structural units contain triazole ring, sulfamoxil-aminomethyl, and sulfamoyl groups. In addition, more compounds have been obtained from the similarity search over the ZINC15 database. All the new compounds have been evaluated using AutoDock Vina and LeDock dockings and our consensus scores, followed by MD simulations and binding free energy calculations. Finally, five new compounds SL7, SL10, SL12, SL14, and SL15 were suggested for further experimental validations. It was expected to find some compounds, but not any herein, from the similarity search over the ZINC database which contains more than ten million chemical compounds. It seems that our method will be helpful to design totally novel compounds that are not in a database. Of course, a huge database such as ZINC should be suitable for virtual screening of other scaffolds for designing non-lactam BlaC inhibitors which can be an alternate way to improve TB treatment [34][35][36][37]. Some β-lactam BlaC inhibitors were not computed in this work since they do not possess the same core structure [38].
The covalent bond is built between the serine hydroxyl and the lactam carbonyl using a new residue SEP with force filed parameters only suitable for these specific carbapenem structures. Opportunely, these parameters gave proper binding free energy values in agreement with the experimental data and were successfully applied to the MD simulations of these carbapenems. For a wide range of applications, these parameters need to be tested for more observable quantities using sophisticated protocols such as those steps used in the development of CGenFF [39] and its extension DGenFF which is developed using a polarizable force field based on classical Drude oscillator model [40].