Molecular Dynamics for the Optimal Design of Functionalized Nanodevices to Target Folate Receptors on Tumor Cells

Atomistic details on the mechanism of targeting activity by biomedical nanodevices of specific receptors are still scarce in the literature, where mostly ligand/receptor pairs are modeled. Here, we use atomistic molecular dynamics (MD) simulations, free energy calculations, and machine learning approaches on the case study of spherical TiO2 nanoparticles (NPs) functionalized with folic acid (FA) as the targeting ligand of the folate receptor (FR). We consider different FA densities on the surface and different anchoring approaches, i.e., direct covalent bonding of FA γ-carboxylate or through polyethylene glycol spacers. By molecular docking, we first identify the lowest energy conformation of one FA inside the FR binding pocket from the X-ray crystal structure, which becomes the starting point of classical MD simulations in a realistic physiological environment. We estimate the binding free energy to be compared with the existing experimental data. Then, we increase complexity and go from the isolated FA to a nanosystem decorated with several FAs. Within the simulation time framework, we confirm the stability of the ligand–receptor interaction, even in the presence of the NP (with or without a spacer), and no significant modification of the protein secondary structure is observed. Our study highlights the crucial role played by the spacer, FA protonation state, and density, which are parameters that can be controlled during the nanodevice preparation step.

Table S1.Breakdown of LIG binding energies in the FR binding pocket from docking calculations in GBIS implicit water solvent for the FR/LIG 0,1-,2-systems.

H-bonds number
Interaction energy (kcal mol -1 ) VdW Electrostatic FR/LIG 0     We observe that for the PEGylated systems (Figure S6a), the PEG chains bend and the NP approaches the protein.Therefore, both TiO2-LIG b0 and TiO2-FR distances decrease until they stabilize after about 100 ns from the beginning of the simulation.In the case of the FR/TiO2-48-FAs b0 - system (Figure S6b), the two distances are very much constant along the MD run, because LIG b0 , as well as all the other FAs b0 , is covalently bonded to the NP, and since it is a short rigid molecule, it has reduced mobility, therefore its COM is found in a narrow range of distances from the S9 NP center.Hence, it is trivial to understand that convergence is reached in a shorter time for this system.

Replicas of MD simulation of FR/TiO2-PEG-20-FAs b1-system
First, we report the results of the analysis on the replica of the 200 ns MD simulation for the FR/TiO2-PEG-20-FAs b1-system, where the same equilibrated starting structure of the original simulation was used, but different initial velocities were assigned to the atoms.
In Figure S14 we report the time evolution along the MD simulation of the distances between the center of mass of the TiO2 NP and that either of the LIG b1-(TiO2-LIG b1-) or of the FR (TiO2-FR).

S23
In conclusion, replicating the MD run of the FR/TiO2-PEG-20-FAs b1-system with different initial velocities has demonstrated that the observations on the dynamics of LIG b1-inside the FR binding pocket, on the interaction of the functionalized NP with the protein and on the FR secondary structure are consistent, making our conclusions more robust.
Then, to further support our conclusions, we also report the results of a second replica of 50 ns, starting from the structure at 150 ns of the original simulation of FR/TiO2-PEG-20-FAs b1-and assigning different initial velocities (Figures S19-S23).The analysis presented in Figures S19-S23

Figure S2 .
Figure S2.Top: PMF profile relative to LIG binding to FR for the FR/LIG 0 system (black line) with the associated error bar (orange).Bottom: umbrella histograms.

Figure S3 .
Figure S3.Top: PMF profile relative to LIG binding to FR for the FR/LIG 1-system (black line) with the associated error bar (orange).Bottom: umbrella histograms.

Figure S4 .
Figure S4.Top: PMF profile relative to LIG binding to FR for the FR/LIG 2-system (black line) with the associated error bar (orange).Bottom: umbrella histograms.

Figure S5 .
Figure S5.(a) Set of ligand and protein atoms (vdW representation) used to calculate intermolecular distances for SOM training.(b) Silhouette profile.The optimal number of clusters (7) was chosen as the one with the highest silhouette score in the range 5-10.

Figure S9 .
Figure S9.Radial distribution function profiles of PEG, FAs b0 , LIG b0 and FR calculated with respect to the Ti atom at the center of the NP, for the FR/TiO2-PEG-10-FAs b0 (a) and FR/TiO2-PEG-20-FAs b0 (b) systems, averaged on the last 50 ns of the 200 ns MD simulation.

Figure S12 .
Figure S12.Radial distribution function profiles of PEG, FAs b1-, LIG b1-and FR calculated with respect to the Ti atom at the center of the NP, for the FR/TiO2-PEG-10-FAs b1-(a) and FR/TiO2-PEG-20-FAs b1-(b) systems, averaged on the last 50 ns of the 200 ns MD simulation.

Figure S15 .
Figure S15.Radial distribution function profiles of PEG, FAs b1-, LIG b1-and FR calculated with respect to the Ti atom at the center of the NP, for the FR/TiO2-PEG-20-FAs b1-system (first replica), averaged on the last 50 ns of the 200 ns MD simulation.

Figure S16 .
Figure S16.Contact surface area between the functionalized NP and the FR along the 200 ns MD simulation for the FR/TiO2-PEG-20-FAs b1-system (first replica).

Figure S17 .
Figure S17.Time evolution of distances between selected atoms of LIG b1-and FR AAs for the 200 ns simulation of the FR/TiO2-PEG-20-FAs b1-system (first replica) in physiological environment.
are consistent with those of Figures 6/S10-13 for the original simulation and Figures S14-18 for the other simulation discussed above and, therefore, further corroborate the conclusions of this work.

Figure S20 .
Figure S20.Radial distribution function profiles of PEG, FAs b1-, LIG b1-and FR calculated with respect to the Ti atom at the center of the NP, for the FR/TiO2-PEG-20-FAs b1-system (second replica), averaged on the last 50 ns of the 200 ns MD simulation.

Figure S21 .
Figure S21.Contact surface area between the functionalized NP and the FR along the 200 ns MD simulations for the FR/TiO2-PEG-20-FAs b1-system (second replica).

Figure S22 .
Figure S22.Time evolution of distances between selected atoms of LIG b1-and FR AAs for the 200 ns simulations of the FR/TiO2-PEG-20-FAs b1-system (second replica) in physiological environment.