Discovery of 5,5′-Methylenedi-2,3-Cresotic Acid as a Potent Inhibitor of the Chemotactic Activity of the HMGB1·CXCL12 Heterocomplex Using Virtual Screening and NMR Validation

HMGB1 is a key molecule that both triggers and sustains inflammation following infection or injury, and is involved in a large number of pathologies, including cancer. HMGB1 participates in the recruitment of inflammatory cells, forming a heterocomplex with the chemokine CXCL12 (HMGB1·CXCL12), thereby activating the G-protein coupled receptor CXCR4. Thus, identification of molecules that disrupt this heterocomplex can offer novel pharmacological opportunities to treat inflammation-related diseases. To identify new HMGB1·CXCL12 inhibitors we have performed a study on the ligandability of the single HMG boxes of HMGB1 followed by a virtual screening campaign on both HMG boxes using Zbc Drugs and three different docking programs (Glide, AutoDock Vina, and AutoDock 4.2.6). The best poses in terms of scoring functions, visual inspection, and predicted ADME properties were further filtered according to a pharmacophore model based on known HMGB1 binders and clustered according to their structures. Eight compounds representative of the clusters were tested for HMGB1 binding by NMR. We identified 5,5′-methylenedi-2,3-cresotic acid (2a) as a binder of both HMGB1 and CXCL12; 2a also targets the HMGB1·CXCL12 heterocomplex. In cell migration assays 2a inhibited the chemotactic activity of HMGB1·CXCL12 with IC50 in the subnanomolar range, the best documented up to now. These results pave the way for future structure activity relationship studies to optimize the pharmacological targeting of HMGB1·CXCL12 for anti-inflammatory purposes.


INTRODUCTION
Cell recruitment is a fundamental event in the establishment of both acute and chronic inflammatory responses (Charo and Ransohoff, 2006). Chemokines and their receptors organize leukocyte trafficking and migration to the tissues both in healthy and pathological conditions (Griffith et al., 2014). Chemokines bind to and activate seven-transmembrane G-protein-coupled receptors (GPCRs), a receptor superfamily involved in many different diseases (Hauser et al., 2017). The different chemokines and their receptors are objects of intense drug-discovery studies (Viola and Luster, 2008). Increasing evidence shows that chemokines usually require interactions with additional players that modulate the inflammation signaling cascade triggered by their receptors (Koenen et al., 2009;Proudfoot and Uguccioni, 2016;Eckardt et al., 2020); hence, the identification of molecules that disrupt the interactions of chemokines with their modulators might offer new and selective pharmacological opportunities against inflammation-related diseases (von Hundelshausen et al., 2017). A prototype of functional synergic heterophilic interaction is represented by the heterocomplex formed by High Mobility Group Box 1 (HMGB1) and the chemokine CXCL12, the 9 kDa ligand of the GPCR chemokine CXCR4 receptor . The CXCL12/CXCR4 axis is crucial for chemotaxis, cell arrest, angiogenesis, cell survival, and homing of hematopoietic progenitor cells in the bone marrow and their mobilization into the periphery both in physiological and pathological conditions (Tachibana et al., 1998;Alsayed et al., 2007;Pawig et al., 2015). The HMGB1·CXCL12 complex has been shown to trigger specific CXCR4 homodimer rearrangements , ERK activation, and calcium fluxes, along with enhanced CXCR4-dependent monocyte migration and tissue regeneration (Bianchi et al., 2017;Lee et al., 2018;Tirone et al., 2018). HMGB1 is a highly conserved, abundant nonhistone nuclear protein (215 amino acids) that behaves as a DNA chaperone. Upon extracellular release, it acts as a DAMP (Damage Associated Molecular Pattern) and triggers cytokinelike pro-inflammatory and chemoattractant effects (Bianchi et al., 2017;Gorgulho et al., 2019). HMGB1 is structurally organized in two independent L-shaped tandem domains (∼80 amino acids each), BoxA and BoxB, followed by a 30 amino acid long acidic C-terminal tail (Knapp et al., 2004;Stott et al., 2010). Its modular organization, flexibility, different oxidation states Bianchi et al., 2017), and its intracellular and extracellular localization confer to HMGB1 the ability to interact with many different partners, including nucleic acids, heparansulphates, lipopolysaccharides, and proteins, thus exerting different functions spanning from architectural chromatin-binding activity (Thomas and  to the modulation of innate immunity (Andersson and Tracey, 2011). The cytokine-like function of HMGB1 upon inflammatory stimuli after infection or injury relies on the interaction of its reduced form with CXCL12, which promotes CXCR4-dependent recruitment of inflammatory cells to injured tissues  and exacerbates the immune response in pathological conditions (Proudfoot and Uguccioni, 2016;D'Agostino et al., 2018). Thus, the disruption of the HMGB1·CXCL12 heterocomplex with small molecules might offer new and selective strategies against inflammation related diseases (Venereau et al., 2016;VanPatten and Al-Abed, 2018).
Interference with this heterophilic interaction is attractive but challenging for several reasons: (i) the three-dimensional structure of the complex is still unknown, and (ii) the interaction surface is expected to be large and dynamic, thus difficult to be targeted by small molecules De Leo et al., 2019;Fassi et al., 2019). Previous work has shown that it is possible to interfere with the pro-inflammatory properties of HMGB1 using small molecules including glycyrrhizin (Mollica et al., 2007), salicylic acid (SA), and its derivative amorfrutin (Choi et al., 2015). Moreover, we have recently shown that diflunisal, an aspirin-like nonsteroidal anti-inflammatory drug, is able to selectively interfere in vitro and in vivo with the HMGB1/CXCL12/CXCR4 inflammatory axis by disrupting the interaction between HMGB1 and CXCL12 (De Leo et al., 2019). This work has introduced the concept that the HMGB1·CXCL12 heterocomplex is a pharmacological target, herewith opening new perspectives for the rational design of novel inhibitors. While CXCL12 alone has been the object of structure-based drug design studies and virtual screening campaigns aiming at inhibiting its direct activation effect on CXCR4 (Veldkamp et al., 2010;Smith et al., 2016), systematic studies on the ligandability of HMGB1 and on the pharmacophoric requirements of HMGB1 ligands are still missing. This is in part ascribable to the Lshaped boxes of HMGB1, which do not appear to have clearly druggable pockets (Laraia et al., 2015;Scott et al., 2016), and to the absence of high resolution structures of HMGB1 in complex with known inhibitors.
In the present work we have explored the ligandability of HMGB1 and performed a structure-based virtual screening. The results were refined according to the pharmacophoric features of known HMGB1 ligands and tested for binding by NMR. This pipeline resulted in the identification of 5,5 ′ -methylenedi-2,3cresotic acid (2a) as a potent inhibitor of the chemotactic activity of HMGB1·CXCL12, with an IC 50 in the subnanomolar range, and one of the best ligands identified up to now.

In silico Assessment of HMGB1 Ligandability
Hot spots identification was performed using the FTMap computational map server (www.ftmap.bu.edu). The structures of BoxA and BoxB were uploaded into the FTMap server and ran according to instructions (Kozakov et al., 2015). The server uses 16 organic molecules as probes and defines as consensus sites (CSs) the regions of the protein where several different probes clusters bind. The results were visually inspected using PyMol (The PyMOL Molecular Graphics System, 2017).
Ligandability of HMGB1 boxes was further assessed using DoGSiteScorer (https://proteins.plus/). BoxA and BoxB structures were uploaded to the DoGSite server and scrutinized for binding sites, and their corresponding DrugScores were calculated according to the published protocol (Volkamer et al., 2012). Pocket Size and DrugScores (the closer the score to 1, the higher the ligandability) were extracted for all identified sites and annotated to pocket numbers.

Library Preparation
The Zbc Drugs subset (http://zinc12.docking.org/browse/ subsets/special) containing 101,746 biogenic drug-like compounds selected according to Lipinski's rule of five was chosen for virtual screening (VS) studies. The Zbc subset was prepared with the Ligand Preparation wizard available in the Schrödinger Suite, by adding hydrogens and minimizing and optimizing the molecules with OPLS_2005 force field. The ligands generated were ready for docking with GLIDE; for docking with AutoDock Vina (Trott and Olson, 2009) and AutoDock 4.2.6 ( Morris et al., 2009) they were converted to pdbqt-formatted file by in-house scripts. Ligands used to generate the training set for the pharmacophoric model (section hits clustering) were prepared with the same procedure.
GLIDE: BoxA and BoxB were prepared with the Protein Preparation Wizard available in the Schrödinger Suite. Water molecules were removed, hydrogen atoms were added, the protonation states were adjusted according to neutral pH, and finally the structures were minimized and optimized using OPLS2005 force field. For BoxA we generated a grid centered at 31.63, 27.3, 33.8 Å in correspondence to R23, the residue at the center of the experimentally validated binding pocket (De Leo et al., 2019). The size of the inner box, that is the ligand diameter midpoint box, was left at the default values of 10 Å edges; the outer box, that is the box within which all the ligand atoms must be contained, was enlarged to 46 Å edges to allow ligands to find unusual or asymmetric binding modes in the active site. The grid for BoxB was centered on 27.22, 12.96, 7.84 Å in correspondence to R109 (the equivalent position of R23), the inner box was set at 10 Å edges, the outer box was set at 29 Å edges and the option to "dock ligands with length ≤19 Å" was flagged. Glide High-Throughput Virtual Screening (HTVS) workflow was used for VS.
AutoDock VINA and AutoDock 4.2.6: BoxA and BoxB, prepared as described for Glide, were converted to pdbqtformatted file. Charges and non-polar hydrogen atoms were added using the prepare_receptor4.py script from MGLTools. The binding site was defined by AutoGrid. The BoxA grid size consisted in 40 points in each direction; the grid point spacing was set at 0.458 Å and was centered on 31.6, 27.3, and 30.70 on x, y, and z axes, respectively. BoxB grid size was 40 Å in each direction, grid point spacing 0.531 Å, and was centered on 20.5, 12.5, and 3.0 on x, y, and z axes, respectively. The max number of energy evaluations was 2,500,000, 27,000 was the max number of generations, and 10 hybrid genetic algorithm-local search GA-LS runs were performed.

Hits Selection
The hit candidates emerging from Glide, AutoDock Vina, and AutoDock 4.2.6 were selected according to energy ranking, distance filtering, visual inspection and absorption, distribution, metabolism, and excretion (ADME/Tox) descriptors through the QikProp module available in the Schrödinger Suite. In total, 34 properties and descriptors as defined in the Qikprop module were used to filter out the compounds outside the 95% of the value range reported for known drugs.

Pharmacophore Model
The pharmacophore model (PhMOD) was generated using the automatic pharmacophore generation protocol in LigandScout3 (Wolber and Langer, 2005) from Inte:Ligand. To generate the pharmacophore model we took advantage of a training set (TS) (Supplementary Figures 1A,B) composed by seven known active HMGB1-interacting molecules from the literature, i.e. glycyrrhizin (Mollica et al., 2007), carbenoxolone (Mollica et al., 2011), salicylic acid (SA) (Choi et al., 2015), amorfrutin (Choi et al., 2015), acetyl-3-aminoethyl-SA (Ac3AESA) (Choi et al., 2015), diflunisal (De Leo et al., 2019), inflachromene (Lee et al., 2014), and 11 decoy HMGB1 ligands with an activity >50 µM as tested in-house in cell migration experiments (acetaminophen, ibuprofen, naproxen, nimesulide, ketoprofen, folic acid, tetrahydro methotrexate, cortisol, cortisone, and prostaglandin E2) (Supplementary Figure 2). We superimposed the negatively charged moiety of the 3D-structures of the seven active members of the TS and the following four pharmacophore features were generated: five hydrophobic moieties, three aromatic rings four 4 H-bonding acceptors, and one H-bonding donor. The model was next refined maintaining only those features that previous structural studies on glycyrrhizin (Mollica et al., 2007) and diflunisal (De Leo et al., 2019) have indicated as relevant (Supplementary Figure 3A). The resulting PhMOD consisted of two H-bonding acceptors (HBA1 and HBA2) and two hydrophobics (H1 and H2). The tolerance sphere of HBA1-2 and H1 was increased by 0.15 and 0.3 Å respectively, and H2 was set as optional feature to decrease the stringency of the PhMOD. The sensitivity of the PhMOD was measured by the ability to select the HMGB1 inhibitors in the first ranking positions among the TS, whereas specificity was measured by the ability of the PhMOD to identify HMGB1 inhibitors only. Sensitivity and specificity were represented in a receiver operating characteristic (ROC) curve to visualize the performance of the PhMOD (Supplementary Figure 3B). The ROC plot was calculated and visualized using LigandScout 3.02 and the AUC value (area under the ROC curve) was used to evaluate the ROC curve, with values between 0 (lowest) and 1 (highest) performance. The PhMOD was able to select active compounds significantly better than random (AUC = 0.89) and with sufficient sensitivity and specificity for screening procedures (Supplementary Figure 3B). The Screening Perspective in LigandScout was used to filter the hit molecules resulting from the three VS programs against the three-dimensional PhMOD. Only molecules that fulfilled at least three out of four pharmacophore features were retained.

Hits Clustering
The screened subsets were pooled together and clustered using Canvas 3.9 software (Schrödinger, LLC, New York, NY, 2020) (Duan et al., 2010). The chemical fingerprint of each compound was first calculated using a binary 2D linear fingerprint. Hierarchical clustering was then performed using the linear fingerprints with the average linkage method, which considers the average distance between all inter-cluster pairs. We obtained 20 clusters and analyzed only those with a compound population >10 molecules. We obtained six main clusters visualized using Hierarchical Clustering Dendrogram panel in Canvas. The clusters were numbered from the highest to the lowest populated: cluster 1 (148 molecules), cluster 2 (65 molecules), cluster 3 (47 molecules), cluster 4 (40 molecules), cluster 5 (25 molecules), and cluster 6 (20 molecules). The chemical scaffold representative of each cluster was calculated using CanvasMCS, the maximum common substructure facility. We set the fraction of molecules that must match the MCS as at least half the population of each cluster.
Recombinant human CXCL12 (Accession code P48061, and residues M1-K69) was cloned into pET30a vector with restriction enzymes NdeI and BamHI. The resulting plasmid was transformed into E.coli BL21(DE3) and cells were grown at 37 • C in Luria-Bertani medium. CXCL12 expression was induced by the addition of isopropyl-beta-D-thiolgalactopyranoside (IPTG) at a final concentration of 1 mM when cultures reached an optical density of 0.6-0.8 at 600 nm. Induced cultures were grown for an additional 4 h at 37 • C, harvested by centrifugation at 6,000 rpm, and stored at −20 • C until further processing. The cell pellets were resuspended in 20 ml buffer containing 50 mM TrisHCl (pH 8), 100 mM NaCl, 1 mM EDTA, 5 mM DTT, 0.1 mg DNAse, 0.1 mg RNAse, and 5 mg lysozyme. The resuspended cells were lysed by sonication pulsed 1 s on and 1 s off for 2 min at 60% power. The inclusion bodies were washed twice with buffer containing 50 mM TrisHCl (pH 8), 100 mM NaCl, 1 mM EDTA, and 5 mM DTT, and finally with 50 mM TrisHCl (pH 8), 100 mM NaCl, 1 mM EDTA, 5 mM DTT, and 0.1% Triton X100, followed by solubilization in buffer supplemented with 6 M Guanidinium HCl and 50 mM HEPES pH 6.5 overnight at room temperature. The solubilized fraction was cleared by centrifugation at 18,000 rpm at 4 • C for 30 min and diluted dropwise into 250 ml refolding buffer containing 50 mM TrisHCl pH 7.4, 50 mM NaCl, 0.1 mM reduced glutathione, and 0.1 mM oxidized glutathione and kept overnight at 4 • C with stirring. Prior to chromatography, the protein solution was centrifuged at 18,000 rpm for 30 min. CXCL12 was purified by cation-exchange chromatography using a SP Sepharose resin (SP Sepharose HP, GE Healthcare Bioscience AB, Uppsala, Sweden). The protein was washed with buffer A, supplemented with 50 mM TrisHCl (pH 7.4) and 50 mM NaCl, and eluted with buffer B, containing 50 mM TrisHCl pH 7.4 and 1 M NaCl. CXCL12 was finally dialyzed against a buffer containing 20 mM phosphate buffer pH 6, 20 mM NaCl.
Protein concentrations were determined considering molar extinction coefficients at 280 nm of 21,430 and 8,700 M −1 cm −1 for HMGB1 and CXCL12, respectively.
Saturation Transfer Difference and Water-Ligand Observed via Gradient Spectroscopy. STD and waterLOGSY experiments have been performed on 0.5 mM compounds (1a, 1b, 2a, 2b, 4a, 4b, 5, 6) in the presence of 0.05 mM HMGB1 in NMR buffer. STD experiments were acquired using a pulse scheme (Bruker pulse sequence: stddiffesgp.3) with excitation sculpting with gradients for water suppression and spin-lock field to suppress protein signals. The spectra were acquired using 128 scans, a spectral width of 9,600 Hz, and 64 K data points for acquisition. For protein saturation, a train of 60 Gaussian-shaped pulses of 50 ms was applied, for a total saturation time of 3 s. Relaxation delay was set to 3 s. On-resonance irradiations were set at 10 ppm for 1a, 1b, 4b, and 5 and at 0 ppm for 2a, 2b, 4a, and 6. Off-resonance was always set at 107 ppm. STD spectra were obtained by internal subtraction of the on-resonance spectrum from the off-resonance spectrum. WaterLOGSY experiments were acquired using a pulse scheme as described (Dalvit et al., 2000) with excitation sculpting and flip-back for water suppression. The spectra were acquired using 128 scans, 32 K data points for acquisition, and mixing time was set to 1 s.
Ligand-based competition experiments were performed comparing STD and waterLOGSY spectra acquired on a sample containing 0.5 mM diflunisal and 0.05 mM CXCL12 with and without 1 mM 2a.
Intermolecular nuclear Overhauser effect (nOes) between 2a and Box A were obtained from 3D 13 C-NOESY-HSQC with no evolution on 13 C dimension (2,048 × 1 × 256 increments) experiments with 15 N/ 13 C filter in F1 (mixing time 200 ms); protein and ligand concentration were 0.8 and 1.6 mM, respectively, in D 2 O.

MST Experiments
MST experiments were performed at 24 • C on a NanoTemper R Monolith NT.115 instrument with red filters, using 40% LED power and 60% MST power. Binding experiments were carried out using 6His-tagged HMGB1 and 6His-tagged CXCL12, noncovalently labeled with the NT647 fluorescence dye (Bartoschik et al., 2018).
For binding assays, 2a and 2c were titrated (16-points) on 6His-tagged HMGB1 and CXCL12 (MST buffer containing 20 mM phosphate buffer pH 7.3, 20 mM NaCl, 0.05% Tween). The ligand dilutions were generated as a 1:2 dilution of the stock solution using MST buffer; a constant amount of labeled proteins (50 nM) was added to all dilution steps. Maximum concentration of 2a and 2c in the titration series was 5 mM. Complex samples were incubated for 15 min before loading into NanoTemper premium capillaries.
Competition experiments were carried out by pre-forming a complex between labeled 6His-tagged HMGB1 (50 nM) and unlabeled CXCL12 (10 µM, i.e., 2 times the estimated K d ) (Jerabek-Willemsen et al., 2011). For 16-points titration series of 2a, serial 1:2 dilutions of the 2a stock solution were made into MST buffer, and a constant amount of pre-formed heterocomplex was added to all dilution steps. All samples were incubated for 15 min and centrifuged at 15,000 g for 10 min before measurements. Maximum concentration of 2a in the titrations series was 5 mM. Addition of 2a induced the recovery of the MST signal of HMGB1 toward the unbound state value.
For all MST experiments, data points were the average of three measurements (error bars correspond to standard deviation). All data analyses were carried out using NanoTemper analysis software using the K d model fitting for the binding assay and Hill model for competition experiment.

Data Driven Docking Models
Molecular docking of 2a on BoxA (residues G3-Y77) and BoxB (A93-G173) (coordinates obtained as described in 2.1) were performed using the data-driven software HADDOCK 2.2 (Dominguez et al., 2003;van Zundert et al., 2016) following the classical three-stage procedure, which includes: (1) randomization of orientations and rigid body minimization, (2) simulated annealing in torsion angle space, and (3) refinement in Cartesian space with explicit water. Ambiguous interaction restraints (AIRs) were defined as residues with CSP > Avg + sd were used to define active residues, whose solvent accessible surface neighbors were set as passive (Supplementary Table 1). In the case of CXCL12, only the residues located around the diflunisal binding site (De Leo et al., 2019) were set as active Optimized parameters for liquid simulation (OPLS) were used for the protein (protein-allhdg5-4 and protein-allhdg5-4-caro). The geometric coordinates and parameters for 2a were calculated and optimized using the PRODRG server (Schüttelkopf and van Aalten, 2004). Calculations generated 1,000, 1,000, and 500 structures for the rigid body docking (it0), the semiflexible refinement (it1), and the explicit solvent refinement (water), respectively. The final 500 structures obtained after water refinement were scored with HADDOCK (HADDOCKscore = 1.0 E vdW + 0.2 E elec + 1.0 E desolv + 0.1 E AIR ) for a weighted combination of van der Waals (vdW) and electrostatic energy terms (Lennard-Jones and Coulomb potentials), empirical desolvation term (Fernández-Recio et al., 2004), and ambiguous interaction restraint energy term, which reflects the accordance of the model to the input restraints.
HADDOCK models were clustered (Daura et al., 1999) based on their interface root mean square deviation (rmsd), setting the cutoff and the minimum number of models in a cluster to 1.8 Å and 10 for the boxes, and 2.5 Å and 10 for CXCL12, respectively. Proteins were aligned and fitted on the backbone of active residues, reported in Supplementary Table 1. The rmsd of 2a was calculated only on the heavy atoms of the entire scaffold.
To remove any bias of the cluster size on the cluster statistics, the final overall score of each cluster was calculated on the four lowest HADDOCK scores models in that cluster. For each protein the cluster with the best fitting relative to the experimentally-driven restraints (lowest number of violations) and the best HADDOCK score (cluster 1 for BoxA and BoxB, cluster 3 for CXCL12) was selected (Supplementary Figure 4).
The analysis of the docking calculations was performed by applying in-house python and tcl scripts. Molecular images were generated by PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC, and 3D Protein Imagining online server (Tomasello et al., 2020).

Cell Migration Experiments
For fibroblast chemotaxis, modified Boyden chambers were used with filters (pore diameter 8 µm; Neuro Probe) coated with 50 µg/mL fibronectin (Roche). Mouse 3T3 cells (50,000 in 200 µL) were added to the upper chamber. Serum-free DMEM as negative control, HMGB1, and/or other molecules were added to the lower chamber at the indicated concentration, and then cells were left to migrate for 3 h at 37 • C. Cells were fixed with ethanol and stained with Giemsa Stain (Sigma), then nonmigrating cells were removed with a cotton swab. All assays were done at least in biological triplicate. The migrated cells were acquired with Zeiss Imager M.2 microscope at 10x magnification, then evaluated with an automated counting program. All assays were done at least in biological triplicate and were repeated at least twice.

In silico Assessment of HMGB1 Ligandability
We first searched for potential druggable hot spots on HMGB1 surfaces using the computational solvent mapping programs FTMap (Kozakov et al., 2015) and Dogsite (Volkamer et al., 2012). FTMap docks in silico small organic probes onto the surface of the protein and samples a large number of probe conformations, whereby contiguous areas of multiple probes clustering together might be considered as ligandable. FTMap identified seven total consensus sites (CSs) on BoxA. CS2 CS3, CS5, and CS7 were within a distance of 8 Å, defining the highest density hot spot (39 probes), located at the interface of the two helices forming the short arm of the L-shaped fold ( Figure 1A). Nine CSs were identified on BoxB; CS2 was the most populated site, but did not have any other CS within 8 Å, thus could not be considered a hot spot. Conversely, CS1, CS8, and CS9, that were within 8 Å distance, created between the short helices of the BoxB the most populated hot spot (30 probes) ( Figure 1B). The two hotspot regions identified on the two HMG boxes are characterized by small patches of partially solvent exposed hydrophobic residues, including A16, F17, A19, and F37 and F101, F102, F104, and C105 in BoxA and BoxB, respectively. The ligandability of HMGB1 was further assessed by Dogsite, a web-based open-access algorithm that interrogates rigid protein structures for binding hotspots. The same regions identified by FTMap were the ones with the highest probability to be liganded, with a probability score DrugScore > 0.5 (0.58 for Box A and 0.73 for Box B) (Figures 1C,D). Collectively, both FTMap and Dogsite suggest that both HMG boxes are ligandable.

Virtual Screening
We next screened the Zbc Drugs subset (101746 compounds) on BoxA and BoxB using the docking programs Glide, AutoDock Vina and AutoDock 4.2.6, and their results were combined (Figure 2). The three programs rely on distinct empirical scoring functions (Chaput and Mouawad, 2017;Li et al., 2019) and differently weighted electrostatic and hydrophobic interactions. Thus, merging of their outputs was expected to mitigate the possible biases toward hydrophobic or electrostatic contributions characteristic for AutoDock Vina, AutoDock 4.2.6, and Glide, respectively.
The Virtual Screening (VS) workflow implemented in Glide (High Throughput Virtual Screening, HTVS; single precision, SP; extra precision XP docking) retrieved 105 compounds for BoxA and 109 for BoxB, ranked according to the scoring function assigned to the poses. For this subset the ADME descriptors were predicted, and molecules with ADME values outside the 95% range calculated for known drugs were discarded, yielding 100 and 105 hits for BoxA and for BoxB, respectively. Filtering according to ADME descriptors excluded molecules with low octanol/water partition coefficient (logP oct/wat ), low predicted binding affinity to human serum albumin, high dipole moment, high electron affinity, and low brain/blood partition coefficient (logBB), all properties in line with the known bias of Glide toward charged and hydrophilic molecules (Wang et al., 2016). For the VS performed with AutoDock Vina and AutoDock 4.2.6, the first 50,000 top ranked poses were filtered, retaining the ones in which the ligand was at a distance ≤6 Å from R23 in BoxA and R109 in BoxB. These arginines are crucial for ligand binding both for glygcyrrhizin (Mollica et al., 2007) and diflunisal (De Leo et al., 2019). Next, visual inspection of the AutoDock Vina/AutoDock 4.2.6 poses yielded 160/243 molecules for BoxA and 70/220 molecules for BoxB (Figure 2). The molecules were next filtered according to ADME descriptors, retaining 98/110 molecules for BoxA and 68/100 molecules for BoxB (Figure 2). The characteristics of the discarded molecules reflected the bias of both AutoDock Vina/AutoDock 4.2.6 toward the selection of hydrophobic hits. In fact they were predicted to be or to generate metabolites with a high aromatic component of the solvent accessible surface area (SASA), low aqueous solubility (log S), and low values of H-bond acceptor moieties, hydrophilic SASA, dipole moment, and water/gas partition coefficient.
Overall, VS using the three different programs yielded three sets for BoxA and three for BoxB for a total of 581 possible hits. This merged set was further interrogated using a pharmacophore model (PhMOD) generated considering known active HMGB1 inhibitors as described in Materials and Methods. The PhMOD consisted in two H-bonding acceptors (HBA1, HBA2) and two hydrophobic moieties (H1, H2) (Figure 2). Pharmacophore screening of the 581 possible hits resulted in 380 structures fulfilling at least 3 out of the four pharmacophoric requirements. These structures were then subjected to fingerprint calculation and hierarchical clustering and were finally grouped in six highly populated structural clusters (i.e., containing more than 10 molecules each) with distinct chemical scaffolds (Figure 3). Finally, we shortlisted for experimental validation by NMR spectroscopy 21 compounds representative of the six clusters (Table 1), based on visual inspection of docking poses, logP, and vendor availability.

NMR Validation of VS Hits Identifies 5,5 ′ -Methylenedi-2,3-Cresotic Acid (2a) as a Ligand of HMGB1
Of the 21 compounds selected for experimental validation, only eight turned out to be soluble (>0.1 mM) in water or in maximum 10% DMSO ( Table 1). The binding to HMGB1 of the soluble molecules was next tested using ligand-based NMR spectroscopy methods. To this aim we prepared samples containing 50 µM HMGB1 and a 10-fold excess of ligand on which we performed saturation transfer difference (STD) (Mayer and Meyer, 1999) and water-ligand observed via gradient spectroscopy (waterLOGSY) experiments (Dalvit et al., 2001) (Figure 4 and Supplementary Figure 5). Of the compounds FIGURE 2 | VS workflow. Scheme of the filtering steps adopted during the docking procedures computed with AutoDock Vina and AutoDock 4.2.6 on BoxA (Left) and BoxB (Right). The filtering criteria and the number of compounds (cmpd) selected at each filtering step are explicitly indicated. The red and the yellow spheres of the PhMOD represent the H-bonding acceptors (HBA1 and HBA2) and the hydrophobic (H1 and H2) moieties, respectively.
FIGURE 3 | Cluster analysis. Dendrogram (obtained with Canvas 3.9) representing the hierarchical clustering based on chemical 2D fingerprint. 20 clusters were obtained, only those with a compound population with more than 10 molecules were considered, yielding 6 main clusters: cluster 1 (148 molecules, violet), cluster 2 (65 molecules magenta), cluster 3 (47 molecules, yellow), cluster 4 (40 molecules, orange), cluster 5 (25 molecules, pink), and cluster 6 (20 molecules, light blue). For each cluster the chemical scaffold representative of the cluster, the ZINC code of the selected molecules for experimental validation, and their corresponding IDs are reported. tested, only 2a (5,5 ′ -methylenedi-2,3-cresotic acid) ( Figure 4A) and 6 (rosmarinic acid) ( Figure 4B) appeared to interact with HMGB1 as shown by STD effects and inversion of the sign in waterLOGSY experiments, whereby the strongest effects were observed for 2a ( Figure 4A). We further validated the interaction of these two ligands by protein-based NMR experiments (Williamson, 2013) using 1 H-15 N labeled HMGB1. Indeed, upon stepwise addition of 6 or 2a, we observed chemical shift perturbations (CSPs) in the HMGB1 Heteronuclear Single Quantum Coherence (HSQC) spectra (Figures 5A,B). For molecule 2a the interaction occurred on the fast-intermediate exchange regime on the NMR time scale with the disappearance of a few peaks on BoxB (Figure 5A). On the other hand, the CSPs induced by 6 were in a fast exchange regime and were far smaller as compared to 2a, and mainly involved BoxA ( Figure 5B). Linear fitting of the chemical shifts as a function of added ligand indicated an apparent K d of 0.9 mM and >10 mM for 2a and 6, respectively ( Figure 5C). Interestingly, 6, a polyphenolic component of the leaves of Perilla frutescens with general antinflammatory properties (Luo et al., 2020), is known to attenuate inflammatory responses elicited by HMGB1/TLR4/NF-κB signaling (Yang et al., 2013;Lv et al., 2019). This axis relies on the oxidized form of HMGB1, where Cys22 and Cys44 located on BoxA form a disulphide bond. On the other hand, the pattern of CSPs ( Figure 5A) clearly indicated that 2a recognized both HMG boxes in full-length HMGB1. Similarly to both glycyrrhizin (Mollica et al., 2007) and diflunisal (De Leo et al., 2019), 2a bound independently to the two domains, and the pattern of CSPs observed for the isolated HMG boxes was comparable to the one in the full length protein (Supplementary Figures 6A,B). When mapped on the structures of the two HMG boxes, HMGB1 residues with significant CSP upon 2a addition (F17, Q20, R23, E24, and S41 on BoxA, F102, F104, C105, R109, K111, I121, D123, V124, A125, K127, L128, G129, E130, M131, and W132 on BoxB) defined a small surface between the first and the second helix (Figure 5D), in agreement with the computational analysis that identified this area as a putative hotspot (Figure 1). Indeed, this region on both boxes is characterized by a small solvent-exposed hydrophobic surface, well suited for favorable van der Waals (vdW) interactions with the aromatic rings of 2a ( Figure 5D).

3D Model of the Interactions of 2a With HMGB1 Boxes
CSPs data were next used to generate HADDOCK (Dominguez et al., 2003) data-driven docking models of 2a in complex with BoxA and BoxB. 2a accommodates between the two short helices of the HMG boxes, establishing favorable vdW interactions with the hydrophobic side chains of V19 and F37 in BoxA and of F102, V124, and A125 in BoxB (Figures 6A,B). For each box, the models also suggest the presence of stabilizing electrostatic interactions between one carboxylate of 2a and the R23 and R109 sidechains of BoxA and BoxB, respectively (Figures 6A,B). Indeed, in NMR titrations performed on a 15 N-HMGB1 mutant, where both R23 and R109 were substituted by alanines (R23A/R109A), we observed reduced chemical shift displacements, indicative of reduced interaction, thus supporting the involvement of these arginines in binding to 2a (Supplementary Figure 7). The second aromatic ring and the associated carboxylic group of 2a establish π-π and polar interactions with F37 and S41, respectively ( Figure 6A). Interestingly, titration on HMGB1 with 2c (3-Methylsalicylic acid), a 2a analog that lacks the second methylsalicylic ring, showed a lack of binding as assessed by both NMR and MST experiments (Supplementary Figure 8), thus supporting the notion that both salicylic moieties contribute to the binding. Conversely, in BoxB the second salicylate moiety does not seem to be involved in stabilizing interactions ( Figure 6B). Collectively, these 3D models indicate that the steric and electronic features of 2a, which fulfill the predicted pharmacophoric model (Figure 6C), are indeed appropriate to interact with both HMG boxes. In particular, the two major interactions between the ligand and the target consist of a salt bridge between the carboxylate of 2a and the guanidinium groups of the conserved R23 and R109, and hydrophobic interactions between the phenyl groups of 2a and the hydrophobic patch at the interface of the two helices, forming the short arm of the L shaped HMG boxes.
2a Breaks the HMGB1·CXCL12 Heterocomplex NMR titrations and HADDOCK calculations indicate that 2a targets in part the region that has been shown to be involved in the interaction with CXCL12 De Leo et al., 2019;Fassi et al., 2019). We therefore asked whether 2a was able to interfere with the HMGB1·CXCL12 heterocomplex ( Figure 7A). To this end we performed MST experiments, whereby increasing concentrations of 2a were added to a preformed HMGB1·CXCL12 complex. During the titration we observed a sharp transition at 0.31 ± 0.04 mM, consistent with the detachment of fluorescently-labeled HMGB1 from the heterocomplex (Figure 7A). These results were further confirmed by NMR-based Antagonist Induced Dissociation Assay (AIDA) (Krajewski et al., 2007). In this experiment, a 15 N HSQC spectrum was first acquired on free 15 N HMGB1 (Figure 7B), then on a preformed complex composed by 15 N labeled HMGB1 (0.1 mM) and unlabeled   (Figure 7D). We interpreted this strong line broadening as an effect associated with multiple equilibria involving free HMGB1, HMGB1 bound to 2a, and HMGB1 bound to CXCL12. Upon addition of ten-fold excess of 2a we observed the narrowing of 1 H-15 N signals, providing direct evidence of complex dissociation ( Figure 7E). Taken together, both MST and AIDA experiments indicate that 2a causes the dissociation of the HMGB1·CXCL12 heterocomplex.

2a Binds to CXCL12
Recent work on the inhibition of HMGB1 chemotactic activity by diflunisal has revealed a peculiar mechanism of action in which the breakage of HMGB1·CXCL12 occurs through the dual binding of the ligand to both proteins (De Leo et al., 2019). In the same study we also demonstrated that glycyrrhizin, the first chemical probe identified to bind and inhibit HMGB1 chemotactic activity (Mollica et al., 2007), was also a binder of CXCL12 and an inhibitor of the HMGB1·CXCL12 heterocomplex. We thus asked whether 2a was also able to bind CXCL12. Indeed, NMR titrations of 2a on 15 N labeled CXCL12 confirmed a direct interaction ( Figure 8A). Binding occurred in the fast exchange regime, in line with the 200 µM affinity, as measured by MST (Supplementary Figure 9A) and by NMR (Supplementary Figure 9B). The highest CSPs upon ligand binding involved residues located on the β1 strand (V23, H25, K27) and on the so-called CXCR4 sulfotyrosine (sulfoY21) binding site (Veldkamp et al., 2008) (N45, Q48, V49) ( Figure 8A). Indeed, STD and waterLOGSY competition experiments performed by adding 1 mM 2a to a preformed CXCL12:diflunisal complex (1:10) indicated that 2a competes with diflunisal for the same binding site, as assessed by the reduction of both STD and waterLOGSY signals of diflunisal upon 2a addition ( Figure 8B).
As for diflunisal, we hypothesized that CSPs affecting residues of the β1 strand were due to allosteric effects induced by ligand binding. A HADDOCK model of the 2a-CXCL12 complex indicated that 2a, in analogy to diflunisal (De Leo et al., 2019), accommodates in the CXCR4 sulfoY21 binding site, through which the two salicylate moieties of 2a establish polar interactions with R47, N44, and N45 sidechains and with the backbone carbonyl of Q48 and the amide of C50, respectively. Both aromatic rings are involved in hydrophobic interactions with V18, V49, and L42 ( Figure 8C). Accordingly, deletion of one methylsalicylic moiety drastically reduced the interaction of 2c and CXCL12 (Supplementary Figure 8).
Taken together, these data show that 2a is also a ligand of CXCL12 and targets the CXCR4 sulfoY21 binding site. 2a Inhibits the Chemotactic Activity Elicited by HMGB1·CXCL12 Heterocomplex Finally, having assessed the ability of 2a to simultaneously bind HMGB1 and CXCL12 and to break the HMGB1·CXCL12 heterocomplex, we asked whether 2a affected the chemotactic activity of the HMGB1·CXCL12 heterocomplex. Indeed, 2a inhibited the chemotaxis of mouse 3T3 fibroblasts in a dosedependent manner, with an IC 50 close to 10 pM, the lowest documented so far for inhibitors of HMGB1-induced cell migration. Importantly, 2a did not influence the general motility of fibroblasts, as it did not affect chemotaxis toward fMLP (Figure 9), indicating that 2a selectively targets the chemotactic activity of the HMGB1·CXCL12 heterocomplex.

DISCUSSION
Recent studies investigating the functional synergism of chemokine-based heterodimers have shown that hampering heterophilic interactions interrupts inflammation (von Hundelshausen et al., 2017). Thus, the targeting of these protein-protein interactions is emerging as a valuable strategy for the development of new selective antagonists suitable for the tailored modulation of specific inflammatory responses. In this sense, a very promising example is represented by CCL5derived peptides that inhibit the atherogenic CCL5·CCL17 interaction and hamper CXCL12-driven platelet activity (von Hundelshausen et al., 2017). However, the development of PPI antagonists, either peptidomimetics or small organic molecules, poses incredible challenges when the three-dimensional structure of the complex is unknown. This is the case for HMGB1·CXCL12 heterocomplex, a master regulator of the recruitment of inflammatory cells via the CXCR4 axis, and therefore a valuable target for the development of selective anti-inflammatory compounds. High-resolution atomistic descriptions of this interaction are still elusive, thus making structure-based drug discovery studies extremely difficult. One reductionist strategy that can be adopted to cope with this problem is to target the isolated components of the heterocomplex. While CXCL12 as direct interactor of CXCR4 is a popular target for drug development, as shown by several structure-based drug design studies (Veldkamp et al., 2010;Smith et al., 2016), no systematic structure-based small-molecule screening has been reported until now for HMGB1. This is probably due to the structural features of both HMG boxes, which do not present obvious targetable pockets. Still, computational analysis of their structures through the mapping servers FTMap (Kozakov et al., 2015) and Dogsite (Volkamer et al., 2012) performed in this work highlighted a binding hot spot at the interface of helix I and helix II of both HMG boxes. Interestingly, the same region is recognized by several known HMGB1 binders, including DNA (Sánchez-Giraldo et al., 2015), small ligands glycyrrhizin (Mollica et al., 2007(Mollica et al., , 2011 and diflunisal (De Leo et al., 2019), peptides (Sgrignani et al., 2020), or proteins such as p53 (Rowell et al., 2012), speaking in favor of virtual screening approaches targeting this area. Prompted by these results, we screened the Zbc library with three different ligand docking programs that rely on different search algorithms and scoring functions, with the aim to compensate for their individual FIGURE 9 | 2a inhibits HMGB1·CXCL12-induced chemotaxis, Mouse 3T3 fibroblasts were subjected to chemotaxis assays in Boyden chambers, 0.25 nM HMGB1, 1.5 nM of CXCL12 or no chemoattractant was added in the lower chamber, together with the indicated concentrations of 2a. 2a does not inhibit fMLP-induced chemotaxis at the highest concentrations tested for HMGB1·CXCL12 induced chemotaxis. Data points with average ± standard deviation (Avg ± SD; n = 3, each point represents a biological replicate) in a representative experiment. Statistics: one-way ANOVA (P < 0.0001), followed by Dunnett's post-tests. ****P < 0.0001 relative to no 2a addition.
weaknesses (Chaput and Mouawad, 2017;Li et al., 2019). The ligands of the best docking poses emerging from these three docking programs were scrutinized and filtered according to a pharmacophoric model based on known HMGB1 binders, then clustered according to structural similarity, and subsequently validated through ligand-based and protein-based NMR methods. Interestingly, at the end of this pipeline, rosmarinic acid (6) emerged as a weak ligand of HMGB1. Indeed, this molecule has been already reported to be an attenuator of the HMGB1/TLR4/NF-κB-dependent inflammatory response (Yang et al., 2013;Lv et al., 2019), supporting the validity of our virtual screening approach. We then focused on the novel hit stemming from our pipeline, 5,5 ′ -methylenedi-2,3-cresotic acid (2a). 2a was previously reported to have modest DNA demethylase inhibition activity (Kuck et al., 2010) and derivatives thereof have been shown to work as weak tyrosine phosphate inhibitors (Shrestha et al., 2007). We demonstrated that 2a interacts with a high micromolar affinity with the predicted hotspot of the single boxes in a way similar to what has been observed with diflunisal, targeting the surface directly involved in the formation of the heterocomplex with CXCL12 (De Leo et al., 2019). Accordingly, NMR and MST competition experiments demonstrated that 2a is able to disrupt the interaction of HMGB1 with CXCL12, with an EC 50 of 0.3 mM. Moreover, cell migration experiments showed that 2a is able to inhibit the chemotactic activity of HMGB1·CXCL12 heterocomplex with an IC 50 of 10 pM, the most potent inhibition activity documented until now, and two orders of magnitude smaller than that of diflunisal. Intriguingly, NMR CSPs and MST data indicate that 2a also binds to CXCL12, as does diflunisal. It targets with high micromolar affinity the CXCR4 sulfotyrosine binding pocket and induces CSPs on CXCL12 β1 strand, likely ascribable to allosteric effects (Ziarek et al., 2013), contributing to the destabilization of the interaction with HMGB1. The dual binding ability to HMGB1 and CXCL12 seems to be a common theme in inhibitors of HMGB1·CXCL12 chemotactic activity, as observed for glycyrrhizin and diflunisal (De Leo et al., 2019). Notably all these molecules, which are active on different inflammation-related proteins, can be thus considered as multitarget drugs, whose promiscuity and deviation from the "one target, one drug" dogma might translate into improved efficacy through the simultaneous engagement of different targets involved in the same pathological process (Bolognesi and Cavalli, 2016). In a polypharmacology perspective (Anighoro et al., 2014), the previously documented weak inhibition activity of 2a toward DNMT1 and tyrosine phosphate, two enzymes involved in detrimental inflammation pathways (Través et al., 2014;Al-Kharashi et al., 2018), might synergistically contribute to the efficacy of 2a as anti-inflammatory modulator. Overall, results emerging from this and previous studies strongly suggest that salicylate derivatives are well suited for the interaction with HMGB1, with the carboxylic and hydroxylic group of the aromatic rings fulfilling the predicted pharmacophoric characteristics required for HMGB1 ligands. Intriguingly, these pharmacophoric features seem in part to also satisfy the chemical requirements for CXCL12 ligands. In fact, previous structurebased design studies aiming at identifying CXCL12 inhibitors have suggested that a carboxylic group, preferentially combined to an aromatic group, is a good mimic of the negatively charged CXCR4 sulfotyrosine (Veldkamp et al., 2010;Smith et al., 2016). Thus, we speculate that some known CXCL12 inhibitors might also be ligands of HMGB1, conceivably exerting a synergic inhibition effect on the disruption of the HMGB1·CXCL12 heterocomplex and as direct inhibitors of the interaction of CXCL12 with CXCR4.
The double targeting action of 2a, diflunisal, and glycyrrhizin only in part reconciles the discrepancy that exists between their high-micromolar ability to disrupt the HMGB1·CXCL12 complex and their nanomolar or picomolar chemotactic inhibition efficacy. Chemotaxis relies on many different signaling cascades at the cell membrane, where the local concentration of both HMGB1 and CXCL12 might increase through their direct interaction with cell surface glycosaminoglycans (GAGs) (Merenmies et al., 1991;Rueda et al., 2012). Additional mechanisms and cooperative binding phenomena among multiple actors at the cell surface might occur (Whitty, 2008) synergistically and/or allosterically, contributing to complex destabilization. We posit that the inhibition activity of these ligands might rely on multistep dynamic processes that, besides the interaction with both CXCL12 and HMGB1, also include the direct binding to a high affinity binding pocket at the heterocomplex interface. Once targeted, this site might cooperate with heterocomplex dissociations and/or induce allosteric changes that affect heterocomplex binding to CXCR4 and the related signaling. As long as the threedimensional structure of the heterocomplex is unknown and the molecular and functional details of its interaction with CXCR4 remain underexplored, the details of the mechanisms dictating the inhibitory potency of 2a and of other ligands of the HMGB1/CXCL12/CXCR4 axis will be unresolved. Nevertheless, our reductionist approach showed that small solvent-exposed clefts of HMGB1 boxes, traditionally viewed as a barrier to the success of in silico ligand screening, can actually serve as legitimate sites for drug discovery. Even more interestingly, cytokine heterocomplexes that appeared undruggable are actually promising drug targets.

CONCLUSIONS
In summary, the combination of VS and experimental validation proved successful in identifying 2a as the best molecule documented until now in terms of inhibition of in vitro chemotaxis activity of the HMGB1·CXCL12 heterocomplex. These results pave the way for future structure activity relationship studies for the optimization of the pharmacological targeting of this heterophilic interaction for anti-inflammatory purposes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FDL performed the NMR and MST experiments, computational studies, analyzed the data, and prepared the manuscript. GQ performed the NMR experiments, while MVM expressed and purified recombinant proteins. FDM performed the cell migration experiments. MEB and GM directed the study, were involved in all aspects of the experimental design, data analysis, and manuscript preparation. All authors critically reviewed the text and figures.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2020.598710/full#supplementary-material Training set migration experiments of decoy molecules building and validation of the Pharmacophore model (PhMOD), HADDOCK restraints, HADDOCK clusters, hits binding to HMGB1 by ligand-based NMR methods, Binding of 2a to BoxA and BoxB, comparison 2a binding to HMGB1 wt and R23A/R109A, MST and NMR binding curve of 2a to HMGB1.