Thermodynamic characterization of a macrocyclic Zika virus NS2B/NS3 protease inhibitor and its acyclic analogs

Cyclization of small molecules is a widely applied strategy in drug design for ligand optimization to improve affinity, as it eliminates the putative need for structural preorganization of the ligand before binding, or to improve pharmacokinetic properties. In this work, we provide a deeper insight into the binding thermodynamics of a macrocyclic Zika virus NS2B/NS3 protease inhibitor and its linear analogs. Characterization of the thermodynamic binding profiles by isothermal titration calorimetry experiments revealed an unfavorable entropy of the macrocycle compared to the open linear reference ligands. Molecular dynamic simulations and X‐ray crystal structure analysis indicated only minor benefits from macrocyclization to fixate a favorable conformation, while linear ligands retained some flexibility even in the protein‐bound complex structure, possibly explaining the initially surprising effect of a higher entropic penalty for the macrocyclic ligand.

parts of the S2 and S3 binding pockets. [18][19][20] In solution, both conformations coexist in equilibrium. [21][22][23] Upon ligand or substrate binding to the active site, exclusively cocrystals of the closed conformation were obtained so far.

| Inhibitor design and optimization
Due to the shallow shape of the NS2B/NS3 active site [24] and its preference for basic P1 and P2 residues, the development of drug candidates targeting this site is considered to be very challenging. [25] Consequently, most active site-directed ZIKV NS2B/NS3 protease inhibitors are based on peptidic or peptidomimetic scaffolds with often limited bioavailability featuring multibasic scaffolds. [12,[26][27][28][29] Although they reached considerably high affinities, these compounds exhibited only poor antiviral potency most likely caused by limited membrane permeabilities. [30][31][32][33] The crystal structure of a linear substrate analog inhibitor in complex with the closely related West Nile virus NS2B/NS3 protease revealed a horseshoe-like backbone conformation with both termini in close proximity to each other (Protein Data Bank [34] [PDB]-ID: 2YOL; Supporting Information: Figure S1). [35] This gave rise to a macrocyclization approach not only intended to optimize the inhibitory potency but also increase membrane permeability resulting in compound 1 (Table 1). [27] In the cocrystal structure of inhibitor 1 in complex with ZIKV NS2B/NS3 (PDB-ID: 6Y3B, Figure 1), the P1 guanidine residue is located at the bottom of the S1 pocket and forms an ionic interaction with Asp129. The carbonyl oxygen of Gly159 forms a hydrogen bond directly to the guanidine as well as a water-mediated hydrogen bond. Other H 2 O-mediated interactions are formed by Asp129 to the backbone nitrogen of the glycine linker. An intramolecular interaction between the P4 carbonyl oxygen and the terminal guanidine nitrogen stabilizes the inhibitor's conformation. This interaction can also be found in the linear reference inhibitor (Supporting Information: Figure S1). The hydroxy group of Tyr161 stabilizes the position of the backbone of inhibitor 1 by forming hydrogen bonds between its hydroxy group to the nitrogen of the P1-P2 amide and the backbone carbonyl oxygen of P3 Lys that further interacts with the Gly153 nitrogen. The side chain of the P2 Lys forms polar contacts to the NS2B Asp83 and NS3 Asn152 side chains and to the Gly82 backbone in the S2 pocket. The side chain of P3 Lys is stabilized by polar interactions with the carbonyl oxygen of Phe84 and the S3 forming side chains of Asp83 and Ser85. It is noteworthy that the NS2B residues Asp83 and Ser85 are found in two conformations based on their occupancies in the crystal structure. Therefore, Asp83 can contribute to the stabilization of both the P2 Lys and the P3 Lys side chain.
The complex formation of an oligo-peptidic ligand with a protein is often accompanied by the fixation of rotatable bonds upon binding. [18] This reduction of conformational degrees of freedom is usually accompanied by an entropic penalty and, thus, lowers binding affinity.
In ligand optimization attempts, a widely accepted strategy to increase affinity is the rigidification of rotatable bonds to preorganize the ligand in the binding conformation via (macro-) cyclization. [18,36] Although it is tempting to assume that mainly entropic effects lead to affinity enhancement by cyclization, also enthalpic benefits have been reported, which has led to some controversy about the underlying physical principles. [37][38][39][40] In the case of peptidomimetics, macrocyclization, T A B L E 1 Affinity and buffer-corrected thermodynamic binding profiles obtained by the fluorometric enzyme activity assay and ITC bridging, and cross-linking represent promising strategies to achieve more potent ligands with several benefits. [27,41,42] Besides higher proteolytic stability, [28] properly restrained molecules are less likely to adopt conformations to fit the binding pockets of off-target proteins; thus, macrocyclization offers a route to achieve more selective compounds. [43] Additionally, by derivatization of polar terminal groups and reduction of the number of rotatable bonds, macrocyclization offers a convincing benefit to membrane permeability. [29,30,44,45] Although the preorganization of ligands is a regularly utilized technique, the driving forces behind the improved affinities remain largely elusive. In rational drug design campaigns, isothermal titration calorimetry (ITC) is an indispensable tool to characterize binding thermodynamics. It allows to quantify the enthalpic (ΔH) and entropic (-TΔS) contributions to the overall binding free energy (ΔG) and to evaluate the impact of structural modifications on these parameters.

| Thermodynamic binding profiles
In ITC experiments, the same affinity trend was observed but with significantly higher K d values compared to K i (Table 1), as occasionally observed when comparing different techniques. [46,47]  potency. [46,47] However, this hypothesis requires further elucidation and is beyond the scope of this manuscript. The buffer-corrected thermodynamic binding profiles (Table 1, buffer correction calculation and signature plots: Supporting Information: Figure S2  effects. [48,49] For inhibitors 1-3, negative ΔC P -values of -3.77, -2.16, and -2.88 kJ·mol -1 ·K -1 were determined (Supporting Information: Figure S3). Differently, for compound 4, no linear correlation between ΔH and temperature was observed. While the negative ΔC P -values of 1-3 indicate the burial of hydrophobic areas, the absolute values are very large for pure protein-ligand interactions. [50][51][52] Most likely, the shifted equilibrium between open and closed conformations upon active site-directed ligand binding is also accompanied by burial of hydrophobic areas. This is also suggested by a decrease in the fraction of hydrophobic accessible surface area from 52% in the open to 48% in the closed ZIKV NS2B/NS3 conformation and seems to be a common feature for flaviviral proteases (Supporting Information: Table S3). In the closed (active) conformation, the NS2B cofactor is tightly wrapped around NS3 while being more solvent exposed in the open conformation. Hence, the largest ΔC P of compound 1 compared to 2 and 3 could be interpreted in different ways, like the stronger burial of hydrophobic patches of the ligand or a more efficient shift of the conformational equilibrium toward the closed conformation. However, ITC alone cannot finally answer that question and several layered effects may contribute to the observed results. [53][54][55]

| Molecular dynamics
To further elucidate this binding behavior, molecular dynamic simulations (MDs) were conducted starting from the ZIKV NS2B/ NS3-1 complex structure (PDB-ID 6Y3B). [27] For inhibitors 2-4, investigated. In analogy to the SPAM method, [56] which is used to estimate binding thermodynamics of water molecules from MD simulations, the distributions of ligand torsion energies were computed ( Figure 2c). One should be aware that the transfer of this method from water to a more complex system neglects many parts of molecular recognition and might not be used as a  Table S4, Figure S4a,b). Subsequently, order parameters (S²) for both backbone and sidechains were investigated (Supporting Information: Table S5, Figure S4c,d). S² are associated with the conformational entropy of residues [57] with S² = 1 indicating low flexibility and S² = 0 indicating high flexibility. [58,59] Backbone S² are fairly in line with the RMSF with only Gly159 and Ser160 being a little more flexible in the NS2B/NS3-2 complex. Sidechain S² showed overall more differences but no obvious trends that discriminate inhibitors 1 and 4 binding from inhibitors 2 and 3.
However, even though differences in S² for certain residues both within the binding site and more distant were observed, the corresponding configurational entropies of the residues seem to approximately cancel each other out. Therefore, the analysis of protein residue dynamics remains inconclusive, but slight differences might not hold the explanation for the observed thermodynamic binding profiles from the ITC experiments.

| X-ray crystallography
These findings were further supported by crystal structure analysis.
By the time, co-crystal structures of inhibitors 2-4 in complex with NS2B/NS3 could be obtained (Supporting Information: Table S6) (Table 1). Finally, 4 is the only ligand containing an acidic, negatively charged moiety reducing its potency as the active site has a highly negative potential (Supporting Information: Figure S5). Furthermore, strong interactions often result in a loss of conformational degrees of freedom for both ligands and proteins and subsequently are accompanied by an enthalpy-entropy compensation. [53,60] Therefore, B-factors as a surrogate for residue flexibility of the NS2B/NS3-inhibitor complexes were further analyzed. While raw B-factors are highly influenced by resolution, crystal packing, or refinement methods used, B-factors were normalized with the BANΔIT webserver for comparability between the different structures. [16] This analysis of normalized B′-factors indicated that some binding site residues in the complexes with inhibitors 2-4 retain higher flexibility (high B′-factors) while binding-site residues are slightly more rigid (lower B′-factor) for the NS2B/NS3-1 complex.
This holds true in general, but especially for NS3 residues 129-132 and Gly159 (S1 pocket), Val154, Val155, Ser160, and Tyr161 (close to the aromatic linker) and residues 83-86 of the S3 pocket formed by NS2B (Figure 3, Supporting Information: Figure S4E, Supporting Information: Table S7). All these residues are located around the termini of the acyclic ligands and the linking position of macrocycle 1. residues 64-68 and NS3 residues 29-32, Supporting Information: Figure S4e). [16] Therefore, it can be presumed that the highly similar binding modes of 1-4 result in similar global effects for protein dynamics.

| Reagents
All reagents and solvents were purchased from Sigma-Aldrich Chemie GmbH, Thermo Fisher Scientific Inc., and Carl Roth GmbH + Co. KG.

| Synthesis
The synthesis and analytical characterization of inhibitors 1-4 and the substrate PhAc-LKKR-AMC are described in previous publications. [27,35] The InChI codes of inhibitors 1-4, together with some biological activity data, are provided as the Supporting Information.

| Recombinant protein expression and purification
The bivalently expressed ZIKV protease construct NS2B/NS3 (bZiPro) for enzyme inhibition assays and ITC experiments was expressed and purified, as described previously. [26,61] Briefly, the vector (pETDuet) containing bZiPro (#86846, www.addgene.com) was transformed into competent Escherichia coli BL21 Gold (DE3) cells (Agilent Technologies) and grown in the LB medium containing 100 mg·L -1 ampicillin at 37°C   Figure S7). Since the K d value for inhibitor 4 binding to NS2B/NS3 was too low to give a sufficient slope at the inflection point, a low-c titration [63] was performed with a 20-fold molar

| MD simulations
MD simulations were performed starting from the ZIKV NS2B/ NS3-1 complex structure (PDB-ID 6Y3B). [27] An apo structure was generated by the removal of the ligand, while complexes with the ligands 2-4 were generated by manipulation (removal or addition of atoms) of the reference ligand 1 within the complex and subsequent energy minimization of the newly introduced ligand atoms within MOE [64] using the Merck Molecular Force Field (MMFF94). [65] The ligands were parameterized using the Generalized Amber Force Field (GAFF2 [66] ) with AM1-BCC [67] charges within antechamber [68] of the AmberTools20. [69] Complex structures were subsequently built with tleap [69] including crystallographic solvent molecules. After initial relaxation over 200 time steps with sander, counter ions (Na + for the protein-ligand complexes, Cl − for ligand simulations) were added and a TIP3P [70] waterbox exceeding the structure by 10.0 Å was built with tleap.
MDs were performed using NAMD2.14 [71] and the AMBER forcefield (ff14SB [72] ). The system was equilibrated over 1 ns by heating from 100 to 300 K over 500 ps and releasing harmonic constraints on the protein and ligand atoms over the following 500 ps in a constant volume box. MD production runs were performed over 10 ns with an NPT ensemble using periodic boundary conditions and a van der Waals cut-off of 14.0 Å with time steps of 2 fs allowing rigid bond lengths. For torsion energy analysis, ligand simulations without protein were conducted under the same conditions. Trajectories were written every ps and concatenated to include every 10th frame with catdcd 4.0 to result in 1000 frames per 10 ns prior analysis in VMD-1.9.3. [73] Order parameters S² were calculated using the isotropic reorientational eigenmode dynamics (iRED) approach within cpptraj [74] of the AmberTools20. [69]  Diffraction intensities were collected at BESSY MX beamline 14.1. Data processing and scaling were performed with the XDS program package. [75] Structure determination was done by molecular replacement with PHASER MR using the bZiPro structure with the PDB code 5GPI as a model. [76] Simulated annealing was performed within the initial refinement with PHENIX to remove potential bias from the search model. [77,78] The structures were subjected to alternating rounds of manual rebuild using Coot [79] to fit amino acids into σ-weighted 2F o -F c and F o -F c electron density maps and the PHENIX refine program (five cycles). For the calculation of R free , 5% of all data were randomly chosen and were not considered during refinement. Isotropic B-factor refinement with TLS parameters was used. Water and inhibitor molecules were located in the electron density and gradually added to the model. Multiple conformations were built if the minor populated conformation showed at least 20% occupancy. Coordinates and restraints for the inhibitors were generated with the Grade Web server. [80] The structures were deposited in the Protein Data Bank with accession codes 7ZLD (2), 7ZLC (3), and 7ZMI (4). The final data collection and refinement statistics are summarized in Supporting Information: Table S7.
B-factors were normalized with the BANΔIT-webserver [16] to analyze differences in protein binding-site flexibility in complex with the ligands 1-4 using the modified z-score method for C α atoms of the protein backbone.