Multi-scale modeling of folic acid-functionalized TiO 2 nanoparticles for targeted photodynamic therapy of tumors

Strategies based on the active targeting of tumor cells are emerging as smart and efficient nanomedical procedures. Folic acid (FA) is a vitamin and a well-established tumor targeting agent because of its strong affinity for the folate receptor (FR), which is an overexpressed protein on the cell membranes of the tumor cells. FA can be successfully anchored to several nanocarriers, including inorganic nanoparticles (NPs) based on transition metal oxides. Among them, TiO 2 is extremely interesting because of its excellent photoabsorption and photocatalytic properties, which can be exploited in photodynamic therapy. However, it is not yet clear in which respects direct anchoring of FA to the NP or the use of spacers, based on PEG chains, are different and whether one approach is better than the other. In this work, we combine Quantum Mechanics (QM) and Classical Molecular Dynamics (MD) to design and optimize the FA functionalization on bare and PEGylated TiO 2 models and to study the dynamical behavior of the resulting nanoconjugates in an aqueous environment. We observe that they are chemically stable, even under the effect of increasing temperature (up to 500 K). Using the results from long MD simulations (100 ns) and from free energy calculations, we determine how the density of FA molecules on the TiO 2 NP and of the presence of PEG spacers impact on the actual exposure of the ligands, especially by affecting the extent of FA-FA intermolecular interactions, which are detrimental for the targeting ability of FA towards the folate receptor. This analysis provides a solid and rational basis for experimentalists to define the optimal FA density and the more appropriate mode of anchoring to the carrier, according to the final purpose of the nanoconjugate.

TiO2/FA nanoconjugates have also shown an excellent ability to reach folate receptors in FR(+) cancer cells. For example, Ai et al. 17 report that TiO2/FA nanoparticles have better capability of inducing reactive oxygen species (ROS) formation with higher cell apoptosis in osteosarcoma cancer cells than bare TiO2 NPs, and that is related to the enhanced cellular uptake of these nanoconjugates mediated by FR in the presence of FA. Furthermore, Liang et al. 18 have combined The TiO2 NP model was designed by our group in previous works 31, 32 and consists of a spherical anatase nanoparticle, carved from the crystalline bulk anatase structure and fully relaxed, first at the DFTB level of theory with a simulated annealing procedure, followed then by a  optimization. The stoichiometry of the NP is (TiO2)223•10H2O and it is characterized by an equivalent diameter of 2.2 nm (700 atoms, Fig. 1b). In the simulations presented in this study we started from the DFTB-optimized NP to which we added as many N-formylglutamic acid (FGA, Fig.   1c) molecules as possible (52 or 48 depending on the COOH group used to anchor the surface) to reach the highest coverage admissible by the surface available Ti sites and the steric hindrance between attached molecules. FGA is a smaller molecule but analogous to FA. The FGAfunctionalized models underwent some thermal treatment through MD simulations and were then fully relaxed at the DFTB level of theory. After that, we converted FGA molecules into FA and only performed partial DFTB geometry optimization (for further details see below).
The FGA or FA molecular adsorption energy on the spherical nanoparticle, , is defined as: where + , and are, respectively, the DFTB total energy of the whole system, of the bare NP and of a single FGA or FA molecule in the gas phase and is the number of adsorbed FGA or FA molecules.
The TiO2/PEG model was obtained by our group in previous studies, 33,34 where the nanoparticle was grafted with 50 PEG molecules and fully optimized through MD simulations and atomic relaxation at the DFTB level of theory with the MATORG parameterization set. On this model we added either 10 or 20 FA molecules, as randomly distributed as possible, by a covalent linkage through an ester bond involving the terminal C of PEG and the FA -carboxylic group. The FAfunctionalized PEGylated TiO2 systems were only partially relaxed at the DFTB level of theory.

Molecular Mechanics-Molecular Dynamics (MM-MD)
The Molecular Mechanics-Molecular Dynamics method is a classical approach, where the electrons are not described explicitly, and the total potential energy of the system is calculated as a function of a set of empirical parameters and potentials that constitute a force field (FF). In this work, all the classical MD simulations were carried out employing the LAMMPS (19 Aug 2019 version) opensource code. 35 The TiO2 NP was described by an improved Matsui-Akaogi FF, reparametrized by Brandt and Lyubartsev, 36 while the CGenFF [37][38][39] was chosen for the adsorbed FA molecules and the TIP3P water molecules. We simulated several TiO2/FA systems, with different FA coverages and binding modes to the NP, both in vacuum and in water. The TiO2/FA topologies were generated by means of the Moltemplate 40 package for LAMMPS and the systems were immersed in a 100 X 100 X 100 Å 3 cubic water box, built up with the Packmol 41 software. During all the atomistic MD simulations, we fixed the position of the FA atoms directly bound to the NP at the DFTB geometry and we treated the nanoparticle as a rigid body, as previously done by some of us, 42 free to translate and rotate as a whole, fixing its internal degrees of freedom at the DFTB-optimized geometry through the RIGID package 43  In the case of the TiO2/PEG/FA systems, we fixed the position of the nanoparticle and of the binding hydroxyl group of the PEG chains at the DFTB geometry and we made use of a 135 X 135 X 135 Å 3 water box. Long-range electrostatic interactions were evaluated by the particle-particle particle-mesh (PPPM) solver, using a real-space cutoff of 12 Å. Short-range Lennard-Jones interactions were smoothly truncated with a 12 Å cutoff by means of a switching function applied beyond 10 Å. A hundred minimization steps ensured that no overlaps between the atoms occurred, then a 1 ns equilibration followed and finally the phase space was explored for a total production time of 100 ns.

Simulation analysis
Atomic radial distribution functions (rdf) were computed with the Radial Distribution Function plugin of VMD, 45 considering all the atoms at a given distance r from the reference and falling in a spherical crown with a thickness of 0.1 Å.
The hydrogen bonds analysis was carried out by means of the Hydrogen Bonds tool 45 provided by VMD. The criteria for the H-bonds classification were set as 3.5 Å for the maximum donor-acceptor distance and 30° for the maximum acceptor-donor-hydrogen angle.
Average interatomic distances and radius of gyration were computed exploiting the interdist tool, included in the open-source software LOOS 3.1. 46 In particular, for all the systems we analyzed: • the distance between the N atom of the amine group of the FA pterin and the six-fold Ti6c atom in the center of the NP, d N FA -NP center • the distance between the FA center of mass (com) and the closest Ti atom of the NP surface, d com FA -NP surface For the TiO2/PEG/FA systems, we also evaluated: • the end-to-end distance of the PEG chains, <h 2 > ½ PEG , i.e. the mean distance between the first and the last heavy atom of each PEG chain (oxygen of the -OH head and carbon of the -CH3 tail) • the end-to-end distance of the PEG-FA chains, <h 2 The calculations of the non-bonding interaction energies (vdW and electrostatic) were performed through the USER-TALLY package, implemented in the LAMMPS code. 35 The root-mean-square deviation (RMSD) of atomic positions was estimated using the VMD RMSD Trajectory Tool package. 45 The diffusion coefficients, D, of the TiO2/FA and TiO2/PEG/FA nanosystems were estimated from the Einstein equation: where MSD is the mean-square-deviation of atomic positions, n is the dimensionality of the diffusion and t is the simulation time. In particular, the diffusion coefficients were obtained by fitting of the MSD at intervals of the 100 ns MD trajectory where the dependency with t was linear, with the ordinary-least-squares (OLS) method.  47 The mean force, 〈 〉 , exerted along the transition coordinate, , is related to the free energy first derivative according to: where |J| is the Jacobian determinant to transform the cartesian coordinates to generalized coordinates and 〈 〉 denotes the cumulative mean of the instantaneous force, , in bins wide. To obtain the free energy profiles, we followed a protocol similar to that of a previous work by some of us: 48 the reaction coordinate, defined as the z-component of the distance between the center of mass of the FA molecule and the center of the NP, was discretized into six consecutive windows, each 5 Å wide, at a z-distance from 10 to 40 Å with respect to NP center. To keep the FA molecule at the desired zdistance from the center of the NP, its center of mass was harmonically restrained in both the window's edges and the xy plane perpendicular to the transition coordinate with a force constant of 10 kcal/mol/Å 2 . For every window, a total simulation time of 150 ns was reached, after a preliminary 10 ns equilibration phase. The nanoparticle and the binding hydroxyl group of the PEG molecules were fixed to the DFTB geometry all along the simulation, with the NP principal axis aligned to the cartesian axis. To obtain the LJ free energy profiles, the atomic partial charges of the considered FA molecule were set to zero; in order to keep the whole system neutral, the charge of the two H atoms of the last -CH2 unit of the PEG chain bonded to the same FA molecule were adjusted.  49 and therefore it is the group that more easily would dissociate onto the NP surface. However, experimental works claim that the interaction of FA with the folate receptor is much stronger when the -COOH group is free and available for the ligandreceptor binding. 50 For this reason and for comparative purposes, in this work, we have considered both the carboxylic groups as potential anchoring groups on the nanoparticle. Therefore, we built up two different binding models: model (TiO2/n- FA-), where the FA molecules bind the nanoparticle through the -COOH and model (TiO2/n- FA-), where the anchoring bond involves the -COOH group. In addition, to speed up the QM calculations, we preliminary considered an analogue of the FA molecule, the N-formylglutamic acid (FGA, see Fig. 1c), which is shorter but, like FA, contains the glutamic acid portion used to anchor to the NP.

Results and discussion
In Fig. 2, we compare four FGA binding modes involving the -COOH group and the equatorial four-fold Ti sites (Ti 9 and Ti 11 in Fig. 1b, which were proven to be the most reactive 51,52 ) of the TiO2 NP: undissociated monodentate (UM), dissociated monodentate (DM), bidentate (B) and chelated (C), using the FGA molecular adsorption energy (eqn (2)) as a criterion for the relative stability of these different adsorption modes. For all the dissociative mechanisms, the dissociated proton from FGA was transferred to the nearest two-fold oxygen atom (O2c) of the NP. The results clearly show that the dissociative adsorption modes are always favored with respect to the non-dissociative ones: in particular, is 0.26 eV more negative for DM with respect to UM.
Moreover, the B and the C coordination modes result, in the order, the most advantageous FGA binding modes with an that is, respectively, 0.24 and 0.33 eV more negative than that of DM.
On the basis of the results on FGA in this section, we conclude that, by extension, FA would preferentially have the anchoring carboxylic group dissociated on the TiO2 surface, and that the ligand would likely form chelate coordinative bonds with Ti4c atoms and/or bidentate ones with Ti4c or Ti5c atoms, since these adsorption modes are energetically more stable. We rationalize these results as due to the greater energy stabilization by formation of two Ti-O bonds instead of one.

Increasing the coverage of the nanoconjugates
As a next step, we have proceeded with increasing the density of FGA molecules on the TiO2 NP, starting with the equatorial Ti4c sites, based on our previous works where we demonstrated that these are the most reactive sites of the NP. 51,52 Hence, also considering the results of the previous section, we first saturated all the 12 equatorial Ti4c atoms with 12 FGA molecules in a chelated configuration.   We observe from Table 1  Moreover, from the analysis of data in Table 1, one can notice that the FGA coordination through the -COOH (binding model ) is energetically more advantageous, reflecting the greater acidity of the -COOH group, being closer to the electronegative N atom that exerts an electronwithdrawing polar effect. However, it is worth pointing out that the TiO2/12-FGA-system is characterized by a notably negative adsorption energy (-2.77 eV), that is the result of an extended network of intermolecular H-bonds among the amide groups of the FGA molecules (see Fig. S1), which would guide the arrangement of the ligands around the NP equator. Finally, comparing the adsorption energies at highest coverages, the FGA average adsorption energy for TiO2/52-FGAresults more negative (by 0.25 eV) with respect to that for TiO2/48-FGA-.
Therefore, to conclude, although the coordination of any further FGA molecule is always energetically less convenient with respect to the adsorption of the first molecule, the intermolecular FGA-FGA H-bonds can make the average FGA adsorption energy more negative and, thus, contribute to the stabilization of the overall structure. Moreover, according to our results, the -COOH group of FGA is found to be a better anchoring group than -COOH.

Testing the thermal stability of the nanoconjugates
In this section, we evaluate the chemical stability of the TiO2/52-FGA-and TiO2/48-FGAnanosystems by means of DFTB simulated annealing calculations. The simulated annealing is a computational technique, which allows to accelerate the exploration of the potential energy surface (PES) of a system, by heating it up to high temperatures and then slowly cooling it down, so that the system, during the process, is guided towards the global minimum of the PES. The protocol we followed for the simulated annealing calculations is the following: we heated up the TiO2/52-FGAand TiO2/48-FGA-from 0 to 500 K for 1 ps, then we performed a 4 ps equilibration phase at constant temperature and, finally, we slowly cooled down the systems to 0 K in 15 ps. In Fig. S2, the temperature profiles that we adopted in the annealing process are shown. After this simulated annealing process, we took the last snapshot from each simulation and performed a further DFTB geometry optimization at 0 K with the Conjugate Gradient algorithm. Table 2 reports the average FGA molecular adsorption energies, at the beginning (OPT/DFTB (1)) and at the end (MD/DFTB) of the 20 ps simulation and also after the final DFTB geometry optimization (OPT/DFTB (2)). In addition, Fig. 4 presents a schematic and visual representation of the TiO2/52-FGA-system along the simulated annealing.  (1)) and at the end (MD/DFTB) of the simulated annealing and after the final DFTB geometry optimization (OPT/DFTB (2)). From the results in this section, we conclude that the O-Ti bonds between FGA and the TiO2 NP are stable, even at high temperatures, such as those simulated in the annealing process, and that some proton transfers from FGA to TiO2 may spontaneously occur.

TiO2/FA high coverage models
Once   We observe that, in vacuum, FA stays closer to the NP (Fig. 6a,b) with the pterin bent to the TiO2 surface, forming H-bonds with the NP in Fig. 6b and intramolecular H-bonds in Fig. 6a. On the contrary, in water, the ligand is more involved in interactions with the solvent (Fig. 6c,d).
In Table S1 we present the average values of some structural parameters we calculated along the 100 ns simulations, including atomic distances, hydrogen bonds and non-bonding interaction energies. In water, both the center of mass of FA and the terminal amine group of the pterin lie at greater distances from the NP surface and center, respectively, with respect to vacuum. The NP-FA vdW interactions are comparable in vacuum and in water, but that is not true for the electrostatic ones.

MD simulations of a TiO2 NP functionalized with 8 FA molecules
In this section, we analyze the TiO2 nanoparticle that was functionalized with 8 folic acid molecules, corresponding to a FA/TiO2 weight ratio of about 0.2. We focus the attention on this coverage density because it was proposed by Lai et al. to produce nanoconjugates with higher cytotoxicity under photoexcitation. 20 We considered two systems where all the equatorial four-fold Ti 10  To quantify all these interactions, in Table S2 we

MD simulations of high coverage TiO2/FA systems
At last, we discuss the MD simulations of the high coverage TiO2/52-FA-and TiO2/48-FAsystems, which have been partially optimized at the DFTB level (see Fig. 5). Fig. 7 shows the last snapshots from the 100 ns MD simulations of the two systems, in vacuum and in aqueous environment at 303 K. Despite the high ligand density, when in vacuum again, the FA molecules bend towards the NP surface in order to maximize the interactions, whereas, in water we observe a competition between the FA-NP and the FA-water interactions.
To quantitatively evaluate this consideration that comes from a first visual inspection, in Table   S3  To compare all the TiO2/FA systems considered so far, in Fig. 8 we report the rdf of the FA center of mass with respect to the six-fold coordinated Ti atom at the center of the NP, for each model, averaged over the 100 ns of the MD simulations. The shape of the rdf curves becomes broader with the increase in the density of adsorbed folic acid molecules, whereas the maximum intensity, located at about [15][16] Å from the NP center, decreases.
In addition, the radial distribution of the ligands around the nanoparticle is always broader in water than in vacuum, because the FA molecules have more potential stable conformations due to the interaction with the solvent molecules.
Finally, we study the dynamical properties of the nanoconjugates in an aqueous environment.
In The diffusion of the nanoconjugates is a crucial aspect to be studied, as the functionalized nanoparticles, after injection, must reach their cellular target, thus the rate of the transport will influence the efficiency (and therefore, the efficacy) of the nanomedical process. To this purpose, in Table S4 we also report the estimated diffusion coefficients of all the TiO2/n-FA systems under investigation in water, computed by the OLS method at MD trajectory intervals showing a linear dependence between the MSD of atomic positions and time (see section 2.3 for more details). We observe that the values for all systems are of the same order of magnitude (i.e. 10 -10 m 2 /s). Moreover, the diffusion coefficients decrease with increasing FA density on the NP, except for the TiO2/1-FAsystem, whose D value is almost five times lower compared to the TiO2/1-FA-system.

Preparation of the models
In this section, we study the dynamics of a TiO2 nanoparticle grafted with PEG chains which are further functionalized with FA molecules. The PEG chain considered in this work is a methoxy-PEG, H3C−[OCH2CH2]n−OH, with n = 11, for a total molecular weight of 516 g/mol. One FA is covalently bonded to one of the PEG ends through an ester bond involving the FA -COOH. In Fig. S8 23 reported that the performance of folic acid on the selective targeting of the cancer cells was maximized at a FA/TiO2 weight ratio of 0.5, which is close to that obtained with the TiO2/PEG/20-FA model. Moreover, we also built up the TiO2/PEG/10-FA system, with a FA/TiO2 weight ratio of 0.25, that is close to the optimal value proposed in another work by Lai et al. 20 After a thousand geometry optimization steps, performed at the DFTB level, we obtain the two models shown in Fig. S9.

MD simulations of TiO2/PEG/10-FA and TiO2/PEG/20-FA systems
In Fig. 9 we present the last snapshots at the end of the 100 ns MD simulations of the TiO2/PEG/ To support our findings, in Table S5 we  In Fig. 10 we also report the rdf of the PEG and FA center of mass with respect to the six-fold coordinated Ti atom at the center of the nanoparticle, averaged over the 100 ns MD simulations, together with the total and LJ free energy profiles to transfer a FA molecule covalently attached to a PEG chain, from the NP surface towards the farthest distance compatible with the PEG chain length (see section 2.3 for more details). To better understand this behavior, in Fig. S10 we report some representative snapshots of the TiO2/PEG/10-FA system from the MD/ABF calculations, with the FA molecule restrained at different z-distances from the NP surface. In particular, we observe that the p-aminobenzoate ring changes its orientation from parallel (z-distance 10-15 Å, Fig. S10a) to perpendicular (z-distance 15-20 Å, Fig.   S10b) with respect to the NP surface: we believe that this second configuration is associated to a higher free energy. For the sake of completeness, in Fig. S11 we also report the snapshots of the TiO2/PEG/20-FA MD simulations. It is noteworthy to point out that the free energy profiles of Fig.   10 do not reach a plateau at increasing z-distance from the NP center: this is because the restrained FA molecule is covalently linked to a PEG chain, therefore pushing FA at high z-distances from the NP causes an increase of the free energy that comes from the constraint imposed by attached PEG chain, which in turn is covalently bound to the NP. To confirm these considerations, in Table S6 we report the estimated diffusion coefficients of the two TiO2/PEG/FA systems in water, computed by the OLS method at MD trajectory intervals showing a linear dependence between MSD and time. Once again, as in section 3.2.3, the diffusion coefficient values are in the order of magnitude of 10 -10 m 2 /s. In addition, the diffusion results slightly slower when the FA density is higher.

Conclusions
In this work, we have simulated the dynamics of titanium dioxide nanoparticles, when these are functionalized with folic acid, which is an efficient targeting agent of several tumor cell lines. First, we have built, through quantum mechanical modeling at the DFTB level of theory, several TiO2/n-FA systems, with increasing folic acid density and considering both the carboxylic groups of the ligand as anchoring groups. Our results reveal that the dissociative chelated and bidentate binding modes are energetically more advantageous and that the -COOH is the best anchoring carboxylic group. Moreover, the DFTB simulated annealing calculations confirmed the chemical stability of the nanoconjugates, when the effect of the temperature (500 K) is considered. Using the DFTB-optimized geometries as a starting point, we have performed long classical molecular dynamics simulations (100 ns) both in vacuum and in water. In general, we have observed a strong interplay in vacuum between folic acid molecules and the NP surface through hydrophobic and electrostatic interactions, whereas in water there is a fair competition between ligand/NP and ligand/solvent interactions. We have also determined the dynamical properties of the nanoconjugates by computing their diffusion coefficients in water and observed a nice correlation with the total mass of the systems.
Due to the fact that the close contact and the large interaction of folic acid molecules with the NP surface could hamper its binding to the folate receptor, we have decided to design models where a spacer is introduced between the NP and the FA molecule. In particular we prepared two models where the NP is grafted by 50 PEG chains, among which some (10 or 20) were selected to be covalently bonded to folic acid molecules at the opposite end. Based on the classical MD simulations of these two models and on some free energy calculations we could prove that the introduction of the molecular spacer reduces the extent of the NP-folic acid interactions. However, we have noticed that the folic acid content is a crucial parameter that highly affects the amount of the FA-FA intermolecular interactions but not much the diffusion coefficient of the nanosystems. According to our results, a low density of FA molecules in the nanoconjugates is expected to boost their performance in targeting tumor cells, due to a more pronounced availability of the ligands for the folate receptor.