The size matters? A computational tool to design bivalent ligands.

Motivation: Bivalent ligands are increasingly important such as for targeting G protein-coupled receptor (GPCR) dimers or proteolysis targeting chimeras (PROTACs). They contain two pharmacophoric units that simultaneously bind in their corresponding binding sites, connected with a spacer chain. Here, we report a molecular modelling tool that links the pharmacophore units via the shortest pathway along the receptors van der Waals surface and then scores the solutions providing prioritization for the design of new bivalent ligands.
Results: Bivalent ligands of known dimers of GPCRs, PROTACs and a model bivalent antibody/antigen system were analysed. The tool could rapidly assess the preferred linker length for the different systems and recapitulated the best reported results. In the case of GPCR dimers the results suggest that in some cases these ligands might bind to a secondary binding site at the extracellular entrance (vestibule or allosteric site) instead of the orthosteric binding site. Availability and implementation: Freely accessible from the Molecular Operating Environment svl exchange server (https://svl.chemcomp.com/). Supplementary information: Supplementary data are available at Bioinformatics online.

1 Introduction

Bivalent ligands have emerged to interrogate dimer function or chemically induce biological proximity (Stanton ). A bivalent ligand is a single chemical entity composed of two covalently linked, by a spacer chain, pharmacophores with potential to interact simultaneously with two protomers (Lane ; Valant ). Homobivalent ligands contain the same two pharmacophore units, whereas heterobivalent ligands link two different pharmacophores. Bivalent ligands are proving valuable in multiple research areas such as in the field of G protein-coupled receptors (GPCRs) (Hiller ; Shonberg ), proteolysis targeting chimeras (PROTACs) (Corson ), or bivalent antibody/antigen interactions (Mack ), among others.

GPCRs are classically described as monomeric transmembrane (TM) receptors; however, they can also form higher-order complexes of the same (homo) or different (hetero) protomers (Ferre ). The study and understanding of these complexes prompt a new dimension for GPCR research (Maurice ). The spacer length of GPCR bivalent ligands is a key factor and depends on the dimer interface, the structure of the pharmacophores, and the position and direction of the attachment points (Hiller ; Shonberg ). Many studies have reported spacers between 15 and 25 atoms such as ligands targeting dimers of the δ- and μ-opioid (Fig. 1, molecule 1) (Daniels ), serotonin (2) (Soulier ), dopamine (3 and 4) (Kuhhorn ; McRobb ), cannabinoid and μ opioid/cannabinoid (5 and 6) (Le Naour ; Zhang ), A2A-D2 (8) (Soriano ) and A1-β2 (9) (Barlow ) (homo/hetero)dimers to name just a few. Recently, highly active ligands with spacers as large as 44 and 66 atoms (7) (Hubner ) were reported.

Fig. 1.

Examples of selected bivalent ligands from the literature. The reported most active molecules for μ-δ (1, MDAN-21), 5HT4-5HT4 (2), D2-D2 (3 and 4), CB1-CB1 (5), μ-CB1 (6), D2-NTS1 (7), A2A-D2 (8), A1-β2 (9), BRD4-VHL (10) and CA-CA (11) dimers

A bivalent PROTAC ligand also covalently links two pharmacophores with a spacer chain (Corson ; Stanton ). The aim of these ligands is to selectively induce degradation of a target protein, by linking a high affinity binder of the target protein to a strong binder of an E3 ubiquitin ligase. This brings the target and ligase in proximity, triggering polyubiquitination and proteasome-dependent degradation of the target. This post-translational approach offers many benefits to silence proteins in a cell. The very first crystal structure of a bivalent PROTAC ligand (Fig. 1, 10) bound at the bromodomain target protein BRD4 and E3 ligase von Hippel-Lindau protein (VHL) was recently reported (Gadd ). The discovery of this ligand involved the exploration of different length polyethylene glycol (PEG) spacers with compound 10 being preferred (Zengerle ).

A further application of bivalent ligands has involved the study of bivalent antibody/antigen interactions. The group of Whitesides used a synthetic dimer of carbonic anhydrase (CA) as a model system to replicate the binding of bivalent antibodies to antigens (Mack ). A series of bivalent benzenesulfonamides connected by different oligosarcosine spacers were synthesized. The preferred ligand contains four sarcosine units (Fig. 1, 11).

Here we present a MOE-based molecular modelling tool that can define and prioritize the preferred spacer length for simultaneous binding at two protomers within an (homo/hetero)dimer. The tool links the pharmacophore units via the shortest pathway along the proteins’ van der Waals surface, identifying extracellular grooves and plausible conformations for a flexible spacer. We have studied with this tool the spacer size of previously reported bivalent ligands of known dimers of CA, GPCRs and BRD4-VHL. This computational tool can speed-up the design and synthesis of future bivalent ligands.

2 Materials and methods

2.1 Computational models of (homo/hetero)dimers

Crystal structures of the BRD4-VHL heterodimer (PDB id 5T35) (Gadd ) and CA homodimer (3PJJ) (Mack ) were used, as well as computational models of GPCR (homo/hetero)dimers built based on experimentally reported TM1/2, TM4/5 and TM5/6 interfaces (Cordomi ). Crystal structures of GPCR protomers were available for all the cases except dopamine D2 receptor, which was modelled based on the D3 receptor crystal structure. Further details are in supporting information.

2.2 Molecular dynamics simulations

The μ-δ opioid GPCR heterodimer was studied with several bivalent ligands. After system construction and equilibration, several replicas of unrestrained MD trajectories were produced for 1 µs. All simulations used GROMACS v5.01 (Pronk ), with force fields described by Cordomi et). Details for MD simulations are available in supporting information.

2.3 A computational tool to design the size of the spacer

The tool was developed in the Molecular Operating Environment (MOE) software (Chemical computing group Inc., Montreal, QC, Canada). An interaction surface is calculated from a 0-potential contour of the van der Waals potential of a probe atom on a lattice around the receptor atoms. The OPLSA-AA van der Waals parameters for a TIP3P water oxygen were used. The vertices for the triangles making up the interaction surface are used as nodes in a graph, connected where within 1.05 Å. On each of the pharmacophore units an attachment point (a specified atom) to the linker is defined, and the closest vertices on the interaction surface are selected. The shortest path between these points is identified through the graph of the vertices (ignoring duplicate vertices) (Fig. 2). A set of pre-built spacer molecules is read as an input, each with two defined attachment points used to make the connection to the corresponding pharmacophore units. The spacer atoms along the shortest path between the attachment points are identified and transformed to positions on the interaction surface using linear interpolation between the interaction surface vertices. Spacers were ignored if the path across the interaction surface was too long when compared with the sum of the van der Waals radii for the atoms from the shortest path between the linker attachment points. Spacers are defined without hydrogens, allowing the bond to be formed with the pharmacophore attachment point prior to energy minimization A stepped energy minimization protocol is used. Initially, only the spacer atoms are free, while all atoms are included in the later stages of the energy minimization. The strength of the van der Waals and electrostatic repulsion terms are initially scaled by a factor of 10−5; a factor that was increased by 10-fold in each iteration until they are applied at normal strengths. The total force field energy (Etot) and the interaction energy between the spacer atoms and the rest of the system are reported for each of the spacers.

Fig. 2.

Graphical representation of how the tool generates the shortest path through the graph of surface vertices (blue spheres) and uses this to define the final minimized ‘shortest path’ spacer (green atoms). Multiple close surface vertices to each connection atom point are used as start and end points

3 Results

3.1 Design of PROTAC ligands

The crystal structure of the BRD4-VHL heterodimer in complex with the PROTAC bivalent ligand 10 MZ1 is known (Gadd ). MZ1 links the BRD4 inhibitor JQ1 to the VHL ligand VH032 via a three-unit PEG spacer. We used this example to validate the accuracy of the computational tool in predicting the conformation of the spacer as the structure of MZ1 is known. Clearly, the experimental and predicted conformations of the spacer are in very good agreement (Fig. 3). As the PROTAC field expands, more structures will reveal the range of distances between pharmacophore binding sites, and this tool can assist design of bivalent ligands (Chessum ; Crew ).

Fig. 3.

Computational assessment of spacer length for BRD4-VHL bivalent ligand. The tool was used to test a database 1- to 7-unit (n) PEG spacers. The energy of the system (E(tot), kcal/mol) shows the best results use two and three PEG repeats. Both correspond to experimentally valid linkers. The computationally derived three-unit PEG structure from the tool (green) superposes excellently to the X-ray structure (grey)

3.2 Design of CA bivalent ligands

Mack ) constructed a synthetic dimer of CA, and a series of bivalent ligands with different lengths, as a model system for the binding of bivalent antibodies to antigens. They obtained the structure of the CA homodimer (although without bivalent ligand) and synthesized a series of benzenesulfonamides connected by oligosarcosine spacers with n = 0, 2, 4, 6, 8 and 10 sarcosine units (LSarnL). The total energy of the systems is reported in Supplementary Figure S1. In agreement with the experimental results (Mack ), our MOE-based tool predicts the LSar4L bivalent ligand (11) as the most favorable.

3.3 Design of bivalent ligands of GPCR (homo/hetero)dimers

Despite enormous effort devoted to their design, rational development has been difficult due to the large variety of GPCR members with different extracellular domains, different ligand-binding modes, and the lack of knowledge concerning their dimeric forms.

3.3.1 Factors influencing the spacer size of bivalent ligands

When considering the design of bivalent ligands, the first step understands the structure of the protein dimer complex. In the case of GPCRs, computational models are usually built (Navarro , 2018a, b; Vinals ) from crystal structures that contain a variety of symmetric dimerization interfaces (summarized in Cordomi ). This mainly includes the TM1/2, TM4/5, TM5 or TM5/6 interfaces (Supplementary Fig. S2). Notably, the TM5/6 dimerization interface reproduces the shortest distance of 27 Å between orthosteric sites, the TM1/2 interface the longer distance of 42 Å, whereas the distances of the TM5 or TM4/5 interfaces are in between and similar to each other (Supplementary Fig. S2). Thus, the optimal length of the spacer of the bivalent ligand depends on the TM interface to target.

The second factor is the influence of ligand orientation. Orthosteric ligand binding typically occurs in a cavity between the extracellular segments of TMs 3, 5–7 or in a minor binding cavity located between the extracellular segments of TMs 1–3, 7 (Rosenkilde ). Ligands can penetrate to different depths within these pockets (Venkatakrishnan ) and leave atoms exposed available for linking at variable locations as shown in Supplementary Figure S3 using crystal structures of the adenosine A2A receptor. Hence, the variation in orientation and buriedness within the same receptor influences spacer length for bivalent ligands.

Clearly protein surfaces are non-uniform and approximating spacer lengths based only on inter-binding site distances is unsatisfactory. This is especially relevant in GPCRs where the extracellular domain is highly variable in sequence, length, and structure (Gonzalez ). For instance, in GPCRs for lipid mediators the extracellular N-terminus and extracellular loop (ECL) 2-folds over the ligand-binding pocket (e.g. CB1R in Supplementary Fig. S4). In contrast, ECL2 of peptide receptors, formed by two β-strands, fully exposes the binding site to the extracellular environment (e.g. opioid receptors in Supplementary Fig. S4). On the other hand, a helical segment forms ECL2 of adrenergic receptors (e.g. β2-adrenergic in Supplementary Fig. S4). Again, this is a further consideration that influences the spacer length of bivalent ligands and prompts the need for a structure based approach.

3.3.2 Design of the optimal spacer size for the μOR–δOR heterodimer

To demonstrate the tool in the more complex GPCR systems, we first study the μ and δ opioid receptor (μ–δ) heterodimer because (i) crystal structures of μ-OR (Manglik ) and δ-OR (Granier ) are known, (ii) the structure of μOR revealed putative dimer interfaces (Manglik ) and (iii) seminal work from the Portoghese group reported a series of μ-δ heterodimer bivalent ligands (Daniels ; Yekkirala ). The structure of μ-OR shows receptor protomers associated into dimers through two different interfaces: via TM1/2 or TM5/6. The bivalent ligand MDAN-21 (1) for μ-δ contains the μ-OR agonist oxymorphone (OXY) and the δ-OR antagonist naltrindole (NTI) (Daniels ) (Fig. 1). As the name indicates, MDAN-21 contains 21-atoms between pharmacophores, including spacer (7 atoms) and linker (2 × 7 atoms). Regarding nomenclature (Shonberg ), linker refers to the chemical group that attaches the spacer chain to each pharmacophoric ligand, and hence the spacer is the chain between. Figure 4A shows the μ-δ dimer interacting via the TM5/6 interface (the TM1/2 interface was not considered because the longer distance between binding pockets cannot accommodate ligands such as MDAN-21) bound to MDAN-21. A database of spacers ranging from 2 heavy atom units (16-atom spacer and linker combined) to 14 units (28-atom spacer and linker) was constructed. The tool constructs the linked bivalent ligands and assesses their fit across the dimer surface. Figure 4B shows the best solution obtained from the linking script after energy minimization, which contains 12 heavy-atom units (26-atom spacer and linker, we name it BL-26).

Fig. 4.

A computational tool to design the spacer size of bivalent ligands. (A) Computational model of the μ–δ heterodimer bound to MDAN-21, the two pharmacophore units (OXY and NTI, respectively) are shown in magenta spheres, carbon spacer is grey. (B) The preferred solution from the computational tool showing the spacer optimally placed across the protein van der Waals surface and connected to each ligand

3.3.3 Computational validation of the calculated spacer size

The tool generates energetically favorable spacers that can be used directly as recommendations for synthesis or for further computational study. In the latter regard we used unbiased molecular dynamics (MDs) simulations in the μs time-scale to evaluate the stability of MDAN-21, the designed BL-26, and an intermediate sized BL-23 ligand at the μ-δ heterodimer interacting via the TM5/6 interface. The stability of the ligands was monitored by the root mean-square deviation (rmsd) of the heavy atoms along the MD trajectories (Supplementary Fig. S5A), whereas the binding of the pharmacophoric units to the orthosteric binding site was monitored by the salt bridge distance between the protonated amine of NTI or OXY and Asp3.32 of δ-OR (Supplementary Fig. S5B) or μ-OR (Supplementary Fig. S5C), respectively. Clearly, the shortest MDAN-21 and BL-23 ligands cannot fulfil the simultaneous binding of the pharmacophoric units to the orthosteric binding site as revealed by both rmsd values and salt bridge distances. In contrast, the pharmacophoric units of the designed BL-26 ligand remain stable at the orthosteric binding cavities through the unbiased 1 μs MD simulation.

3.3.4 The extracellular vestibule

During the MD simulation of MDAN-21 and BL-23, the NTI pharmacophoric element remained bound to the orthosteric site of δ-OR, whereas the OXY pharmacophore unit left the orthosteric site of μ-OR However, before reaching solvent the OXY pharmacophore ligand was retained by electrostatic (Glu2315.35 and Lys2355.39 in TM5) and van der Waals (Trp2285.32, Trp3207.35) interactions in the extracellular domain of μ-OR, as shown by the plateau of the rmsd values (Supplementary Fig. S5A). The presence of a small cavity at the entrance of the orthosteric binding site [named extracellular (Dror ) or membrane (Stanley ) vestibule or secondary (Gonzalez ) or metastable (Fronik ) binding site] has been described. To further evaluate the proposed binding of MDAN-21 (NTI in the orthosteric binding site of δ-OR and OXY in the extracellular vestibule or secondary binding site of μ-OR), we performed two additional MD replicas starting at this conformation. The stability of this binding mode was monitored and confirmed by the low rmsd values of ligand heavy atoms (Supplementary Fig. S6). Clearly, MDAN-21 can simultaneously bind the orthosteric site of δ-OR and the vestibule site of μ-OR, providing a plausible alternative explanation for MDAN-21’s interesting pharmacology (Portoghese ).

3.3.5 Calculation of the optimal spacer size in reported GPCR bivalent ligands

Having shown how our MOE-based tool can help define preferred spacer lengths of bivalent ligands for the μ-δ heterodimer, we have further applied the tool to the D2–D2 and CB1–CB1 homodimers and μ–CB1, D2–NTS1 and A1–β2 heterodimers (Table 1) for which bivalent ligands are also reported in the literature (Fig. 1). As described above, the TM1/2, TM4/5 and TM5/6 interfaces were considered, in each case permitting different possible linker and spacer lengths. Images of the preferred calculated result for each case are provided in supporting information Supplementary Table S1. The differences between experimentally reported linker lengths versus the corresponding calculated values are immediately apparent. A fundamental result of our systematic approach is the proposition that many reported ligands are too short to link orthosteric sites of adjacent monomers via a receptor-surface bound spacer. Most experimental reports show spacers between 15 and 25 atoms, but even in the shortest possible path such as in the TM5/6 interface for opioid receptors, with little extracellular hindrance (Supplementary Fig. S4), the spacer still requires 26 heavy atoms. Looking in more detail, several trends can be observed. TM1/2 dimers typically prefer a heavy atom linker plus spacer of around 36–40 heavy atoms, the TM4/5 interface is in the range of 50–55, while the TM5/6 interface presents the shortest number of atoms. However, there are exceptions. The TM1/2 interface of the CB1 and D2 homodimers prefer fewer heavy atoms than the other interfaces, because the linker is directed towards TM1 (these differences can be seen by inspecting the images in Supplementary Table S1). Hence, as we have discussed above, specific details of each system impact the spacer length, meaning it is recommended to model each new case in order to optimally apply the approach.

Table 1.

Comparing the preferred experimentally reported linker plus spacer size versus the preferred calculated size for the GPCR (homo/hetero)dimer formed via the TM1/2, TM4/5 or TM5/6 interfaces

Experimental results from literature			Calculated results, preferred number of heavy atoms for each respective dimera
Systemsb	Bivalent ligandc	Number of heavy atomsd	TM1/2	TM4/5	TM5/6
μ–δ	1	21	40	42	26
D2–D2	4	16	36	56	47
CB1–CB1	5	15	33	54	53
μ–CB1	6	20	41	55	42
D2–NTS1	7	44 (66)e	36	52	27
A1–β2	9	25	40	55	33

Note: The comparison is performed in terms of number of heavy atoms.

Models were built according to experimental details, the script was applied in each case to determine the preferred linker size.

For references see Section 1 of main text.

Ligand number refers to Figure 1, and represent the most active ligand identified in each associated experimental study.

Based on the preferred ligand.

This study reported that molecule 7 (m = 2, Fig. 1) and another example with m = 3 (66 atom spacer) were both preferred bivalent ligands.

4 Discussion

We have presented a computational tool to assist the design of the optimal spacer and linker size connecting two pharmacophore units bound in a dimer. It provides a scoring based on energy minimization of the resulting ligand receptor complex. Our results suggest inconsistencies between claimed bivalent ligands of GPCRs and what seems plausible assuming an extracellular surface binding mode for the spacer, as has also been pointed by others (Glass ; Lane ; Shonberg ). The alternative is for the spacer to bind between TM helices and disrupt interactions of several helical turns, which would incur a high energetic cost. Thus, these ligands with short spacers may adopt alternative binding modes (schematized in Fig. 5). Some ligands with very short spacers might bind the orthosteric site as well as a secondary binding site at the extracellular entrance (vestibule or allosteric site) of the same protomer (Fig. 5A). However, this mode of binding to a single protomer can also allosterically modulate the binding of a ligand to the orthosteric site of a second protomer (Lane ). Second, ligands with larger spacers might bind the orthosteric site of Protomer A and the secondary binding site at the extracellular entrance of Protomer B (Fig. 5B). It has been suggested that the conformations of orthosteric and secondary binding sites of the same protomer are coupled such as that the presence of a ligand in one affects the shape of the other (Bokoch ; Dror ). And finally, actual bivalent ligands, with the optimal linker size, can simultaneously bind the orthosteric binding sites of Protomers A and B (Fig. 5C).

Fig. 5.

Binding modes of bivalent ligands to GPCR (homo/hetero)dimers. (A) Ligands can bind the orthosteric site as well as a secondary binding site at the extracellular entrance of the same protomer. (B) Ligands with larger spacers might bind the orthosteric site of protomer A and the secondary binding site at the extracellular entrance of Protomer B. (C) Bivalent ligands, with the optimal spacer size, can simultaneously bind the orthosteric binding sites of Protomers A and B

Theoretically, the affinity of a bivalent ligand could approach the sum of the binding free energies for the individual components along with entropic corrections (Foreman, 2017; Mammen ). With linkers that are sufficiently long, impressive potency gains have been seen (Mack ) but various aspects contribute to a less than expected binding affinity. Negative cooperativity between two protomers can result in the binding of a ligand to the first protomer decreasing the affinity of the ligand for the second protomer (Ferre ). Secondly, as mentioned above, the spacer may not permit the two pharmacophores to bind simultaneously in their respective sites. Meanwhile, entropy, is expected to be beneficial because the second pharmacophore unit is present at a much higher local concentration once tethered in a close radius above the dimer (Numata ). Nevertheless, binding of the spacer and linker will have an entropic cost once fixed on the receptor surface and recent work to address this has prompted the study of less flexible spacers (Tanaka ). Kramer and Karpen (1998) and Mack have shown that flexible spacers much longer than the distance between sites still bind tightly. In many cases however, GPCR bivalent ligands lack the potency gain and our results suggest they are too short to span the high affinity binding sites. Therefore, does the size matter? It becomes clear that a minimum size is required and evaluating this using the presented tool is expected to be a pragmatic and valuable approach. Future reports will show the prospective use of this tool for the design of highly potent bivalent ligands of GPCR homodimers.

Click here for additional data file.