MD Modeling of Pharmacological Vector-Receptor Pairs for Specific Drug Delivery to the Tumor: Atomic/Molecular Mechanisms of RGD-Peptide Embedding in the αvβ3-Integrin Receptor ()
1. Introduction
The work investigates promising and, in our opinion, the most feasible mechanisms of specific binding and delivery of drugs (molecular reagents) to the tumor in the “VECTOR-RECEPTOR” format [1]-[6]. It also describes the existing mechanism of amino acid transfer through the LAT1 receptor used in boron-neutron capture therapy (BNCT) [7]-[10]. The method of computational molecular dynamics (MD) is used to calculate the interactions of pharmacological vector-receptor pairs. It should be noted that the experimental study of the above issues is difficult, and therefore it is advisable to use computational and simulation methods of analysis. Undoubtedly, effective MD methods are among the most powerful approaches applied to the tasks of modeling interactions and specific effects of drugs and receptors. In this paper, we aim to study the dynamic and structural features of peptide-protein complexes and pharmacological vector-receptor pairs. The following vectors are proposed as vectors for targeted delivery of boron-containing nanoparticles:
1) A peptide containing the amino acid sequence RGD (L-arginine, glycine, L-aspartic acid). This sequence is able to specifically bind to αvβ3 integrin, which is overexpressed in the tumor node.
2) Folic acid, which binds to the folic acid receptor FR-α, overexpressed by many tumor cells.
Thus, the purpose of these computational experiments is to simulate the interaction processes of the pharmacological pair “VECTOR-RECEPTOR” (Table 1):
Table 1. Pair RGD-integrin αvβ3.
Vector |
Receptor |
RGD (peptide containing amino acid sequence—L-arginine, glycine, L-aspartic acid) |
Integrin αvβ3 |
Here, αvβ3 is an integrin, a vitronectin receptor, a marker of tumor neoangiogenesis, which consists of two parts: integrin alpha V and integrin beta 3. It should be noted that RGD peptides are able to bind to each of the subunits individually or both simultaneously. It is currently unknown whether the specificity of RGD peptides differs for murine and human αvβ3 integrins. It is proposed to calculate the binding of RGD to both mouse and human αvβ3 integrins. The same applies to other possible routes of specific boron delivery to the tumor. It is also proposed for comparison to evaluate the specificity of the only currently used mechanism of boron delivery, amino acid transport using LAT1 formulations (see Table 2):
Table 2. RGD-integrin pair αvβ3.
Vector |
Receptor |
RGD (peptide containing amino acid sequence—L-arginine, glycine, L-aspartic acid) |
Integrin αvβ3 |
folic acid |
folate receptors FR-α (the tumor marker) FR-β |
boron phenylalanine |
amino acid transporters: LAT1 (mostly) LAT2 ATB (0, +) |
Computational methods offer important ways to describe and understand atomic and molecular phenomena in many physical and biochemical systems. These methods are widely used in physics, chemistry, and biological sciences, as well as in the computer-aided design of new materials and drugs [11]-[15]. The source file of the PDB (Protein Data Bank) is 3ZE2 (the full module of the built-in memory alja-beta-3 and the additional module RGD; PDB DOI: https://doi.org/10.2210/pdb3ZE2/pdb) [16]. The above molecular structure consists of several chains (MOL_ID: 1; CHAIN: A, C; FRAGMENT: RESIDUE 32-488; MOL_ID: 2; CHAIN: B, D; FRAGMENT: RESIDUES 27 - 498; …) containing peptides (MOL_ID: 5; MOLECULE: RGD PEPTIDE; CHAIN: I, J; … ORGANISM_SCIENTIFIC: HOMO SAPIENS; HUMAN; CHINESE HAMSTER) (Figure 1 and Figure 2).
For the numerical MD simulation implemented in this study, the multi-purpose AMBER-18 package is used with a fast implementation of the “pmemd.cuda” module on a computer with a CPU/GPU cluster. The numerical experiment significantly accelerates the operation of the AMBER-18 computer code (“pmemd.cuda”) and is aimed at determining the specificity or binding of the RGD peptide system to the αvβ3 integrin receptor. MD analysis data provide localization and conformation of RGD-peptide in the region of binding to the αvβ3 integrin receptor at the atomic and molecular level. It is worth noting that in modern life science, the issues of studying peptide-protein and protein-protein interactions, interactions of membranes, receptors and enzymes with coenzymes, peptides and substrates are relevant, taking into account various chemical interactions with solvent molecules and specific surface fragments. These are the most intriguing problems in connection with very important applications in nanobiotechnology, biochemistry, etc. These interactions are determined by hydrogen, hydrophobic, ionic,
![]()
Figure 1. The general view of the RGD-integrin pair αvβ3. The αvβ3 receptor is an integrin (on the right) and a small peptide containing the amino acid sequence RGD (L-arginine, glycine, L-aspartic acid) (on the left).
Figure 2. The general view of the RGD-integrin pair αvβ3. Structural loops A, B, C, and D and RGD, a peptide containing an amino acid sequence (L-arginine, glycine, and L-aspartic acid), are shown in the αvβ3 integrin receptor.
electrostatic interactions, as well as their combinations.
2. Experimental Part
In this section, we present a description of the main parameters and algorithms used in MD computer modeling. In this study, we used the CPU and GPU-based computing environments to perform MD simulations using the AMBER package. We used several computing environments, such as a 16-core cluster, a Geforce graphics processor (GPU = (GTX 1080Ti)) to implement MD modeling in Amber18 (pmemd; MD/AMBER) with GPU acceleration (pmemd.cuda). We refer to the program code and the Amber 2018 reference manual [11]-[15]. The calculations were performed on the servers of the heterogeneous HybriLIT platform of the Multifunctional Information and Computing Complex (MIC) of the M.G. Meshcheryakov Laboratory of Information Technologies of the Joint Institute for Nuclear Research (MLIT JINR, Dubna). The heterogeneous platform consists of the Govorun supercomputer and the HybridIT training and testing center, which are a two-component system that includes a central processor based on the latest Intel architectures (Intel Xeon Phi and Intel Skylake processors) and a graphics processor based on NVIDIA DGX-1 Volta. Some calculations were performed on the local server of the FLNP (Frank Laboratory of Neutron Physics) at JINR using two 4-core 64-bit Intel Xeon E5-2640 processors with a clock speed of 2.4 GHz and 8 GB of RAM running the Linux CentOS operating system version 8. NVIDIA Corporation GP104 graphics card is installed on the local server [GeForce GTX 1070] with Intel Xeon E7 v4/Xeon E5 v4/Xeon E3 v4/Xeon processors. We have implemented the main MD production models (CPU/GPU) (also common to many other types of simulations) for the PDB ID: 3ze2; DOI: https://doi.org/10.2210/pdb3ZE2/pdb). In this case, the MD simulation was performed using the Amber 18 code (CPU/GPU environment).
The MD modeling on the RGD-peptide + αvβ3-integrin receptor + water molecular system (see Figure 3(a) and Figure 3(b) below) was carried out in three stages: energy minimization, NVT and NPT pulse relaxation procedures. The preparatory steps (minimizing energy consumption, balancing NVT and NPT) were performed for 10 ns during dynamic relaxation changes. It is worth noting that the Amber graphics code is used, and for performance reasons, the graphics code does not recalculate unrelated list cells. Once the system is balanced, fluctuations in the size of the box should be small and should not create problems with productive calculations. We symmetrically surrounded the RGD + αvβ3-integrin + water system in all spatial directions using pressure regulation and periodic boundary conditions. Thus, the entire system, RGD-peptide + αvβ3-integrin, solvated with TIP3P water, was modeled as a very large thermostatically controlled cubic system (see visualizations below, snapshots of the original models in Figure 3. The MD calculations with the included Shake algorithm for hydrogen bonds without estimating the strength of hydrogen bonds and the Langevin thermostat were performed at T = 303 K. Periodic boundary conditions of constant pressure, anisotropic pressure coupling, and an unrelated 10 Å boundary were used. CPU/GPU accelerated Amber18 simulation (pmemd/pmemd.cuda) is performed in three main configuration steps, which contain the main MD production models:
1) Reduce the load on the system to loosen damaged contacts;
2) Slowly heat the system to the set temperature;
3) Bring the system to equilibrium at the set temperature.
Figure 3. The initial system, RGD-peptide + αvβ3-integrin receptor, solvated with water (TIP3P model) for subsequent MD calculations.
At the first stage, we minimized the system by imposing restrictions on the main atoms and on some other atoms. The initial charge distribution of titrated residues corresponds to the jumpprc.constph and jumpprc.conste states in tleap by default. Minimization was performed using a “sander” instead of “pmemd” to ensure and evaluate the energy behavior during the minimization phase. We then performed constant volume heating simulations on minimized structures to slowly heat the system from 10 to 303 K during several initial simulation steps (2.0 ns) with a target temperature of 303 K. We have retained the restrictions on the core atoms, but with a weaker force constant than that used in the minimization. These calculations were performed using pmemd, pmemd.MPI, Sander, or Sander.MPI instead of using GPU-accelerated code with pmemd.cuda. In the next step, we balance the system (RGD-peptide + αvβ3 integrin receptor + water) at a target temperature of 303 K. Then, an explicit simulation of the solvent’s behavior at constant pressure is performed to ensure that the density of the system is stabilized. The simulation is performed for 50 - 100 ns in implicit and explicit solvents.
3. Results and Discussions
Atomic/Molecular Mechanisms of RGD Peptide Insertion into
the αvβ3 Integrin Receptor
The results of the MD calculations are shown in Figures 4-6. Comparison of the spatial positions of the RGD-peptide + αvβ3-integrin receptor system, solvated with TIP3P water, according to MD calculations, in the initial (a) state of the system (t = 0) and at t = 30 ns, t = 60 ns and t = 90 ns relaxed states. In this case, two RGD peptides were modeled, located outside and inside the αvβ3 integrin receptor. One of the two RGD peptides is a peptide from the original PDB file localized inside the αvβ3 integrin receptor. The other RGD peptide in its initial position is located outside the region inside the receptor, freely diffusing throughout the entire region of the imitation cell and naturally contacting and binding to the αvβ3 integrin.
In the process of dynamic changes over a long time interval from 0 to 100 nanoseconds (50 million steps of integration of the equations of motion, for several of the models described above), comparative images of multiple spatial configurations of the RGD-peptide + αvβ3-receptor integrin system solvated TIP3P with water were obtained. It should be noted that MD calculations indicate a slow, long-term relaxation process of the entire RGD-peptide + αvβ3-integrin receptor system, the results of which are the determination of the spatial position of RGD-peptides and their specific binding to the αvβ3-integrin receptor.
In the search for optimal configurations of RGD-peptide the MD results of the atomic and molecular level indicate on a spontaneous (natural) binding of the RGD-peptide-vector to the αvβ3 integrin-receptor. In this aspect, multiple MD model calculations were performed, in which several configurations of the RGD-peptide + αvβ3 integrin receptor system were reproduced, each of which consisted of equilibrium states as a result of dynamic relaxation processes in the time range from zero to 100 nanoseconds. Based on a set of scenarios and analysis of MD statistical data, we discovered and observed the process of natural incorporation of RGD-peptide into the protein-binding pocket of the αvβ3 integrin receptor (Figures 4-6).
As noted above, one of the two RGD peptides, RGD-1, is a peptide in the original PDB file located inside the αvβ3 integrin receptor, which remains localized throughout 50 million steps of dynamic changes and integration of the equations of motion. Another RGD peptide located outside the αvβ3 integrin receptor, diffusing freely throughout the entire region of the modeling cell, naturally contacts and binds to αvβ3 integrin.
Peptide RGD-1 (res.1354-res.1359):
---GLY1354-ARG1355-GLY1356-ASP1357-SER1358-PRO1359---
Peptide RGD-2 (res.1360-res.1365):
---GLY1360-ARG1361-GLY1362-ASP1363-SER1364-PRO1365---
Some design features of RGD peptides should be noted, namely, the “pharmacological” expediency of choosing the structure of peptides in order to effectively bind them to the αvβ3 integrin receptor. The beginning and end of the RGD peptide, respectively, consist of (GLY, PRO) amino acids with hydrophobic uncharged side radicals.; (ARG) is a side-chain amino acid that is positively charged in aqueous solutions with neutral or near-neutral pH values, and the isoelectric point of these amino acids is significantly higher than 7.0, and they have an excessive positive charge at neutral pH, as well as proteins containing a predominant amount of such amino acids. acid residues are called alkaline because their isoelectric point is also greater than 7.0; (ASP) is an amino acid with a negative charge in the side chain of the radical, and the isoelectric point of these amino acids is less than 7.0; (SER) is an amino acid with a hydrophilic uncharged radical (protein regions containing such an amino acid residue are capable of hydration and interaction through hydrogen bonds with other similar residues).
Next, Figure 4 shows the final configuration positions of the RGD-1,2 peptide system with αvβ3 integrin in the final (t = 100 ns) relaxed state. The final stages of the binding of RGD peptides (1 and 2) to the αvβ3 integrin receptor as a result of the natural (spontaneous) incorporation of RGD peptide into the αvβ3 integrin receptor. Relaxed states are shown, (a) configurations with a broken line and (b) configurations with an atomic structure in which two RGD peptides in the αvβ3 integrin receptor are inserted parallel to each other inside a protein-binding pocket.
As noted above, each of the RGD-1,2 peptides consists of six amino acid residues with specific hydrophobic and hydrophilic properties. The first RGD-1 peptide, located inside the αß3 integrin receptor, remains localized throughout all dynamic changes. Another peptide, RGD-2, is located outside the αvβ3 integrin receptor, can freely spread over the entire area of the modeling cell and, thus, naturally contacts and binds to αvβ3 integrin in the same place as RGD-1 – inside the αvβ3 integrin receptor. It is noteworthy, however, that the externally introduced RGD-2 peptide is oriented perpendicular to the RGD-1 peptide when approaching it, but then, due to the spatial structure, feature, or nature of the “protein binding pocket”, it will be oriented parallel to the RGD-1 peptide (see Figures 4-6). This fact may be of great importance from the point of view of “pharmacological” perspectives, namely for understanding the atomic and molecular mechanisms of the specificity and binding of RGD peptides to the αvβ3 integrin receptor, predicting the functional aspects of the interaction of the VECTOR-RECEPTOR pharmacological pair (RGD-peptide + αvβ3 integrin receptor) for the delivery of a specific drug to tumors based on computer-aided design and modeling of protein structure.
Next, Figure 5 shows four consecutive snapshots (configuration states) of the RGD-1 and RGD-2 peptide vectors in the dynamic process of their interaction with the αvβ3 integrin receptor. Figure 5 shows the molecular configurations of RGD peptides as a result of the natural (spontaneous) dynamic diffusion unification of RGD-2-peptide from the outside with RGD-1-peptide inside the region of the αβ3 integrin receptor for temporary states (a) t = 0, (b) t = 10 ns, (c) t = 80 nc and (d) t = 100 ns, respectively. In this case, the spatial positions and transformations of structures from the initial (GLY) to the final (PRO) amino acid residues in the peptide chains are indicated. The configurations snapshots are shown in Figure 5. The configuration images demonstrate how two RGD peptides in the αvβ3 integrin receptor eventually integrate parallel to each other into the protein-binding pocket shown above in Figure 4.
![]()
![]()
Figure 4. (a)-(b) The final stages of the binding of RGD peptides (1 and 2) to the αvβ3-integrin receptor as a result of the natural (spontaneous) incorporation of RGD-peptide into the αvβ3-integrin receptor. Relaxed states are shown, (a) configurations with a broken structure, and (b) configurations with an atomic structure in which two RGD peptides in the αvβ3 integrin receptor are inserted parallel to each other into a protein-binding pocket (several amino acids are shown).
It is worth noting that each of the two peptides-biovectors RGD-1 and 2 consists of six amino acid residues, with the RGD-1 peptide located inside the αvβ3 integrin receptor remaining localized throughout all dynamic changes. Another peptide, RGD-2, is located outside the αvβ3 integrin receptor, freely diffuses over the entire area of the model cell and, thus, naturally contacts and binds to αvβ3 integrin in the same place as RGD-1, inside the αvβ3 integrin receptor. It is noteworthy that the RGD-2 peptide introduced from the outside is oriented perpendicular to the RGD-1 peptide when approaching it, but then, due to the peculiarity of the spatial structure (or the nature of the “protein binding pocket”), it will be oriented parallel to the RGD-1 peptide (see Figures 4-6 and images of RGD-1 and RGD peptides-2 in the diagram below).
![]()
Figure 5. (a)-(d) Sequential snapshots of the interaction of RGD peptides (1 and 2) in the αvβ3 integrin receptor as a result of natural (spontaneous) dynamic diffusive administration of RGD peptide from outside the αvβ3 integrin receptor region. The configuration states are shown (a) t = 0, (b) t = 10 ns, (c) t = 80 ns, and (d) t = 100 ns, where two RGD peptides in the αvβ3 integrin receptor are eventually inserted parallel to each other into the protein-binding pocket of amino acid residues shown in Figure 6.
According to the above scheme 1, the distribution of all possible distances between the peptides Δ [RGD-1 (GLY, PRO) − RGD-2 (GLY, PRO) was calculated], the results of which are presented using distance diagrams (Figure 6 and Figure 7). The MD results shown in Figure 6 and Figure 7 also demonstrate the weakened spatial orientation of the RGD-1 and RGD-2 peptide vectors in the structure of the αvβ3 integrin receptor. In this case, the RGD-2 peptide, located outside the αvβ3 integrin receptor, eventually integrates almost parallel to the
Figure 6. The binding specificity and orientation features of the RGD-1 and RGD-2 peptide vectors in the spatial structure of the αvβ3-integrin receptor in the initial t = 0, intermediate and final t = 100 ns relaxed states (top view, side view) and the dynamics of the distances between the peptides Δ [RGD-1 (GLY, PRO) − RGD-2 (GLY, PRO)] (at the bottom).
RGD-1 peptide during prolonged 100-ns interactions and dynamic changes and spontaneous diffusion movement. The diagrams of the spatial distribution of distances in Figures 6-8 indicate a “dense” introduction of the external peptide RGD-2 with the peptide RGD-1 localized inside the αvβ3 integrin receptor.
Figure 7. The binding specificity and orientation features of the RGD-1 and RGD-2 peptide vectors in the spatial structure of the αvβ3-integrin receptor in the initial t = 0, intermediate and final t = 100 ns relaxed states (top view) and the dynamics of the distances between the peptides Δ [RGD-1 (GLY, PRO) − RGD-2 (GLY, PRO)] (at the bottom).
Figure 8 shows the comparative molecular configurations of the RGD-1 and RGD-2 peptide vectors in the structure of the αvβ3 integrin receptor in the initial t = 0 and final t = 100 ns relaxed states (side and top view). Figure 8 illustrates the specifics of the binding and orientation of the RGD-1 and RGD-2 peptide vectors within the spatial structure of the αvβ3 integrin receptor.
Figure 8. The binding specificity and orientation features of the RGD-1 and RGD-2 peptide vectors in the spatial structure of the αvβ3-integrin receptor in the initial t = 0 and final t = 100 ns relaxed states (side and top view).
Thus, the MD simulation time of 100 ns has seen long enough, specifically the binding of the peptide to the receptor occurs at around 30 - 40 ns due to the high mobility of the RGD peptide, which is initially located outside the receptor region. So far, important conformational changes and stable binding of the peptide to the receptor have occurred much faster, in less time than 100 ns, due to the specific or diffusive nature of the molecular system modeled by the MD method. For the experimental validation we relate the peptides binding to integrin receptors, such as RGD peptides, where in the experimental aspect, the peptides-receptors were conjugated and specifically targeted tumor cells expressing integrin receptors. The RGD-peptide was chosen as the target molecule because of its conveniently small size, specificity, and the very important role that RGD’s main target, αvβ3 integrin, plays in tumor angiogenesis, proliferation, and metastasis.
4. Conclusions
In this work, computer MD calculations of the interaction processes of the pharmacological pair “VECTOR-RECEPTOR” were performed to simulate promising mechanisms and processes of delivering a specific drug to a tumor. As a result of long 100-second MD calculations, the configuration states were obtained and the spatial positions of the RGD-peptide + αvβ3-integrin receptor system, solvated with water (TIP3P model), were determined. Based on a comparative MD analysis of the configuration states of the RGD-peptide + integrin αvβ3 system, the spatial positions of two RGD peptides located outside and inside the αvβ3 integrin receptor were identified. Distance diagrams between two RGD (RGD-1), which is a peptide from the original PDB file (localized inside the αvβ3 integrin receptor), and RGD-2, a peptide outside the receptor that freely diffuses over the entire area of the model cell, were calculated and compared. The configuration images and diagrams of the interaction distance of RGD-1,2 peptides with αvβ3 integrin demonstrate the features of natural contact and specific binding at the atomic-molecular level.
In conclusion, in current study the MD simulations were used to evaluate the RGD peptide binding to αvβ3 integrin receptor for drug delivery, thereby modelling two RGD peptides and their dynamic interactions over 100-ns, showing diffusional processes with a spontaneous binding of the external peptide to the receptor. Peptides binding to integrin receptors, such as RGD peptides, were experimentally conjugated and specifically targeted tumor cells expressing integrin receptors. In this relation, the RGD peptide was chosen as the target molecule because of its conveniently small size, specificity, and the very important role that RGD’s main target, αvβ3 integrin, plays in tumor angiogenesis, proliferation, and metastasis. With regard to the experimental validation, the MD simulated peptide-conjugated boron-containing compounds for BNRT are being developed as a new strategy for targeted delivery of anti-cancer therapeutics. It is worth noting that the peptide-binding receptors play a crucial role in cancer therapy, because they offer a specific and selective target for the delivery of anti-cancer agents. The peptide component is MD simulated and designed to bind to a specific peptide receptor on the surface of cancer cells, which makes it possible to selectively deliver a drug load to cancer cells without affecting normal cells.
Acknowledgements
The work was performed within the framework of the state assignment of the Ministry of Science and Higher Education of the Russian Federation (No. 1024011000011-7-1.4.2; 3.5.2 Conjugates of boron-containing quantum dots with biovectors for the diagnosis and boron-neutron capture therapy of superficial malignant tumors (FEEM-2024-0011)). The work was carried out within the framework of a joint JINR (Joint Institute for Nuclear Research) and Dubna State University research program. Numerical calculations were performed on the servers of the heterogeneous Hybridity platform of the Multifunctional Information and Computing Complex (MIC) of the M.G. Meshcheryakov Laboratory of Information Technologies (MLIT, JINR, Dubna, Moscow Region, Russian Federation).