Dispersion state phase diagram of citrate-coated metallic nanoparticles in saline solutions

The fundamental interactions underlying citrate-mediated chemical stability of metal nanoparticles, and their surface characteristics dictating particle dispersion/aggregation in aqueous solutions, are largely unclear. Here, we developed a theoretical model to estimate the stoichiometry of small, charged ligands (like citrate) chemisorbed onto spherical metallic nanoparticles and coupled it with atomistic molecular dynamics simulations to define the uncovered solvent-accessible surface area of the nanoparticle. Then, we integrated coarse-grained molecular dynamics simulations and two-body free energy calculations to define dispersion state phase diagrams for charged metal nanoparticles in a range of medium’s ionic strength, a known trigger for aggregation. Ultraviolet-visible spectroscopy experiments of citrate-capped nanocolloids validated our predictions and extended our results to nanoparticles up to 35 nm. Altogether, our results disclose a complex interplay between the particle size, its surface charge density, and the ionic strength of the medium, which ultimately clarifies how these variables impact colloidal stability.


Model development
As described in the main text, by minimizing the change of free energy ΔG cap we recover the value of N of the highest likelihood. The first contribution, ΔG desolv , can be expressed in terms of the desolvation energy of a single ligand (ΔG desolv lig ) multiplied by N and that of the metallic NP (ΔG desolv M ). The linear trend followed by ΔG desolv lig is plotted as a function of N in Supplementary   Fig. 1. Note, moreover, that ΔG desolv M cancels out with the apolar contribution from the solvation energy of the capped nanoparticle ΔG solv apolar .

a. Solvation free energy calculationsthermodynamic integration protocol
In order to calculate the solvation free energy ΔG desolv lig , we adopted the thermodynamic integration methodology. In our case, we first immersed an explicit coarse-grained citrate molecule in a box of polarizable water molecules parametrized with the refPol force field. 1-3 The free energy calculation was stratified into 31 runs with increasing values of the coupling parameter, . The spacing between windows was of 0.05 for < 0.6 and 0.02 elsewhere. A value of 0.5 was used as the alpha parameter in the soft core potential. 4,5 For each window in the calculation, two minimizations were performed for a maximum of 5000 steps each. The first of them used the steepest descent method, and the second one followed the l-BFGS algorithm. Once the systems were minimized, we proceeded to thermalize and pressurize to 310K and 1 bar with a 10 ns long simulation in the NPT statistical ensemble. These employed the Berendsen thermostat (τT = 2.0 ps) and the isotropic Berendsen barostat (τP = 5.0 ps and κ = 4.5 × 10 -4 bar -1 ). After the equilibration, a production run was performed for each window for 25 ns using a time-step of 20 fs and switching the Berendsen barostat for its Parrinello-Rahman counterpart (τP = 14.0 ps and κ = 4.5 × 10 -4 bar -1 ). 6 The variations in the Hamiltonian and its derivative respect to were saved every 10 steps (0.2 ps). Finally, the solvation free energy was estimated with Bennett's acceptance ratio method. 7 The aforementioned procedure was performed twice to sequentially switch off the electrostatic and van der Waals interactions, respectively.

S4
Supplementary Fig. 1 The different contributions accounted by the main terms included in the theoretical model as a function of the number of bound citrate molecules. The plots shown are built for α = 3. a. Breakdown of ΔG desolv . b. Breakdown of ΔG bind . c. Breakdown of ΔG solv .

b. Citrate binding electronic energy
In order to determine the number of equivalent binding modes of citrate onto gold surfaces, we first generated a realistic, all-atom representation of a gold NP. For this, we used the NanoCrystal web server. 8 The circumradius parsed to the server was 1.5 nm, and the surface energies utilized for the gold facets were those derived by Skriver and co-workers. 9 The result was an NP formed by six (100) planes and eight (111)  Here, is the electronic binding energy associated with the -th binding mode, its multiplicity, and − the corresponding Boltzmann weighting factor. Following this scheme, we estimated the average electronic binding energy between citrate and gold to be -40.9 kcal mol -1 .

S5
Supplementary Fig. 2 Graphical representation of all the possible binding modes for acetate, succinate, and glutarate as described by Al-Johani, et al 10 and using their same nomenclature. Each binding mode is accompanied by its multiplicity, that is, equivalent binding orientations in the grid considering the interchangeability of the carboxylate groups and their respective oxygen atoms. Orange spheres are gold atoms from a particular lattice facet.

c. Entropic contribution to binding
In order to account for the entropic contribution of binding a ligand (e.g., citrate) onto a metallic surface, three terms were considered. These terms account for the decrease in translational and rotational degrees of freedom as well as the insurgence of vibrational modes upon binding. The total change in entropy can thus be written as in Supplementary Equation (2).
The translational loss of entropy is calculated from the entropy of an ideal gas employing the Where is the molar concentration of citrate, is the ideal gas constant, is Avogadro's number, and Λ is the thermal de Broglie wavelength given by Supplementary Equation (4): Where is the mass of one ligand molecule, is Boltzmann's constant, and the temperature which in this case was set to 310K. Similarly, the loss of molecular rotational entropy was calculated from Supplementary Equation (5).
Where is the geometric mean of the molecule's moments of inertia, as stated by Supplementary Equation (6).
Where 1 , 2 , and 3 are the three principal moments of inertia calculated from the diagonalization of the tensor of inertia. In our case, we derived the moments of inertia from an all-atom molecular dynamics simulation of a single citrate molecule immersed in water. The setup of this simulation is described in Supplementary Discussion 1f. Finally, the decrease in vibrational entropy arises from the formation of metal-ligand bonds. Their contribution was estimated from Supplementary Equation (7).
Where is the stretching frequency of the formed bond. For the specific case of citrate, the frequency for the Au-O bonds was taken as 1638 cm -1 , as measured by Wulandari et al. 11 It is noteworthy that, as evidenced in Supplementary Fig. 1, the entropic correction for citrate binding onto gold is negligible compared to the process's enthalpy.

d. Derivation of /
The electric potential energy of a charge density distribution ( ) is given by Supplementary Equation (8).
In which ( ) is the electric potential generated by the charge distribution, ( ), and the integral covers the entire space. For the case of a hollow sphere with uniform surface charge density, the charge distribution function can be written as in Supplementary Equation (9).

S7
( , ) = ( ) 4 2 ( − R NP ) (9) Where R NP is the radius of the NP, ( ) is the total charge, and is the Dirac delta function. In order to generalize our model, we consider a generic, polyprotic ligand that adopts a mean deprotonation state α when bound to the NP. For a ligand with three acid protons like citrate, the mean deprotonation state lies between 0 and 3, that is α ⊆ [0,3]. Then, the total charge may be written as in Supplementary Equation (10).
On the other hand, plugging (9) into Gauss's Law, it is possible to obtain the electric potential generated by the charge distribution in a medium of dielectric constant . This is given by Supplementary Equation (11).
Substituting (9) and (11) into (8), and solving the integral, we obtain an expression for the electric potential energy as a function of the number of bound ligands (Supplementary Equation (12)).
Note that this expression assumes each incoming charge as independent from the rest; however, each ligand constraints the incorporation of α fundamental charges to the system. Thus, we can write ΔE lig/lig as in Supplementary Equations (13) and (14).
S8 Supplementary Equations (14) and (16) imply that the citrate-citrate and citrate-solvent interactions is mainly electrostatic. Thus, the theoretical framework presented here remains valid for small, charged ligands that do not affect greatly the NP's hydrophobicity.

f. Explicit citrate coarse-grained model for solvation free energy calculations
For the calculation of the solvation free energy of one citrate molecule (i.e., ΔG desolv lig in Supplementary Discussion 1a), we adopted an explicit coarse-grained representation. The model for citrate corresponded to three consecutive beads of type Qda as implemented in the standard Martini force field. Each bead had a charge of -1 to reproduce the most populated protonation state of fully solvated citrate molecules. The parameters used for the chemical bonds were kb = 5000 kJ mol -1 nm -2 and x0 = 0.52 nm, the typical inter-bead distance in the Martini force field. The parameters used for the angle between the beads were kθ = 5000 kJ mol -1 rad -2 and θ0 = 111.11°.
The explicit-citrate coarse-grained model was parametrized against atomistic MD simulation. In detail, the equilibrium angle θ0 was obtained from a 100 ns-long, all-atom MD simulations of citrate (parametrized with the GAFF force field 13 ) in TIP3P water. 16 This simulation was also used to calculate the three main moments of inertia for a citrate molecule (i.e., Ii in Supplementary Discussion 1c). In addition, the same atomistic simulation was used to obtain the projected area per citrate molecule onto our AuNPs (i.e., a in Supplementary Discussion 2).

Estimation of citrate surface coverage
The fraction of the surface area covered by citrate molecules, φ, was calculated by means of Supplementary Equation (17).
In which is the total surface area of the NP, is the number of citrate molecules bound to the metallic core, and is the area covered by one citrate molecule. for a series of MD-accessed conformers and calculating the projected area. Specifically, the structure of citrate was withdrawn every 10 ps from the trajectory file. For each of these states, 100 random rotations were generated. Then, the molecular surface of each of these configurations was projected onto the XY plane, and the area was estimated. The distribution of all the calculated areas is shown in Supplementary Fig. 3. Supplementary Fig. 3 Illustration of the method employed for the calculation of the projected surface area of citrate. For every frame in the trajectory, a series of orientations were sampled. The van der Waals surface was then projected onto the XY plane, and its area was calculated. The figure also shows the probability distribution function of the projected area, which has a maximal probability at 0.16 nm 2 . Carbon atoms are shown in cyan, oxygen atoms in red, and hydrogen atoms in white. The XY plane is colored green.

Counterion residence time analysis and sampling of slow degrees of freedom
For this study, we had to ensure that the rotational degrees of freedom of the implicit-citrate coarsegrained NPs were exhaustively sampled at short values of CV1. To understand the effectiveness of our PMF calculations, we evaluated four distinct sampling protocols. The first protocol (referred to as P1) consisted of equal, short production runs. In P1, all the windows were simulated for 25 ns, each. In the second protocol, P2, the setup was similar. The frames where CV1 > 1.2 nm were simulated as in P1, but the windows with CV1 ≤ 1.2 nm were simulated for 250 ns. P3, the third protocol, implied more brusque modifications onto the coordinate files of frames with CV1 ≤ 1.2 nm. From these frames on P2, the counterions residing 0.7 nm from any metallic bead for at least 50 ns were exchanged with water molecules in bulk. Once the ions were displaced, the same frames were simulated for 250 ns. The last protocol, P4, included the simulated annealing described in the Methods section of the main text. The described protocols were tested on two systems. Both systems involved NPs with a total charge of -40 e (charge density 1.9 e nm -2 ), and they had respective ionic strengths of 30 mM and 70 mM.

S10
As can be seen from Supplementary Fig. 4a, the estimated free energy of dimerization varies significantly depending on the assumed sampling protocol. It is important to note that, even though the charged beads follow a uniform distribution on the NPs, their placement is not perfectly symmetrical. This leads to the distinguishability of states depending on the relative orientation of both NPs. When P1 and P2 are contrasted, it becomes evident that the energy profiles are smoothened as more time was simulated in each window. This indicates that the NP explored more widely their rotational degrees of freedom. The local energy minima visited by the particles is a consequence of the formation of a suboptimal network of salt bridges between metallic, charged beads, and sodium counterions at different NP orientations.
To understand better the dynamics of the surrounding ions, we computed the residence time of each binding event onto the metallic surfaces. Supplementary Fig. 4b shows the residence time for P2 (i.e., sampling for 250 ns when CV1 ≤ 1.2 nm) in the two test cases, that is, an NP with total charge -40 e (charge density 1.9 e nm -2 ) at an ionic strength of 30 and 70 mM. Even though most of the binding events persisted for less than 50 ns, there was a non-negligible amount of ions permanently bound to the particles. In order to ensure that these ions were not kept in their position by steric traps, we calculated the same residence times as a function of the inter-particle distance (CV1, Supplementary Fig. 4c). Interestingly, as the NPs came closer together, the number of highresidence binding events increased. This trend suggests that the binding of the ions was mediated by electrostatic interactions, which in turn were enhanced as the two rigid bodies approached each other.
Supplementary Fig. 4 Assessment of the best sampling protocol. a. PMF profiles obtained for NPs with a net charge of -40 e (charge density 1.9 e nm -2 ) in a solution of ionic strength 30 mM (purple) and 70 mM (orange). b. Residence times of sodium ions binding onto the surface of the NPs for the protocol P2. The solid lines show the average and standard deviation calculated from the windows with CV1 ≤ 1.2 nm. c. The residence time of sodium ions as a function of the interparticle distance (CV1) for protocol P2.
Based on these results, we adopted two approaches that could enhance the sampling of the diffusion of the sodium counterions. The first of these approaches, P3, consisted of the manual replacement of the persistent ions with bulk water molecules. Upon this replacement, a new set of ions appeared in the same halo-like structure (Fig. 3). Since the relative position of the hydrophobic NPs and the respective salt bridges differ from P1 and P2, the free energy profiles varied slightly ( Supplementary Fig. 4a). Nonetheless, P3 introduced an abrupt alteration to the free energy landscape leading to an instantaneous push of the system away from its thermodynamic equilibrium. In light of these observations, a final strategy, P4, was considered. In the protocol P4, the conformational space of the counterions was extensively sampled by performing a simulated annealing prior to the data collection run. As can be seen from the two test cases exhibited in Supplementary Fig. 4a, the proper phase space sampling leads to a significant decrease in the free energy of aggregation. P4 was the chosen sampling protocol for the simulations discussed in the main text of this work.

S12
It is important to note that the sodium ions' residence times computed from our coarse-grained models offer a qualitative description of the lability of the counterions present, rather than structural insights on their exact binding mode onto metallic NPs. Thus, in order to verify the results of our coarse-grained models, we performed MD simulations with atomistic representations of equivalent citrate-capped AuNPs in electrolytic solutions. These simulations enabled us to calculate, at a finer resolution, the residence time of the sodium counterions. Indeed, our atomistic simulations of citrate-capped AuNPs displayed the same trend, that is, as the charge density of the NP increased, the sodium ions became more affine for the NP's surface, thus bonding for longer times ( Supplementary Fig. 5). and 30 monohydrogencitrate molecules, i.e., surface charge density of 0.9, 1.9, and 2.9 e nm -2 , respectively. The atomistic metallic core was constructed with the Wolff method implemented in the NanoCrystal web server, 8 as described in Supplementary Discussion 1, and its nonbonded parameters were taken from Heinz and co-workers. 12 The monohydrogencitrate molecules were parametrized with the GAFF force field 13 (charges derived with the RESP method 14,15 ). The coating ligands, which were frozen throughout the simulation, were placed with their two terminal S13 carboxylate groups pointing towards the gold core, as elucidated with high-level quantum mechanical calculations. 10 The system was solvated in electrolytic solutions (I = 0, 30, 70, and 170 mM) using the TIP3P water model 16 and AMBER14SB ions. 17 Each combination of surface charge and ionic strength was simulated for 250 ns.

Correlation between well-depth ε and free energy of dimerization
Supplementary Fig. 7 Correlation between the hydrophobicity of the metallic NPs (as described by the well-depth ε) and the free energy of NP dimerization.

Coarse-grained models for metal NPs a. The building of coarse-grained hydrophobic NPs
The metallic nanoparticles (NPs) employed in this study consisted of 126 beads arranged in rings stacked on top of each other as described elsewhere (Supplementary Fig. 9a). [18][19][20][21] The beads were assigned the C1 type as incorporated in the Martini v2.2 force field, 22,23 and which corresponds to a purely hydrophobic moiety. Moreover, the mass of an NP represented at the atomic level was distributed on all the available beads of the coarse-grained model. This mass corresponded to a value of 556 u.m.a per bead. The spherical shape of the NPs was retained by imposing an elastic network with a force constant of 15000 kJ mol -1 nm -2 . The elastic network united each bead with its six nearest neighbors as well as their farthest neighbor ( Supplementary Fig. 9a).

b. The building of charged NPs for the building of dispersion state phase diagrams
In this study, we used implicit-ligand models for our charged metallic NPs. In detail, our capped NPs were represented as hydrophobic cores with a partial charge in certain surface beads. For the construction of these bodies, we employed the same algorithm as for pristine NPs (Supplementary Discussion 7a) modifying the charge of selected hydrophobic beads according to the target charge S15 of the NP (Supplementary Fig. 9b). The sites were selected by uniformly shuffling all the beads according to their spherical coordinates. The charge of each bead was chosen as -2 e to mimic the most populated deprotonation state of citrate molecules onto spherical gold NPs. 24 For example, an NP with a total charge of -40 e (charge density 1.9 e nm -2 ), would have 20 surface beads with a charge of -2 e, each. Importantly, our charged-sphere models are valid for gold NPs coated with small, dianionic capping agents. Consequently, the phase diagrams derived from these models NPs describe the dispersion state of gold NPs in the presence of ligands like citrate, cyclic oxocarbons, and dicarboxylic acids.
Supplementary Fig. 9 Graphical representation of the models employed for the various studied components. a. The pristine metallic NPs (orange) were constructed from a series of concentric rings stacked over one another and held together by an elastic network. b. An explicit coarsegrained model for citrate (cyan) was built as three charged beads, and it was used only to calculate ΔG desolv lig . The citrate-capped NPs used to derive the dispersion state phase diagrams were modeled as hollow spheres conformed by hydrophobic beads, some of which carried a formal charge. The number of charged sites was proportional to the net charge of the NP. c. Initial placement of two NPs immersed in a saline solution for the potential of mean force calculations. The collective variable CV1 was defined as the closest distance between the two bodies' surface.