Molecular dynamics simulations of cRGD-conjugated PEGylated TiO 2 nanoparticles for targeted photodynamic therapy

The conjugation of high-affinity cRGD-containing peptides is a promising approach in nanomedicine to efficiently reduce off-targeting effects and enhance the cellular uptake by integrin-overexpressing tumor cells. Herein we utilize atomistic molecular dynamics simulations to evaluate key structural-functional parameters of these targeting ligands for an effective binding activity towards α V β 3 integrins. An increasing number of cRGD ligands is conjugated to PEG chains grafted to highly curved TiO 2 nanoparticles to unveil the impact of cRGD density on the ligand’s presentation, stability, and conformation in an explicit aqueous environment. We find that a low density leads to an optimal spatial presentation of cRGD ligands out of the “stealth” PEGylated layer around the nanosystem, favoring a straight upward orientation and spaced distribution of the targeting ligands in the bulk-water phase. On the contrary, high densities favor clustering and internalization of cRGD ligands in the inner region of the PEGylated layer, driven by a concerted mechanism of enhanced ligand-ligand interactions and reduced water accessibility over the ligand’s molecular surface. These findings strongly suggest that the ligand density modulation is a key factor in the design of cRGD-targeting nanodevices to maximize their binding efficiency into over-expressed α V β 3 integrin receptors.


Introduction
Active targeting strategies, exploiting the biological interaction between ligands on the surface of nanoparticles (NPs) and the cell target, are receiving increasing attention because, in several cases, they have been found to increase the therapeutic efficacy with respect to passive targeting strategies based on physicochemical properties of the NPs and on the enhanced permeability and retention (EPR) effect of tumor cells [1].
Among cell targets, integrins are a promising class of receptors for two main reasons: first, they are overexpressed in many types of cancer cells [2] and tumor microvasculature and, secondly, they regulate the communication between cells and their microenvironment. Therefore, the surface functionalization of nanomaterials with integrin-specific ligands not only increases their binding affinity to cancer cells but also enhances the cellular uptake of the NPs thanks to the intracellular trafficking of integrins [3,4] or through more complex routes, such as phagocyte hitchhiking [5].
Recent advances in integrin-targeted therapeutics have proved that fine-tuning the density of high-affinity integrin ligands conjugated to the surface of nano-delivery systems, e.g., metal or metal oxide NPs, is a viable way to improve the clinical outcome [6]. On the one hand, it is intuitive to assume that increasing the availability and, therefore, the spatial density of targeting ligands on the NP might lead to a higher probability of binding events to targeted integrins, enhancing their internalization by tumor cells [7]. On the other hand, there is experimental evidence that increasing ligand density as a tunable parameter of a nano-delivery system might cause "offtargeting" effects impairing their targeting activity, e.g. off-targeting effect due to impaired stealth properties [8], over-crowding effects due to ligand-ligand interactions impairing their proper orientation for binding, over binding to cell receptors decreasing the uptake efficiency by target cells, among others [9].
From an empirical perspective, Abstiens et al. [10] have studied the interplay between the ligand density and the integrin-mediated cellular uptake of polymeric NPs decorated with cyclic tripeptide Arg-Gly-Asp molecules (cRGDfK). The authors found that NPs with short poly(ethylene glycol) (PEG) chains and high cRGDfK-conjugation density have superior efficiency in targeting αvβ3 integrins than NPs with longer PEG chains and lower cRGDfK-conjugation densities. In a more recent paper [8], the same authors have introduced a new concept of "steric shielding" for cRGDfK-conjugated polymeric NPs to avoid their off-targeting to endothelial integrins. They argued that, by fine-tuning the cRGDfK surface density and the length of PEG chains, undesired offtarget accumulation may be prevented, while a satisfactory on-target efficacy may be maintained. In the same vein, Valencia et. al. [11] have found that copolymeric NPs composed of cRGDyK-conjugated PLGA(poly lactic-co-glycolic acid)-PEG chains show optimal ligand densities up to 50% for enhanced tumor cell uptake and minimum phagocytic clearance.
Theoretical predictions based on mathematical tumor models have shown that controlling the density of RGD-targeting ligands is a workable way to enhance their tumor-penetration efficiency with a stout effect on the integrin binding kinetics [12].
The authors have demonstrated that the addition of cyclic RGD targeting ligands to poly(amidoamine) (PAMAM) dendrimers, by impairing the ligand's binding affinity to integrins, enhances the tumor penetration by PAMAM-RGD conjugates. Recently, Biscaglia et al. [13] have investigated, experimentally and by means of atomistic MD simulations, the role of a cationic spacer on the presentation and orientation of cRGDyK ligands conjugated to PEGylated gold NPs. They found that cRGDyK ligands linked to PEGylated Au NPs through an oligolysine spacer have enhanced targeting activity towards integrins, whereas, in the absence of this spacer (when the cyclic RGD is directly bonded to PEG chains), the targeting efficacy is impaired.
Titanium dioxide (TiO2) NPs are attracting increasing attention in nanomedicine because of their unique photocatalytic properties, high biocompatibility, low toxicity [14], low cost, and high chemical stability [15]. In particular, they are used in photodynamic therapy for cancer treatment [16][17][18], being excellent ROS (reactive oxygen species) generators under light irradiation, as drug delivery systems [19,20], and for cell imaging [21][22][23]. Coating TiO2 NPs with PEG is a common strategy to elongate their circulation time [24,25]. On top of that, active targeting with RGD ligand could improve the selectivity of photodynamic therapy towards tumor cells.
In this work, we aim at elucidating, by means of fully atomistic molecular dynamics simulations, the mechanism of how the density of cRGDyK (from now on cRGD) conjugation of PEGylated TiO2 nanoparticles affects the ligand's conformation and presentation in an aqueous environment and, consequently, the ligand's targeting ability.
First, in Section 2, we present the computational methods and the details of the analysis performed in this work. Then, in Section 3.1, we present the separated studies of the single components of the nanodevice (cRGD, PEG chain, PEGylated-NP, etc.), which are useful not only to understand the actors involved but also to assess the reliability of the chosen setup and models, through comparison against experimental data and previous studies. Section 3.2 is devoted to the investigation of the whole targeting nanodevice: first, we examine the effect of cRGD density on the cRGD spatial distribution, also through a systematic comparison against available experimental data, and, secondly, we shed light on the role played by the cRGD density on the nanodevice's stabilization and diffusion in water. By analyzing the total interaction energy and its electrostatic and vdW components, we discuss the driving forces controlling the extent of exposure of cRGD ligands in the nanosystem and elucidate the underlying molecular mechanisms. Finally, we compare the resemblance of the various cRGD structures found in our simulations, at different ligand densities, to the X-ray atomic structure of a similar ligand in complex with a αvβ3 integrin for a factbased assessment of the potential targeting ability. In the Conclusions section, we wrap up the paper addressing the implications of ligand density in the development of new nano-delivery systems for biomedical applications.

Single cRGD, PEG500, and PEG500-cRGD in water
The Ligand Reader & Modeler module of CHARMM-GUI [26,27] was employed to generate the three-dimensional structures of PEG500 chains and their conjugation with cRGD ligands was done using the GaussView program [28]. Bonded and non-bonded parameters for cRGD, PEG500 and PEG500-cRGD were assigned by analogy using the CHARMM36m [29][30][31] and CGenFF [32][33][34] FFs. The partial atomic charges of cRGD and PEG500 are reported in the Supporting Material (Fig. S1, Table S1, Fig. S2, Table  S2). The cRGD, PEG500, and PEG500-cRGD molecules were solvated by cubic boxes (56x56x56 Å 3 , 60x60x60 Å 3 , and 76x76x76 Å 3 , respectively) filled up with mTIP3P [35][36][37] water molecules at the experimental density of about 0.99 g/cm 3 using the CHARMM-GUI Solvation Builder [26,38]. One Clcounter-ion was added to the systems containing the cRGD molecule to maintain the system's overall electroneutrality due to the positive net charge acquired by this ligand at physiological pH conditions. First, a minimization phase with 10000 steps using the conjugate gradient algorithm was performed, followed by a subsequent equilibration phase carried out in an NPT ensemble up to 50 ns and extended for more 50 ns of MD production for each system. The long-range solver Particle Mesh Ewald (PME) scheme [39] handled the electrostatic interactions with a cut-off distance of 12.0 Å in the real space with forces smoothly switched to zero after 10.0 Å. A Langevin thermostat with a damping coefficient at 1.0 ps -1 and a semi-isotropic Nosé-Hoover Langevin-piston scheme with a piston period of 50 fs and a piston decay of 25 fs were coupled to the systems to keep the ensemble temperature and pressure constant at 303.15 K and 1 atm, respectively. Newton's equations of motion were solved in time using a 2 fs time step by the impulse-based Verlet-I/r-RESPA integrator [40]. A combination of SHAKE/RATTLE [41,42] and SETTLE [43] algorithms were employed to impose holonomic constraints on all covalent bonds involving hydrogen atoms. To validate the CHARMM-FF assignment of topologies and parameters for cRGD, PEG500 and PEG500-cRGD, we carried out, separately, MD simulations for each of these molecules using the NAMD2.13 code [44].

PEGylated and cRGD-conjugated PEGylated TiO2 nanoparticle in water
To keep the PEGylated TiO2 NP model consistent with our previous work [45], we used a bare 2.2 nm TiO2 NP [46,47] as the substrate for the PEGylated TiO2 NP model.
In short, we grafted the bare 2.2 nm TiO2 NP with 50 methyl-terminated PEG500 chains, whose -OH terminal groups bind to 4-coordinate Ti atoms, which are the most reactive sites, and 5-coordinate Ti atoms on the TiO2 NP surface. The resulting grafting density (σ) is equal to 225.2 chains Å -2 , which was previously demonstrated to be near the transition from mushroom to brush polymer conformation [45]. Further details on the protocol employed to the TiO2 PEGylation can be found elsewhere [45]. The above-cited PEGylated TiO2 NP model was then conjugated with different cRGD ligand densities. Three cRGD ligand densities were studied: a lower one, with 5 out of 50 PEG chains conjugated with cRGD and corresponding to a degree of conjugation of 10 %; an intermediate one, with 10 out of 50 PEG chains conjugated with cRGD ligands reaching a degree of conjugation of 20 %; a higher one, with 25 out of 50 PEG chains conjugated with cRGD ligands reaching a degree of conjugation of 50 %. Here, the degree of conjugation is defined as the percentage of cRGD-PEG500 chains out of the total PEG500 chains, whether conjugated or not with cRGD. PACKMOL [48] was employed to center the systems in cubic boxes (100x100x100 Å 3 for the PEGylated TiO2 NP and for the cRGD-conjugated PEGylated TiO2 NPs with a degree of conjugation of 10% and 20% and 120x120x120 Å 3 for the NP with a degree of conjugation of 50%) which were filled up with mTIP3P [35][36][37] water molecules at the experimental density of about 0.99 g/cm 3 . The CGenFF force field [32][33][34] was used to assign the partial atomic charges and LJ parameters for the PEG500 and PEG500-cRGD chains. The partial atomic charges and Lennard-Jones (12,6) parameters for the TiO2 NP were assigned according to the coordination number of titanium and oxygen atoms using an optimized version of the original Matsui-Akaogi FF [49] further refined by Brandt et al. [50] This FF has been tested and validated by one of us for a TiO2 NP model tethered with small organic molecules [51]. Lorentz-Berthelot rules were used to obtain the cross-term parameters for LJ interactions between unlike atoms. The long-range solver Particle-Particle Particle-Mesh (PPPM) handled the electrostatic interactions with a real-space cut-off of 10 Å and a threshold of 10 -5 for the error tolerance in forces [52]. The Lennard-Jones (12,6) interactions were truncated with a 10 Å cut-off with a switching function applied beyond 8 Å. A minimization phase was carried out to minimize the total potential energy of the system, followed by an equilibration phase 200 ns long, keeping the TiO2 NP core fixed at its QM-optimized geometry. In the MD production phase, the TiO2 core was treated as an independent rigid body able to rotate and move in the system, as done in a previous work by some of us [53]. This approach holds the DFTB-optimized geometry and avoids any misshaping of the TiO2 core during the MD simulation. The MD production phase explored 100 ns of the phase space in the NVT ensemble at 303.15 K held constant by the Nosé-Hoover thermostat with a damping coefficient at 0.1 ps -1 . Only the last 50 ns of production phase were considered for the analysis.
Non-bonded interactions were treated according to the CHARMM potential energy functions [31,54,55], with a force switching distance of 10 Å and a real-space cut-off distance of 12 Å under Periodic Boundary Conditions. Newton's equations of motion were integrated in time using the Velocity-Verlet integrator with a timestep of 2.0 fs.
The SHAKE algorithm [41] was used to impose holonomic constraints on all covalent bonds including hydrogen atoms. All MD simulations including the nanoparticle were carried out using the CHARMM implementation in the LAMMPS code (version 29 Oct 2020, http://www.lammps.org) [56].

Cluster Analysis
To identify the ensemble of relevant conformations acquired by cRGD, either as a single molecule free to diffuse in water or conjugated to the PEG500 chains, we carried out a cluster analysis using the GROMOS method [57]. To determine the cluster membership, we based the analysis only on heavy atoms of cRGD with a cutoff for RMSD of 2.0 Å.

Hydrogen-bonding
The hydrogen-bond formation through the MD trajectory was tracked using the Hydrogen Bonds tool implemented in VMD [58]. We counted as an H-bond formation when the following geometrical criteria were satisfied: (1) the distance between the Hdonor and the H-acceptor heavy atoms was less than 3.5 Å; (2) the angle formed between the H-donor and H-acceptor heavy atoms, having the hydrogen atom as the vertex, was less than 30 °.

Non-bonding interaction energy
To evaluate the relevant aspects of the intermolecular interactions established by the nanodevice components and their molecular surrounding during the MD simulation, we utilized the USER-TALLY package implemented in LAMMPS. The vdW and electrostatic components of intermolecular interactions are reported.

Structural parameters of polymers
To validate the CHARMM-based PEGylated TiO2 NP model against existing literature data and further characterize the polymeric chain behavior in water, we analyzed some polymer structural descriptors, as detailed in the following.
To assess the average length of PEG chains, we measured the end-to-end distance (⟨h 2 ⟩ 1/2 ), defined as the distance between the first and the last heavy atom of a polymeric chain. To estimate the approximated spherical volume occupied by the PEG chains tethered to the TiO2 NP surface, we used the radius of gyration (Rg) analysis defined as follows: where 9 : is the position of each : th heavy atom of the PEG500 chain, 9 ;<=> is the position of the center of mass of the chain and ( is the total number of heavy atoms.

Self-diffusion coefficient
To evaluate the diffusional properties of cRGD molecules, either free to diffuse in water or conjugated to PEG500 chains grafted on the TiO2 NP surface, we estimated the selfdiffusion coefficient based on the Einstein relation that relates the MSD to the selfdiffusion coefficient D as follows: where n is the desired dimensionality and D is estimated by fitting the MSD at intervals in which a linear dependency with time t is observed.

Building up the targeting nanodevice model
In the subsequent sub-sections, we will provide all the relevant details for each modeling step followed in the building-up process of the cRGD-conjugated PEGylated

cRGD model
In the following, we discuss the relevant chemical aspects to be considered for the accurate modeling of cRGD and we validate the CHARMM-based cRGD model through an extensive conformational and structural comparison against available experimental data [59]. In particular, we underline three major features of the cRGD molecule, which could be underestimated at a first sight: Even when the starting-point structures of cRGD contained Gly-Asp and Lys-Arg cis bonds, we observe that, after the annealing cycle, the most stable cRGD conformation is all-trans. Thus, the latter conformation is adopted throughout this work.
To assess the structural convergence and the most recurrent conformational states acquired by cRGD in water, we analyze the RMSD distribution for the whole cRGD structure and its amino acid side-chains (Fig. 2b). Fig. 2a illustrates cRGD structure, the definition for the Ψ (N-Cα-C-N) and Φ (C-N-Cα-C) torsional angles in cRGD, and the color code assigned for its side-chains (the same color codes are consistently adopted throughout the text). Fig. 2c shows the Ramachandran plots for the backbone torsional angles of cRGD analyzed over the MD production phase.
Upon analysis of Fig. 2b, we see that within the converged RMSD region (last 50 ns of the MD simulation), i.e. the production phase, the cyclic portion of cRGD (green) alternates between two main structural conformations, whereas the most mobile sidechains are those of Arg and Lys. Cilengitide (CGT), which is a cyclic peptide (Arg-Gly-Asp-D-Phe-NMe-Val), whose chemical composition resembles the one of cRGD studied in this work and thus provides a fair reference. Black crosses and black diamonds refer to the crystal structure of CGT in complex with the αvβ3 integrin (PDB code: 1L5G [59]) and to MD predictions made by Marelli et al. [62] for CGT in water, respectively. The comparative analysis shows that the single cRGD ligand in the bulk-water phase is able to sample the backbone torsional dynamics in close agreement with the experimental crystal data. In addition, we also identify the sampling of a second dihedral angle distribution for Arg and Gly, which is not identified in the reference X-ray structure [59], but it is in the MD reference [62] in the case of Arg.
To identify the most populated conformations sampled during the MD production phase by the cRGD sidechains, we carry out the cluster analysis based on the RMSD of all heavy atoms, as detailed in Section 2.2.1. We find 8 clusters and report the cluster occupancy of the first three most populated clusters in Fig. 3a, since the other 5 clusters represent less than 1% of the conformations adopted during the production phase.  Previous experimental and computational findings [63][64][65][66] have proven a direct link between the biological activity of RGD peptides and some of their intrinsic structural parameters. Three criteria were proposed for the assessment of biological activity of RGD-containing peptides [63,64] and they are illustrated in Fig. 3c: (i) the distance d1 between the positively (Asp) and negatively (Arg) charged groups, i.e.
between the C ζ Arg and the C γ Asp; (ii) the distance d2 between the C β of Arg and Asp; (iii) the angle θ formed between the C β Arg and the C β Asp with the C α Arg as the vertex. Overall, we find a fair agreement for d1 and d2 between our MD predictions (cRGD structure in Cluster 1, 2 and 3) and those found for CGT in the crystal structure [59] and by DFT calculations [66]. On the contrary, we find a considerable deviation between the cRGD structure in Clusters 1, 2 and 3 and the CGT structure in complex with αvβ3 integrin for θ, as seen in Table 1.
These results above confirm the reliability of the CHARMM-based cRGD model to adequately sample the backbone conformational transitions and relevant geometrical parameters that may drive a successful binding to αvβ3 integrin. In Section 3.2.3, the discussion is resumed to investigate how the cRGD attachment impacts the above-cited geometrical criteria for a pre-evaluation of its bioactivity.

PEG500 chain model
Herein, we will check the capability of the CHARMM-FF to accurately describe relevant experimental polymeric parameters of PEG500 in an aqueous solution against previous predictions using the AMBER-FF [45,67], which were extensively validated against experimental data [27,28].
The polymeric chain model considered in this study contains 11 monomeric units with an overall molecular weight of 500 Da, thus PEG500, and it is terminated with a hydroxyl group on the one side, which will be used to graft the chain to the NP surface, and a methyl group on the other side. To assess the correctness of the empirical FF parameters for the classical description of PEG500 chains in water, we estimate two polymeric properties that are often quantified experimentally, namely end-to-end distance (⟨h 2 ⟩ 1/2 ) and radius of gyration (Rg), comparing them against available theoretical and experimental literature data [45,[67][68][69].
On the basis of the measurement of ⟨h 2 ⟩ 1/2 distance, the PEG500 molecule described by the CHARMM-FF is, on average, slightly more extended (16.1 (±5.5) Å) in water than the one simulated using AMBER FF (14.3 (±4.0) Å) [67]. Moreover, we find a fair agreement for the Rg between the CHARMM-based PEG500 model and available numerical MD simulation [45] and experimental measurements [67][68][69] for an identical PEG500 chain (Table 2). In addition, we also validate the PEG500 model in the bulk-water phase for the CHARMM-FF through a comparative analysis of the polymer dihedral distribution against the same MD descriptor estimated using the AMBER-FF. We observe two maximum peaks for the O-CH2-CH2-O dihedral angle: one positioned at -69.9 ° and a second one at +72.1 ° (see Fig. S4 in Supporting Material). These estimations using the CHARMM-FF agree with those reported in Ref. [45] for the same system and in Ref. [70] and [71]. Based on the results above, one can infer that CHARMM-FF is a reasonable FF choice to describe polymeric linkers based on polyethylene glycol units under explicit water solvation. In the next section, we will proceed with the validation of a PEG500 molecule conjugated with cRGD.

PEG500-cRGD chain model
After implementing and validating the cRGD and PEG500 components separately in Sections 3.1.1 and 3.1.2, we now conjugate both cRGD and PEG500 models into a single molecule through the formation of an amide bond, as seen in Fig. 4a, mimicking the experimental structure found in [25]. The resulting zwitterionic cRGD-conjugated PEG500 chain model (PEG500-cRGD) will be used here to functionalize the TiO2 NP following the protocol reported in Section 2.1.2.
To find out whether the conjugation of cRGD to the PEG500 chain affects the ligand conformation or not, we repeat the same conformational analysis as previously done for a single cRGD in Section 3.1.1. Fig. 4b displays the Ramachandran plot analysis for Φ and Ψ torsional angles of a cRGD ligand conjugated to the PEG500 chain. Upon analysis of the Ramachandran plot in Fig. 4b, we can observe a resemblance between the conformational states acquired by the cRGD ligand, whether conjugated to PEG500 chains or not (Fig. 4b vs. Fig. 2c). Thus, we can infer that cRGD conjugation to PEG500 linkers has no impact on the ligand's torsional dynamics. Overall, we confirm the reliability of the CHARMM-FF parameters to describe the torsional dynamics of cRGD conjugated to PEG500 chains in explicit bulkwater solvation and their suitability to compose more complex systems, as will be done in Section 3.1.6.

TiO2 nanoparticle model
The spherical TiO2 NP model (diameter 2.2nm, (TiO2)223·10H2O) used in this study has been developed by some of us in previous works using high-level quantum mechanical (QM) methods [46,47]. This model was proven to be chemically stable upon water adsorption [72] and solvation using QM and QM/MM approaches [73].
Moreover, the validity and accuracy of the chosen force field parameters [50] for the description of bioinorganic nanohybrids was also assessed in previous works by some of us against QM ad QM/MM calculations [45,51]. number of H-bonds descriptors, whose definitions can be found elsewhere [45].

PEGylated TiO2 nanoparticle model
In Table 3, we compile the structural parameter data estimated using the CHARMM-based PEGylated TiO2 NP model and the available reference data. For all the polymeric parameters analyzed here, we observe a fair agreement between our MD results and the simulation predictions in the literature.

N. of Hbonds
This work 225. Three different cRGD densities are modeled: 1) 10% of cRGD-conjugated PEG500 chains (from now on NP-PEG-cRGD10 model in Fig. 5a); 2) 20% of cRGD-conjugated PEG500 chains (from now on NP-PEG-cRGD20 model in Fig. 5b); 3) 50% of cRGDconjugated PEG500 chains (from now on NP-PEG-cRGD50 model in Fig. 5c). Likewise, we have decided to use short PEG500 chains to avoid the well-known effects of ligand entangling and clocking due to long polymeric chains as linkers [10]. Recent experimental studies have also found that shorter PEG chains with low molecular weight (i.e., PEG350) exhibit similar circulation half-life as the longer and heavier ones do (i.e., PEG2000) [74] and may enhance the cellular uptake in breast cancer and myeloma cells [75].

Analyzing the nanodevice: impact of cRGD ligands density
In the following subsections, we investigate the role of ligand density on some biophysical properties of relevance to the cRGD ability to target the αvβ3 integrin. In

On the cRGD presentation
Cell adhesion is of utmost importance to improve the target-specific delivery of cRGDconjugated nanomaterials, which directly depends on the degree of targeting-ligands exposure and their orientation out of the "stealth" polymeric coating layer towards the biological milieu. Experimental findings have suggested that the spatial organization of ligands plays a crucial role in the T-cell immune response [76]. The authors discovered that ligand's presentation and disposition, at optimal nanometric distances from each other on DNA flat sheets, have the ability to suppress T-cell signaling in vivo [76]. Moreover, a recent combined experimental and MD simulation study has shown the implications of cRGD presentation and orientation on the targeting efficiency of cRGD-conjugated PEGylated Au NPs towards αvβ3 integrin. They found that the enlargement of ligand-conjugated chains by a short molecular spacer is sufficient to enhance the cRGD moiety presentation out of the polymeric stealth layer, which is corroborated by a substantial increase in the water accessibility towards the cRGD surface [13].
Here, we investigate how different degrees of ligand conjugation impact the cRGD presentation and modulate the ligand's molecular environment in the system.
To quantify the extent of cRGD exposure from the stealth PEG region towards the bulk-water phase, we first compute the number density profiles for the relevant moieties of the NP-PEG-cRGD10, NP-PEG-cRGD20 and NP-PEG-cRGD50 nanosystems. Fig. 6 shows the number density profiles for water molecules, PEG500  Interestingly, we identify a direct correlation between the degree of PEG500 chains conjugated with cRGD ligands and the distribution of cRGD ligands in the system (Fig. 6). At 10% of cRGD conjugation, the cRGD density has a log-normal distribution with the longest tail towards the bulk-water phase among all systems. The latter observation implies that cRGD ligands can reach farther regions beyond the PEG500 layer with considerable exposure to the bulk-water phase. By increasing the cRGD conjugation density from 10% to 20%, we identify an increased cRGD density in the outer region of the PEG500 shell, confirmed by a slight shift in the density maximum peak towards the bulk-water phase in the latter system. Further raising of the cRGD density from 20% to 50% shows no improvement in the exposure of cRGD towards the bulk-water phase. On the contrary, we see an increasing number of cRGD ligands located in the inner region of the PEG500 layer, with the resulting number density profile mostly resembling that estimated at 20% of cRGD density.
Although the number density profiles can accurately predict the average location of cRGD in the system, we still have no information on the orientation and entangling tendencies of cRGD at different ligand densities. To this end, we analyze the tilt angle acquired by cRGD against their distance from the NP center (Fig. 7). For the angular analysis, two vectors are considered: one vector is defined from the NP center to the C α Lys and the second vector from the C α Lys (vertex) to the N atom of Asp in the cRGD cycle (Fig. 7a).  Fig. 7 reveals that the ligand density modulation affects the degree of cRGD exposure and orientation. At 10% of cRGD conjugation, we find that stand-up configurations (tilt angles higher than 120°) are denser populated at farther distances from the NP surface (Fig. 7b) compared to the systems at higher ligand densities ( Fig.   7d and f), while at 20% of cRGD conjugation, we identify a clear tendency of cRGD in populating a straighter stand-up configuration at shorter distances, as seen in Fig. 7d.
Further increase in the cRGD density leads to a substantial modification in the interplay between the cRGD molecules' orientation and their distance from the NP surface. Fig.   7f reveals multiple hotspots at closed-up configurations (tilt angles lower than 120°) being densely populated at 50% of cRGD conjugation, which does not occur at lower cRGD densities (Fig. 7b and 7d).
To unveil the molecular reasons behind this cRGD behavior noticed at different densities, we conduct a volumetric analysis to quantify the surface area of contact between cRGD and other nanosystem components (PEG500 chains, solvent molecules, and other cRGD ligands) with the aid of Voronoi tessellation [77]. In this way, we can quantify the surface area exposure of ligands to distinct parts of the system and find their preferential molecular neighborhood along the MD trajectory. Table 4 shows the surface area of contact between the cRGD ligands at 10%, 20% and 50% of cRGD conjugation normalized by the total number of ligands.  Table 4, we see that the molecular surface of cRGD ligands is mainly accessed by water molecules, although the solvent surrounding the ligand surface decreases with the increase of cRGD conjugation. The largest decrease of water contact with the cRGD surface occurs going from 20% to 50% of cRGD conjugation, corresponding to a reduction of about 100 Å 2 in the water accessibility towards the ligand surface. We may notice that, in opposite trend but as a compensating effect, the surface area of contact between cRGD ligands themselves gets larger to a similar extent (of about 120 Å 2 ).
At first glance, the volumetric data (Table 4) seem to contradict the number density data reported in Fig. 6, where a higher cRGD density leads to an enhanced number of cRGD ligands populating the outer region of PEG500. However, further volumetric analysis reveals a substantial enlargement in the surface area of contact between cRGD-cRGD ligands at higher cRGD densities conjugation (Table 4), reducing the solvent accessibility to their surface. Hence, the observation that cRGD ligands populate a more stand-up conformation at higher ligand density is likely linked to undesirable ligand-ligand interactions than their exposure to the bulk-water phase, which is rather evident when the cRGD density is the highest among all.
From the above results, we can infer that cRGD ligands are mainly located at the PEG500-water interface, and their exposure to the bulk-water phase in an upward mode is preferred at the lowest ligand density. In addition, we also notice that the cRGD density directly impacts the ligand-ligand interactions, which may impair the water accessibility towards the ligand molecular surface.

On the nanodevice stability and diffusion in water
Recent experimental investigations have hypothesized a direct link between the targeting ligand solubility and their availability on the nanosystem surface. A recent experimental investigation has verified that targeting ligands with high water solubility tend to escape the hydrophobic polymeric region of PEGylated nanoparticles [11]. In the same vein, experimental findings have speculated that enhanced water solubility of polymeric chains may increase the RGD ligand diffusion and, therefore, impair the binding affinity of RGD to integrins [60,78]. A recent experimental work has also unveiled the role of ligand diffusion on specific α5β1 and αvβ3 integrin activation and signaling. These findings suggested an enhanced cell adhesion and differentiation by αvβ3 integrin on covalently-immobilized cRGD ligands compared to higher ligand diffusion regimes [79].
The interplay between the nanosystem's solubility and the ligand diffusion may have significant implications on the ability of cRGD to bind integrins, and therefore, we investigate the nanodevices' stability and the ligand diffusion upon increasing of cRGD density on its surface. To do so, we first estimate the intermolecular interaction energy between the NP-PEG-cRGD10, NP-PEG-cRGD20 and NP-PEG-cRGD50 nanosystems and water molecules. After that, we decompose the total energy of interaction into its respective electrostatic and vdW components, as detailed in Section 2.2.3. Table 5 compiles the nonbonded interaction energy and its vdW and electrostatic contributions between each system and water and between different nanodevice components for the NP-PEG-cRGD10, NP-PEG-cRGD20 and NP-PEG-cRGD50 systems.

NP-PEG-cRGD50
System -water -9.4·10 3 (±2·10 2 ) -1.87·10 3 (±6·10) cRGD -water -190 (±5) Analysis of Table 5 reveals that the main energetic stabilization comes from intermolecular interactions between the nanosystem and solvent. Interestingly, we identify that increasing the cRGD content considerably intensifies the energetic stabilization in water compared to the reference model (NP-PEG), being the electrostatic contribution the most relevant to it. On the other hand, the role of vdW interactions seems to be limited to a stabilization of the non-electrostatic forces between the nanosystem components, particularly between cRGD and PEG500 chains.
Noteworthy, the decreasing intensity of the electrostatic interactions between cRGD and water in Table 5 well agrees with the volumetric analysis in the previous section (Table 4), where we noticed a reduction of the surface area of contact between cRGD and the solvent molecules at increasing cRGD ligand density.
As seen in Table 5, the energetic stabilization coming from intermolecular interactions between the nanodevice's components are way less intense than those with the solvent. However, we identify that ligand-ligand interactions become relevant at 50% of cRGD conjugation mainly strengthened by electrostatic interactions, while they remain negligible at lower ligand density. Further discussion on the implications of ligand-ligand interactions on the optimal ligand presentation for integrin binding can be found in the previous section.
Finally, we investigate the effect of anchoring the cRGD ligands to the polymeric chains on their diffusional properties, which we expect to be reduced. To this end, we calculate the diffusion coefficient of cRGD either free to diffuse in bulk-water phase (reference) or bonded to a PEG500 chain in the NP-PEG-cRGD10, NP-PEG-cRGD20, and NP-PEG-cRGD50 nanosystems (Fig. 8).  Table S3. Comparing the self-diffusion coefficients in Fig. 8, we identify a decrease of one order of magnitude in the self-diffusion coefficient of cRGD molecules due to their conjugation to PEG500 chains in the NP-PEG-cRGD10 system. Increasing the cRGD conjugation to 20% and 50% further slows down the cRGD self-diffusion coefficients, reaching values about two orders of magnitude smaller than those observed for the free cRGD molecule diffusing in water. This latter remark is intimately associated with the increasing of the cRGD-cRGD interactions identified in the previous volumetric and energetic analyses (Tables 4 and 5). Hence, we can infer that not only the ligand restraints due to their conjugation to polymeric chains but also the ligand density parameter affects the cRGD diffusion.
Based on the above results, we can infer that increasing the cRGD molecules' density on the nanosystem surface makes them more stable in water than their nonconjugated counterpart. Also, the physical restrains imposed on the cRGD ligands upon their conjugation to the polymeric chains grafted on the NP surface has a stout effect on their diffusional properties, which might affect their biological activity in vivo.

On the optimal conformation of cRGD for integrin targeting
The biological activity of cRGD-containing nanodevices is intimately related to the spatial disposition of cRGD ligands and their moieties prior to and after the binding event [60]. Likewise, a successful binding of cRGD to αvβ3 integrin depends on the proper conformation acquired by this ligand in solution, and it has been demonstrated that molecular simulations are well-suited to reproduce the pre-organized structure of cyclic RGD peptides [80]. For instance, recent findings have shown that slight changes in the molecular structure of cRGD-based peptides considerably impact their conformational states leading to substantial implications on the binding affinity to αvβ3 integrin [81].
To investigate how the ensemble of conformations acquired by the cRGD ligands conjugated to the PEGylated TiO2-based nanodevice is affected by different ligand densities, we analyze the Φ and ψ dihedral angles involved in the torsional transitions of the cRGD cycle, as defined in Fig. 2a. Fig. 9 displays the Ramachandran plot for the Φ and ψ dihedral angles of cRGD during the MD production phase. To the best of our knowledge, no crystallographic data is available for the cRGD-αvβ3 integrin complex to date. Thereby, we use as reference the X-ray structure values of Φ and ψ dihedral angles of CGT bounded to αvβ3 [59], being CGT (Arg-Gly-Asp-D-Phe-NMe-Val) a cyclic peptide based on the same RGD sequence and that possesses a similar structure to the one of the cyclic peptide used in this work. As seen in Fig. 9a, b and c, we observe a fair agreement between the Φ and ψ dihedral angles sampled for Arg, Gly, Asp, and Lys side chains and their equivalent in the Xray αvβ3 integrin-CGT complex. Besides, for the D-tyrosine side chain in cRGD, we notice that the Φ and ψ values found in the crystal structure for D-Phe are only marginally sampled during the MD simulation.
Next, we analyze the main cRGD conformations sampled during the production phase by means of cluster analysis (details in Section 2.2.1). We perform the analysis on the most exposed cRGD of each system, that is the one with the longest average distance from the center of the NP. Interestingly, the clusters of cRGD conformation in all the three nanosystems resemble the ones found for the free cRGD in water (Section 3.1.1). The molecular structures for the most relevant conformations are reported in Fig. S5 in Supporting Material.
As previously stated in Section 3.1.1, several experimental works suggest the evaluation of a particular set of geometrical criteria (see Fig. 3c for further details) to assess the potential biological activity of RGD-containing peptides towards integrins [63,64]. Here, we retake the same analysis done for the single cRGD molecule in water (Table 1), further evaluating this set of geometrical parameters for the cRGDconjugated ligands in the NP-PEG-cRGD10 NP-PEG-cRGD20 and NP-PEG-cRGD50 nanodevices. Table 6 compiles the structural parameters evaluated for CGT X-ray structure, the single cRGD molecule solvated in water, and the cRGD ligands. Analysis of Table 6 reveals that the C ζ Arg-C γ Asp d1 and C β Arg-C β Asp d2 distance values are in fair agreement between the X-ray crystal structure and the cRGD ligands conjugated on the PEGylated TiO2 NP surface. Thereafter, we see no significant changes on these two geometrical parameters for cRGD ligands either as a single molecule free to diffuse in water or conjugated to PEG500 polymeric linkers. However, we observe a slight deviation of the C β Arg-C α Gly-C β Asp θ angle of cRGD compared to the CGT crystal structure in complex with αvβ3. This deviation is likely linked to the dynamic nature of cRGD in explicit water solvation in contrast to the molecular confinement of CGT into the binding pocket of αvβ3 integrin.
Although the biological activity requires further investigation, we can infer from the analysis above that the cRGD conjugation to PEG linkers has little effect on the conformational states acquired by this ligand. Hence, this result suggests that cRGD is expected to retain its biological activity upon PEGylated conjugation to TiO2 NPs.

Conclusions
Enhancing active targeting strategies by rigging PEGylated nanodevices with highaffinity cRGD ligands has been proven as an efficient way to reduce off-targeting and improve the uptake of NPs into overexpressing-integrin tumor cells. On the one hand, much has been discussed in the experimental literature on the effects of ligand density and linker lengths on the efficacy of targeted polymeric NPs. On the other hand, little is known about the molecular mechanisms that drive the ligand presentation, orientation, and conformation when conjugated to flexible polymeric chains linking them to NPs.
In this study, we evaluate the effect of different degrees of cRGD conjugation to PEG chains grafted on highly curved TiO2 NPs on the ligand's availability, conformation, and presentation, which is believed to play a pivotal role in the translation of these nanodevices as targeting agents in vivo. Our findings confirm that cRGD conjugation to polymeric chains slows down the ligand diffusion and has negligible effect on the cRGD ability to sample optimal conformations, meeting two important requirements for an effective binding activity towards αvβ3 integrin. However, we found that the chosen cRGD density is a critical parameter and profoundly impacts the availability of cRGD targeting ligands on the nanodevice surface. Although the cRGD exposure out of the polymeric layer is little affected by the increasing number of cRGD-conjugated polymeric chains in the nanodevice, further volumetric analysis reveals a considerable increase of ligand-ligand interactions impairing their proper solvation, which might affect the optimal presentation of cRGD ligands for a successful binding to the αvβ3 integrin receptor.
These findings represent a starting point for further investigations on cRGDconjugated nanosystems for integrin-targeting and provide helpful insights for the design of effective targeted nanoparticles.

Conflicts of interest
The authors declare no conflict of interest.  Table S1. Partial atomic charges for c(RGDyK) atoms: atom names, CHARMM atom types, partial atomic charges.

Simulated annealing of cRGD
The same simulation protocol described in the main text was used both for the all-trans conformation Fig. S3 (a) and for the cis-trans one in Fig. S3 (b), where G-D and K-R bonds are cis and the other ones trans. After a production phase of 100 ns was run, the two c(RGDyK) conformations were subjected to the following simulated annealing procedure, which allows for a minimum global search: first, the temperature is heated up to 700 K at a rate of 2 K/step; second, the system is let evolve in time for 50 ns; then it is cooled down to the original temperature at a rate of 2 K/step and finally allowed to evolve for another 50 ns.
For both the all-trans conformation and the cis-trans one the temperature of 700 K allowed us to observe amide bond cis-trans isomerization. In particular, Gly-Asp amide bond isomerized more often than other ones, but eventually the structures ended up in an alltrans conformation during the cooling down phase. Thus, the all-trans conformation was adopted for this work.  Table S3. Self-diffusion coefficients (D) of single cRGD in water and bonded to PEGylated TiO2 nanodevices at 10%, 20% and 50% of cRGD conjugation with their standard deviation.