{'fraction_non_alphanumeric': 0.03128061422094639, 'fraction_numerical': 0.034333492812065707, 'mean_word_length': 5.150203166596609, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery of 5,5’-methylenedi-2,3-cresotic acid as a potent inhibitor of 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 to the recruitment of inflammatory cells forming a heterocomplex with the chemokine CXCL12 (HMGB1•CXCL12), herewith 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, 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 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 responses1. Chemokines and their receptors organize leukocyte trafficking and migration to the tissues both in healthy and pathological conditions2. Chemokines bind to and activate seven-transmembrane G-protein-coupled receptors (GPCRs), a receptor superfamily involved in many different diseases3. The different chemokines and their receptor are object of intense drug-discovery studies4. Increasing evidences show that chemokines usually require interactions with additional players that modulate the inflammation signaling cascade triggered by their receptors5-7; hence, the identification of molecules that disrupt the interactions of chemokines with their modulators might offer new and selective pharmacological opportunities against inflammation related diseases8. 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 receptor9. 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 conditions10-12. The HMGB1•CXCL12 complex has been shown to trigger specific CXCR4 homodimer rearrangements9, ERK activation and calcium fluxes, along with enhanced CXCR4-dependent monocyte migration and tissue regeneration13,14. HMGB1 is a highly conserved, abundant non-histone 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 cytokine-like proinflammatory and chemoattractant effects 13,15. 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 tail16,17. Its modular organization, flexibility, different oxidation states13,18 along with 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 activity19 to the modulation of innate immunity20. 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 tissues9 and exacerbates the immune response in pathological conditions5,21. Thus, the disruption of the HMGB1•CXCL12 heterocomplex with small molecules might offer new and selective strategies against inflammation related diseases22, 23. Interference with this heterophilic interaction is attractive but challenging for several reasons: i. the three-dimensional structure of the complex is still unknown, ii. the interaction surface is expected to be large and dynamic, thus difficult to be targeted by small molecules 9,24,25. Previous work has shown that it is possible to interfere with the pro-inflammatory properties of HMGB1 using small molecules including glycyrrhizin26, salicylic acid (SA) and its derivative amorfrutin27. 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 CXCL1224. 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 object of structure-based drug design studies and virtual screening campaigns aiming at inhibiting its direct activation effect on CXCR428,29, 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 L-shaped boxes of HMGB1, which do not appear to have clearly druggable pockets30,31, 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 structurebased 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,3-cresotic acid (2a) as a potent inhibitor of the chemotactic activity of HMGB1•CXCL12, with an IC50 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 instructions32. 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 PyMol33. 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 protocol34. Pocket Size and DrugScores (the closest the score to 1, the higher the ligandability) were extracted for all identified sites and annotated to pocket numbers.

Structure Based Virtual Screening Library preparation
The Zbc Drugs subset (http://zinc12.docking.org/browse/subsets/special) containing 101746 biogenic drug-like compounds selected according to the 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, minimizing and optimizing the molecules with OPLS_2005 force field. The ligands generated were ready for docking with GLIDE; for docking with AutoDock Vina35 and AutoDock 4.2.636 they were converted to pdbqt-formatted file by in-house scripts. Ligands used to generate the training set for the pharmacophoric model (section 2.3.5) were prepared with the same procedure.

VS Docking studies
Protein-based VS studies were performed with three different programs: Glide (Schrödinger), AutoDock Vina35 and AutoDock 4.2.636. Ligands were docked onto BoxA and BoxB. 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, 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 of R23, the residue at the center of the experimentally validated binding pocket 24. 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 of 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 pdbqt-formatted 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, 30.70 on x, y, z axes, respectively. BoxB grid size was 40 Å in each direction, grid point spacing 0.531 Å and was centered on 20.5, 12.5, 3.0 on x, y, z axes, respectively. 2500000 max number of energy evaluations, 27000 max number of generations, 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 LigandScout337 from Inte:Ligand. To generate the pharmacophore model we took advantage of a training set (TS) (Supplementary Figure S1A-B) composed by 7 known active HMGB1-interacting molecules from the literature, i.e. glycyrrhizin26, carbenoxolone38, salicylic acid (SA)27, amorfrutin27, acetyl-3-aminoethyl-SA Ac3AESA27, diflunisal24, inflachromene39 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, prostaglandin E2) (Supplementary Figure S2). We superimposed the negatively charged moiety of the 3D-structures of the 7 active members of the TS and the following 4 pharmacophore features were generated: 5 hydrophobic moieties, 3 aromatic ring, 4 H-bonding acceptor and 1 H-bonding donor. The model was next refined maintaining only those features that previous structural studies on glycyrrhizin26 and diflunisal24 have indicated as relevant (Supplementary Figure S3A). The resulting PhMOD consisted of two H-bonding acceptor (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 among the TS the HMGB1 inhibitors in the first ranking positions, 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 S3B). 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  S3B). 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)40. The chemical fingerprint of each compound was first calculated using 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 6 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.

Protein expression and purifications
Recombinant HMGB1 constructs, including BoxA (residues 1-89), BoxB (residues 90-175), wildtype protein (Accession code P63158, residues 1-214), and mutant HMGB1 (R23A/R109A) in labelled and unlabeled forms were produced as described24. After expression and purification, proteins were dialyzed against NMR buffer, containing 20 mM phosphate buffer pH 7.3, 150 mM NaCl, 1 mM DTT. 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 additional 4 hours at 37°C, harvested by centrifugation at 6000 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 second on and 1 second off for 2 minutes at 60% power. The inclusion bodies were washed twice with buffer containing 50 mM TrisHCl (pH 8), 100 mM NaCl, 1 mM EDTA, 5 mM DTT and finally with 50 mM TrisHCl (pH 8), 100 mM NaCl, 1 mM EDTA, 5 mM DTT, 0.1% Triton X100, followed by solubilization in buffer supplemented with 6M Guanidinium HCl and 50 mM HEPES pH 6.5 overnight at room temperature. The solubilized fraction was cleared by centrifugation at 18000 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, 0.1 mM oxidized glutathione and kept overnight at 4°C with stirring. Prior to chromatography, the protein solution was centrifuged at 18000 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 21430 and 8700 M-1 cm-1 for HMGB1 and CXCL12, respectively. Proteins used for cell-based assays were provided by HMGBiotech (Milan).
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 9600 Hz, 64K 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, 5 and at 0 ppm for 2a, 2b, 4a, 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 described43 with excitation sculpting and flip-back for water suppression. The spectra were acquired using 128 scans, 32K data points for acquisition, 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.
Titrations: For NMR titrations, at each titration point a 2D water-flip-back 1H-15N-edited HSQC spectrum was acquired with 2048 (160) complex points for 1H (15N), respectively, apodised by 90° shifted squared (sine) window functions and zero filled to 256 points for indirect dimension. Titrations of 2a and 6 into 15N HMGB1 were performed adding 0.5, 1, 2, 3, 5, 10 equivalents of ligand into the labelled protein. In order to minimize dilution and NMR signal loss, titrations were carried out by adding small aliquots of concentrated ligands (10 mM in 20 mM phosphate buffer, pH 7.3, 150 mM NaCl) to the 15N

MST experiments
MST experiments were performed at 24°C on a NanoTemper® 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, non-covalently labelled with the NT647 fluorescence dye 47. For binding assays 2a was titrated (16-points) on 6His-tagged 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 labelled proteins (50 nM) was added to all dilution steps. Maximum concentration of 2a in the titration series was 5 mM. Complex samples were incubated for 15 minutes before loading into NanoTemper premium capillaries. Competition experiments were carried out pre-forming a complex between labelled 6His-tagged HMGB1 (50 nM) and unlabeled CXCL12 (10 μM, i.e. 2 times the estimated Kd)48. 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 minutes and centrifuged at 15,000 g for 10 minutes 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 towards 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 Kd 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), BoxB (A93-G173) (coordinates obtained as described in 2.1) were performed using the data-driven software HADDOCK 2.249,50 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 follows: residues with CSP > Avg + sd were used to define active residues, whose solvent accessible surface neighbors were set as passive (Supporting Table S1). In the case of CXCL12, only the residues located around the diflunisal binding site24 were set as active (Supporting Table S1), as STD competition experiments of 2a in the presence of diflunisal demonstrated that they both compete for the same site. In the case of BoxA, intermolecular nOes were included as unambiguous restraints in the calculations only in the semi-flexible refinement stage, setting the maximum distance of the nOe H pairs to 5 Å (Supporting Table S1). 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 server51. Calculations generated 1000, 1000, 500 structures for the rigid body docking (it0), the semi-flexible refinement (it1) and the explicit solvent refinement (water), respectively. The final 500 structures obtained after water refinement were scored with HADDOCK (HADDOCKscore = 1.0 EvdW + 0.2 Eelec + 1.0 Edesolv + 0.1 EAIR) for a weighted combination of van der Waals (vdW) and electrostatic energy terms (Lennard-Jones and Coulomb potentials), empirical desolvation term52 and ambiguous interaction restraint energy term, which reflects the accordance of the model to the input restraints. HADDOCK models were clustered53 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 Supplementary Table S1. 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 S4A-C). The analysis of the docking calculations was performed 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 server54.

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 hours at 37°C. Cells were fixed with ethanol and stained with Giemsa Stain (Sigma), then non-migrating 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 FTMap32 and Dogsite34. 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 7 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 highest 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, F37 and F101, F102, F104, 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 (DS = 0.58 for Box A and 0.73 for Box B) ( Figure 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 three different docking programs, i.e. Glide, AutoDock Vina and AutoDock 4.2.6, and their results were combined (Figure 2). The three programs rely on distinct empirical scoring functions55,56 and differently weighted electrostatic and hydrophobic interactions. Thus, merging of their outputs was expected to mitigate the possible biases towards 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 (logPoct/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 towards charged and hydrophilic molecules57. For the VS performed with AutoDock Vina and AutoDock 4.2.6 the first 50000 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 glygcyrrhizin26 and diflunisal24. 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 towards 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, 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 4 pharmacophoric requirements. These structures were then subjected to fingerprint calculation and hierarchical clustering and were finally grouped in 6 highly populated structural clusters (i.e. containing more than 10 molecules each) with distinct chemical scaffolds ( Figure 3). Finally, out of these clusters, based on visual inspection of docking poses, logP and vendor availability, we shortlisted for experimental validation by NMR spectroscopy 21 compounds representative of the 6 clusters (Table 1).

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 8 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)58 and water-ligand observed via gradient spectroscopy (waterLOGSY) experiments59 (Figure 4 and Supplementary Figure S5). Of the compounds 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 experiments60 using 1H-15N labelled HMGB1. Indeed, upon stepwise addition of 6 or 2a we observed chemical shift perturbations (CSPs) in the HMGB1 Heteronuclear Single Quantum Coherence (HSQC) spectra ( Figure 5A,B).  Figure S6A-B). When mapped on the structures of the two HMG boxes, HMGB1 residues with significant CSP upon 2a addition (F17, Q20, R23, E24, 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 HADDOCK49 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 ( Figure 6A,B). For each box, the models suggest also the presence of stabilizing electrostatic interactions between one carboxylate of 2a and the R23 and R109 sidechains of BoxA and BoxB, respectively ( Figure 6A,B). Indeed, in NMR titrations performed on a 15N-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 S7). The second aromatic ring and the associated carboxylic group of 2a establish π-π and polar interactions with F37 and S41, respectively ( Figure 6A). 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 in i. a salt-bridge between the carboxylate of 2a and the guanidinium groups of the conserved R23 and R109, and ii. 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 CXCL1224,25,64. 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 labelled HMGB1 from the heterocomplex ( Figure 7A). These results were further confirmed by NMR-based Antagonist Induced Dissociation Assay (AIDA) 45. In this experiment a 15N HSQC spectrum was first acquired on free 15N HMGB1 ( Figure 7B), then on a preformed complex composed by 15N

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 interaction occurs through the dual binding of the ligand to both proteins24. In the same study we also demonstrated that glycyrrhizin, the first chemical probe identified to bind and inhibit HMGB1 chemotactic activity26, was 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 15N labelled 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 S8A) and by NMR (Supplementary Figure S8B). 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 site65 (N45, Q48, V49) ( Figure 8A). Indeed STD and waterLOGSY competition experiments performed 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 diflunisal24, accommodates in the CXCR4 sulfoY21 binding site, whereby 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, L42 ( Figure 8C). 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 dose-dependent manner, with an IC50 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 targets selectively 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 inflammation8. Thus, the targeting of these PPIs 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 CCL5-derived peptides that inhibited the atherogenic CCL5•CCL17 interaction and hampered CXCL12-driven platelet activity8. 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 antiinflammatory compounds. However, 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 designs studies29,66, no systematic structure-based small-molecule screening has been reported until now for HMGB1. This is probably due to the structural features of both the HMG boxes, that do not present clearly targetable pockets. Still, computational analysis of their structures through the mapping servers FTMap32 and Dogsite34 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 DNA67, small ligands glycyrrhizin26,38 and diflunisal24), peptides68 or proteins such as p5369, 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 weaknesses55,56. The ligands of the best docking poses emerging from these three docking programs were i. scrutinized and filtered according to a pharmacophoric model based on known HMGB1 binders, ii. clustered according to structural similarity and iii. 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 62,63, 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 activity70. We demonstrated that 2a interacts with high micromolar affinity with the predicted hotspot of the single boxes in a way similar to what observed with diflunisal, targeting the surface directly involved in the formation of the heterocomplex with CXCL1224. Accordingly, NMR and MST competition experiments demonstrated that 2a is able to disrupt the interaction of HMGB1 with CXCL12, with an EC50 of 0.3 mM. Moreover, cell migration experiments showed that 2a is able to inhibit the chemotactic activity of HMGB1•CXCL12 heterocomplex, with an IC50 of 10 pM, the most potent inhibition activity documented until now, and 2 orders of magnitude larger than that of diflunisal. Intriguingly, NMR CSP and MST data indicate that 2a binds also 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 effects71 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 diflunisal24. 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 satisfy also the chemical requirements for CXCL12 ligands. In fact, previous structure-based 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 sulfotyrosine29,66. Thus, we speculate that some known CXCL12 inhibitors might be also ligands of HMGB1, conceivably exerting a synergic inhibition effect on the disruption of the HMGB1•CXCL12 heterocomplex and as direct inhibitors the interaction of CXCL12 with CXCR4.
The double targeting action of 2a, diflunisal and glycyrrhizin only in part reconciles the discrepancy existing 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)72,73. Additional mechanisms and cooperative binding phenomena among multiple actors at the cell surface might occur74 herewith 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 to heterocomplex dissociations and/or induce allosteric changes that affect heterocomplex binding to CXCR4 and the related signaling. As long as the three-dimensional 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.