Molecular dynamics simulation of the Pb(II) coordination in biological media via cationic dummy atom models

The coordination of Pb(II) in aqueous solutions containing thiols is a pivotal topic to the understanding of the pollutant potential of this cation. Based on its hard/soft borderline nature, Pb(II) forms stable hydrated ions as well as stable complexes with the thiol groups of proteins. In this paper, the modeling of Pb(II) coordination via classical molecular dynamics simulations was investigated to assess the possible use of non-bonded potentials for the description of the metal–ligand interaction. In particular, this study aimed at testing the capability of cationic dummy atom schemes—in which part of the mass and charge of the Pb(II) is fractioned in three or four sites anchored to the metal center—in reproducing the correct coordination geometry and, also, in describing the hard/soft borderline character of this cation. Preliminary DFT calculations were used to design two topological schemes, PB3 and PB4, that were subsequently implemented in the Amber force field and employed in molecular dynamics simulation of either pure water or thiol/thiolate-containing aqueous solutions. The PB3 scheme was then tested to model the binding of Pb(II) to the lead-sensing protein pbrR. The potential use of CDA topological schemes in the modeling of Pb(II) coordination was here critically discussed.


Introduction
The coordinative bond uncommonly combines high strength with reversibility, probably because of its dual, electrostatic and covalent nature. The HSAB theory introduced a paradigm that allows to assess the stability of a coordinative bond based on the metal and ligand properties, mostly reflecting their propensity to form mostly electrostatic (hard) or mostly covalent (soft) coordination [1][2][3][4]. In some instance, however, a metal ion may manifest either hard or soft behavior depending on the ligands available in a certain molecular system. The biological media are fair examples of molecular systems containing both hard and soft ligands. The extracellular matrix is generally approximated by an aqueous solution, thus being assumed as a strong hard environment. On the other hand, protein and/or other biomacromolecules may contain hydrophobic pockets confined out of the water bulk, where soft ligand groups may be present and generate compartmentalized soft environments. The Pb(II) ion is a typical borderline metal ion [5] for which either hydration or coordination at soft groups of proteins or other biological macromolecules might be predicted to occur. These properties may help to comprehend the high pollutant potential of lead, able to reach high concentrations in water or other hard ligand sites and, concomitantly, saturate also the available soft ligand sites.
Nowadays lead is omnipresent in batteries, herbicides, plumbing materials, gasoline and paints [6][7][8]. Unfortunately, its ubiquitous use ushers the widespread contamination of soil and water, causing the contamination of the whole food chain which leads to severe intoxication in humans, especially in children [9,10].
Primary targets for the Pb 2+ ions in bacteria as well as in human cells and intracellular liquids are thiols in neutral or deprotonated states (thiolates) [11][12][13][14][15][16][17][18][19][20][21][22]. The cysteinebased metal-binding domain is the basis of the lead defence mechanism in Cupriavidus metallidurans CH34, where the three-thiol motif of the pbr-resistant operon (pbrR) allows a flexible modulation of the lead affinity in response to the required metal-chaperone functionality-changing the protonation states of thiols, the chelating protein is able to 1 3 20 Page 2 of 12 induce the lead sequestration high affinity or release low affinity [23].
From a chemical point of view, the selective coordination of the Pb 2+ ion, operated by chelating proteins, unveils several intriguing aspects.
The presence of the 6s 2 electron pair in the valence configuration of this cation is known to affect both coordination number and coordination geometry, thus leading to either holo-directed or hemi-directed complexes [11-15, 17, 18, 20, 24]: in holo-directed complexes, ligands are evenly distributed around the metal center by occupying all coordination sites, whereas in hemi-directed complexes, some coordination sites are unoccupied and the ligand distribution around the metal center is uneven.
The sequestration of Pb(II) ion from the bulk, operated by protein thiols, can be described as a sequence of ligand exchange occurrences that gradually replace one or more waters in the first hydration shell with S-donating thiol groups. The computational modeling of these processes must be able to adequately describe this sort of competition between hard and soft coordination of Pb(II), i.e., between hydration and thiol coordination, respectively. As known, quantum chemistry (QM) approaches such as density functional theory (DFT) can accurately describe the structure and the stability of metal complexes; however, their employment in the investigation of metal coordination at biomolecules presents some limitations: (i) the application of standard QM methods to systems with more than one hundred atoms is limited by the high computational burden; (ii) generally, QM methods are addressed to either molecular or pseudomolecular models where bulk effects are often poorly incorporated; (iii) the exploration of the conformational space of biomacromolecules, often necessary to make the binding sites accessible to the metal ions, cannot be realistically afforded at QM level of theory.
Classical molecular dynamics (MD) might be considered an alternative computational tool by overcoming most of the above mentioned limitations of QM approaches; however, the description of metal ion coordination via empirical force fields is particularly challenging. On the one hand, the use of both bonded and non-bonded potentials may result in an adequate description of the metal-ligand bonding; on the other hand, the reversibility featuring the coordinative bond would be completely untreated.
In this study, the modeling of Pb(II) coordination in aqueous solutions containing thiol ligands was carried on by the employment of cationic dummy atom schemes. Preliminary density functional theory studies allowed to identify the most important structural parameters featuring the hydrated Pb(II) ion and both Pb(II)-thiol and Pb(II)-thiolate complexes; the latter being representative of the most common protein binding sites accessible to this cation. Two cationic dummy atom schemes, with three and four dummy sites, were designed to fraction part of the mass and charge of Pb(II). The corresponding PB3 and PB4 models of Pb(II) were incorporated in the Amber force field and employed to simulate the coordination of this cation in either pure water, thiol and/or thiolate-containing solutions, or in the three-cysteines binding site of pbrR, a lead-sensing protein characterized by a high affinity for Pb(II). Our results evidenced how the use of CDA schemes may effectively lead to some improvement in the description of the Pb(II) coordination and allowed to see how this methodology would be ameliorated in perspective.

DFT calculations
All DFT calculations were performed with the Gaussian 09 A.02 [25] quantum chemistry package. The DFT functionals are known to give a good description of geometries and reaction profiles for transition metal-containing compounds [26][27][28][29][30]. Optimizations, frequency calculations, electronic and solvation energy evaluations were carried out in solvated phase (C-PCM) [31,32] and by using the range-corrected density functional wB97X [33] which is known to give a good description of geometries and reaction profiles for transition metal-containing compounds [34][35][36]. The LANL08 (d) basis set [37,38] was used for lead both for optimizations and single-point energy calculations. All other atoms were described by the double-zeta basis set 6-31+G* [39][40][41][42] during optimizations and frequency analyses as well as by the triple-zeta basis set 6-311++G** [39,41,43,44] during single-point electronic and solvation energy calculations.
Frequency calculations were performed to verify the correct nature of the stationary points and to estimate zero-point energy (ZPE) and thermal corrections to thermodynamic properties.
Single-point electronic energy calculations were performed on the solvated-phase geometries. The C-PCM continuum solvent method was used to describe the solvation [32]. It has recently been shown to give considerably smaller errors than those for other continuum models for aqueous free energies of solvation for neutrals, cations and anions, and to be particularly effective for the computations of solution properties requiring a high accuracy of solution free energies [45]. Solvation free energies, taken as the difference between the solution energies and the gas phase energies, were added to the gas phase enthalpies and free energy values to obtain the corresponding values in the aqueous solution.

Molecular dynamics simulations
The molecular dynamics simulations were performed in the AMBER03 force field [46] by using the Gromacs package [47]. The 3-point water model SPC [48] was used with the typical density of liquid water at 300 K and 1.0 atm. The LINCS algorithm was adopted to constrain all bond lengths [49], and the long-range electrostatics was computed by the particle mesh Ewald method [50]. Trajectory analyses were carried out by using suitable Gromacs utilities with the support of VMD graphical interface [51].
The simulation of models PB1, PB3, and PB4 with water were performed via the following computational scheme: (1) local energy minimization, (2) equilibration NVT run of 5 ns at 300 K, (3) production run of 50 ns at 300 K in an isothermal/isobaric ensemble, using the velocity rescaling scheme for temperature and the isotropic Berendsen coupling scheme for pressure [52]. The investigated model was placed in a cubic box with dimensions of about 30 3 Å 3 and solvated with 880 water molecules. Being + 2 the charge of the considered models, electrical neutrality was gained by adding eight Cl − counterions. In the systems with (a) added 8 thiols, (b) added 8 thiolated, and (c) added 4 thiols and 4 thiolates, the number of water molecules was: (a) 880, (b) 866, (c) 864, whereas the added counterions were: (a) two Cl − , (b) six Na + , (c) two Na + .
The starting structure of pbrR691 (5GPE.pdb), a homologue of lead(II)-binding, protein pbrR was taken from a previous study [23]. The simulation of the protein with water and with models PB1 and PB3 added into the metal-binding locus was performed as follows: (1) local energy minimization, (2) equilibration NVT run of 50 ns at 300 K, (3) production run of 250 ns at 300 K in an isothermal/isobaric ensemble, using the velocity rescaling scheme for temperature and the isotropic Berendsen coupling scheme for pressure [52]. The protein system was placed in a rectangular box with dimensions about 95 × 72 × 64 Å 3 and solvated with 12,879 water molecules. Being + 8 the charge of the considered protein, electrical neutrality was gained by adding eight Clcounterions.
Model PB1 is a single Pb atom with the mass 207.2 a.u., charge + 2, van der Waals parameters sigma of 0.2245 nm and epsilon 0.17018074 kJ/mol. Model PB3 is constructed of a Pb atom with the mass 171.2 a.u., 0 charge, and 3 dummy atoms (charge + 2/3 on each) in a tetrahedral orientation situated at distances of 1.0 Å from the Pb atom. The dummy-Pb-dummy angles are 73.1°. Model PB4 is constructed of a Pb atom with the mass 159.2 a.u., 0 charge, and 4 dummy atoms (charged + 0.5) in a distorted trigonal bipyramidal orientation-one lacking position-situated at distances of 1.0 Å from the Pb atom. The dummy-Pb-dummy angles are 72.5°. In models PB3 and PB4, the van der Waals sigma and epsilon are 0.3 nm and 0.8 kJ/mol, whereas these parameters for the dummy atoms are 0.05 nm and 0 kJ/mol.
Parameter data files of PB3 and PB4 topological schemes are available (Note S1).

General considerations and results outline
The formation of metal complexes is the result of the ligand exchange equilibria occurring in solution and involving both water molecules and coordinating solutes.
We can ideally sketch the complexation of a generic M n+ metal ion in a water solution containing the coordinating ligand L at a certain concentration [L].
It is clear that, beside their relative concentrations, the chemical nature of M n+ and L plays a pivotal role in determining the most representative species, i.e., the possible values of m and k parameters in the Scheme 1.
The HSAB theory allows to qualitatively predict the stability of metal complexes by taking into account the electronic structure of metal and ligands involved [1][2][3][4]: Hard metals (high positive charge, low ionic radius) form stable complexes with hard ligands (highly electronegative, low polarizable atoms), while soft metals (low positive charge, high ionic radius) form stable complexes with soft ligands (low electronegative, highly polarizable ions).
The coordination of the Pb 2+ ion in water solution containing either thiol or thiolate ligands, being a focus of the present study, can be considered a case limit. Indeed, while the hard and soft nature of water and thiol/thiolate ligands, respectively, can be easily identified, the Pb 2+ ion has been classified within a group of borderline metals [53,54]. Hence, we can say that Pb 2+ could form stable complexes with either hard or soft ligands, and that the coordination of both water and thiol/thiolate molecules is expected to be facile.
Thus, the computational simulation of either Pb(II)water-thiol or Pb(II)-water-thiolate systems can be only performed whether the trade-off between Pb 2+ -O and Pb 2+ -S coordination may be adequately assessed.
In this study, the classical molecular dynamics (MD) simulation of water solutions of Pb 2+ containing or not either Scheme 1. Hydration and ligand exchange equilibria for an aqueous solution of the metal ion M n+ and neutral ligand L 20 Page 4 of 12 methanethiol (CH 3 SH) or methanethiolate (CH 3 S -)-hereafter these two species are often called thiol or thiolate, respectively-was carried out by testing several topological schemes based on the use of non-bonded potentials for the treatment of coordinative interactions. This choice was motivated by the requirement of an adequate modeling of the trade-off between O coordination and S coordination and by the fact that ligand exchange events can be detected in MD trajectories only via non-bonded potentials describing the Pb-O and Pb-S coordinative bonds. A brief outline of the presented results: a preliminary assessment of the structure of Pb-water and Pb-thiol/thiolate complexes via density functional theory (DFT) calculations and a description of the topological schemes employed in the MD simulations (Sect. 3.2); analyses of the MD trajectories and comparison with DFT and/or experimental data for the Pb(II)-water system (Sect. 3.3); MD simulations of Pb(II)-water-CH 3 SH and Pb(II)-water-CH 3 Ssystems (Sect. 3.4); MD simulations of the Pb 2+ coordination at the pbrR binding site (Sect. 3.5).

The Pb(II) coordination in pure water or thiol/ thiolate aqueous systems
The coordination of Pb(II) ion in pure water was initially assessed. The geometries of [Pb(OH 2 ) n ] 2+ complexes with n = 1-9 were optimized at DFT level of theory by the use of the wB97X exchange-correlation functional as described in the Methods section. All obtained structures of Pb-water clusters were characterized by a minimum Pb-O distance in the 2.27-2.43 Å range (Table 1), although, by increasing the number of water molecules (n > 4), we also detected slightly longer Pb-O distance, typically in the 2.47-2.87 Å range (Table 1), consistently with available computational data [55,56]. The calculated values of O-Pb-O angles for the n < 7 species were mostly detected in the range of either 64-75 (proximal) or 132-135 (distal) degrees, in agreement with octahedral-like spatial orientations of the ligands bound at the metal center as also reported by others [55][56][57][58][59][60][61][62].
Interestingly, all water molecules were detected in the first coordination sphere of both n = 7 and n = 9 species, whereas the Pb(II) center of the n = 8 complex was found to be coordinated by only six water molecules ( Table 1). The DFT energy values for the incremental addition of one water molecule in the considered Pb(II)-water complexes evidenced a stabilization of the coordination number with n = 4 or 6; indeed, for n > 6, the addition of water molecules was found to destabilize the metal complex. Although our DFT calculations cannot definitively disentangle the stability of [Pb(OH 2 ) n ] 2+ clusters with n = 4-6, [Pb(OH 2 ) 4 ] 2+ and [Pb(OH 2 ) 6 ] 2+ were assumed as the most stable species characterized by hemi-directed distorted tetrahedral and holo-directed octahedral coordination geometries, respectively (Fig. 1).
The DFT investigation of the str uctures of [Pb(SHCH 3 ) n ] 2+ and [Pb(SCH 3 ) n ] +2−n through the incremental binding of either thiol or thiolate ligands unveiled that a lower coordination number 3 and mostly hemi-directed coordination geometries affect these complexes (Fig. 1). The most representative geometrical parameters of Pb-thiol and Pb-thiolate complexes are reported in Table 1. The influence exerted by the 6s lone pair on the metal coordination produces a pronounced effect in Pb-thiol and Pb-thiolate complexes characterized by trigonal pyramidal and distorted tetrahedral coordination geometries, respectively. As shown, the average Pb-S distances differ by a 1.17 factor from the average Pb-O water distance which is very close to the 1.16-1.20 ratio between the S and O van der Waals radii [5, Table 9]; thus, stating that Pb-O and Pb-S bond lengths mostly reflect the different size of the ligand atom.

Topological schemes to model the Pb(II) coordination in a bulk of water molecules
The preliminary DFT assessment of the structure of either Pb(II)-water or Pb(II)-thiol/thiolate binary systems provided a body of information to initiate the parametrization of these species into the classical Amber force field [46]. Indeed, DFT calculations indicate that either hemi-or holodirected coordination of water at Pb(II) may be observed with cn of 4 or 6, respectively. Thus, the coordination of Pb(II) ions in water solution was simulated by using three different topological schemes named PB1, PB3, and PB4 (Fig. 2). In the simplest PB1 scheme, Pb(II) is essentially a spherical, bivalent cation, and its interactions with other atoms are modelled by only non-bonded potential terms of the Amber force field, i.e.,  electrostatic and Lennard-Jones potentials. Because of the high positive charge and spherical symmetry, the PB1 scheme controls the coordination of Pb(II) mostly through strong electrostatic interactions with negatively charged atoms and induces symmetrical arrangements of ligands around the metal center, i.e., consistently with the holodirected complexation. Alternatively, the simulation of the hydrated complex was performed by anchoring three of four dummy centers around the Pb(II) ion, giving rise to the PB3 and PB4 topological schemes, respectively. These schemes were designed to promote a hemi-directed coordination of Pb(II) by delocalizing part of the Pb mass and a fraction of the positive charge on the dummy atoms (Du) (Fig. 3). The PB3 and PB4 representations of the Pb(II) cation were mainly extrapolated from the CDA topologies of Ca 2+ and Mg 2+ reported by Sept et al. [63], which provided the criteria to assign the Pb-Du equilibrium distances and both charge and mass fractionation (See Methods section), whereas the equilibrium Du-Pb-Du angles were assigned their corresponding DFT calculated values (Fig. 2).
Molecular dynamics simulations of a Pb(II) ion in a bulk of TIP3 water molecules was then performed by using either schemes PB1, PB3, or PB4.
Radial distribution function (rdf) analyses of the Pb-O atomic pairs evidenced the most relevant differences in the modeling of the Pb(II)-water coordination by the different topological approaches. The rdf profile calculated on the PB1 trajectory showed a single peak at about 2.6 Å that is consistent with a symmetrical arrangement of water molecules around the metal center as expected (Fig. 3). The corresponding cumulative curve showed that 6-8 water molecules are on average coordinated to the Pb(II), by assuming the distance range 2.6-2.75 Å (Fig. 3). On the other hand, the rdf analyses of both PB3 and PB4 evidenced the presence of two peaks at about 2.4 and 2.7 Å that indicate a hemi-directed arrangement of water molecules in the Pb(II) first coordination sphere (Fig. 3). This outcome is in excellent agreement with our DFT data that also indicated two Pb-O coordination distances in models with 6-8 ligated waters, e.g., 2.43 and 2.66 Å detected in the [Pb(H 2 O) 6 ] 2+ model ( Table 1). Further accreditation of the CDA topological schemes was found by the calculation of either Pb(II) or water diffusion coefficients. All three topologies yielded good estimates of diffusion coefficients which were found to be consistent with available computational or experimental data, whereas the 3-point model PB3 was found to be slightly more accurate in the description of this phenomenon ( Table 2).

MD simulations of the Pb(II) coordination by either CH 3 SH or CH 3 Sin water at different thiol/thiolate concentrations
By testing the PB1, PB3, and PB4 schemes in the modeling of the Pb(II)-water system, we found that a higher consistency with both DFT results and available experimental data is achieved by the use of CDA topologies. Therefore, the preliminary DFT assessment of the structures of Pb(II)-thiol/ thiolate complexes evidenced a stability gain at a coordination number of 3. Thus, PB3 and PB4 were subsequently employed in the simulation of a quaternary system formed by Pb(II)-thiol-thiolate-water in the 1:4:4:864 ratio. This quaternary system was assumed to model water solutions characterized by a high molar concentration of both lead and sulfur ligand, where either thiol or thiolate coordination is expectedly favored. In fact, our calculations showed that Pb(II)-thiolate is the only possible sulfur coordination detectable within either PB3 or PB4 topological schemes. Moreover, the simulation of a ternary Pb(II)-thiol-water with a 1:8:880 ratio, with enhanced excess of thiol in the water solution, also showed no evidence of sulfur coordination. To better interpret the lack of thiol coordination in the simulated systems, we repeated the simulation of the Pb(II)-thiol-thiolate-water system within the topological scheme PB1. Also in this case, the only sulfur coordination to lead detected was ascribed to thiolate, with no evidence of thiol coordination. The rdf profiles of the Pb-S distances for the PB3 and PB4 simulations showed a single peak at 2.5 Å for both (Fig. 4), suggesting a substantially similar Pb(II) coordination sphere within either PB3 or PB4 molecular systems. Notably, the average Pb-S distance of 2.5 Å for the thiolate coordination to lead is very close to the average distance of 2.52 Å calculated at DFT level of theory (Table 1). Therefore, both CDA schemes were characterized by a very similar curve of Pb-S cumulative number that clearly states the presence of two thiolates in the first coordination sphere of Pb(II) center. Conversely, the rdf analysis of PB1 trajectory showed a single peak at approximately 2.7 Å together with a cumulative number of 2.5 thiolates per Pb(II) center (Fig. 4).

The MD simulation of the Pb(II) coordination on pbrR
The PB1 and PB3 topological schemes were then employed to simulate the binding of a Pb(II) ion at the corresponding binding site of pbrR, a Pb(II)-responsive bacterial protein [11,14]. The apo-structure of this lead-binding protein has been previously investigated by us [23], so that we report here only a brief description of this system. The structure of Pb-pbrR obtained by X-ray diffraction analyses consists of eight protomers divided into four functional homodimers; each dimer has two Pb(II) binding sites formed by the antiparallel interactions of the two so-called dimerization helix domains and composed of Cys113 and Cys122 (on one unit) and Cys78′ (on the other unit) (Fig. 5). Either PB1 or PB3 topological schemes were then applied to simulate the Pb(II)-pbrR adduct obtained by placing the metal ion at the geometrical centroid of the positions of Cys113, Cys122, and Cys78′ sulfur atoms.
In the former topological scheme, the Pb(II) metal ion is released within only 2 ps of the preliminary NVT simulation into the water bulk where it stably settled. The same result was also obtained by inducing a more gradual equilibration by the temporary freezing of the lead ion translation, thus testifying that within the PB1 scheme the lead ion is much more stable in its hydrated form rather than in its pbrR-bound form. Notably, such an outcome was found to be unaffected by the ionization status of the lead-binding cysteines; indeed, even by assuming the thiolate form of either Cys113, 112, or 78′, the lead(II) ion is released in a short time into the water bulk.
The use of PB3 topological scheme allowed to gain a better trade-off between hydration and pbrR binding of the Pb(II) ion. In the all-thiolate status of the binding site, the PB3-Pb(II) metal ion was found to remain stably coordinated at the side chain of Cys122, while the side chains Cys113 and Cys78′ detached in a few simulation steps to be replaced by O-coordinating ligands (Fig. 5). The Pb(II)-pbrR adduct obtained within the PB3 topological scheme is characterized by the stronger coordination of Cys122 thiolate group (2.55 Å), the carboxylate of Glu75 (2.26 Å), and the oxygen of a water molecule (2.43 Å) beside the weaker coordination of other O-coordinating groups (> 2.71 Å) as depicted in Fig. 5.

Discussion
The presence of proteins specialized in the sensing, trafficking, and elimination of Pb(II) ion in the bacterial proteome testifies the fact that evolution has impacted the elusive; dual nature of this species that can be considered a hard/soft borderline metal ion. The computational simulation of Pb(II) in biological media is thus issued by assessing the correct trade-off between hydration and protein coordination, or, in other words, between hard and soft coordination.
The use of an empirical force field for the simulation of metal-protein coordination in water allows the calculation of long time trajectories together with an adequate exploration of the protein conformational events, although the treatment of metal coordination in a classical force field is challenged by the high degree of covalency in such bonds.
In this study, we assessed the use of cationic dummy atom (CDA) topological schemes for the classical MD simulation of either hydration or protein coordination of the Pb(II) ion. Preliminary DFT calculations were carried out to investigate in detail the composition and the structure of both Pb(II)-water and Pb(II)-thiol/thiolate binary systems, the latter being small model of the lead coordination occurring in proteins. The investigation of the Pb(II) hydration at high level of theory reported by others [55,56,[58][59][60][61] predicted the holo-directed coordination of up to nine water molecules. In our DFT study, the incremental coordination of water molecules to the Pb(II) showed that, although optimized structures with up to nine H 2 O in the first coordination sphere could be obtained, the most stable [Pb(H 2 O) n ] 2+ complexes were characterized by n = 4 or 6 and, again, holo-directed coordination geometry ( Table 1). The NBO analysis of these Pb(II) complexes evidenced a particularly high s-character of the lone pair residual on the metal center, consistently with an almost uniform electron density on the holo-coordinated lead (Table 3). On the other hand, DFT calculations showed the presence of two different Pb-O distances of about 2.4 and 2.6 Å in the [Pb(H 2 O) 6 ] 2+ complex: These structural data are consistent with a small, but still present, degree of covalency affecting the Pb-O bond, somehow highlight the slightly stronger coordination of three out to six water in hydrated Pb(II). The DFT study of either [Pb(CH 3 SH) n ] 2+ or [PbCH 3 S) n ] +2−n clusters was also carried out through the incremental coordination up to n = 4 of either thiols or thiolates, respectively. Interestingly, our calculations indicated that n = 3 is the highest limit of the coordination number for both Pb(II)-thiol and Pb(II)-thiolate clusters (Table 1). We interpreted this outcome as a corroboration to what we already detected in the Pb(II)-water clusters: There are up to three stable coordination sites on the Pb(II) metal center, tailored by a high degree of covalency in the metal-ligand bond. Noticeably, both proximal and distal X-Pb-X angles were found in the 70°-95° range in either water or thiol/ thiolate complexes with cn = 3 (Table 1): These data can be easily explained by assuming the donation of either O or S lone pair to the three, mutually perpendicular 6p orbitals of lead; the presence of these orbital combinations was confirmed by NBO analyses.
These DFT calculations have already disclosed the different natures of Pb-O and Pb-S bonds: Higher coordination numbers (cn = 4-6) and holo-directed coordination geometries were detected in the Pb(II)-water binary systems, whereas lower coordination numbers (typically cn = 3) and hemi-directed geometries in Pb(II)-thiol/thiolate systems. The lower directionality detected in the coordination geometries of hydrated Pb(II) ion reflects the more pronounced electrostatic nature of the Pb-O bond compared to the Pb-S bond. Such an assumption was also corroborated by the electronic charge distribution in the most stable Pb-water and Pb-thiol clusters, i.e., [Pb(H 2 O) 6 ] 2+ and [Pb(CH 3 SH) 3 ] 2+ , respectively (Table 3). Indeed, the S coordination induced a lower (less positive) charge on the lead ion compared to the O coordination because of the more pronounced ligand-tometal charge transfer tailoring the Pb-S bond.
Our DFT results confirm that the Pb(II) ion may be involved in either less directional, more electrostatic or more directional, more covalent coordination, but, in either cases, at least three sites of stronger coordination may be identified.
CDA schemes were then sketched and tested to enhance the description of the Pb(II) ion within the Amber force field. Such an approach has been introduced for the investigation of mostly hard cations like Mg 2+ and Ca 2+ [63,65] characterized by symmetrical coordination and mostly electrostatic metal-ligand interactions. Here, two topological schemes composed of three and four cationic dummy atoms, i.e., PB3 and PB4, were taken into account and compared with a more classical scheme, i.e., PB1. In the PB3 and PB4 schemes, the 207.19 amu mass and the + 2 charge of Pb(II) ion were redistributed in a four-and five-body systems, respectively, as depicted in Fig. 2. The main idea behind the use of CDA anchored on the Pb(II) ion comes from the need to impart the hemi-directed coordination to this metal ion, especially when it is bound to soft ligands, and to decrease the electric field around Pb(II) through the delocalization of the + 2 charge on the dummy sites. Molecular dynamics (MD) simulations were thus carried out in either pure water, or, CH 3 SH or CH 3 Saqueous solutions in order to test the predictive capabilities of the PB1-4 schemes. Above all, the three topological schemes provided for an adequate description of Pb(II) hydration. On the other hand, the PB3 and PB4 schemes allowed to better describe the composition and geometry of the first hydration shell, showing two values of the Pb-O distance, approximately 2.4 and 2.7 Å (Fig. 3), in close agreement with the results of DFT calculations ( Table 1). The single value of Pb-O distance, gained with PB1, testifies that this topological scheme is less accurate in the modeling the Pb(II) hydration shell. A further corroboration to the enhancement gained with CDA schemes is represented by the predicted values of water and Pb(II) diffusion constants; the PB3 scheme showed the closest agreement with available experimental data (Table 2). Notably, all these MD outcomes are consistent with our assumption that three sites, accessible to a stronger coordination, exist on the Pb(II) ion center.
In presence of thiol or thiolate ligands, the hydrated Pb(II) ion is expected to undergo to ligand exchanges where one or more water molecules are replaced by S-coordinating ligands. The implementation of Amber parameters describing the Pb(II) coordination allowed to simulate the ligand exchange events taking place in water solution of either CH 3 SH or CH 3 S − , assumed as reduced models of the leadbinding protein sites. Again, the CDA schemes PB3 and PB4 provided for a better description of the first coordination sphere, yielding the structures with average Pb-S distances in close agreement with the corresponding DFT data. On the other hand, a limitation was seemingly evidenced in all analyzed topological schemes: No evidence of thiol coordination at Pb(II) was detected even at high concentrations of CH 3 SH, while the coordination of only two thiolates was found in all systems containing the CH 3 S -. These results can be interpreted by the different nature of the Pb-S coordinative bond in Pb(II)-thiol and Pb(II)-thiolate complexes; indeed, a higher electrostatic contribution is expected to characterize the coordination of thiolate. In all PB1-4 topological schemes, the lead coordination is treated by the use of non-bonded potentials, i.e., Coulomb and Lennard-Jones potentials, that are expected to adequately describe only the electrostatic and exchange contributions to the Pb-ligand interaction, respectively, whereas the orbital contributions, expected to hold most of the covalency in the coordinative bond, are only indirectly taken into account in the PB3 and PB4 schemes. The Pb(II)-thiolate complexes, where a higher weight of electrostatics is expected, are thus better described by both CDA schemes, whereas the description of Pb(II)-thiol complexes is lacking the important contribution of orbital interactions. Although these outcomes might represent a limitation to the employment of non-bonded topological schemes for the description of Pb(II) coordination by sulfur ligands, it should be also taken into account that the binding of thiol at metal centers may often be followed by deprotonation, leading to the corresponding thiolate complex. Several theoretical and experimental evidences have confirmed that a marked decrease of the thiol pKa occurs upon coordination at even soft metal centers such as Au(I) [66][67][68][69], Au(III) [70] and also Pb(II) [23]. Thus, the use of empirical force field implemented with non-bonded potentials describing the Pb(II) coordination at sulfur ligands may still be considered valuable by assuming the thiolate, and not thiol, as the metal-binding form. For example, the classical MD simulation of the binding of Pb(II) at the thiol side chains of cysteines may be afforded with the caveat that the presumed binding residues are assumed in their thiolate forms. The pbrR protein is a Pb(II) sensing protein, structured as a homodimer, and characterized by two lead-binding sites each formed by a cysteine triad. As previously described by us [23], each cysteine triad composed of Cys113, Cys122, and Cys78′ (Fig. 5) is predicted to assume the thiol-thiolate-thiolate, respectively, form in the coordination of Pb(II).
In the present study, the simulation of the coordination of one Pb(II) at pbrR was performed by the employment of either PB1 or PB3 topological schemes and by assuming the all-thiolate protonation state of the cysteine triad. As expected, the accentuated hard character induced by the PB1 scheme caused the release of the Pb(II) ion in the water bulk, i.e., a fairly hard environment, in only 2 ps, whereas a stable confinement of the Pb(II) ion in the pbrR binding site can be detected only with the PB3 scheme (Fig. 5). In fact, the Pb(II) coordination detected in the equilibrated PB3 trajectory is not consistent with the X-ray structure of lead-bound pbrR [14]: we found that only the thiolate group of Cys122 is kept bound at Pb(II), while both Cys113 and Cys78′ are replaced in a few steps by oxygen ligands (Fig. 5). Again, these outcomes showed that although CDA topological schemes may contribute to better describe the coordination of borderline metal ions such as Pb(II), in terms of both directionality and attenuation of electrostatics in the involved coordinative bonds, the lack of explicit terms describing the ligandto-metal orbital interaction is a limit to their predictive capabilities.

Conclusions
In this paper, the employment of cationic dummy atom schemes to describe the Pb(II) coordination in water and in presence of thiol/thiolate ligands has been assessed with the ultimate aim of modeling the binding of this pollutant metal ion at its protein biomolecular target via classical MD simulations. Preliminary DFT investigations indicated that the Pb(II) coordination by water gives rise to holo-directed structures with a maximum of six ligated waters, although two different lengths of the Pb-O bond were detected. On the other hand, both thiol and thiolate complexes of Pb(II) resulted to be hemi-directed with at most three S-ligands coordinated to the metal center. The use of three dummy sites delocalizing the + 2 charge and part of the Pb mass, in the topological scheme named PB3, resulted to be the most effective in the modeling of both hydration and thiolate coordination of the Pb(II). The PB3 scheme provided a good prediction of the composition and geometry of the Pb(II) coordination sphere, and, also, a good estimate of either water and Pb(II) diffusion constants. Although the use of this CDA scheme represents a clear breakthrough compared to the standard, onebody description of Pb(II)-here named PB1 scheme, the lack of explicit terms accounting for the covalency of the Pb-ligand interaction was found to be a potential limit.
Nevertheless, the use of the PB3 topological scheme may be taken into consideration as a starting point: Within this scheme, the electrostatic and exchange contributions to the Pb(II) coordinative bonds, as well as their inherent directionality, can be adequately described. In perspective, we envision that the PB3 scheme could be improved by the introduction of pairwise potentials, e.g., Pb-S(thiol) or Pb-S(thiolate), that are able to account for the orbital interaction contributing to each bond.
Funding Open Access funding provided by Università degli Studi G. D'Annunzio Chieti Pescara within the CRUI-CARE Agreement. The authors gratefully acknowledge the University G. D'Annunzio (FAR 2019) for the financial support.
Data availability All data are available on request to the authors.