Core motif substructure searches were performed using SMILES filters applied to the MolPort library and the FKKTlib. The core motifs of cyclobutenedione (with an additional nitrogen atom) and thiourea were defined by SMILES expressions O=C1C=C(N)C1=O and NC(N)=S, respectively. First, MolPort library filtering was performed using the RDKit substructure filter node (RDKit: Open-source cheminformatics, 2021) as implemented in the KNIME analytics platform (Berthold et al., 2007). 2,452 cyclobutenediones were extracted from the MolPort database and docked to the prepared CCR7 receptor using Glide XP (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020) (Friesner et al., 2006) as described below. Of the 100 highest scoring compounds, 16 diverse compounds were selected for purchase. Second, the MolPort database was searched for the thiourea core motif, which yielded more than 90k available compounds. Duplicates, PAINS (Dahlin et al., 2015), and compounds with reactive functional groups were removed to yield 63k compounds. Docking with Glide SP (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020) was performed as described below, and after clustering, from the top 500 hits, 13 diverse compounds were purchased. Second, the FKKTlib was filtered and all 9 and 13 compounds available in solid form with cyclobutenedione and thiourea core motifs, respectively, were experimentally evaluated.

Ligand-based virtual screening (LBVS) of the ZINC database of purchasable compounds using LiSiCA software (Lešnik et al., 2015) was performed with the bioactive 3D conformation of cmp2105 as the reference compound (PDB ID: 6QZH, ligand JLW) (Jaeger et al., 2019). The double bonds of the thiadiazole-dioxide of the JLW ligand were correctly assigned, since they are missing in the PDB structure. Both 2D and 3D options of the LiSiCA were used with all other settings set to default values. From the 200 compounds most similar to the reference cmp2105 according to the Tanimoto score, 27 diverse compounds were purchased—12 of them arising from the 2D method and 15 were discovered with the 3D method.

The MolPort library was screened using ROCS 3.3.2.2 software (OpenEye Scientific Software, Inc, Santa Fe, NM, United States; www.eyesopen.com) (Hawkins et al., 2007). For model A, a 3D pose for the first query, navarixin, was obtained by docking with Glide XP to the prepared CCR7 receptor, as described below. The bioactive 3D conformation of cmp2105 (PDB ID: 6QZH, ligand JLW) (Jaeger et al., 2019) was used to create models B and C. All models were validated with the set of active compounds and generated decoys. The default settings of ROCS were used for virtual screening of all three queries. Virtual hits were prioritized based on the ComboScore, which considers similarity of 3D shape (“ShapeTanimoto”) and chemical pattern (“ColorScore”). For each query, 27–30 top scoring compounds were purchased from the clustered list of top 100 scoring hits.

Before the release of the CCR7 X-ray crystal, a homology model of CCR7 was built using the structure of human CCR9 (PDB ID: 5LWE; B chain). The template was identified by running 10 PSI Blast iterations on the starting CCR7 (UniProt ID: P32248) sequence to identify five top scoring templates (PDB IDs: 5LWE, 5UIW, 4YAY, 5WB2 and 5UNF) (Müller et al., 1999). The alignment and template was used to build the homology model using YASARA Twinset software (Krieger et al., 2002; Krieger and Vriend, 2015) using the following parameters: speed: Slow, EValue Max: 0.5, Templates Total: 5, Templates SameSeq: OligoState: 4, Alignments: 5, LoopSamples: 50 and TermExtension: 10.18 models were built and each model subjected to an unrestrained energy minimization with explicit water molecules by simulated annealing employing the YASARA2 force field (Krieger et al., 2009). The models were rated according to a quality Z-score and the best scoring model was used. The latter contained 276 of 378 target residues (73.0%) aligned to template residues. The sequence identity was 46.0% and the sequence similarity 68.1% (BLOSUM62 > 0). The monomer homology model after full unrestrained simulated annealing minimization was rated as optimal by YASARA with internal quality Z-score of 0.110, comprised amino acids 47-352, and was further checked with WHAT-IF test set.

The CCR7 receptor was prepared from the X-ray crystal structure (PDB ID: 6QZH) (Jaeger et al., 2019) using Protein Preparation Wizard (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020) (Madhavi Sastry et al., 2013). Briefly, missing side chains and missing loop 255-261 were modelled with Prime (Jacobson et al., 2004), hydrogen atoms were added, residues were protonated at pH 7.0, the hydrogen bonding network was refined, waters beyond 3.0 Å from other heteroatoms were removed, and restrained minimization was performed. The double bonds of the thiadiazole-dioxide core in the co-crystalized ligand cmp2105 (PDB ID: 6QZH, ligand JLW) were assigned (Jaeger et al., 2019). Only the allosteric binding site was considered for pharmacophore-based screening and molecular docking.

The prepared CCR7 receptor was used to generate a structure-based pharmacophore model using LigandScout 4.4 (Inte:Ligand GmbH) (Wolber and Langer, 2005). Exclusion volumes defining regions based on the shape of the binding site residues were generated, and all features were converted to vectors. One hydrogen bond donor and one hydrophobic feature were marked as optional. This model was validated with the set of active compounds and generated decoys. Default settings in LigandScout were used. Virtual screening of the diversity set of compounds available from trusted commercial vendors yielded 78 virtual screening hits, which were then visually inspected, clustered according to Morgan fingerprints, and 23 diverse compounds were purchased.

Molecular docking was performed sequentially with FRED 3.4.0.2 (OpenEye Scientific Software, Inc, Santa Fe, NM, United States; www.eyesopen.com) and Glide software (Schrödinger Suite 2020-2: Glide, Schrödinger, LLC, New York, NY, 2020). First, Make Receptor 3.4.0.2 (OpenEye Scientific Software, Inc, Santa Fe, NM, United States; www.eyesopen.com) was used to define grid box of the allosteric binding site of the prepared receptor. The volume of the box was 6,725 Å³ (17.75 Å × 21.15 Å × 17.92 Å) and size of the outer contour was reduced to 1,139 Å³. Re-docking of the co-crystalized ligand cmp2105 using FRED and Glide SP resulted in a root-mean-square deviation (RMSD) of 0.77 Å and 0.34 Å, respectively, confirming the validity of the pose prediction during docking. In the same manner, docking with Glide XP (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020) (Friesner et al., 2006) was used to obtain the bioactive 3D conformation of navarixin. The MolPort library was docked with FRED and 100,000 highest scoring hits were used for sequential docking with Glide. A 3D structure of one stereoisomer was generated using LigPrep (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020). The prepared receptor’s grid box was centered on the co-crystallized ligand and docking was performed using Glide SP (Schrödinger Suite 2020-2, Schrödinger, LLC, New York, NY, 2020) (Friesner et al., 2004). The 100 highest scoring virtual hits were clustered according to Morgan fingerprints, and 46 diverse compounds were purchased.

Molecular docking with ProBiS-Dock algorithm (Konc et al., 2022) was performed with the prepared CCR7 receptor. Similar receptors with allosterically bound ligands that were known at the time of screening, i.e., CCR9 (PDB ID: 5LWE, ligand 79K [vercirnon]) (Oswald et al., 2016) and CCR2 (PDB ID: 5T1A, ligand VT5) (Zheng et al., 2016), were aligned to the prepared CCR7 receptor. All three ligands, namely cmp2105, vercirnon and VT5, were extracted and used as template ligands that are required for molecular docking with ProBiS-Dock. Re-docking of the co-crystalized ligand cmp2105 with an RMSD of 0.77 Å was performed for validation. From the MolPort library, one million compounds were randomly selected and used for virtual screening. The top 450 virtual hits were clustered according to Morgan fingerprints, and 40 diverse compounds were purchased.

For the LBVS coupled to SBVS approach, we examined similar protein complexes published in the PDB database using ProBiS server (https://probis.nih.gov/) to identify all possible binding sites and small-molecule binding modes on the CCR7. Therefore, the built CCR7 homology model was used as an input for ProBiS calculation and one binding site identified (binding site one in ProBiS; proximity of ligand vercirnon from CCR9 PDB ID: 5LWE) (Konc and Janežič, 2012). The postulated binding site comprised of two pockets was defined by residues: Thr93, Leu147, Ile150, Val264, Ile265, Val268 and Val79, Val80, Thr82, Tyr83, Phe86, Asp94, Thr95, Leu97, Leu98, Leu100, Asp110, Asp336. Model was later validated with an all-atom RMSD of 1.942 Å towards PDB 6QZH crystal with allosteric site correctly identified relative to crystal JLW ligand. With the binding site defined, the receptor structure was generated using OEDocking 3.2.0.2 software package (OpenEye Scientific Software, Inc, Santa Fe, NM, United States; www.eyesopen.com) with MakeReceptor. A box with the volume of 5,805 Å³ (18,33 × 19,00 × 16,67 Å) was defined around reference ligand vercirnon. A balanced site shape potential was calculated where docking volume was 587 Å³. No constraints were used. Docking of the chemokine receptor targeted library (539,814 compounds) to the prepared receptor was performed using FRED from OpenEye as described above with dock_resolution parameter set to High. Top 100 scoring compounds according to Chemgauss4 score were collected, clustered according to Morgan fingerprints (20 clusters, average RMSD linkage) and best scoring representatives selected for purchase and testing.

Before we initiated an extensive in silico-based screening program, several compounds described in the literature were either resynthesized (see Supplementary material) or purchased from commercial vendors to confirm their CCR7 antagonism in various pharmacological assays. Compounds previously shown to be CCR7 antagonists (i.e., cmp2105, navarixin, CS-2, and CS-3), as well as inactive control compounds (maraviroc and vercirnon), were included in this study (Figure 1; Table 1). A competition binding assay was established based on the specific interaction of fluorescently labeled CCL19 with CCR7 overexpressed on whole living cells (Figure 2A). Using CCL19^AF647 as a tracer, the binding affinity for cmp2105, navarixin, CS-2 and CS-3 was evaluated (Figure 2B; Table 1). Whereas dose-dependent inhibition of CCL19^AF647 binding was confirmed for cmp2105 and navarixin, CS-2 and CS-3 did not show any CCR7 affinity. The previously observed stabilizing effect of CS-2 and CS-3 in thermal stability experiments was much smaller than for cmp2105 and navarixin (Jaeger et al., 2019), suggesting a very low binding affinity for CCR7. It should further be noted that for cmp2105 an IC₅₀ value of 35 nM was previously reported when this compound was assessed in membrane-based competition experiments using radioactively labeled CCL19 (Jaeger et al., 2019). The fact that in our assay whole cells are used instead of membrane preparation, which requires the compound to first enter the cell before reaching its intracellular binding pocket, may therefore partly explain the increased apparent IC₅₀ value observed here.

Reference compounds were also evaluated in a CCR7 kinetic, fluorescence-based calcium mobilization assay. Cmp2105 and navarixin showed IC₅₀ values in the 5–15 µM range for antagonizing the CCL19-induced calcium response (Figure 3; Table 1) in line with their CCR7 antagonistic activity previously determined in a β-arrestin recruitment assay (Jaeger et al., 2019). In agreement with the lack of observed binding affinity, CS-2 and CS-3 were also inactive in this CCR7 calcium mobilization assay (Table 1). Furthermore, the absence of activity of vercirnon and maraviroc in the calcium mobilization assay is in agreement with their lack of activity in the thermal shift assay (Jaeger et al., 2019).

To expand the current set of potent CCR7 modulators, we launched a virtual screening campaign. Based on known CCR7 ligands, in particular cmp2105 and navarixin, an LBVS was performed. Furthermore, a recently published crystal structure of the receptor (Jaeger et al., 2019) was used for structure-based virtual screening (SBVS). In addition to libraries of commercially available compounds, we also used the FKKTlib academic library for screening (Figure 4).

LBVS is commonly used in drug discovery and is based on the assumption that structurally similar compounds have similar biological properties. Various metrics are used to express similarity between compounds (Maggiora et al., 2014). We started with a simple substructure search for core motifs that are typical for small molecule ligands targeting the intracellular binding sites of various chemokine receptors. Thioureas bind to an intracellular binding pocket of CXCR2 (Nicholls et al., 2008) and cyclobutenediones have been shown to bind intracellularly to CXCR2 (Liu et al., 2021, 2) and CCR7 (Jaeger et al., 2019). This approach was applied to both the FKKTlib academic library and a library of commercially available compounds.

Second, we used ligand-based virtual screening software LiSiCA (Lešnik et al., 2015) to find compounds with different scaffolds and core motifs than those of the reference compound cmp2105. LiSiCA is based on a graph-theoretical representation of molecules and uses a fast maximum clique algorithm (Konc and Janežič, 2007) to search for 2D or 3D similarities between a reference compound and a database of target compounds. The similarities found are expressed by the Tanimoto coefficients. A library of commercially available compounds was compared to the reference based on both 2D and 3D molecular representations.

Third, 3D shape-based virtual screening was performed by rapid overlay of chemical structures (ROCS) (Hawkins et al., 2007). This method is based on the concept that compounds have a similar shape if their volumes, described by a Gaussian function, overlap well. In addition to molecular volume, a color force field is used to describe other molecular features, such as hydrogen bond donors and acceptors, anions, cations, hydrophobes and rings (Kirchmair et al., 2009). As a starting point for the ROCS search, we modeled a 3D conformation of navarixin by docking with Glide XP and used it to generate ROCS model A (Figure 5A). The bioactive conformation of cmp2105 was extracted from the co-crystal structure (PDB ID: 6QZH) and used directly to create models B and C, which differed in the selection of color features. Only relevant hydrogen bond donors and acceptors based on the distances in the crystal structure were used for the model B (Figure 5B). For model C, only color features in the inner part of the binding pocket were selected, leaving more degrees of freedom for the part of the molecule that extends toward the solvent (Figure 5C). All three models performed well in screening a set of active compounds and generated decoys. The results are presented in the form of receiver-operating characteristic (ROC) curves (Figure 5). Subsequently, the models were used to screen a library of commercially available compounds.

A structure-based pharmacophore model was constructed from the crystal structure of cmp2105 (Figure 6A). The model consisted of a hydrogen bond acceptor for the sulfonyl moiety, two hydrogen bond donors for the secondary amines (one labeled as optional), a hydrogen bond donor and acceptor for the phenol moiety, four hydrophobic features, and exclusion spheres. A hydrophobic feature for the methyl moiety on furan ring was also labeled as optional. The results of screening a set of active compounds and generated decoys were visualized by ROC plot, with the rate of active compounds on the y-axis and the rate of decoys on the x-axis (Figure 6B). Three of eight active compounds were detected by this model, which was in turn used to screen a library of commercially available compounds.

In the next SBVS approach, molecular docking into the allosteric binding site of CCR7 was employed (Figure 7A). First, we used FRED software (McGann, 2012), which was capable of high-throughput docking of a prepared library containing more than 3.5 million commercially available compounds. Then, only 100,000 highest-scoring compounds were used for subsequent docking with Glide software (Friesner et al., 2004), which was expected to be more successful in enriching a virtual hit list but is also more computationally intensive (Kellenberger et al., 2004; McGaughey et al., 2007). Besides, Glide performed better in a re-docking experiment with cmp2105, achieving an RMSD of 0.34 Å, compared to FRED with an RMSD of 0.77 Å.

A SBVS approach was also explored using the ProBiS-Dock algorithm (Konc et al., 2022). Allosteric binding sites of other chemokine receptor crystal structures available at the time of our study were aligned and compared. Accordingly, three template ligands were selected: cmp2105 (CCR7), vercirnon (CCR9), and VT5 (CCR2). The template ligands were used together with the CCR7 crystal structure (PDB ID: 6QZH) (Jaeger et al., 2019) as input to the ProBiS-Dock algorithm (Figure 7B). When docking the library of commercially available compounds, both the docked compound and the receptor were treated as fully flexible to account for the induced fit of ligand binding. The obtained poses were scored using a combination of a site-specific and a generalized statistical scoring function. A site-specific scoring function scores the docked compounds based on their overlap with the template ligands, while a generalized statistical scoring function scores the compounds based on their interactions with the receptor.

Finally, a coupled virtual screening approach was explored, in which a chemokine receptor targeted compound library containing 539,814 compounds was docked using FRED software. This library covered similar compounds to the literature actives on all chemokine receptors. The allosteric binding site of a CCR7 homology model, 5 Å around the vercirnon ligand from the template structure of CCR9, was used for docking.

In total, 287 virtual screening hits were selected from the various approaches and experimentally evaluated as potential CCR7 antagonists in the calcium mobilization assay (Supplementary Excel File). One compound (C040, Figure 8A) showed dose-dependent inhibition of the CCR7-mediated calcium response, affording an EC₅₀ value of 13.14 µM (Figure 3A). Given the conserved nature of the intracellular binding pocket targeted by cmp2105 and navarixin (Jaeger et al., 2019), the inhibiting effect of C040, alongside the reference compounds cmp2105 and navarixin, on the intracellular calcium mobilization mediated by several other chemokine receptors was evaluated (Figure 3B–D). Cmp2105 and navarixin inhibited the CCR7 and CXCR2 mediated calcium mobilization, in line with literature data (Gonsiorek et al., 2007; Jaeger et al., 2019), but had no (or only very limited) effect on CXCR4 and CCR5 mediated responses. In contrast, C040 completely lacked receptor specificity as it inhibited the calcium mobilization downstream all tested chemokine receptors, with similar potencies (IC₅₀ values in the 3–18 µM range). These data suggested that the chemokine receptor antagonistic activity of C040 may, at least partially, be due to interference with the fluorescent assay readout. To further explore this hypothesis, the ability of C040 to inhibit calcium responses not mediated by human chemokine receptors was investigated. It is known that stimulation of CHO-K1 cells with adenosine triphosphate (ATP) leads to a rapid release of Ca²⁺ from intracellular stores (Iredale and Hill, 1993). Using CHO-K1 cells, exactly the same experimental set-up as for the chemokine receptor expressing U87 cells was applied, essentially including the same fluorescent calcium dye (Fluo-2) for cell loading. Also in this experimental setup C040 was able to dose-dependently inhibit the measured calcium response induced by ATP (10 µM final concentration) (Figure 8B), confirming its interference with this particular fluorescent readout. Furthermore, when C040 was assessed in the CCR7 competition binding assay described above, it was inactive at the highest concentration tested (25 µM). Altogether, these data indicate that C040, despite showing activity in the CCR7 calcium assay, should not be selected as a hit compound, for a medicinal chemistry-based optimization campaign.