Atomistic modelling and structural characterisation of coated gold nanoparticles for biomedical applications

This study presents the results of atomistic structural characterisation of 3.7 nm diameter gold nanoparticles (NP) coated with polymer polyethylene glycol (PEG)-based ligands of different lengths (containing $2-14$ monomers) and solvated in water. The system size and composition are selected in connection to several experimental studies of radiosensitisation mechanisms of gold NPs. The coating structure and water distribution near the NP surface are characterised on the atomistic level by means of molecular dynamics simulations. The results of simulations carried out in this study, combined with the results of our recent study [J. Phys. Chem. A 126 (2022) 2170] and those from the field of polymer physics, are used to calculate key structural parameters of the coatings of radiosensitising gold NPs. On this basis, connections between the coating structure and distribution of water are established for different NP sizes as well as lengths and surface densities of coating molecules. The quantitative analysis of water distribution in the vicinity of coated metal NPs can be used to evaluate the radiosensitising effectiveness of a particular NP system based on the proximity of water to the NP metal core, which should impact the production of hydroxyl radicals and reactive oxygen species in the vicinity of metal NPs exposed to ionising radiation.


I. INTRODUCTION
Metal nanoparticles (NPs) have garnered much attention in recent years due to their wide range of biotechnology and biomedical applications, including medical imaging [1], drug delivery [2], and radiotherapy of cancers [3,4]. When injected into the body and localised within a tumour region, metal NPs can act as radiosensitising agents through irradiation with x-rays and ion beams [3,5,6]. The interaction of metal (typically gold) NPs with ionising radiation results in the strong emission of secondary electrons [7][8][9], which interact with molecules of the surrounding medium (that is, to a large extent, water). As a result of such interactions, hydroxyl radicals and other reactive oxygen species are produced, which may deliver damage to tumour cells [10,11].
At the same time, the use of metal NPs for biomedical applications can lead to different cytotoxic effects, including immune responses [12], inflammation [13], and damage to cell membranes [14]. In order to reduce toxicity, metal NPs proposed for biomedical applications are usually synthesised with organic coatings, which also help to increase the uptake and longevity of NPs within cells [15]. One of the most popular coatings used in experiments with radiosensitising NPs is poly(ethylene glycol) (PEG), often chosen as it increases the stability and reduces the toxicity of NPs injected into biological systems [16].
Structural and energetic properties of coated (also commonly referred to as functionalised or ligandprotected) metal NPs have been widely studied computationally by means of density functional theory (DFT), * M.D.Dickers@kent.ac.uk † verkhovtsev@mbnexplorer.com classical molecular dynamics (MD), and Monte Carlo methods [17][18][19][20][21][22][23][24][25][26]. In particular, several studies [21,25] were devoted to the atomistic MD-based modelling and structural characterisation of coated metal NPs in view of their radiosensitising properties. In the cited studies, gold NPs with a diameter of 1.6 nm were coated with thiol-PEG n -amine (S-(CH 2 ) 2 -PEG n -NH 2 ) molecules and solvated in water. The structural properties of the coatings and the distribution of water in the vicinity of the NPs were analysed at the atomistic level. In Ref. [27], the transport of low-energy secondary electrons and the production of hydroxyl radicals in the vicinity of 1.6-nm gold NPs coated with thiol-PEG-amine ligands has been studied theoretically and computationally. It was demonstrated that the production of secondary electrons and OH • radicals in the vicinity of coated metal NPs is highly dependent on water distribution near the NP surface. A greater amount of water close to the NP surface results in a larger yield of hydroxyl radicals formed in that spatial region. As such, determining water distribution in the vicinity of coated metal NPs is essential to quantify the production of secondary electrons and free radicals and, therefore, to evaluate the effectiveness of such NPs as radiosensitising agents for a particular NP size and elemental composition. The distribution of water surrounding the gold NPs is highly dependent on the structure of the coating, with higher-density coatings preventing water molecules from penetrating close to the NP surface [11,25]. It is, therefore, essential to understand, on the atomistic level, the structure of the NP coating to predict the distribution of water molecules.
Gold NPs studied in experiments on NP radiosensitisation have a wide range of sizes from ∼2 to 100 nm [28]. For NPs with a size of several (∼1-10) nanometers, the structural characterisation of the metal core, coating, and surrounding molecular environment can be performed us-ing the all-atom MD-based approach [21,25].
However, as the system size increases, this approach becomes computationally expensive and even infeasible for NPs with a size of ∼10 nm and above, considered by some experimental studies [10,11,28]. A more viable approach is to develop a model for describing and predicting the structural properties of coated metal NPs, such that the coating structure and the distribution of water molecules near the metal NP surface can be analysed for any NP size as well as for any length and surface density of ligand molecules.
This paper aims to take a further step in this direction. In the present paper, the structure of an experimentally relevant system -a gold NP with a diameter of 3.7 nm functionalised with thiol-PEG n -amine (n = 2 − 14) molecules and solvated in water is characterised on the atomistic level using the MBN Explorer [29] and MBN Studio [30] software packages. Based on the results of structural characterisation carried out in the present study and those from Ref. [25], key characteristics and common features of the coating structure and distribution of water are analyzed for a broad range of coating densities.
The interconnection of structural and chemical parameters of the coating (composition, length, and surface density of ligand molecules) and NP size with the structural properties of water near the NP surface is an interesting and open research question, as these relations cannot be readily determined experimentally. Computer simulations can provide necessary insights into these characteristics at the atomistic level, but they are limited to relatively small NP and ligand sizes.
A 3.7 nm gold NP has been selected as a case study in this paper. NPs of that size, functionalised with PEG, were used in radiosensitisation experiments [28,31]. In those experiments, the particular system size was chosen for the relatively large surface area of gold NPs for efficient attachment of PEG molecules, low toxicity, and rapid cell penetration. The chosen system's size also falls within the range of "ideal" NP sizes of ∼1-50 nm in diameter to pass through the cell membrane without causing cellular damage [32]. The 3.7-nm NPs are also large enough to be visible using transmission electron microscopy with a typical image resolution of ∼1 nm [28].
Polymer physics can also provide insight into the structure of organic coatings of radiosensitising metal NPs. There is a significant body of research into the physics of polymers attached to flat and curved surfaces, stemming from the work of Flory [33] and, in particular, de Gennes [34]. The structural properties of polymer chains attached to flat surfaces is a well-understood problem, studied for the first time by Alexander [35] and de Gennes [36]. The extension of these models to account for polymers attached to curved surfaces was first made by Daoud and Cotton [37] and further developed by Wijmans and Zhulina [38], providing a theoretical description of the height of a polymer chain attached to a curved surface in the limit of both large and small cur-vatures (that is, metal cores with small and large radii, respectively). The results of MD simulations presented in this study, in combination with the results of our recent study [25] and those from the field of polymer physics [26,38,39], are used to calculate key structural parameters of the coatings of radiosensitising gold NPs and establish connections between the coating structure and distribution of water for different NP sizes as well as lengths and surface densities of coating molecules. The results of the quantitative analysis of water distribution in the vicinity of coated metal NPs can subsequently be used in future studies to evaluate the radiosensitising effectiveness of a particular NP system based on the proximity of water to the NP metal core.

II. METHODOLOGY
This section briefly outlines the computational methodology for the creation and atomistic-level characterisation of 3.7-nm gold NPs coated with thiol-PEGamine molecules and for complementary simulations of thiol-PEG-amine molecules attached to a flat gold surface. The computational protocol used to simulate coated 3.7-nm gold NPs solvated in water follows that outlined in the earlier study [25], including the creation of a metal core, coating of the metal core with thiol-PEG-amine molecules of different lengths, annealing of the coated systems in a vacuum, and, finally, solvation of systems in water. Modifications from the previously reported procedure are described in greater detail in the following sections.
In this study, classical MD simulations have been performed using the MBN Explorer software package [29] for advanced multiscale modelling of complex molecular and nanoscale systems. MBN Studio [30], a multitask toolkit and a dedicated graphical user interface for MBN Explorer, was used to construct the systems, prepare input files, and analyse simulation outputs.
A. Thiol-PEG-amine coated gold nanoparticles Three PEG n molecules with the number of ethylene glycol monomers n equal to 2, 8 and 14 have been considered in this study. The initial geometries of these molecules were taken from the LigandExpo database [40]. Each PEG n molecule is functionalised with a thiol (-SH) group at one end and an amino (-NH 2 ) group at the other end, forming thiol-PEG n -amine (n = 2, 8, 14) molecules. Chemical formulas and characteristics of the considered ligand molecules are listed in Table I. A spherical gold NP with a diameter of 3.7 nm, containing 1553 atoms, was created using the modeller tool of MBN Studio [30]. In experiments involving radiosensitizing gold NPs, the size of the NPs cannot be controlled precisely, and the typical size distribution of NPs is characterized by the standard deviation of about 1 nm. The TABLE I. Specification of three thiol-PEGn-amine (n = 2, 8,14) ligand molecules considered in this study. Natom is the number of atoms in each molecule; their chemical formulas and molecular weights are listed in the last two columns. The corresponding ID for the original molecules taken from the LigandExpo [40]  constructed NPs thus represent a prototype of a statistically probable system synthesized in experiments involving radiosensitizing NPs. For large NPs the crystallographic planes within the gold core may impact the final coating morphology [41], however for the NP size considered in this work these effects are insignificant. Thiol-PEG n -amine molecules were then attached nonsite-specifically to the NP surface as described in Ref. [25]. The number of attached molecules varied from 44 to 300, corresponding to the range of ligand surface densities σ = 1.0 − 7.0 molecules per nm 2 . This range of ligand densities has typically been considered in experimental studies of coated gold NPs for biotechnology applications [42]. The surface density of ligands σ is defined as the number of ligands per unit surface area and is also commonly referred to as grafting density [26,43,44]. Seven values of ligand surface density have been considered for each ligand length. The total number of atoms in free gold NPs coated with thiol-PEG n -amine (n = 2 − 14) varied from ∼2600 atoms for the shortest ligand molecule and lowest surface density to ∼33,500 atoms for the longest ligand molecule and highest density cases, respectively.
In the computational protocol utilized in the present paper, the interaction between all gold atoms of the metal core (including atoms located on the NP surface and in its interior) has been modelled by means of the manybody Gupta potential [45]. The covalent interaction between the gold atoms located on the metal core surface and sulphur atoms of the ligands is modeled using the Morse potential. Parameters of this interaction and noncovalent interactions between gold atoms and atoms of the ligands were previously validated through DFT calculations; see Ref. [25] for details. All interatomic interactions within the thiol-PEG-amine ligands and their interactions with water molecules after solvation were described using the CHARMM molecular mechanics force field [46]. CHARMM parameters for the bonded and non-bonded interactions between atoms of the ligands and gold atoms are listed in Ref. [25].
After the ligand attachment, the coated metal NPs were annealed (by heating the system from 0 to 600 K at a rate of 1 K/ps and then cooling the system back to 0 K at the same rate) before being introduced to a solvent. The process of solvation of coated metal NPs in water follows the protocol outlined in Ref. [25]; the TIP3P water model [47] was used to describe the solvent. The number density of water molecules was set to 33.4 nm −3 , corresponding to the accepted number density of water molecules in a liquid medium at ambient conditions. All simulations were conducted for the canonical (NVT) ensemble of particles at 300 K and at constant pressure. The geometry of solvated NPs was first optimised using the velocity quenching algorithm over 20,000 optimisation steps. Then the systems were equilibrated at 300 K using the Langevin thermostat for 5 ns. The last 2 ns of each simulated trajectory were used for the structural characterisation analysis described in Section III A.
In the case of the gold NPs coated with thiol-PEG 2amine ligands, a simulation box of (10 × 10 × 10) nm 3 was used, and the total number of atoms in the system (including solvent) was equal to ∼106,000. For the case of thiol-PEG 8 -amine ligands, the simulation box size was set to (12 × 12 × 12) nm 3 , with a characteristic system size of ∼182,000 atoms. Finally, for thiol-PEG 14 -amine ligands, a simulation box of (15 × 15 × 15) nm 3 was used, with the characteristic number of atoms in such systems equal to ∼360,000 atoms. 21 systems (7 ligand surface densities at 3 ligand lengths) have been modelled and characterised in this study. The total simulation time was about 100,000 CPU hours.

B. Thiol-PEG-amine coated gold surfaces
Complementary simulations have also been performed to characterise the structure of thiol-PEG n -amine (n = 2 − 14) ligands attached to a flat gold substrate. A gold slab with a size of (4.08 × 4.08 × 1.5) nm 3 containing five layers of gold atoms was created from an ideal gold crystal by means of the modeller tool of MBN Studio [30]; the created substrate contained 1000 gold atoms. The simulations were carried out using periodic boundary conditions, such that the system is periodically translated in xy-plane.
As the next step, the gold substrate was coated with thiol-PEG n -amine (n = 2, 8, 14) molecules. Ligand molecules were evenly distributed above the substrate surface. Three surface densities of ligands were chosen, namely 1.0, 2.3 and 7.0 molecules per nm 2 . These particular values were chosen to describe the case of low, intermediate, and high surface densities corresponding to the range of values used for thiol-PEG n -amine coated gold NPs. The number and spacing of thiol-PEG n -amine molecules on top of the flat gold surface were dictated by the chosen surface density.
The process of attaching thiol-PEG n -amine molecules to the substrate followed the same protocol outlined in Section II A and in Ref. [25]. To avoid the translational motion of the whole system in the z-direction, the bottom layer of gold atoms was fixed during the process of ligand attachment. Following this, a series of constanttemperature MD simulations have been conducted to equilibrate the systems at 300 K for a period of 0.6 ns.
The equilibrated systems were solvated in water following the protocol outlined above. The simulation box size was set equal to (4.08 × 4.08 × 6.0) nm 3 for thiol-PEG 2amine and thiol-PEG 8 -amine and (4.08 × 4.08 × 7.0) nm 3 for thiol-PEG 14 -amine. In the case of the highest ligand surface density of 7.0 nm −2 , the simulation box size in the z direction was increased to 7.5 nm due to the elongation of the ligands as a result of electrostatic repulsion between ligand atoms. The solvated systems underwent energy minimisation to avoid any possible overlap between water and ligand molecules.
In order to create physically correct system states thermalised at a given temperature, geometry optimisation calculations were followed with two subsequent MD simulations. The first was a short 20 ps simulation at 300 K to pre-equilibrate the systems, followed by a longer 5 ns long MD simulation to thermalise the systems at 300 K.

III. RESULTS AND DISCUSSION
Structural analysis of the coatings of nanometer-sized gold NPs can provide valuable insight into the effect of coating structure on the distribution of water in the vicinity of the solvated NPs. Such an analysis is performed for NPs of two sizes (1.6 and 3.7 nm in diameter), different ligand lengths, and ligand surface densities. Based on the outcomes of this analysis, the relationship between the radial distribution of atoms of the coating and the distribution of water in the vicinity of the NPs is explored.
A. Structural characterisation of thiol-PEG-amine coatings and their conformation states Figure 1(a) shows the atomic number density of a coating made of 44 to 300 thiol-PEG 8 -amine molecules attached to 3.7 nm gold NPs (corresponding to the ligand surface density in the range from 1.0 to 7.0 nm −2 ). The number density is plotted as a function of the average distance from the NP surface. The atomic number density has been calculated as the number of atoms of the coating in spherical bins of 0.1 nm thickness, divided by the volume of each bin. The curves presented in Figure 1 and further below in this section were averaged over 50 independent frames taken from the final 2 ns long segments of the simulated trajectories. Figure 1(b) shows the results of a similar analysis [25] conducted for 1.6 nm diameter gold NPs coated with thiol-PEG 8 -amine molecules for the same range of ligand surface densities (the number of molecules varies from 8 to 56). Figure 1 illustrates that the profiles of the atomic number densities of the coating for the NPs of two different sizes are qualitatively similar. The maximum atomic number density of the coating is located at ∼0.35 nm from the NP surface in the case of a 3.7 nm NP and at ∼0.3 nm from the surface for a 1.6 nm gold NP. The maximum atomic number density for the coating of the 3.7 nm NP is generally higher than that for the 1.6 nm NP, implying a dependence of the maximum number density of the coating on the NP radius.
In both cases, a shoulder in the number density distributions is formed at distances ∼ 1.9 − 2.4 nm from the surface of a 3.7 nm NP and at distances ∼ 1.6 − 2.2 nm from the surface of a 1.6 nm NP. This shoulder appears at intermediate ligand surface densities of σ ∼ 4.0 nm −2 (172 and 32 ligands; see the light-green curves) but is most pronounced at densities larger than 5.0 nm −2 (216 and 40 ligands; see the dark-green curves). The appearance of the shoulder is associated with a change between conformation states of the coating molecules from "mushroom"-like structures to "brush"-like structures [25]. In polymer physics, conformation states describe the structures of polymer coating attached to surfaces. There are two main conformation states, namely "mushroom" and "brush" conformations [36,48,49]. The "mushroom" conformation occurs when ligands are loosely packed together and do not interact with each other. In this case, the ligands are strongly bent. "Mushroom" confirmations typically occur at lower ligand surface densities. The "brush" conformation occurs at higher ligand densities when the ligands begin to interact with each other and repel each other. This interaction causes the ligands to straighten [34,36,38,39,48]. Figure 2 shows the dependence of the maximum number density values of thiol-PEG n -amine ligands on ligand surface density. In all the cases considered, a decrease in the derivative of this function is observed at characteristic surface densities, which approximately correspond to the surface densities of the transition between the "mushroom"-like and "brush"-like conformations. The surface density value, at which the derivative decreases, varies with ligand length, occurring at lower surface densities for longer ligands. The details of this transition are elaborated below in Section III A 3. This transition is less pronounced for the smaller 1.6 nm diameter NP (dashed curves in Fig. 2), where a smooth increase in the peak number density is observed, while for the larger 3.7 nm diameter NP (solid curves), the change in the derivative is more evident.
The atomic number density distributions plotted in Figure 1 can be understood further by analysing the transition from the "mushroom"-like to "brush"-like structures. Figure 3 shows the probability distribution of finding nitrogen atoms of the terminal NH 2 groups of thiol-PEG n -amine molecules (n = 2, 8, 14) at a given distance from the center of mass (CoM) of the NP. This is determined by summing the number of nitrogen atoms in 0.1 nm thick spherical bins and dividing it by the number of ligand molecules in the coating. The probability distribution is plotted as a function of ligand surface density and the corresponding number of ligand molecules.
As the nitrogen atoms are present in the terminal group of the ligands, their position is directly related to the conformation state seen [25]. For lower density (fewer ligand) cases, an island of increased probability is seen close to the NP surface (located at the distance of ∼1.85 nm from CoM), corresponding to the "mushroom"like structures. For increasing surface density, the probability distribution shifts to favour "brush"-like conformations. In all cases, a transition region is also seen, with a mixture of conformations being present for ligands surface densities up to ∼ 3.0 − 4.0 nm −2 . At higher ligand densities, the conformation state of the coating is exclusively "brush"-like. In Figure 1(a), the shoulder begins to appear for surface densities of between 4.0 nm −2 and 5.0 nm −2 (corresponding to 172 and 216 ligands, respectively). As discussed below, the origin of the shoulder can be attributed to the position of the terminal NH 2 group. Thus, the presence of a mixture of "mushroom"like and "brush"-like conformation states may suppress the formation of a shoulder at lower coating densities.

Coating thickness
Different conformation states in thiol-PEG n -amine ligands can also be identified through the calculation of the average thickness of the coating. Figure 4(a) shows the thickness of the coatings made of thiol-PEG n -amine (n = 2, 8, 14) ligands attached to the surface of a 3.7 nm gold NP. Dashed curves with circles show the average distance between sulphur and nitrogen atoms,r S−N , and the solid curves with squares show the distance between the average position of sulphur atoms and the distance from CoM that encloses 95% of all coating atoms, r 95% −r S .
At small numbers of attached ligands (i.e., at low surface densities), the distance (r 95% −r S ) exceeds the distancer S−N for all the cases considered, indicating that at lower surface densities, ligands are strongly bent, and thus they attain a "mushroom"-like conformation. For the cases of thiol-PEG 8 -amine and thiol-PEG 14 -amine (blue and green curves, respectively) ther S−N distance ex- These are plotted as a function of surface density of ligands / the corresponding number of ligand molecules, and distance from NP center of mass (CoM). The color scale shows the probability of finding a nitrogen atom at a given distance from CoM. Note the different color scale in each panel.
ceeds the (r 95% −r S ) distance at higher ligand densities. The S-N distancer S−N becomes larger than (r 95% −r S ) for the coatings made of ∼172 to 216 ligands. Such shows the calculated thickness L of thiol-PEGn-amine coating as a function of the number of monomers n for the 3.7 nm NP considered in the present study (solid symbols) and the 1.6 nm NP studied in Ref. [25] (open symbols). The calculated values are fitted to the L ∝ n 0.6 and L ∝ n 0.8 dependencies known from the polymer physics for the "mushroom" and "brush" polymer conformations, respectively [39]. a transition from "mushroom"-like to "brush"-like structures is also seen in the nitrogen probability distributions shown in Figure 3, indicating that at higher ligand densities the coating takes on a "brush"-like conformation. In the case of thiol-PEG 2 -amine, the 95% distance exceeds the S-N distance for all considered values of ligand densities, indicating that a purely "brush"-like regime is not observed in this case. Rather, at lower ligand densities, the thiol-PEG 2 -amine ligands have the "mushroom"-like shapes, while at higher densities the coating represents a mixture of "mushroom" and "brush"-like conformations.
The calculated values of coating thickness obtained from the MD simulations have been compared to the values predicted by polymer theory [39]. According to the latter, the thickness L of a coating made of polymer molecules (containing n repeating units) grafted to a curved surface scales as L ∝ n 0.6 in the "mushroom" regime and L ∝ n 0.8 in the "brush" regime [39]. Figure 4(b) shows a fit of the calculated thiol-PEG n -amine coating thickness as a function of the number of ethylene glycol units n to these theoretical scaling laws. Solid symbols correspond to the case of a 3.7 nm gold NP considered in the present study, while open symbols denote the results from Ref. [25] for a 1.6 nm gold NP. Figure 4(b) demonstrates that the calculated coating thicknesses agree with the theoretical values, indicating that the results of present MD simulations resemble those expected from polymer physics.

Radius of gyration
Further insight into the structural properties and conformation of ligands in thiol-PEG-amine coatings can be obtained from the analysis of the radius of gyration of the ligands. The radius of gyration, R g , is commonly used to characterise the size of a coiled polymer chain; it is defined as the root-mean-square distance from all the atoms of the molecule to its center of mass [50]: Here M is the mass of the ligand molecule and r CoM is the position vector of its center of mass; m i and r i are the mass and the position vector of a constituent atom i, respectively.
In the present study, the radius of gyration for thiol-PEG n -amine (n = 2 − 14) molecules has been determined through complementary MD simulations of single thiol-PEG n -amine molecules solvated in water and thermalised at 300 K. For each molecule, the radius of gyration was calculated according to Eq. (1) by averaging over 50 frames of the simulated trajectory. The calculated values of R g for thiol-PEG n -amine molecules of different sizes are shown in Figure 5 by red symbols.
The calculated values of R g for thiol-PEG n -amine (n = 2 − 14) molecules have been compared to the radii of gyration for PEG n molecules (also known as polyethylene oxide, PEO) reported in the literature. Figure 5 presents the comparison of our data with those obtained from other MD simulations [26,[51][52][53] and experimental data obtained from static light scattering measurements [54]. In the latter case, the following scaling law between the radius of gyration (given in nanometers) and the molecular weight M w of a polymer molecule was proposed: R g = 0.0202M 0.55 w [54]. The green crosses in Figure 5 describe the values of R g calculated using the indicated scaling law for the molecular weights of thiol-PEG n -amine ligands listed in Table I and taken from Ref. [25]. Overall, the R g values determined in the present study agree nicely with the results reported in the literature. Open blue symbols in Figure 5 show the results of a recent study [26] where the R g values were determined for different PEO n (n = 12, 20, 36, 45) molecules on the basis of MD simulations. The figure shows the dependencies of R g on the number of monomers n, obtained in the present study in the range of n = 2 − 14 and in Ref. [26] for n = 12 − 45, follow a similar trend. Fitting of the combination of the present R g values and those from Ref. [26] yields a power law dependence R g ∝ n ν , where ν ≈ 0.594 (see the dashed curve in Figure 5). This result is in agreement with the power law dependence of R g ∝ n 0.549 determined through MD simulations for polyacrylic acid (PAA) polymers [55]. The exponent in these equations is known as the Flory Exponent ν and defines the type of polymer chain described. In the case of a real polymer chain, ν = 3/5 [56], which is close to the value ν = 0.594 determined in the present study.

Ligand density for the mushroom to brush transition
Using the values of R g determined in the previous section, one can evaluate the ligand surface density of thiol-PEG n -amine molecules at which the transition from "mushroom"-like to "brush"-like structures occurs for the ligands of different sizes. As follows from the polymer theory [39], at low ligand surface density, σ < ∼ 1/R 2 g , the radius of gyration of the ligands does not exceed the spacing between the neighboring ligands attached to a surface. In this case, individual ligands do not interact with adjacent ligands, thus acquiring "mushroom"-like shapes. Consequently, the coating thickness in this regime is equal to L ≈ 2R g [36,39,57]. When the "mushroom"like structures begin to overlap (which happens at surface densities σ > ∼ 1/R 2 g ), a "mushroom"-to-"brush" transition occurs, and the ligands acquire "brush"-like shapes.
The transition ligand density between the "mushroom" and "brush" regimes has been estimated from the results of performed MD simulations by analysing the probability distributions of the terminal nitrogen atoms within the NP coatings, similar to the results plotted in Fig. 3. In this analysis, we have summed up the probabilities of finding nitrogen atoms at distances r < 2R g and r > 2R g from the NP metal surface for various surface densities of ligands. In the "mushroom" regime, the height of the polymer coating is approximately equal to twice the radius of gyration, h ≈ 2R g [57]. Hence, by taking the sum of the nitrogen atom probability in the range of radial distances R c ≤ r ≤ R c + 2R g (where R c is the radius of the metal NP core), the transition ligand density has been evaluated for each particular case. It has been assumed that the transition density corresponds to the density at which less than 10% of nitrogen atoms are located at distances r < 2R g and more than 90% of nitrogen atoms are located at r > 2R g from the metal core surface. Figure 6(a) illustrates this method for the case of the 1.6 nm gold NP coated with thiol-PEG 8 -amine at different surface densities. The transition ligand density calculated this way for 1.6 and 3.7 nm NPs coated with thiol-PEG n -amine ligands of different lengths is plotted in Figure 6(b). The discrepancy between the data points of the estimated transition surface density and the 1/R 2 g dependence (see the dashed curve) can be attributed to the nature of this transition. As shown in Fig. 3, at lower surface densities, the coating structure is not exclusively "mushroom" but rather a mixture of both "mushroom" and "brush"-like conformations.

B. Model for the number density of thiol-PEG-amine coating
The results of the present MD simulations for a 3.7 nm gold NP and the earlier results for a 1.6 nm gold NP [25] indicate that radial dependence of the atomic number density of the thiol-PEG n -amine coatings has several common features for different NP sizes, ligand lengths, and surface densities of ligands. As such, one can develop a model that describes the number density profile of the coating for different parameters of the coating. The considerations presented below in this section concern the range of intermediate and high ligand surface densities (σ > ∼ 4.0 nm −2 ) when the ligand molecules attain "brush"-like shapes.
The atomic number density profiles for different thiol-PEG n -amine coatings shown in Figure 1 can be split into three regions, which are illustrated in Figure 7 for the case of 216 thiol-PEG 8 -amine ligands (surface density σ = 5.0 nm −2 ) attached to a 3.7 nm gold NP.
The first region proximal to the NP surface (denoted as Region I; see the dashed red curve in Fig. 7) is characterised by a rapid increase in the atomic number density of coating atoms up to the maximum value. Figure 1 shows that for a particular NP size, the gradient of this rise is approximately the same for different values of the ligand surface density. As such, the radial distribution of number density in this region is determined by the distribution of ligands over the gold NP surface. Figure 8 shows the radial distribution of the number density of sulphur atoms in the thiol group of thiol-PEG-amine ligands and the number density of gold atoms in the metal core. The distribution of sulphur atoms is centred around the average radius of the NP, R c = 1.85 nm, while the num- ber density of gold atoms starts to decrease at smaller distances (of ∼1.6 nm) from the CoM. This indicates that the metal NP core, after the annealing in a vacuum followed by solvation in water and thermalisation at 300 K, becomes non-spherical. Indeed, the NP surface has many facets, and sulphur ligands are bound across the gold surface based on the NP shape. Therefore, the distribution of ligands over the NP surface is directly linked to the distribution of the surface gold atoms of the NP core. As such, it is expected that atomic number density profiles for Region I will depend on the NP size. The maximum number density for atoms of the coating can be determined using the following geometric model outlined in Ref. [25]. The volume occupied by each repeating structural unit in thiol-PEG-amine ligands can be approximated as a cylinder with height a (that is a characteristic monomer length) and base radius r 0 . Considering the NP core as a sphere of radius R c , the number of ligands that can be distributed across the surface, N lig , is determined through the ratio of the NP core surface area, A core , to the cross-sectional area of a single ligand, A lig , The number density of atoms in the coating is given by the number of atoms in spherical bins of constant thickness ∆r divided by the volume of each bin. There are, on average, eight atoms per monomer in a thiol-PEG-amine ligand. Taking the bin thickness equal to the characteristic linear size of a single monomer, ∆r = a, one obtains the following expression for the maximum number density of atoms in the coating: Using a monomer length of a ≈ 0.25 nm and crosssectional radius r 0 ≈ 0.26 nm, one obtains the maximum number density n max ≈ 132 nm −3 for a 3.7 nm NP and n max ≈ 112 nm −3 for a 1.6 nm NP, which agree with the data plotted in Fig. 1. For Region I, the dependence of atomic number density on the radial distance from the CoM can be described as part of a sigmoid function, where L 1 defines the maximum value of the sigmoid, equal to the maximum number density for a given ligand surface density. k 1 defines the gradient of the sigmoid and is, therefore, a function of the distribution of sulphur and gold atoms (see Fig. 8) that depend on the NP radius but not on the number of monomers. r ′ 1 defines the midpoint of the sigmoid. The analysis performed in this study suggests that all parameters that define Region I are dependent on both the NP size and surface density of ligands, but not the number of monomers in ligand molecules.
In Region II (see the yellow dashed curve in Fig. 7), the coating number density has a power law dependence on the radial distance from CoM, Here L 2 is a scaling parameter defined by the NP size and ligand surface density, and α ≈ 2. The n 2 (r) ∝ 1/r 2 dependence stems from the distribution of the carbon and oxygen atoms that make up the backbone of the thiol-PEG-amine ligands [25]. Figure 9 shows the number of carbon, oxygen, and nitrogen atoms in spherical bins of 0.1 nm thickness in the coating made of 216 thiol-PEG 8amine ligands (surface density σ = 5.0 nm −2 ) attached to a 3.7 nm NP. The distribution of carbon and oxygen atoms in such a coating is approximately uniform (with the relative variation of ∼20%) at distances from ∼2.0 to 4.0 nm from the CoM, while the surface area of each spherical bin increases as r 2 . The approximately constant number of atoms in each bin results in the 1/r 2 decrease in the number density. The atomic number density for all ligand lengths and surface densities considered in the present study and Ref. [25] was fitted with the power law function, Eq. (5), in the range of distances from CoM where the radial distribution of carbon and oxygen atoms is approximately uniform. The resulting values of the exponent α are 1.97 ± 0.25 for the 1.6 nm gold NP and α = 2.03 ± 0.35 for the 3.7 nm NP. This result indicates that the n 2 (r) ∝ 1/r 2 dependence is a general structural feature of the studied polymer coatings, independent of the NP size and length of ligand molecules.
Region III (see the purple dashed curve in Fig. 7) is characterised by the appearance of a shoulder in the radial distribution of the coating number density. This shoulder is located at the position of the terminal NH 2 group of the thiol-PEG-amine ligands (see the distribution of nitrogen atoms in a thiol-PEG 8 -amine coating made of 216 ligands in Figure 9). The distribution of nitrogen atoms shown in Figure 9 is centred around a distance from the CoM of ∼4.0 nm, which corresponds to the position of the shoulder in the number density profile shown in Fig. 1(a).
This shoulder appears at ligand surface densities σ > ∼ 4.0 nm −2 and, therefore, is also related to the conformation state of ligand molecules within the coating. Higher density coating results in a "brush" conformation due to electrostatic repulsion between densely packed ligands, resulting in the position of the terminal NH 2 groups being extended to larger distances from the CoM (see Fig. 3). As demonstrated in Figure 1, the position of the shoulder is shifted to larger radial distances from CoM as the ligand surface density increases. Note that no such shoulder has been observed in the radial distribution of coating number density at lower ligand densities σ < ∼ 4.0 nm −2 (see Fig. 1). This behaviour is explained by the "mushroom"-like conformation of ligand molecules at low ligand surface densities so that the terminal NH 2 groups are scattered more uniformly over a broader range of distances from CoM.
The formation of a shoulder in the number density distribution can be understood from the random walk of the free end of a ligand molecule. For ligands attached to a curved metal NP surface, each ligand molecule can be described as a flexible polymer with one of its ends being fixed at the tip of a cone-shaped area.
The problem of a polymer constrained inside a cone was studied theoretically using the random walk theory [58]. In the cited study, the probability density for a polymer chain end to be located at a given point in space after a given number of random steps has been determined. It was found that the probability density as a function of the radial coordinate is described by a Gaussian distribution [58]. Based on these considerations, the radial dependence of atomic number density for Region III can also be characterised by a part of a sigmoid function: Here L 3 is the maximum value of the sigmoid function (equal to the lowest number density for Region II); k 3 depends on the width of the distribution of nitrogen atoms and, therefore, it becomes steeper with increasing surface density; r ′ 3 is the sigmoid's midpoint. In contrast to Regions I and II, the parameters characterising Region III depend on the NP size, ligand surface density, and the number of monomers in a ligand molecule.
The three regions of the number density profile introduced above can be combined into a piecewise function to describe the radial dependence of the atomic number density of ligands in terms of several geometrical parameters of the system: Here R 1 and R 2 are the radial distances from CoM defining the borders of the regions (see black dashed lines in Fig. 7), which depend on the NP radius; apart from that, R 2 is dependent on the number of ligand monomers. The continuity of the resulting function can be achieved by combining the piecewise function, Eq. (7), with a set of linking functions f L which take the form [59] f L (x) = 1 − 1 Here a is the point where the two functions intersect, and b is a parameter that makes the joining of the linked functions smoother. The resultant continuous function is then given by where f 1 and f 2 are the two parts of the piecewise function, Eq. (7), that are being linked. To link a third part of the function, one uses another linking function along with the output of the first two linked functions. Linking all parts of the piecewise function, Eq. (7), yields the following continuous function for the radial dependence of the atomic number of density of the coating: where n 1 (r), n 2 (r) and n 3 (r) are given by Eqs. (4), (5) and (6), respectively.

C. Structural characterisation of thiol-PEG-amine ligands grafted to a flat gold surface
It is also interesting to investigate how the structural properties of thiol-PEG-amine coatings on top of a gold surface depend on the system's geometry. For that, a similar analysis to the one presented above in Section III A has been performed for the thiol-PEG-amine coatings attached to a flat gold surface. The process for calculating the atomic number density of the coating on a flat substrate is much the same as that for coated gold NPs. The number of coating atoms in slabs of 0.1 nm thickness, starting from the gold surface, was determined and then divided by the volume of each slab. Figure 10 illustrates the number density of thiol-PEG 8 -amine and thiol-PEG 14 -amine coatings as a function of the distance from the substrate for three values of ligand surface density, σ = 1.0, 2.3 and 7.0 nm −2 . Figure 10 displays several notable features. First is the presence of a series of peaks in the number density distributions for ligand surface densities of 7.0 nm −2 (blue curves). Each of these peaks corresponds to an individual ethylene glycol monomer unit, implying that no water has penetrated between the ligands. In the lower density cases, σ = 1.0 and 2.3 nm −2 , the number density profiles resemble those shown in Figure 1 for similar ligand surface densities.
It should be noted that the in the case of a flat surface, there is no power law dependence of n(r) on r, as described by Region II for the case of a metal NP. As such, the power α characterising the power law behaviour in Region II should be related to the geometry of the metal surface, α ≈ 2 for the spherical geometry and α ≈ 0 for a flat geometry.

D. Distribution of water around coated NPs
Analysis of water distribution near the surface of metal NPs is one of the important aspects for evaluating the ability of coated metal NPs to act as radiosensitising agents. Figure 11 presents the number density of water molecules near the surface of a 3.7 nm NP coated with thiol-PEG 8 -amine ligands at ligand surface densities σ = 1.0 − 7.0 nm −2 . In the region close to the NP surface, the number density distribution of water molecules differs significantly for the different values of ligand surface density. At the lowest ligand density considered, the number density distribution of water molecules has a maximum, and the amplitude of this peak decreases with increasing ligand surface density. Such behaviour has a simple geometrical explanation: As ligand density increases, there is no space left for water molecules to penetrate the coating region toward the NP surface. The amount of water in proximity to the NP surface is essential for the production of hydroxyl radicals [10,27], which may have a strong impact on the radiosensitising effectiveness of metal NPs.
For the systems considered in this study, the coating region around the metal NP comprises only thiol-PEGamine and water molecules. Therefore, by determining the radial dependence of the number density of thiol-PEG-amine, one can obtain the number density of water based on the excluded volume principle [60]: n water (r) = 1 − n PEG (r) n max PEG n a water . Here n PEG (r) is the atomic number density of thiol-PEGamine ligands as a function of the radial distance from the CoM, n max PEG is the maximum number density of thiol-PEG-amine given by Eq. (3), and n a water = 33.4 nm −3 is the number density of water molecules at ambient conditions. Figure 12 compares the number density profile of water molecules determined from MD simulations (symbols) with that evaluated using Eq. (11) (dashed curves). The evaluated water density profile depends on the maximum number density of ligands, n max PEG , estimated using Eq. (3). The latter, in turn, depends on r 0 , the cross-sectional radius of a cylinder surrounding a ligand molecule, which should approximately be equal to a characteristic van der Waals radius of atoms of the ligands. Setting this value to r 0 = 0.26 nm provides a good correspondence to the simulation results. Therefore, the main features of the water number density profile obtained from the atomistic MD simulations can be evaluated with good accuracy by means of Eq. (11).

IV. CONCLUSION
This study has been devoted to atomistic modelling and structural characterisation of a 3.7 nm gold NP coated with thiol-PEG n -amine molecules of different lengths (n = 2 − 14) and solvated in water. The system size and composition have been selected in connection to several radiosensitisation experiments [28,31]. The coating structure and water distribution near the NP surface have been characterised on the atomistic level using the MBN Explorer [29] and MBN Studio [30] software packages. It has been found that the radial distribution of the number density for coating molecules of different lengths and surface densities is general between different NP sizes and can be divided into three spatial regions. The first region is characterised by the distribution of ligands over the NP metal surface and, therefore, the surface geometry and size of the NP. The second region is characterised by the distribution of the carbon and oxygen atoms of the PEG backbone and is characterised by inverse square dependence on the distance from the NP center of mass. The third region is characterised by the specific conformation state of the thiol-PEG-amine coating, with ligands acquiring either "mushroom"-like or "brush"-like shapes depending on the surface density of ligands attached to the NP.
The results of MD simulations presented in this study, in combination with the results of our recent study [25] and those from the field of polymer physics, have been used to calculate key structural parameters of the coatings of radiosensitising gold NPs and establish connections between the coating structure and distribution of water for different NP sizes as well as lengths and surface densities of coating molecules. The results of the quantitative analysis of water distribution in the vicinity of coated metal NPs can be used in follow-up studies to evaluate the radiosensitising effectiveness of a particular NP system based on the proximity of water to the NP metal core, which should impact the production hydroxyl radicals and reactive oxygen species in the vicinity of metal NPs.