Interaction of Radiopharmaceuticals with Somatostatin Receptor 2 Revealed by Molecular Dynamics Simulations

The development of drugs targeting somatostatin receptor 2 (SSTR2), generally overexpressed in neuroendocrine tumors, is focus of intense research. A few molecules in conjugation with radionuclides are in clinical use for both diagnostic and therapeutic purposes. These radiopharmaceuticals are composed of a somatostatin analogue biovector conjugated to a chelator moiety bearing the radionuclide. To date, despite valuable efforts, a detailed molecular-level description of the interaction of radiopharmaceuticals in complex with SSTR2 has not yet been accomplished. Therefore, in this work, we carefully analyzed the key dynamical features and detailed molecular interactions of SSTR2 in complex with six radiopharmaceutical compounds selected among the few already in use (64Cu/68Ga-DOTATATE, 68Ga-DOTATOC, 64Cu-SARTATE) and some in clinical development (68Ga-DOTANOC, 64Cu-TETATATE). Through molecular dynamics simulations and exploiting recently available structures of SSTR2, we explored the influence of the different portions of the compounds (peptide, radionuclide, and chelator) in the interaction with the receptor. We identified the most stable binding modes and found distinct interaction patterns characterizing the six compounds. We thus unveiled detailed molecular interactions crucial for the recognition of this class of radiopharmaceuticals. The microscopically well-founded analysis presented in this study provides guidelines for the design of new potent ligands targeting SSTR2.


■ INTRODUCTION
Most efforts of modern medicine are addressed toward personalized medicine, in which each patient is treated according to the molecular features of the disease of interest. 1,2 In this contest, radiopharmaceuticals have been extensively used to specifically target unhealthy tissues. 3,4 According to the decay properties of the radionuclide, compounds can be employed for diagnostic or therapeutic purposes, or both (theranostics). 5 Radionuclides emitting γ or β+ (e.g., 111 In and 68 Ga) are exploited for imaging with single-photon emission computed tomography (SPECT) and positron emission tomography (PET), respectively, while those emitting βor α (e.g., 177 Lu and 211 At) are used for therapeutic treatments. 6 In this last case, after the binding of a radiopharmaceutical to the given target and its subsequent internalization, a cytotoxic dose of radiation is delivered to the cancer cell. 7 In some cases, the radionuclide emits both β+ and β-( 64 Cu), or γ and β-( 177 Lu), making their use suitable for theranostics. 4,8 However, the inability to precisely quantify the radiation doses supplied to tumors and normal tissues has been one of the main drawbacks of radionuclide-based treatments. For example, 111 In decays by electron capture emitting relatively high-energy γ photons with a half-life (t 1/2 ) of 67.2 h, resulting in suboptimal imaging resolution and high radiation exposure in patients, which is even more pronounced when using shortlived isotopes, such as 68 Ga (t 1/2 = 1.13 h). Therefore, to alleviate this problem, it is possible to use a longer-lived radionuclide, allowing a more accurate assessment of biodistribution and tissue clearance. An example of alternative diagnostic agent is represented by the positron-emitting isotope 64 Cu (t 1/2 = 12.7 h, β+ = 17.4%, E max β+ = 653 keV). 9,10 Both 68 Ga and 64 Cu are widely used in peptide receptor radionuclide therapy, a prominent example of which is represented by the treatment of neuroendocrine tumors (NET), 11 where cancer cells are detected by exploiting a high concentration of somatostatin receptors on their surface. 12 These receptors are class A G-protein-coupled receptors (GPCRs) and include the five distinct isoforms SSTR1−5. 13 The isoform 2 (SSTR2), belonging to the SRIF1 sub-class together with the isoforms 3 and 5, 14 is the most expressed in these types of tumors 15,16 and, as a result, several drugs have been developed to specifically bind this receptor. 17 To date, eight peptide-based radiopharmaceutical compounds targeting somatostatin receptors have been approved by FDA, and are routinely used in clinics for different applications (Table 1). 6 Radiopharmaceuticals targeting SSTR2 share a similar three components structure made by (1) a biovector mimicking the structure of the endogenous ligand somatostatin, that is conjugated with (2) a chelator moiety carrying (3) a radionuclide ( Figure 1). 19 In the last year, different structures of SSTR2 in multiple conformational states have been published. 15,20,21,13,22,23 Noteworthy, most of these structures are in complex with agonist ligands, which is often somatostatin or its analogous, like the octa-peptide octreotide. Both experimental and computational studies explored the conformational features of SSTR2 that are common to those of other class A GPCRs, 24 and the key elements characterizing the binding with different types of ligands (see Figure S1 for an overview of the three-dimensional structure of SSTR2 and its main domains). 25 The availability of such structural data can burst the development of new SSTR2 ligands able to bind the receptor with high affinity. 26 However, a detailed molecular and atomistic-level description of the interaction of the radiopharmaceutical/SSTR2 complex is missing, thus hampering the rational design of new effective ligands of this family. Therefore, exploiting the available structural knowledge, in this work, we carefully analyzed the key dynamical features and detailed interactions of SSTR2 in complex with six radiopharmaceuticals. We focused on compounds loaded with either 68 Ga or 64 Cu: the former is the leading β+ emitting radiometal for PET imaging and is contained in two approved drugs ( 68 Ga-DOTATATE and 68 Ga-DOTATOC, Table 1), the latter is used for theranostic purposes and is contained in two approved drugs as well ( 64 Cu-SARTATE and 64 Cu-DOTATATE, Table 1). The two radionuclides were simulated in complex with three different chelators: 1,4,7,10-tetraazacyclododecane-N,N′,N″,N‴-tetraacetic acid (DOTA), 1,4,8,11-tetraazacyclotetradecane-1,4,8,11-tetraacetic acid (TETA), and 3,6,10,13,16,19hexazabicyclo[6.6.6]icosane (SAR) (Figure 1). For the peptide portion, we considered three derivatives of octreotide, namely, TOC, TATE, and NOC: the first one is octreotide with the replacement of F3 into Y3, the second differs from TOC at the last residue (threonine T8 instead of threoninol T-ol8), and the third is octreotide with the replacement of F3 with naphthalene (Nal3) (Figure 1).
The choice of the radiopharmaceuticals was driven by the aim of exploring the influence of the different portions of the ligands in the interaction with the receptor by (1) considering the same radionuclide-chelator ( 68 Ga-DOTA) and changing the peptide (TOC, TATE, and NOC), (2) considering the same chelator-peptide (DOTA-TATE) and changing the radionuclide ( 68 Ga and 64 Cu), and (3) considering the same radionuclide-peptide ( 64 Cu-TATE) and changing the chelator (DOTA, TETA, SAR). Through multicopy μs-long molecular  dynamics (MD) simulations based on a previous investigation on SSTR2 in different states, 25 here we found analogies and differences in the interaction patterns characterizing the binding of the six compounds with SSTR2, and we discovered how each moiety can influence the dynamical behavior of the complexes. The detailed molecular-level analysis presented in this study, thoroughly mapping the SSTR2/ligand interactions, revealed previously unknown structural and mechanistic insights into molecular recognition processes of radiopharmaceuticals at SSTR2.

■ RESULTS AND DISCUSSION
We performed multicopy all-atom MD simulations of six metal-based radiopharmaceutical compounds in complex with SSTR2 (total simulation time of 15 μs per system). We focused on the influence that each component exerts on the dynamic properties of the complexes and the resulting interaction pattern. In the following, we analyze the role of the peptide moiety, the radionuclide, and the chelator by changing only one component at a time and comparing the MD results in terms of dynamics and detailed molecular interactions. Following this strategy, the role of the different components in the interaction could be evaluated more accurately. The Ballesteros-Weinstein numbering scheme for class A GPCRs is adopted throughout the paper. 27 For better clarity, SSTR2 and ligand residues are indicated using the three-and one-letter nomenclature, respectively. Generally, according to root-mean-square deviation (RMSD) of the ligand heavy atoms with respect to the initial frame of the production run, all compounds were highly stable inside the binding pocket, following the order 64 Cu- Figure S2). As expected, most fluctuations were found at the terminal portions of the ligands (i.e., the chelator moiety and the last residue of peptide T8 or T-ol8, Figure S3). Overall, for all compounds, we found the known conserved interactions involving residues located in the bottom part of the binding pocket of SSTR2 (i.e., Asp122 3.32 , Gln126 3.36 ) and the D W4 and K5 motif of the ligands ( Figure  S4), as well as other residues already reported in previous works. 13,20,25 In the following, we focus only on the comparative analysis of protein−ligand interactions characterizing the selected radiopharmaceuticals under investigation.
Small Changes in the Peptide Structure Strongly Affect the Dynamics of the Complex. In the MD simulations of 68 Ga-DOTATOC/TATE/NOC, we did not change the chelator-radionuclide portion ( 68 Ga-DOTA), but only the peptide biovector (TOC, TATE, NOC). As a result, we were able to focus on the influence that small variations in the peptide structure (T8, T-ol8, Y3, or Nal3) exert on the interaction with SSTR2. In all cases, cluster analysis of MD trajectories (see the Computational Details section and Table  S1) reveals a dominating cluster (population in the range of 60−80%) that does not differ significantly from the other two (RMSD in the range of 0.9−4.7 Å), confirming the overall stability of the binding modes (Table S1). Inspecting how the population of the dominant binding mode changes with time in all replicas, we found that it is the most populated one along the whole μs-long time-scale simulation or starting from a few hundreds of ns (see Figure S5). The representatives of the most populated clusters for the three cases are shown in Figure  2. Interestingly, the peptide portion overlaps neatly with the cryo-EM conformation of octreotide in complex with SSTR2 (RMSD octreotide vs TOC/TATE/NOC portions: 1.3/1.4/ 1.5 Å). 20 Looking at the structures, we found in all cases that the 68 Ga-DOTA moiety was placed between TM5 and TM6, with  (Table S1). The extracellular loop ECL2 is colored in red; the gallium ion is represented as a pink sphere; the receptor is represented as a transparent white cartoon; the chelator DOTA is colored in dark cyan; the peptide portions TOC, TATE, and NOC are colored in magenta, violet, and pink, respectively. The top panels report the superimposition of the representative structural clusters with octreotide (green transparent sticks), taken from the PDB 7T11. In the bottom panels, the main interactions are shown as black dotted lines. some of the contacts involving ECL2 and ECL3 as well, while T8/T-ol8 interacted only with ECL2 ( Figures 2 and S6). By combining the clustering of MD trajectories with the interaction fingerprint analysis (see the Computational Details section), we could identify detailed protein−ligand interactions stabilizing the complexes ( Figure 3).
First, we discovered the prominent role of residue Tyr205 5.35 in interacting with the 68 Ga-DOTA moiety for 61, 63, and 42% of the total simulation time in 68 Ga-DOTATOC/TATE/ NOC, respectively. In addition, the binding of ligands is stabilized by hydrophobic interactions of Y3/Nal3 with Tyr205 5.35 (persistence of 74, 86, and 78% for the three compounds, respectively). In turn, the ligand Y3/Nal3 residue interacts also with Ile195 (50, 51, and 45%) and Val280 6.59 (71, 82, and 74%) ( Figure 3). Interestingly, previous works reported the key role of Tyr205 5.35 and Ile195 in the interaction of SSTR2 with SST14 and octreotide (through F7 and F3, respectively), and it was also pointed out that these two residues can contribute to the selectivity of the different SSTR isoforms. 20,21,13,23,28 Furthermore, Phe294 7.35 and Ser279 6.58 (belonging to the hydrophobic sub-pocket constituted by TM6−7 and ECL3 13 ) seem to be involved in isoform selectivity as well, 23 and we found consistently their interaction with the disulfide bridge featured by all compounds (75/70% for 68 Ga-DOTATOC, 86/81% for 68 Ga-DOTA-TATE, 77/41% for 68 Ga-DOTANOC). Although the aforementioned interactions are conserved, we noticed some differences in the persistence of the DOTA-Tyr205 5.35 interaction. In particular, the replacement of tyrosine with naphthalene at position 3 in 68 Ga-DOTANOC results in a higher steric hindrance that destabilizes the interaction between the chelator and the protein, letting the 68 Ga-DOTA fluctuate more than that for the other compounds (see Figure S3). Conversely, the difference between 68 Ga-DOTATOC and 68 Ga-DOTATATE, which share the tyrosine residue at position 3, should be searched in the terminal residue T-ol8/T8. Indeed, the presence of a carboxylic negative charge on the T8 of 68 Ga-DOTATATE allows the peptide to interact with the basic Arg190, located at ECL2, for 59% of the simulation time ( Figures 2B and 3), while this interaction was found neither in 68 Ga-DOTATOC nor in 68 Ga-DOTANOC. At the same time, this polar interaction appears also to stabilize a second one between 68 Ga-DOTA and Arg184 (belonging to ECL2 as well) that in 68 Ga-DOTATATE was found in 64% of the simulation time, compared to 36 and 19% of 68 Ga-DOTATOC and 68 Ga-DOTANOC, respectively. Interestingly, previous studies reported the interaction between Arg184 and somatostatin 20 and the Arg184Ala mutation was found to decrease the potency of somatostatin, but not that of octreotide. 22 This last finding supports the crucial role of the deprotonated C-terminus (T8) in the interaction with the receptor. This can possibly explain also the higher selectivity of 68 Ga-DOTATATE toward the isoform 2, which is characterized by the presence of two arginine residues at the ECL2

Journal of Chemical Information and Modeling
pubs.acs.org/jcim Article (Arg184 and Arg190), whereas 68 Ga-DOTATOC binds the isoform 5 (belonging to the same sub-class SRIF1) that features on the ECL2 an acid residue (Glu182) instead of a basic one. 20,29 Further simulations of this class of radiopharmaceuticals interacting with both SSTR2 and SSTR5 are needed to confirm this hypothesis. The small differences in the peptide structure reflect not only on single-protein residue interaction but also on the overall dynamics of the receptor, especially of the very mobile ECL2. This loop is known to play a key role in the interaction with ligands 30,31 and it is characterized by opening and closing movements. 20,25 For this reason, we computed the percentage of MD frames in which the loop was found closed, according to the threshold values established in our previous work. 25 These thresholds refer to geometric parameters, namely, distances and angles, characterizing the movements of this loop. As a result, ECL2 was closed in about 50, 7, and 18% of the simulation time in 68 Ga-DOTATOC/TATE/NOC, respectively ( Figure S8). The marked differences can trace back to the characteristic behavior of DOTA and T-ol8/T8 moieties in the three compounds described above: in 68 Ga-DOTATOC, the 68 Ga-DOTA portion stably interacts with Tyr205 5.35 (thanks also to the presence of residue Y3 of the peptide), moving this group away from ECL2, and allowing its closure. In 68 Ga-DOTATATE the interaction with Tyr205 5.35 is still present, but the 68 Ga-DOTA moiety also strongly interacts with Arg184, mediated by the T8-Arg190 interaction. Since both arginine residues are located at the ECL2, their involvement in the interaction with the ligand very likely impairs its closure (see below). Differently from the other two compounds, in 68 Ga-DOTANOC the chelator loosely interacts with Tyr205 5.35 , leading to a higher oscillation that prevents the closure of ECL2.
Substitution of 68 Ga 3+ with 64 Cu 2+ Affects the Persistence of Ligand/SSTR2 Interactions. After assessing the role of the peptide moiety, we focused on the influence of the radionuclides by comparing 64 Cu-and 68 Ga-DOTATATE. Both gallium and copper ions are hexa-coordinated when in complex with DOTA (by four nitrogen and two oxygen atoms), showing a pseudo-octahedral geometry. 8 Due to the intrinsic properties of the two radionuclides (e.g., electric charge, van der Waals radius, Jahn−Teller distortion 8 ) the coordination geometries differ, showing a more elongated one in 64 Cu 2+ , compared to 68 Ga 3+ (Table S2). Keeping in mind the limitations associated with classical/force field-based MD simulations when describing such challenging types of atoms, 32 these differences reflected in the conformation assumed by the DOTA group during the MD simulations, where the free/  Figure S9). Focusing on the whole ligands, inspection of the dynamical behavior of 64 Cu-DOTATATE and 68 Ga-DOTATATE reveals that these compounds share overall the same pattern of interactions with SSTR2. The only relevant differences are found for residue T8 that interacts with Arg184 and Ser192 in 68 Ga-DOTATATE, and for residue Y3 that interacts with Asn276 6.55 in 64 Cu-DOTATATE (Figures 3 and 4). However, these differences do not significantly affect the dynamics of ECL2 (which was found to be closed in the 4 and 7% of the simulation time for 64 Cu-DOTATATE and 68 Ga-DOTATATE, respectively) ( Figure S8).

Changes of the Chelator Moiety Influence the Interactions at the Peptide C-Terminal.
In the third part of this work, we considered the same radionuclide ( 64 Cu 2+ ) and the same peptide (TATE) while considering three different chelators (DOTA, TETA, SAR). TETA differs from DOTA just for six atoms (60 vs 54 atoms, respectively), and both coordinate the copper ion through four nitrogen and two oxygen atoms (from the carboxylic groups). Differently from DOTA, in TETA the carboxylic groups are located one above, and one below the plane formed by the nitrogen atoms, conferring a slightly higher steric hindrance (average Connolly surface area 33 computed on the MD trajectories: 307 ± 3 Å 2 vs 329 ± 2 Å 2 , respectively). SAR has the lowest surface area (299 ± 2.0 Å 2 ) but, differently from the other chelators, it is associated with a butanediamide linker that increases its effective steric hindrance (405 ± 4 Å 2 ) as well as its flexibility ( Figure S3F). Besides the presence of a linker, another important difference between SAR and the other two chelators is the absence of negatively charged groups (i.e., carboxylic acid moieties), as the chelator coordinates the copper ion through its six nitrogen atoms.
Focusing on the MD simulations, Figure 5 shows the representatives of the most populated cluster of the 64 Cu-DOTA/TETA/SAR-TATE compounds. Consistently with what has been reported above for 68 Ga-based systems, also in this case, we found that all compounds interact with Phe294 7.35 and Ser279 6.58 via their disulfide bridge and with Tyr205 5.35 , Ile195, and Val280 6.59 through Y3 and the chelator moiety ( Figure 5). Noteworthy, in contrast to the other radiopharmaceuticals, in 64 Cu-SARTATE residue Tyr205 5. 35 interacts with the linker portion and not with the chelator (SAR) that remains thus free to oscillate during the MD trajectories.
Combining the clustering of MD trajectories with the interaction fingerprint analysis, we found that DOTA in 64 Cu-DOTATATE interacts with SSTR2 with the overall higher persistence compared to the other two chelators (TETA and SAR) ( Figure 4). As expected, this suggests that the presence of a hindering chelator destabilizes the interaction between the peptide and SSTR2; nonetheless, the improved stability of copper inside such chelator is known to yield high-quality images. 34 Interestingly, although the TATE peptide was common to all Cu-labeled compounds, the change of the chelator affected the interactions involving the C-terminus and the residues located at the ECL2. In detail, 64 Cu-DOTA/TETA-TATE interacted with Arg190 through the terminal T8, while the chelator moieties were involved with Arg184. Conversely, in 64 Cu-SARTATE, the terminal T8 was found to interact mainly with Arg184 (54% of the simulation time) and to a lesser extent with Arg190 (23%), whereas the SAR portion interacted poorly with the ECL2 (only 26% with Glu200) compared to .3%, respectively. ECL2 is colored in red, the copper ion is represented as a green sphere, the receptor is represented as a transparent white cartoon, the peptide TATE is colored in lilac, the chelators DOTA in dark cyan, TETA in purple, SAR in gold, and its linker in dark gray. The top panels report the superimposition with octreotide (green transparent sticks) taken from PDB7T11. In the bottom panels, the main interactions are shown as black dotted lines. the other chelators (Figures 4 and 5). This different behavior can be traced back to the total +2 net charge of 64 Cu-SAR (compared to −1 of 64 Cu-DOTA/TETA) that, despite the presence of a negative C-terminus, penalizes the interactions with basic residues.
Interestingly, when simulating SSTR2 in complex with 64 Cu-DOTA/TETA-TATE, the ECL2 was found to be closed in 2.7 and 3.5% of the simulation time, respectively, which is consistent with what was registered for 68 Ga-DOTATATE (6.6%). On the contrary, in 64 Cu-SARTATE, the ECL2 was able to close upon the binding pocket in 19.1% of the simulation time, similarly to 68 Ga-DOTANOC (18.1%) ( Figure S8). As mentioned above, both 64 Cu-DOTA/TETA-TATE were able to strongly interact with both Arg184 and Arg190 (such as 68 Ga-DOTATATE), while 64 Cu-SARTATE interacted only with Arg184 (like 68 Ga-DOTATOC/NOC). Therefore, these results suggest that the closure of the ECL2 loop is mostly impaired by the presence of strong polar interactions with the ligands, but also by the high fluctuations of the chelator.

■ CONCLUSIONS
In this study, we investigated the interaction of six radiopharmaceuticals with SSTR2, a key drug target for neuroendocrine tumors. We predicted the binding modes of these compounds and rationalized the role of the three different moieties characterizing this class (i.e., radionuclidechelator-peptide). Starting from the experimental structure of the receptor in complex with the somatostatin analogous octreotide, we generated the protein−ligand complexes that underwent to overall 15 μs of MD simulation time each. The analysis of the MD trajectories revealed that the substitution of the radionuclide ( 68 Ga 3+ with 64 Cu 2+ ) did not influence the dynamics and the main interactions established by the ligand, while the pattern of interaction of the C-terminus is strongly affected by changes of the chelator moiety (DOTA, TETA, SAR). The radionuclide-chelator portion is stabilized by crossinteractions between Tyr205 5.35 , Ile195, and the third residue of the peptide (Y3 for TOC and TATE, Nal3 for NOC). Furthermore, we found that upon small changes in the peptide structure (at the C-terminal T8/T-ol8 and at the third residue Y3/Nal3), the dynamics of both the chelator portion and SSTR2 strongly differ, possibly paving the way to a molecular rationalization of the differences in SSTR isoform selectivity.
The detailed molecular-level analysis presented in this study and the overall computational platform can be extended to other radiopharmaceuticals of this class, thus contributing to the rational design of new potent ligands targeting SSTR2. ■ COMPUTATIONAL DETAILS System Setup. The starting 3D structure of SSTR2 was retrieved from PDB ID7T11, 20 in which the receptor was solved in complex with the synthetic agonist octreotide and the G-protein. Missing atoms were added by structure refinement using Modeller10.2. 35 Given the close similarity between the peptide portion of the six radiopharmaceuticals and octreotide, it was reasonable to assume as the initial position of the ligands in the binding pocket those obtained by direct superimposition between the peptide portion and octreotide. The stability of the initial binding modes was thoroughly tested by monitoring the MD trajectories through analysis of RMSD/F values (see below). To reduce the computational cost, we did not include the G-protein in the structures. The ionization state of the residue side chains, the tautomeric states of histidine residues, and the Asn/Gln flipping were checked by the H++ server. 36 The CHARMM-GUI server 37 was used to embed the protein into a double layer of phosphatidylcholine (POPC, 70%) and cholesterol (30%). 38 The system was inserted in an OPC water box 39 and neutralized by adding K + and Cl − ions, reaching a 0.15 M concentration. The AmberTools20 software 40 was used to assign the force field lipid17 to POPC and cholesterol 41 and ff19SB to the protein. 42 The peptide portion of the ligands was obtained by manually changing the experimental structure of octreotide, solved in complex with SSTR2 (PDB ID7T11 20 ). The chelator structures were retrieved from the Cambridge Structural Database, 43 choosing those entries solved in complex with a radionuclide (DOTA ID 1136299, 44 TETA ID 624742, 45 SAR ID 915824 34 ). Finally, the chelator portions carrying the radionuclide were manually bound to the peptide N-terminal. For the generation of the ligand force field parameters, we combined two approaches: (1) one for the peptide and (2) one for the chelator-radionuclide portion. (1) The force field ff19SB was assigned to the peptide, and nonstandard residues (i.e., D-phenylalanine, naphthalene, D-tryptophan, threoninol) were parametrized as described previously. 25 (2) Given the peculiarity of the metal coordination bond involving 68 Ga 3+ and 64 Cu 2+ and their challenging parametrization, we used the Metal Center Parameter Builder (MCPB.py) procedure 46 implemented in Amber20. In detail, the 3D structures of the chelator-radionuclide and the first residue of the ligands ( D F1), obtained as described above, underwent quantum mechanics (QM) calculations at the Density Functional Theory (DFT) level 47 with the B3LYP functional using the Gaussian16 package (Revision A.03). 48 We performed geometry optimization on 68 Ga-DOTA using different basis sets in order to identify the best one for our system (Table S2). We compared the coordinaion distances between the experimental and QM optimized structure, and we computed the mean absolute error results (MAE). We observed that a large basis set does not lead to big differences in geometry optimization. Therefore, we employed the hybrid B3LYP functional, 49 in conjunction with the split-valence 6-31G(d,p) Gaussian basis set 50 to save computational time. For each compound, the ground-state structure was optimized, and then a full vibrational analysis was performed. In the case of 64 Cu-SAR, solvation effects were calculated using the integral equation formalism of the Polarized Continuum Model (IEF-PCM), 51 with water as the solvent 52 to avoid the collapse of the chelator. In all cases, the DFT-based structural parameters are in good agreement with the available experimental data (Table S2). The vibrational analysis results were used by MCPB.py to generate the bonded terms of the force fields. Then, on the optimized geometry we performed B3LYP/6-31G(d,p) single-point energy calculations to generate the atomic partial charges fitting the molecular electrostatic potential. We used the Merz−Kollman scheme 53 to construct a grid of points around the molecule under the constraint of reproducing the overall electric dipole moment of the molecule. Atomic partial charges were then generated through the two-step restrained electrostatic potential (RESP) method 54 implemented in the Antechamber package. 55 These steps enabled the generation of the force field of the chelator-radionuclide moieties using the General Amber Force Field 2 (GAFF2). 56 Journal of Chemical Information and Modeling pubs.acs.org/jcim Article MD Simulations. Each system underwent an energy minimization combining the steepest-descent and the conjugated gradient algorithms (2500 steps each) and applying positional restraints on the protein−ligand complex (10.0 kcal mol −1 Å −2 ) and on cholesterol and phosphate groups of phosphatidylcholine molecules (2.5 kcal mol −1 Å −2 ).
NVT and NPT equilibrations followed minimization, in which the positional restraints were incrementally reduced. The NVT equilibration was divided into two steps: first 125 ps with the same positional restraints of the minimization, then further 125 ps decreasing the restraint strength to 5.0 kcal mol −1 Å −2 for the protein and the ligands, and keeping 2.5 kcal mol −1 Å −2 for cholesterol and the phosphate groups of phosphatidylcholine molecules (overall NVT equilibration time = 250 ps). The following NPT equilibration was divided into four steps: (1) 125 ps using positional restraints of 2.5 kcal mol −1 Å −2 for the protein−ligand, and 1.0 kcal mol −1 Å −2 for the membrane components, (2) 500 ps using 2.5 kcal mol −1 Å −2 for protein−ligand and 1.0 kcal mol −1 Å −2 for the membrane, (3) 500 ps using 0.5 kcal mol −1 Å −2 for the protein−ligand and 0.1 kcal mol −1 Å −2 for the membrane, and finally (4) 500 ps using 0.1 kcal mol −1 Å −2 for the protein− ligand and leaving the membrane completely free to move (overall NPT equilibration time = 1.625 ns). We used the Langevin thermostat (1 ps −1 as collision frequency, 310 K) and the Berendsen barostat (1 atm), a cutoff of 9 Å, the time step was incremented from 1 to 2 fs with the SHAKE algorithm, 57 the Particle Mesh Ewald method for long-range electrostatics was employed. 58 The production runs were carried on for 3 μs, using the NPT ensemble and 4 fs as a time step, adopting the hydrogen mass repartition scheme. 59 Five replicas were generated for each system, resulting in a 15 μs simulation time. The MD simulations were conducted using the PMEMD module of Amber20, 40 and the trajectory frames were written every 100 ps.
MD Trajectory Analyses. MD replicas were first concatenated and the CPPTRAJ software 60 was used to perform a cluster analysis. A hierarchical algorithm 61 was used to group all frames into conformational clusters, according to the compound RMSD. Considering all cases, we found that three clusters sample significantly different conformations of the ligands which are, at the same time, reasonably populated (see Table S1). In all cases, the RMSD values of the ligands were computed on all of the heavy atoms, after aligning the backbone of the receptor in the MD trajectories, with respect to the first frame of the production run. Interaction fingerprints were computed using the ProLIF Python library 62 on all of the frames of the MD trajectories. The numbers of interactions were combined for all replicas and converted into persistence of interactions (%).