Simulation of uranyl-biomolecule interaction using a cationic dummy atom model

We report a multisite cationic dummy atom model of the uranyl cation, suitable for atomistic simulation of the interaction of this ion with water and biomolecules. This reproduces geometry of model complexes with typical ligands, interactions and diffusion coefficient in water, and the structure of an adduct with a uranyl binding protein. A protocol for prediction of binding mode and strength of complexes formed with small peptides is proposed.


Introduction
Uranium is an actinide element which occurs in sea water and mineral ores, which poses serious radiological and toxic risk for all living organisms because even after millions of years of evolution, it was not incorporated in any form of life [1].In humans, its targets are kidney and bones, yet no consensus exists on the specific interactions that account for toxicity.Uranium possesses several redox states, the most stable of which in aqueous solutions is the linear uranyl [U(VI)O 2 ] 2+ , which coordinates 4-6 ligands in equatorial plane with affinity for hard oxygen donors; it is present in < 10 ppb concentrations in typical environments [2].
EXAFS spectra of uranyl in aqueous solutions demonstrated that the first solvation shell of the uranyl cations contains five waters located in equatorial plane at 2.4 to 2.5 Å [3,4].The coordination to five water molecules was corroborated by X-ray scattering on uranyl in aqueous solution, albeit with 12% of uranyl ions bound to four waters [5].More recently, neutron scattering was used to determine an average water coordination number of 4.6 in 1 M UO 2 Cl 2 [6].
Design of uranyl sequestering agents is complex due to need for affinity and selectivity over ions such as Ca 2+ , requiring femtomolar (10 − 15 M) affinity.Amidoxime, [7] phosphine and hydroxy-phosphino, [8] and carboxylate-based ligands [2] have been tested.Peptides robustly chelate metals, making them ideal candidates: [9] the metalbinding motifs inspired by osteopontin, [10] calmodulin, [11] and albumin [12] were employed for uranyl chelation, demonstrating excellent affinities, with carboxylates of glutamate, aspartate and C-termini as well as phosphorylated residues heavily involved [13].The engineered "super uranyl-binding protein" SUP demonstrated exceptional affinity and selectivity for uranyl with a K d of 7.4 femtomolar (fM) and > 10000-fold selectivity over other metal ions [14].A uranyl-selective protein designed on the basis of a nickel(II)-responsive protein binds uranyl with a dissociation constant K d = 53 nM [15].
Computational methods are a valuable complement to experiment for study of toxic, radioactive uranyl species.DFT studies are numerous: simulation of aqueous solvation corroborates coordination of uranyl by 5 waters, although relative energies of four-and six-coordinated complexes differ between methods [16,17].In uranyl-protein studies, DFT was used for the reduction mechanism of uranyl by G. sulfurreducens and D. acetoxidans, [18] where coordination to carboxylates of glutamate and carboxylate agreed with experiment.Uranyl complexes were investigated in the super uranyl binding protein, [19] and with calmodulin via a combination of spectroscopy and DFT [20].Uranyl coordination to human serum apotransferrin protein was predicted to occur through tyrosine, aspartate, and carbonate ion, in agreement with experimental data [21].
The computational cost of DFT limits application to small systems and/or few conformations.Classical molecular mechanics allied to molecular dynamics (MD) allows consideration of whole, flexible protein structures.Nevertheless, it is a demanding task to describe the coordination of metal via empirical force fields.The first force field for UO 2 2+ was developed in the 1990s, [22] and non-additive potentials allowed development of polarizable force fields [23,24].These reproduced the distance between uranyl and oxygen of water, yet the polarizable potential was computationally costly.One breakthrough was a force field [25] which was able to reproduce uranyl's hydration free energy and uranium-oxygen distance.This force field facilitated the emergence of hydrogen bonds between water molecules and the axial oxygen atoms of uranyl.Another series of force fields for actinyls AnO 2 n+ (n = 1, 2) ranging from U to Am [26] were parametrized to fit the potential energy surface obtained via B3LYP.The hydration free energies for AnO 2 2+ were less exothermic than elsewhere [27] and the actinide (VI)-oxygen distances were overestimated with respect to experiment.MD simulation of uranyl-protein systems are scarce.Recently, molecular structure of the surface-immobilized super uranyl binding protein was studied by discrete molecular dynamics simulation, providing an in-depth understanding of molecular interactions between SUP and the surface and the effect of uranyl ion binding on the SUP interfacial structures [28].MD studies of SUP demonstrated that uranyl chelation is determined by amino acids in the binding site, waters in the first coordination sphere, and the sturdiness of the second-sphere hydrogen bond network [19].MD simulations gave structural models for uranyl coordination in the native and phosphorylated calmodulin siteI, the EF-hand motif of which has a 1000-fold affinity toward uranyl with respect to calcium [29].The interaction between uranyl and ubiquitin was characterized via MD simulations, disclosing its structural transformations at atomistic level [30].Uranyl in this case also was found to coordinate to glutamate and aspartate.
In this study, the modelling of uranyl-biomolecules coordination was performed via use of non-bonded potentials for the description of the uranyl-ligand interaction.We developed a cationic dummy atom scheme for uranyl cation, characterized by parts of the mass and charge of the uranyl being fractioned in five sites rigidly anchored to uranium.Cationic dummy model is a robust concept in bioinorganic chemistry, whose use has been shown to address shortcomings in description of polyvalent ions using conventional, atom-centred potentials.Cationic dummy models for have been developed for a very wide range of maingroup and transition metal ions [31], [38]and successfully utilized for atomistic description of biological problems, mostly in metallo-and metal-binding proteins, but to date no application to uranyl has been reported.After confirming that our scheme reproduces the correct coordination geometry of model uranyl complexes, we employed the molecular dynamics simulation of uranyl with a set of model biomolecules.

Methods
DFT calculations were performed using the Orca package v5.0.3 [32] using the BP86-D3BJ functional [33][34][35] and basis set consisting of either small-core (60 electron) ECP [36] on U and def2-TZVP [35] on remaining atoms for smaller models, or large-core (78 electron) ECP on U and def2-SVP on light atoms for larger systems.We note that structural data are not strongly dependent on functional choice (Table S1).All such calculations included CPCM implicit model of aqueous solvation [37].Parameters for molecular mechanics description of uranyl, including equilibrium bond length and angle and force constants, were extracted from harmonic frequency calculation on model systems A cationic dummy atom model was constructed by adapting the nonbonded parameters reported in reference [38].5 dummy atoms were placed in a regular pentagon in the equatorial plane of [UO 2 ] 2+ , each bearing a charge of + 0.5 e. Charges of − 0.25 on each uranyl O are used to balance the positive charge on U, and reflect the relatively low Lewis acidity of these atoms.Lennard Jones well-depths for U and O were taken from ref. [26], and each dummy atom assigned a small LJ radius and mass of around 5% that of U. Fig. 1 shows the general structure of the CDA.
Molecular mechanics simulations were performed using the AMBER16 package [39].Library and forcefield modification files encoding the CDA parameters noted above were constructed, and used together with standard amino acid and solvent libraries to construct systems of interest in the leap utility of AmberTools.Initial tests used steepest descent and conjugate gradient minimisation in GBSA implicit solvent model [40] for comparison with DFT-optimised or experimental data.Subsequent dynamics simulation used explicit solvent in TIP3P model of water [41].Analysis was performed with cpptraj utility of AmberTools, along with Chimera and VMD.We also include a prototypical peptide complex, [UO 2 (Asp) 3 ], for which coordination mode around equatorial region of uranyl is comparable to DFT but peptide conformation varies.

Initial
Further validation was sought from explicit solvent simulation of the uranyl cation, placed within a cubic box of TIP3P water with at least Å from ion to box edge.After heating in NVT ensemble from 100 to K, followed by equilibration of box size with NPT simulation at 298 K and 1 bar, both over 1 ns, production MD in NPT ensemble at 298 K for 150 ns was used to extract radial distribution function of water molecules' solvation of uranyl.RDF of water molecules to central U(VI) centre (Fig. 3) shows clear first and second solvation shells, the former between 2.5 and 2.9 Å, the latter between 3.0 and 3.7 Å.These integrate to 5 and 11 water molecules, respectively, the former in excellent agreement with previous reports of equatorial coordination numbers.RDF for water distribution (Fig. S1) around uranyl O atoms indicates first solvation shell between 2.0 and 3.5 Å, reflecting hydrogen bonding to these atoms, and integrates to 15 molecules in total for both oxygen atoms.
The same simulation was used to extract diffusion coefficient of uranyl, taking only the final 100 ns that shows approximately linear trend in mean square displacement with time.This yields a value of 1.1 × 10 -5 cm 2 s − 1 , identical to that reported by Perez-Conesa et al. [42] To account for systematic errors resulting from viscosity effects in periodic MD simulation, the ratio of diffusion coefficients of uranyl and water from the same simulation is found to be 0.254, close to the experimental value of 0.3 ± 0.01 but notably smaller than Perez-Conesa et al's value of 0.4.
As a further test, we extracted coordinates of a single chain of the uranyl binding protein (PDB entry 4FZP), a protein that exhibits femtomolar affinity by rational design of suitable aspartic and glutamic acid residues [19].A CDA model of this protein was constructed from the experimental data, retaining two water molecules that are in close proximity to the uranium centre, and orienting the pentagon of dummy atoms in the frame of reference defined by experimental uranyl coordinates.The resulting structure was minimised in TIP3P water over 2000 steepest descent then conjugate gradient steps.Fig. 4 shows two views of the resulting structure, the first showing the equatorial coordination of Glu17, Asp68 and two water molecules in the equatorial plane of uranyl, as well as the Arg residue that lies close to both acidic groups.Essentially the same pattern of interactions remains when the same procedure was performed in implicit GBSA water model.
Following minimisation, molecular dynamics simulation was performed, first heating the structure from 100 to 298 K over 1 ns, then a further 20 ns of NVT dynamics, in TIP3P water.Within 0.5 ns of MD, Glu65 enters the coordination sphere of uranyl, along with Glu18, Asp69   and 2 H 2 O noted above, leading to 5-coordinate equatorial binding of uranyl.However, we see no evidence of Asp14 engaging in coordination, although this residue does form a persistent hydrogen bond with Glu18 as part of the secondary coordination sphere noted in ref. [19].
Following this equilibrium MD, we estimated the free energy of binding of uranyl, using steered MD to increase the distance between U and the centroid of the Cγ/Cδ of Glu18, Glu65 and Asp69 by 10 Å over the course of 20 ns.This results in a free energy of binding of − 34.9 kcal/ mol, which is greater than but comparable to the experimental value of − 19.3 kcal/mol derived from the reported binding constant of 7.4 fM.We also carried out MM-PBSA estimation of free energy of binding, arriving at ΔHbind = -64.64± 0.47 kcal/mol and TΔSbind = -17.85± 1.99 kcal/mol, giving an overall free energy of binding of − 46.8 kcal/ mol.We hope to report improved attempts to calculate the free energy of binding, for example using thermodynamic integration or umbrella sampling, in future work.
A particular goal for development of a non-bonded model of uranyl coordination was to determine whether binding modes to peptides could be predicted a priori.As a test, we used two short peptides derived from bovine milk proteins, Val-Glu-Ser-Lys (VESL) and its serinephosphorylated analogue VES P L [43].Each peptide was built in extended geometry and combined with CDA model of uranyl, and solvated in implicit GBSA water.Minimisation and conventional MD of the separated systems did not result in any specific uranyl-peptide contacts.We therefore used short, steered MD to push the uranyl cation into contact with the peptides, using distance between U and Cγ of Glu as a general coordinate.For VESL, this showed a barrier of ca. 4 kcal/mol for approach, before a sharp drop to a stable orientation at − 3.8 kcal/mol corresponding to coordination of uranyl by Glu (Fig. S2).We then used self-guided Langevin dynamics (SGLD) to explore the conformational flexibility of the complex formed, which quickly located a more stable form in which C-terminal carboxylate and backbone O also bind to U (Fig. 5).In this, U⋯O distances were found to be 2.54 Å to Glu, 2.64 Å to Ser backbone, and 2.51 Å to C-terminus.Conventional MD over 20 ns retained this coordination mode, although backbone O from Glu and Ser interchanged in first coordination sphere.Following this equilibration, a second, slower steered MD run was used to estimate the binding free energy of − 20.2 kcal/mol (Fig. 5).Here, the distance of U from Cα of Ser was chosen as the coordinate, to give general measure of separation without bias toward any particular binding residue.Fig. 5b shows a well-defined minimum in which all contacts shown in Fig. 5b are present, as well as a broad, shallow plateau in which contact with C-terminus is lost but coordination through Glu is retained.
Applying the same protocol to the phosphorylated peptide VES P L yields similar results, although here the initial steered MD already locates a chelated complex in which both Glu and Ser P bind in the equatorial plane of uranyl.SGLD conformational search and equilibration MD do little to change this picture, and no evidence of backbone or Cterminal coordination is found (Fig. 6).We ascribe this to the relative stability and rigidity of the Glu-U-Ser P chelate formed, which has U⋯O distances of 2.62 Å to Glu and 2.58 Å to Ser P .Steered MD to escape from this low energy chelate gave a binding energy of − 21.9 kcal/mol (Fig. 6), in qualitative agreement with Zänker et al's report of increased strength of binding on phosphorylation.Here again a well-defined minimum is evident, along with a broad plateau in which coordination to Ser P is lost but Glu coordination is retained.
In summary, we have developed and tested a cationic dummy approach to describing equatorial coordination of the uranyl cation.By distributing 5 dummy atoms around the equator to carry positive charge, this model correctly reproduces the structure of model complexes as well as the adduct formed with a uranyl-binding protein.It also gives reasonable prediction of binding energy as well as diffusion coefficient.We then show that prediction of binding modes to small peptides can be achieved with a simple protocol of conventional and steered MD validation concentrated in comparison of DFT and CDA minimised structures of model uranyl complexes [UO 2 (H 2 O) 5 ] 2+ , [UO 2 (OC(H)NH 2 ) 5 ] 2+ , and [UO 2 (O 2 CCH 3 ) 3 ].Superposition of structures, as shown in Fig. 2, indicates that key features are reproduced, most notably the classic equatorial coordination of ligands around the axial uranyl group.Key distances are reproduced well, including uranyl = 1.80 Å (DFT value 1.76 Å), equatorial U-OH 2 = 2.56 Å (DFT value 2.49 Å).U-O = C(H)NH 2 = 2.58 Å (DFT 2.40 Å), U-OC(O)CH 3 = 2.54 Å (DFT 2.51 Å).

Fig. 1 .
Fig. 1.Structure of cationic dummy atom, with U shown as teal sphere, O as red spheres, and dummy atoms in cyan.

Fig. 4 .
Fig. 4. Minimised structure of 4FZP uranyl binding site, with comparison with experimental coordinates (minimised geometry shown as RGB, experimental as cyan).

Fig. 5 .
Fig. 5. (a) low energy conformation of vesl-uranyl identified by sgld; b) free energy profile for removal of uranyl from low energy conformation.

Fig. 6 .
Fig. 6.(a) low energy conformation of ves P L-uranyl identified by SGLD; b) Free energy profile for removal of uranyl from low energy conformation.