**1. Introduction**

Biomolecular materials containing the arginine-glycine-aspartate (RGD) sequence are always at the center of biophysics research in their application such as in the bone scaffold, synthesis, and regeneration of tissue and cartilage [1,2], in imaging as radiotracers [3–6], for cancer therapy [7–9], and in targeted drug delivery [10–15]. The RGD motif peptide serves as a primary integrin recognition site in extracellular matrix proteins since it has a strong binding affinity to integrins, which are heterodimeric cell surface receptors and mediate cell-extracellular matrix adhesion [16–19]. Out of 24 known integrins, one third of them bind to the RGD motif as the primary recognition sequence in their ligands, which makes RGD an attractive target in numerous drug delivery systems [20–24]. The increased application of integrins in drug development and their functions in the physiological processes requires a complete understanding of the structure and properties of the RGD peptides at body temperature (310K), crucial in the design of selective inhibitors. Therefore, a detailed study of the structure and properties of 1FUV in RGD peptides at 310 K is of particular

**Citation:** Baral, K.; Adhikari, P.; Jawad, B.; Podgornik, R.; Ching, W.-Y. Solvent Effect on the Structure and Properties of RGD Peptide (1FUV) at Body Temperature (310 K) Using Ab Initio Molecular Dynamics. *Polymers* **2021**, *13*, 3434. https://doi.org/ 10.3390/polym13193434

Academic Editor: Célio Bruno Pinto Fernandes

Received: 17 September 2021 Accepted: 29 September 2021 Published: 7 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

interest in delineating the modifications that occur with a rise or fall of body temperature. Moreover, the short- and long-range interactions are crucial for the molecular recognition and self-assembly of biological macromolecules, and consequently determining the microscopic parameters such as partial charge and frequency-dependent dielectric functions based on more accurate quantum mechanical calculation will facilitate the fundamental understanding of electrostatic, polar, and van der Waals-London dispersion interaction in biological (macro)molecules.

The peptide-water interfacial reactions are of paramount importance in biological systems as solvation is ubiquitous in biological interaction with water as an inevitable component in blood and multicomponent body fluids. Understanding the peptide-water interaction is also essential regarding the RGD-targeted drug delivery mechanisms. It is, therefore, crucial to investigate the consequences of peptide-water interaction and aqueous solvation and its effects on internal bonding, charge distribution, and dielectric response. The solubility of 1FUV in water is generally attributed to the hydrogen bonding between water molecules and oxygen from the peptide backbone. To explore the details of peptide solvation properties, molecular dynamics (MD), molecular mechanics, and Monte Carlo methods, all based on specific force field parameters, are widely used. Nevertheless, the information on the role played by the hydrogen bond (HB) within RGD and at the RGD-water interface is in general still missing, as are the relevant ab initio MD (AIMD) studies [12,25–27]. Classical MD studies cannot provide a quantitative description of essentially quantum HBs (O–H and N–H) since they depend on the details of the force field parameterization. Consequently, the need for AIMD simulation in the study of dry and aqueous solvated peptide/protein over classical MD has been repeatedly recognized [28–31] and a more accurate AIMD calculation with a sufficiently large number of water molecules is highly desirable for a complete understanding of 1FUV aqueous solvation.

Although experimental techniques such as nuclear magnetic resonance (NMR), surfaceenhanced Raman spectroscopy, X-ray photoelectron spectroscopy, etc., are used to investigate the structure and properties of RGD containing peptides [16,32–35], the experimental studies specifically involving 1FUV are scant. Moreover, although experimental studies can probe the consequent changes of the peptide structure due to water contact, they cannot provide critical information and details of microscopic properties and interfacial reaction mechanisms at the atomistic scale. Computer-assisted atomistic simulation thus seems the best alternative to investigate the molecular behavior and provide new insights to complement the meager experimental data [36,37], in many cases guiding the experiments as well. Nevertheless, only a few computational studies have so far been reported on the 1FUV peptide [38].

To answer the issues raised above and provide new insights into the study of the 1FUV peptide, this study is designed to simulate four models of dry and solvated 1FUV at 0 K and body temperature, namely 1FUV0, 1FUV310, 1FUVS0, and 1FUVS310, respectively. Our study provides a deeper understanding of the temperature-dependent structure, fundamental properties, peptide-water interfacial reaction mechanism, and solvation effect. To the best of our knowledge, our effort is the first example of a study dedicated to the 1FUV peptide using AIMD that ensures the reliability of the simulated structures and calculated properties. The rest of the paper is structured as follows. In the next section, a brief description of the modeling technique followed by the simulation workflow is described. The main part is the result and discussion section, in which we discuss our findings and articulate the future prospects for applying AIMD to other complex biomolecular systems.

#### **2. Materials and Methods**

#### *2.1. Modeling 1FUV Peptide*

The RGD peptides that specifically bind to *αvβ<sup>3</sup>* and *αvβ<sup>5</sup>* integrins have medical significance in designing the inhibitors of tumor and retinal angiogenesis [39,40], in osteoporosis [41], and in targeted drugs for tumor vasculature [42]. Considering the importance

of *αvβ<sup>3</sup>* and *αvβ<sup>5</sup>* integrins [16] and as 1FUV binds only to these integrins, we choose to focus on 1FUV in this work.

The initial structure of 1FUV was taken from the Research Collaboratory for Structural Bioinformatics protein data bank (PDB) generated from the NMR measurements [43]. The stored PDB file contains about 19 such models with the same composition, the number of atoms, molecular weight, and the chemical formula, but slightly different molecular volumes. In this study, we chose the first model out of 19, which has 135 atoms comprised of 11 amino acids. The initial atomic positions of 1FUV were enclosed in a sufficiently large supercell of size 37.12 Å × 34.41Å × 28.95 Å to avoid any unintended interaction due to periodic boundary conditions. This structure was then fully relaxed by applying a density functional theory (DFT)-based Vienna ab initio simulation package (VASP) [44,45] known for its efficiency in structure optimization and energy minimization to obtain 1FUV0, the structure at 0 K. We used a projector-augmented wave method [46,47] with Perdew-Burke-Ernzerhof potential for the exchange correlation functional within the generalized gradient approximation [48]. We employed a relatively high cutoff energy of 500 eV and set the electronic convergence criterion at 10−<sup>5</sup> eV. The force convergence criterion for ionic relaxation was set at 10−<sup>3</sup> eV/Å. A single K-point calculation was used as the simulated models are sufficiently large supercells.

The fully relaxed structure at 0 K was then heated slowly to 310 K, equivalent to human body temperature, using NVT ensemble within 4 ps. A Nose thermostat was used to control the temperature of the heat bath [49,50]. A time step of 0.5 fs was used for the ionic motion during the simulation to ensure accuracy for integration of the equation of motion, especially for the light atoms such as hydrogen. The peptide model at this elevated temperature was then equilibrated for 5 ps to obtain the equilibrium structure at 310 K, viz. 1FUV310. Out of 5 ps run, the fluctuations of temperature and pressure were minimized by the initial 4 ps run, and the averaged output structure is taken from the last 1 ps run, thus used as a production run. The solvated model at 0 K (1FUVS0) was constructed based on 1FUV0 by surrounding the peptide structure with 80 water molecules. The 1FUVS0 was then again fully relaxed by using VASP. Similar to the dry model, this solvated structure was then heated to 310 K to obtain the 1FUVS310 model using NVT ensemble within 4 ps. The same protocol as used in the dry model was adopted to obtain the equilibrated structure of the solvated model at 310 K.

To reveal the structural stability of simulated models, we present the root-meansquare deviation (RMSD) of atomic positions of the peptide structure during equilibration at 310 K in Figure 1a. The time-RMSD profile of global structure shows a smooth increase followed by a leveling off at ~ 1.90 Å in the dry model and ~ 1.15 Å in the solvated model, respectively, indicating that the equilibrium has been reached. The constant RMSD shown in the inset for the last 1 ps run, used as a production run, ensures the sufficient equilibration of simulated structures. Figure 1b shows the velocity autocorrelation function (VACF), including exponential fitting, during the equilibration of dry and solvated structures at 310 K which confirms the loss of initial velocity from the previous configuration and attest to sufficient equilibration of the systems. It also implies that the amplitude of VACF is smaller in the solvated model and dies out quicker than that in the dry model. In addition, the exponential fitting shows the solvated model reaches equilibrium faster than the dry model.

#### *2.2. Calculation of Properties*

The AIMD simulated and VASP relaxed structures discussed above are then used as input for the electronic structures and dielectric response calculations using the inhouse developed package of the orthogonalized linear combination of atomic orbitals (OLCAO) methodology [51]. The combination of these two DFT methods, VASP and OLCAO, is found to be robust and highly effective in studying the electronic structure and bonding properties, especially for complex biomolecules [52–55]. In recent years, we have further demonstrated that the combination of these two methods is extremely effective

in dealing with new and emerging problems for complex protein structures related to SARS-CoV-2 [56–60]. *Polymers* **2021**, *13*, x FOR PEER REVIEW 4 of 16

**Figure 1.** Calculated (**a**) RMSD and (**b**) VACF including exponential fitting, of simulated dry and solvated 1FUV peptide at 310 K. **Figure 1.** Calculated (**a**) RMSD and (**b**) VACF including exponential fitting, of simulated dry and solvated 1FUV peptide at 310 K.

[56–60].

*2.2. Calculation of Properties* The AIMD simulated and VASP relaxed structures discussed above are then used as input for the electronic structures and dielectric response calculations using the in-house developed package of the orthogonalized linear combination of atomic orbitals (OLCAO) methodology [51]. The combination of these two DFT methods, VASP and OLCAO, is found to be robust and highly effective in studying the electronic structure and bonding properties, especially for complex biomolecules [52–55]. In recent years, we have further The OLCAO is an all-electron method using local density approximation of DFT in which the Gaussian-type orbitals are employed for the expansion of the atomic basis set. There are many advantages of using OLCAO for electronic structure calculation, such as the flexibility of the basis sets choice, lower computational cost, and ease of bonding and charge transfer analysis using the Mulliken scheme [51]. The implementation of localized atomic orbitals in the basis expansion enables us to quantify the charge transfer and interatomic bonding via effective charge (*Q*\*) on each atom and bond order (BO) value between pairs of atoms (*ραβ*) using the Mulliken scheme [61,62] as

$$\mathbf{Q}\_{\mathbf{a}}^{\*} = \sum\_{i} \sum\_{n, \text{occ}} \sum\_{j, \beta} \mathbf{C}\_{\text{i}n}^{\*n} \mathbf{C}\_{j\beta}^{n} \mathbf{S}\_{\text{i}n, j\beta} \tag{1}$$

$$\rho\_{a\beta} = \sum\_{n,\alpha c} \sum\_{i,j} \mathbf{C}\_{ia}^{\*\eta} \mathbf{C}\_{j\beta}^{\mathbf{u}} S\_{ia,j\beta} \tag{2}$$

There are many advantages of using OLCAO for electronic structure calculation, such as the flexibility of the basis sets choice, lower computational cost, and ease of bonding and charge transfer analysis using the Mulliken scheme [51]. The implementation of localized where *Siα,j<sup>β</sup>* is the overlap matrix between the basis Bloch sums of the orbital index (*i*, *j*) and atomic specification (*α, β*). N is the band index, *i*, *j* are the orbital quantum numbers, and *Cj<sup>β</sup>* is the eigenvector coefficient.

atomic orbitals in the basis expansion enables us to quantify the charge transfer and interatomic bonding via effective charge (*Q\**) on each atom and bond order (BO) value between pairs of atoms ( ) using the Mulliken scheme [61,62] as \* \* , , . *n n i j i j i n occ j <sup>Q</sup> C C S* <sup>=</sup> (1) \* , , , *n n i j i j n occ i j C C S* <sup>=</sup> (2) where *Siα,jβ* is the overlap matrix between the basis Bloch sums of the orbital index (*i*, *j*) and atomic specification (*α, β*). N is the band index, *i*, *j* are the orbital quantum numbers, From the calculated value of *Q*\* the charge transfer between the ions due to atomic interaction can be quantified in terms of the partial charge (PC). It is defined as the deviation of *<sup>Q</sup>*\* from the charge of neutral atom *<sup>Q</sup>*<sup>0</sup> and is given by <sup>∆</sup>*<sup>Q</sup>* <sup>=</sup> *<sup>Q</sup>*<sup>0</sup> <sup>−</sup> *<sup>Q</sup>*\*. The BO value calculated from the above equation gives the direct quantitative measure of bond strength between a pair of bonded atoms. The sum of all BO values in the system is the total bond order (TBO), which is a quantum mechanically calculated parameter. It is a single quantum mechanical metric helpful in assessing the internal cohesion and strength in a material [63]. The use of TBO to characterize the internal strength of a material and correlate it with the calculated physical properties is a novel and highly appealing concept.

and *Cjβ* is the eigenvector coefficient. From the calculated value of *Q*\* the charge transfer between the ions due to atomic interaction can be quantified in terms of the partial charge (PC). It is defined as the devi-OLCAO is equally suitable to calculate the imaginary part, ε<sup>2</sup> (¯*h*ω), of the frequencydependent complex dielectric function within the random phase approximation of interband optical transition theory according to the following equation:

$$\begin{split} \varepsilon\_{2}(\hbar\omega) &= \frac{\epsilon^{2}}{\pi m\omega^{2}} \int\_{BZ} dk^{3} \sum\_{\begin{subarray}{c} \mathbf{n}, \mathbf{l} \\ \mathbf{n}, \mathbf{l} \end{subarray}} \left| \left\langle \Psi\_{\mathbf{n}}(\stackrel{\rightarrow}{\mathbf{k}}, \stackrel{\rightarrow}{\mathbf{r}}) \right| - i\hbar \stackrel{\rightarrow}{\nabla} \left| \Psi\_{\mathbf{l}}(\stackrel{\rightarrow}{\mathbf{k}}, \stackrel{\rightarrow}{\mathbf{r}}) \right\rangle \right|^{2} \times f\_{\mathbf{l}}(\stackrel{\rightarrow}{\mathbf{k}}) \left[ 1 - f\_{\mathbf{n}}(\stackrel{\rightarrow}{\mathbf{k}}) \right] \\ &\quad \delta \left[ E\_{\mathbf{n}}(\stackrel{\rightarrow}{\mathbf{k}}) - E\_{\mathbf{l}}(\stackrel{\rightarrow}{\mathbf{k}}) - \hbar\omega \right] \end{split} \tag{3}$$

rial [63]. The use of TBO to characterize the internal strength of a material and correlate it

with the calculated physical properties is a novel and highly appealing concept.

#### **3. Results**

#### *3.1. Analysis of 1FUV Structures*

The RGD peptide, 1FUV, studied in this work contains 11 amino acid sequence ACDCRGDCFCG, namely alanine (Ala, A), cysteine (Cys, C), asparagine (Asp, D), arginine (Arg, R), glycine (Gly, G), and phenylalanine (Phe, F). As this amino acid sequence contains four C, this peptide is also known as RGD-4C, in which the cysteine position is at 2, 4, 8, and 10th place of the sequence, named as Cys1, Cys2, Cys3, and Cys4, respectively. Each Cys amino acid contains the S atom and the arrangement of two Cys, and hence the S–S bond plays an important role in the configuration and dipole moment of this peptide [64]. Although the presence of four Cys groups in RGD-4C allows three possible combinations of S-S bonds, only two natural configurations of this peptide exist based on the S–S bonding arrangement [16]. The peptide with Cys1–Cys4 and Cys2–Cys3 disulfide bonding arrangement is known as RGD-A isomer, while the one with Cys1–Cys3, Cys2– Cys4 bonding arrangement is known as RGD-B isomer [16]. As RGD-A is a far stronger binder to integrin *αvβ<sup>3</sup>* than RGD-B [16], it is reasonable to explore details of the structure and properties of RGD-A.

The relaxed structures for RGD-A peptide (1FUV) at 0 K and 310 K in the dry and aqueous environment are shown in Figure 2. 1FUV0 shows only a minor change in its structure due to VASP relaxation compared with the original NMR structure [16] taken from the PDB. However, there are noticeable changes in its conformation after AIMD simulation at 310 K. It shows a rotation of amino acid groups and changes in the bond length (BL) and BO value, which will be discussed in detail later. In the case of the solvated model, there again appear notable changes in the peptide structure. More water molecules interact with the peptide and there is a rotation of amino acid groups at 310 K as compared with 0 K. In both dry and solvated models, the S–S bond containing amino acids come closer at 310 K as compared with those at 0 K. However, the average S–S BL increases and hence its strength decreases at 310 K. It is noticed that the configurations of HBs, O–H, change at 310 K as compared with 0 K, resulting in the overall decrease in its strength in the dry model. However, in the solvated model, a large contribution of O–H bonding strength comes from peptide-water interaction due to which its value increases at 310 K. These differences are, although small, important on issues related to human health. *Polymers* **2021**, *13*, x FOR PEER REVIEW 6 of 16

**Figure 2.** Structures of dry (**a**, **b**) and solvated (**c**, **d**) 1FUV peptide at 0 K and 310 K, respectively. (N = Blue, C = Brown, O = Red, H = Black, and S = Yellow). **Figure 2.** Structures of dry (**a**,**b**) and solvated (**c**,**d**) 1FUV peptide at 0 K and 310 K, respectively. (N = Blue, C = Brown, O = Red, H = Black, and S = Yellow).

simulated models is shown in Figure 3. The top of the valence band is set at 0 eV. The overall features of TDOS at 0 K and 310 K in dry and solvated models look similar as they contain the same amino acid groups; however, there is a decrease in the highest occupied molecular orbital (HOMO)‐lowest unoccupied molecular orbital (LUMO) gap with the increase in temperature. The position of the energy difference between HOMO and LUMO plays the main role in various chemical reactions. The number of peaks decreases in solvated models as compared with the dry model, but the overall intensity increases. The HOMO‐LUMO gap for 1FUV0 is ~ 3.40 eV, while its value in 1FUV310 decreases to ~ 2.10 eV due to a temperature rise, indicating that the molecule is chemically more reactive at a higher temperature. The HOMO‐LUMO gap decreases largely in solvated models as compared with dry models. In 1FUVS0, this gap is ~ 0.3 eV, but it disappears in 1FUVS310.

*3.2. Electronic Structure and Interatomic Bonding*

#### *3.2. Electronic Structure and Interatomic Bonding*

Electronic structures are the fundamental properties of a material that help to understand many other physical properties. The calculated total density of states (TDOS) for the simulated models is shown in Figure 3. The top of the valence band is set at 0 eV. The overall features of TDOS at 0 K and 310 K in dry and solvated models look similar as they contain the same amino acid groups; however, there is a decrease in the highest occupied molecular orbital (HOMO)-lowest unoccupied molecular orbital (LUMO) gap with the increase in temperature. The position of the energy difference between HOMO and LUMO plays the main role in various chemical reactions. The number of peaks decreases in solvated models as compared with the dry model, but the overall intensity increases. The HOMO-LUMO gap for 1FUV0 is ~ 3.40 eV, while its value in 1FUV310 decreases to ~2.10 eV due to a temperature rise, indicating that the molecule is chemically more reactive at a higher temperature. The HOMO-LUMO gap decreases largely in solvated models as compared with dry models. In 1FUVS0, this gap is ~0.3 eV, but it disappears in 1FUVS310. *Polymers* **2021**, *13*, x FOR PEER REVIEW 7 of 16

**Figure 3.** Calculated TDOS for dry and solvated 1FUV peptide at 0 K and 310 K. **Figure 3.** Calculated TDOS for dry and solvated 1FUV peptide at 0 K and 310 K.

The nature of bonding and its quantitative assessment is important in the electronic structure study of the material. While most of the MD studies interpret the bonding analysis based exclusively on the geometric distance configuration, OLCAO can quantify the strength of each bonded pair by providing the BO value. We calculate the TBO values in all 4 models and resolve them into partial pair‐resolved BO (PBO) components. This is an unorthodox approach that generalizes the bonding between individual bonded pairs simply and effectively to reveal the hidden details which cannot be extricated by other methods, either experimental or computational. The nature of bonding and its quantitative assessment is important in the electronic structure study of the material. While most of the MD studies interpret the bonding analysis based exclusively on the geometric distance configuration, OLCAO can quantify the strength of each bonded pair by providing the BO value. We calculate the TBO values in all 4 models and resolve them into partial pair-resolved BO (PBO) components. This is an unorthodox approach that generalizes the bonding between individual bonded pairs simply and effectively to reveal the hidden details which cannot be extricated by other methods, either experimental or computational.

In general, the BO value scales inversely with BL; however, the actual BO value of a bonded pair depends on its bonding environment and the atoms surrounding it. The calculated values of TBO and PBO are listed in Table 1, showing that at 310 K the TBO value decreases in the dry model, resulting from the elongation of some bonds. The TBO value of the solvated model increases at 310 K, opposite to the case of the dry model, which is due to the formation of more bonds between the peptide backbone and water at this elevated temperature. Figure 4 shows the comparison of PBO values for each bonded pair in simulated models at 0 K and 310 K in the form of a vertical bar diagram. The scale breaks In general, the BO value scales inversely with BL; however, the actual BO value of a bonded pair depends on its bonding environment and the atoms surrounding it. The calculated values of TBO and PBO are listed in Table 1, showing that at 310 K the TBO value decreases in the dry model, resulting from the elongation of some bonds. The TBO value of the solvated model increases at 310 K, opposite to the case of the dry model, which is due to the formation of more bonds between the peptide backbone and water at this elevated temperature. Figure 4 shows the comparison of PBO values for each bonded pair in simulated models at 0 K and 310 K in the form of a vertical bar diagram. The scale breaks

are used along the vertical axis to delineate the minor differences between PBO values

model. The PBO for the O–H pair increases sharply in the solvated model at 310 K as compared with the model at 0 K, which arises due to the interaction of more water molecules with the peptide, and an increase in the HBs population at 310 K. It should be emphasized again that the small changes in these numbers at different temperature can be

We provide further information on the strength of individual bonded pairs presenting a complex distribution of BO versus BL in Figure 5. It clearly shows a wide variety of bonding pairs ranging from a strong covalent bond originating from N–H, N–C, C–H, C– C, C–O, and O–H pairs to weaker covalent bonds at longer BL. HBs are extremely weak but not entirely negligible because of their sheer numbers. The first group of strong covalent bonds arises from N–H, C–H, and O–H pairs with BL less than 1.2 Å . The next group of covalent BO pairs arises from C–C, C–O, and N–C bonds, whose BL lies in the range

significant in the context of the human body.

are used along the vertical axis to delineate the minor differences between PBO values and show their magnified image. Figure 4 shows that there is a minor change in the PBO values due to the rise of temperature, except in the case of the O–H pair in the solvated model. The PBO for the O–H pair increases sharply in the solvated model at 310 K as compared with the model at 0 K, which arises due to the interaction of more water molecules with the peptide, and an increase in the HBs population at 310 K. It should be emphasized again that the small changes in these numbers at different temperature can be significant in the context of the human body. *Polymers* **2021**, *13*, x FOR PEER REVIEW 8 of 16 1.2–1.7 Å . The C–C bonds show two clusterings, one having a high BO greater than 0.5 e from C–C double bonds, whereas the other group has a BO value less than 0.4 e from single C–C bonds. The BO values for N–C pairs also exhibit a group separation, having BO values greater than 0.4 e and less than 0.4 e.


**Table 1.** Calculated TBO and PBO values for each bonded pair in simulated models. **Table 1.** Calculated TBO and PBO values for each bonded pair in simulated models.

**Figure 4.** Comparison of PBO values for each bonded pair in dry and solvated 1FUV at 0 K and 310 K. **Figure 4.** Comparison of PBO values for each bonded pair in dry and solvated 1FUV at 0 K and 310 K.

The HBs are distributed in the range of 1.5–2.5 Å , with some O–H pairs having relatively higher BO values up to 0.1 e. In solvated models, the covalent O–H bonds within the water molecule have BO values higher than 0.2 e, while those bonds originating from the next neighbor or weak HBs beyond 2.5 Å have almost zero BO values. The bonding distributions due to C–S and S–S pairs have BO values comparable or slightly higher than those of strong HBs. The S–S bonding configurations are important in peptide structure, and they show significant changes in BO values due to the rise of temperature and solvation. Especially in the solvated model at 310 K, the BO of S–S bonds has a large difference. Wang et al. studied the influence of disulfide bond on properties of RGD peptide using MD and found that this bond causes the restriction on peptide mobility and affects its We provide further information on the strength of individual bonded pairs presenting a complex distribution of BO versus BL in Figure 5. It clearly shows a wide variety of bonding pairs ranging from a strong covalent bond originating from N–H, N–C, C–H, C–C, C–O, and O–H pairs to weaker covalent bonds at longer BL. HBs are extremely weak but not entirely negligible because of their sheer numbers. The first group of strong covalent bonds arises from N–H, C–H, and O–H pairs with BL less than 1.2 Å. The next group of covalent BO pairs arises from C–C, C–O, and N–C bonds, whose BL lies in the range 1.2–1.7 Å. The C–C bonds show two clusterings, one having a high BO greater than 0.5 e from C–C double bonds, whereas the other group has a BO value less than 0.4 e from single C–C bonds. The BO values for N–C pairs also exhibit a group separation, having BO values greater than 0.4 e and less than 0.4 e.

dipole moment [64]. The bonded pairs H–H and H–S with BL larger than 2.5 Å have BO values close to zero. The HBs are distributed in the range of 1.5–2.5 Å, with some O–H pairs having relatively higher BO values up to 0.1 e. In solvated models, the covalent O–H bonds within the water molecule have BO values higher than 0.2 e, while those bonds originating from the next neighbor or weak HBs beyond 2.5 Å have almost zero BO values. The bonding distributions due to C–S and S–S pairs have BO values comparable or slightly higher than those of strong HBs. The S–S bonding configurations are important in peptide structure, and they show significant changes in BO values due to the rise of temperature and solvation. Especially in the solvated model at 310 K, the BO of S–S bonds has a large difference. Wang et al. studied the influence of disulfide bond on properties of RGD peptide using MD and found that this bond causes the restriction on peptide mobility and affects its dipole

moment [64]. The bonded pairs H–H and H–S with BL larger than 2.5 Å have BO values close to zero. *Polymers* **2021**, *13*, x FOR PEER REVIEW 9 of 16

**Figure 5.** BO versus BL plot for dry and solvated 1FUV at 0 K and 310 K. **Figure 5.** BO versus BL plot for dry and solvated 1FUV at 0 K and 310 K.

#### *3.3. Hydrogen Bonding Analysis 3.3. Hydrogen Bonding Analysis*

Hydrogen bonding is especially important in a biomolecular system since it provides the key information to understand many intriguing phenomena, making the detailed analysis of HBs in the solvated models crucial. Unfortunately, in most MD studies, the analysis of HBs is carried out based exclusively on structural data, such as the BL and the location of HBs, failing to provide a quantitative assessment of HB strength and their local bonding arrangements. On the contrary, we bring forth quantitative information and deeper analysis of HBs using ab initio calculation in terms of BO‐BL plot, as shown in Figure 5. It shows there exists some strong HB with BO value ~ 0.1 e, while most of the HBs have BO less than 0.1 e. At the body temperature, the overall strength of HB from the Hydrogen bonding is especially important in a biomolecular system since it provides the key information to understand many intriguing phenomena, making the detailed analysis of HBs in the solvated models crucial. Unfortunately, in most MD studies, the analysis of HBs is carried out based exclusively on structural data, such as the BL and the location of HBs, failing to provide a quantitative assessment of HB strength and their local bonding arrangements. On the contrary, we bring forth quantitative information and deeper analysis of HBs using ab initio calculation in terms of BO-BL plot, as shown in Figure 5. It shows there exists some strong HB with BO value ~ 0.1 e, while most of the HBs have BO less than 0.1 e. At the body temperature, the overall strength of HB from the O–H bonding pair slightly decreases in the dry model while it increases in the solvated model.

O–H bonding pair slightly decreases in the dry model while it increases in the solvated model. In the case of solvated models, the change in the HB pattern is more complicated due to the presence of H2O molecules which increase the population and strength of O–H bonds at 310 K. The more HBs imply the stronger binding force of biomolecule with solvent. There is one HB originating from the N–H pair within the peptide backbone whose BO value at 310 K increases in the dry model and decreases in the solvated model because of the rotation and changes in the configuration of the peptide structure. The results on In the case of solvated models, the change in the HB pattern is more complicated due to the presence of H2O molecules which increase the population and strength of O–H bonds at 310 K. The more HBs imply the stronger binding force of biomolecule with solvent. There is one HB originating from the N–H pair within the peptide backbone whose BO value at 310 K increases in the dry model and decreases in the solvated model because of the rotation and changes in the configuration of the peptide structure. The results on HB analysis in solvated models indicate a significant difference in the human body as compared with the vicinal environment of a dry peptide at zero temperature.

#### HB analysis in solvated models indicate a significant difference in the human body as *3.4. Partial Charge Distribution*

compared with the vicinal environment of a dry peptide at zero temperature. *3.4. Partial Charge Distribution* The PC distribution of biomolecules is important for determining the long-range electrostatic and polar interaction involving the molecule. It entails a deep qualitative understanding of the structure and reactivity of a molecule. Rather than just making a rudimen-The PC distribution of biomolecules is important for determining the long-range electrostatic and polar interaction involving the molecule. It entails a deep qualitative understanding of the structure and reactivity of a molecule. Rather than just making a rudimentary judgment by dividing PC into positive, negative, and neutral as adopted in many simplified calculations, we provide a quantitative assessment of PC of every atom in

tary judgment by dividing PC into positive, negative, and neutral as adopted in many simplified calculations, we provide a quantitative assessment of PC of every atom in the

the simulated models. As the PC gives the quantity of charge transfer in atomic interaction, a positive value implies a loss, and a negative value means the gain of electronic charge. *Polymers* **2021**, *13*, x FOR PEER REVIEW 10 of 16

> The distribution of calculated PC for each atom in dry and solvated models is shown in Figures 6 and 7, respectively, which provide a wealth of information that corroborates the electronic structure results discussed above. The solid symbols represent PC of ions at 0 K, while the hollow symbols represent PC at the body temperature. O (P) and H (P) denote O and H atoms from the peptide backbone while O (W) and H (W) denote those from the water molecule. In both dry and solvated models at 310 K, small changes appear in the PC of ions compared with their values in corresponding models at 0 K. Figures 6 and 7 show that H and S lose charges, N gains charge, while C is both losing or gaining charges. In the dry model, O gains charge, while in the solvated model (Figure 7), few O atoms from the peptide backbone lose charge, and the remaining gain charge. The calculated PC values for each atom are scattered in a certain range rather than clustering and overlapping to each other, which implies the PC value is not the same even for the same atom in the same model, as their values depend on the local bonding environment. The PC of O and H in water and peptide backbone is not the same, implying that their local bonding characteristics are completely different (Figure 7). The dispersed distribution of calculated PC for each atom implies that the quantum mechanical calculation is essential to deal with such type of problem rather than using MD studies in which the PC values are predefined and kept constant in the simulation where the local environment changes. simulated models. As the PC gives the quantity of charge transfer in atomic interaction, a positive value implies a loss, and a negative value means the gain of electronic charge. The distribution of calculated PC for each atom in dry and solvated models is shown in Figures 6 and 7, respectively, which provide a wealth of information that corroborates the electronic structure results discussed above. The solid symbols represent PC of ions at 0 K, while the hollow symbols represent PC at the body temperature. O (P) and H (P) denote O and H atoms from the peptide backbone while O (W) and H (W) denote those from the water molecule. In both dry and solvated models at 310 K, small changes appear in the PC of ions compared with their values in corresponding models at 0 K. Figures 6 and 7 show that H and S lose charges, N gains charge, while C is both losing or gaining charges. In the dry model, O gains charge, while in the solvated model (Figure 7), few O atoms from the peptide backbone lose charge, and the remaining gain charge. The calculated PC values for each atom are scattered in a certain range rather than clustering and overlapping to each other, which implies the PC value is not the same even for the same atom in the same model, as their values depend on the local bonding environment. The PC of O and H in water and peptide backbone is not the same, implying that their local bonding characteristics are completely different (Figure 7). The dispersed distribution of calculated PC for each atom implies that the quantum mechanical calculation is essential to deal with such type of problem rather than using MD studies in which the PC values are predefined and kept constant in the simulation where the local environment changes.

**Figure 6.** Calculated partial charge for each atom in dry 1FUV at 0 K and 310 K. (The subscripts 0 and 310 denote the temperature 0 K and 310 K, respectively.). **Figure 6.** Calculated partial charge for each atom in dry 1FUV at 0 K and 310 K. (The subscripts 0 and 310 denote the temperature 0 K and 310 K, respectively). *Polymers* **2021**, *13*, x FOR PEER REVIEW 11 of 16

*3.5. Dielectric Response*

crystals and glasses [69–72].

**Figure 7.** Calculated partial charge for each atom in solvated models of 1FUV at 0 K and 310 K. (The subscripts 0 and 310 denote the temperature 0 K and 310 K, respectively.). **Figure 7.** Calculated partial charge for each atom in solvated models of 1FUV at 0 K and 310 K. (The subscripts 0 and 310 denote the temperature 0 K and 310 K, respectively).

The dielectric response of biomolecules in an aqueous solution has been studied for

Figure 8 (a) shows the spectrum for ε<sup>2</sup> (ℏω) in 1FUV0 calculated using OLCAO, displaying a sharp peak at ~ 5.70 eV and one broad peak at ~ 14.30 eV. A similar spectrum is observed in the dry model at 310 K; however, the broad peak position shifts slightly to the lower energy side. In solvated models, in addition to the two peaks observed in dry models, there appears another sharp peak at ~ 0.1 eV, but with slightly less intensity (Figure 8 (b)). One of the striking features observed is that at the body temperature, the absorption end of ε<sup>2</sup> (ℏω) shifts to lower energy which is consistent with the decrease in the HOMO‐ LUMO gap at 310 K. We want to make it clear that our aim in this work is to compare the full dielectric spectra at 0 K and body temperature in case of a dry and solvated environment and the actual calculation of the dielectric constant is out of the scope of this work.

solving the Poisson–Boltzmann equation [68]. However, it is interesting to note that the dielectric properties of the 1FUV are seldom reported. We present a detailed analysis of the imaginary part, ε<sup>2</sup> (ℏω), of the frequency-dependent dielectric function, based on the ab initio calculation under random phase approximation such as those used in inorganic

#### *3.5. Dielectric Response*

The dielectric response of biomolecules in an aqueous solution has been studied for decades [65–67] and remains a subject of active research. The static dielectric constant of peptide/protein plays a crucial role in the calculation of the electrostatic field obtained by solving the Poisson–Boltzmann equation [68]. However, it is interesting to note that the dielectric properties of the 1FUV are seldom reported. We present a detailed analysis of the imaginary part, ε<sup>2</sup> (*h*¯ ω), of the frequency-dependent dielectric function, based on the ab initio calculation under random phase approximation such as those used in inorganic crystals and glasses [69–72].

Figure 8a shows the spectrum for ε<sup>2</sup> (*h*¯ ω) in 1FUV0 calculated using OLCAO, displaying a sharp peak at ~5.70 eV and one broad peak at ~14.30 eV. A similar spectrum is observed in the dry model at 310 K; however, the broad peak position shifts slightly to the lower energy side. In solvated models, in addition to the two peaks observed in dry models, there appears another sharp peak at ~0.1 eV, but with slightly less intensity (Figure 8b). One of the striking features observed is that at the body temperature, the absorption end of ε<sup>2</sup> (*h*¯ ω) shifts to lower energy which is consistent with the decrease in the HOMO-LUMO gap at 310 K. We want to make it clear that our aim in this work is to compare the full dielectric spectra at 0 K and body temperature in case of a dry and solvated environment and the actual calculation of the dielectric constant is out of the scope of this work. *Polymers* **2021**, *13*, x FOR PEER REVIEW 12 of 16

**Figure 8.** Calculated imaginary part of the complex dielectric function for (**a**) dry and (**b**) solvated **Figure 8.** Calculated imaginary part of the complex dielectric function for (**a**) dry and (**b**) solvated 1FUV at 0 K and 310 K.

#### **4. Discussions and Summary**

1FUV at 0 K and 310 K.

**4. Discussions and Summary** In this section, we discuss the results presented above and summarize the main findings to answer the queries raised in the introduction. More accurate AIMD and DFT calculations adopted in this work provide us a wealth of quantitative information on the fundamental properties of 1FUV, which are not yet available or cannot be attained by an experimental protocol. Our findings provide a comparative analysis of the structure and properties of 1FUV and a better understanding of its solvation properties at 0 K and body temperature, 310 K. The key message is that there are discernible changes in the peptide structure and properties at 310 K, which could point to the potential implications for hu-In this section, we discuss the results presented above and summarize the main findings to answer the queries raised in the introduction. More accurate AIMD and DFT calculations adopted in this work provide us a wealth of quantitative information on the fundamental properties of 1FUV, which are not yet available or cannot be attained by an experimental protocol. Our findings provide a comparative analysis of the structure and properties of 1FUV and a better understanding of its solvation properties at 0 K and body temperature, 310 K. The key message is that there are discernible changes in the peptide structure and properties at 310 K, which could point to the potential implications for human health issues if such calculations can eventually be applied to large biomolecular systems.

man health issues if such calculations can eventually be applied to large biomolecular systems. Our results show noticeable changes in the conformation of amino acids at body temperature than at 0 K. The major structural differences observed in dry and solvated peptides at body temperature help in designing the specificity of integrins. The disulfide Our results show noticeable changes in the conformation of amino acids at body temperature than at 0 K. The major structural differences observed in dry and solvated peptides at body temperature help in designing the specificity of integrins. The disulfide bridge connectivity greatly determines the orientation of the RGD motif and peptide structure which may dictate selectivity towards different integrins. The HOMO-LUMO gap

bridge connectivity greatly determines the orientation of the RGD motif and peptide structure which may dictate selectivity towards different integrins. The HOMO‐LUMO gap of 1FUV decreases due to the rise of temperature and solvation effect. Such infor-

interaction within 1FUV and other biomolecules. The TBO value decreases in the dry model while it increases in the solvated model at 310 K. The dry model at 0 K is more cohesive than that at 310 K, characterized by a higher TBO value, while it is just the opposite in the solvated models. The BO‐BL plot analysis provides a detailed picture of interatomic bonding within the peptide structure and its explicit solvent environment. The solvated model shows a considerable peptide‐water interactions occur at 310 K, increasing the HBs population than at 0 K. The detailed analysis of HBs and quantification of their strength is crucial in understanding the biological interactions inside the human body and for the design and delivery of RGD peptide-targeted drugs. The partial charge analysis provides another set of significant parameters to define the reactivity of biomolecules and their interactions with the integrins. It shows that the PC value strongly depends on the temperature as well as the solvated environment. Furthermore, the calculated PC values scatter in a certain range, even for the same atom in the same model. Quantum mechanical calculation is, therefore, absolutely necessary to capture all these crucial differences and

to provide more accurate results inaccessible to classical MD studies.

of 1FUV decreases due to the rise of temperature and solvation effect. Such information is quite helpful to understand the chemical reactions involving biomolecules.

The bonding analysis is a significant result that helps to understand the interatomic interaction within 1FUV and other biomolecules. The TBO value decreases in the dry model while it increases in the solvated model at 310 K. The dry model at 0 K is more cohesive than that at 310 K, characterized by a higher TBO value, while it is just the opposite in the solvated models. The BO-BL plot analysis provides a detailed picture of interatomic bonding within the peptide structure and its explicit solvent environment. The solvated model shows a considerable peptide-water interactions occur at 310 K, increasing the HBs population than at 0 K. The detailed analysis of HBs and quantification of their strength is crucial in understanding the biological interactions inside the human body and for the design and delivery of RGD peptide-targeted drugs. The partial charge analysis provides another set of significant parameters to define the reactivity of biomolecules and their interactions with the integrins. It shows that the PC value strongly depends on the temperature as well as the solvated environment. Furthermore, the calculated PC values scatter in a certain range, even for the same atom in the same model. Quantum mechanical calculation is, therefore, absolutely necessary to capture all these crucial differences and to provide more accurate results inaccessible to classical MD studies.

Another important result is the calculation of the dielectric response function and its relation to the dielectric constant. Although this topic has attracted researchers from various backgrounds for a long time, unified and reliable ab initio calculations are still out of reach. As an initial assessment, we analyze the imaginary part of the dielectric function of the 1FUV peptide and study the differences in the dielectric spectra due to the rise in temperature and the solvation effect. Such optical spectra are helpful to estimate the long-range van der Waals–London interaction in biomolecules, as shown in the study of carbon nanotubes [73]. Overall, our results show apparent differences in bonding, PC, and dielectric function between the 1FUV at 0 K and the body temperature in a dry and solvated environment. These results can be quite helpful to develop and refine the force field parameters used in MD simulation. Unfortunately, as there is a lack of theoretical studies in 1FUV and the focus of the existing experimental study [16] is different, we are unable to compare our results with the previous results. However, we anticipate that our results will inspire future experimental and/or theoretical works to investigate further the 1FUV and related peptides, including the effect of salts and different numbers of water molecules, which significantly influence peptide structure and properties.

In the past, we have successfully demonstrated the applicability of AIMD in the study of oxide glassy materials and fluoride salts [74–77]. Our accomplishment in the present work is applying AIMD to a relatively small peptide structure that looks quite promising. It ensures the feasibility and encourages a further study of other possibly more complex biomolecules by using AIMD, of course assuming sufficient computing resources to be available.

**Author Contributions:** W.-Y.C., K.B. and P.A. conceived the project. K.B. and W.-Y.C. performed the calculations and made figures. K.B. and W.-Y.C. drafted the paper with inputs from P.A., R.P. and B.J. All authors participated in the discussion and interpretation of the results. All authors edited and proofread the final manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project is funded by the National Science Foundation of USA: RAPID DMR/ CMMT-2028803.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data that supports the results in this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This research used the resources of the National Energy Research Scientific Computing Center supported by DOE under Contract No. DE-AC03-76SF00098 and the Research Computing Support Services (RCSS) of the University of Missouri System. R.P. would like to acknowledge the support from the 1000-Talents Program of the Chinese Foreign Experts Bureau, the School of Physics, the University of the Chinese Academy of Sciences, Beijing, and the Institute of Physics of the Chinese Academy of Sciences, Beijing.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

