Atomistic simulations of gold surface functionalization for nanoscale biosensors applications

A wide class of biosensors can be built via functionalization of gold surface with proper bio conjugation element capable of interacting with the analyte in solution, and the detection can be performed either optically, mechanically or electrically. Any change in physico-chemical environment or any slight variation in mass localization near the surface of the sensor can cause differences in nature of the transduction mechanism. The optimization of such sensors may require multiple experiments to determine suitable experimental conditions for the immobilization and detection of the analyte. Here, we employ molecular modeling techniques to assist the optimization of a gold-surface biosensor. The gold surface of a quartz-crystal-microbalance sensor is functionalized using polymeric chains of poly(ethylene glycol) (PEG) of 2 KDa molecular weight, which is an inert long chain amphiphilic molecule, supporting biotin molecules (bPEG) as the ligand molecules for streptavidin analyte. The PEG linkers are immobilized onto the gold surface through sulphur chemistry. Four gold surfaces with different PEG linker density and different biotinylation ratio between bPEG and PEG, are investigated by means of state-of-the art atomistic simulations and compared with available experimental data. Results suggest that the amount of biotin molecules accessible for the binding with the protein increases upon increasing the linkers density. At the high density a 1:1 ratio of bPEG/PEG can further improve the accessibility of the biotin ligand due to a strong repulsion between linker chains and different degree of hydrophobicity between bPEG and PEG linkers. The study provides a computaional protocol to model sensors at the level of single molecular interactions, and for optimizing the physical properties of surface conjugated ligand which is crucial to enhance output of the sensor.


Introduction
Nanoscale biosensors [1] have emerged as one of the promising techniques in clinical analysis associated with bio-Nanotechnology Nanotechnology 32 (2021) 095702 (9pp) https://doi.org/10.1088/1361-6528/abc6dc Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. molecular recognition; particularly because of its sensitive yet label free detection. Apart from effective spatial and chemical specificity of the substrate coverage [2,3] followed by suitable immobilization of the bio-receptor [4][5][6], the underlying interfacial phenomena can also yield guidelines to optimize performance of a biosensor. Due to the experimental complexity to investigate the sensor [7,8] at the atomistic level, computational studies can be considered as an alternative treatment to understand potential role of surface chemistry in biosensors. Nevertheless, owing to complex architecture of a biosensor, extensive atomistic simulation incorporating explicit contribution from both the surrounding and the 'substrate-receptor-analyte' conjugate is limited till date. In a recent study, Lowe et al [9] provided proof-of-concept that molecular dynamics (MD) simulation can be used as a complementary tool to experiments for potentiometric biosensor response prediction, by performing extensive MDs of silica-electrolyte-biomolecule interface. On the other hand, Pang et al [10] provided guideline for the performance optimization of the sensors, using first-principle calculations and MDs to disclose the role of temperature fluctuation, along with material thickness on a graphene/Au based surface plasmon resonance biosensor. Similarly, Cholko et al [11] used MDs to explain the implications of the behavior of ssDNA targets near differently functionalized gold surface (hydrophobic, hydrophilic, and oppositely charged surfaces) for nucleic acid biosensing. The potential role of modification of different surface properties of carbon nitride materials in accordance with tuning the structural properties of a genetically engineered metalloenzymes is being highlighted using atomistic simulations by Puente-Santiago et al [12], leading to improvement in their electron-transfer efficiency for the fabrication of highly sensitive amperometric biosensors. Ghoshdastider et al [13] reported atomistic simulations to unravel the nature of the enzyme-graphene interface for the design of graphene-based biosensors. The design of a biosensor by means of atomistic simulations is suggested by Ogorzalek et al [14], identifying biomaterials with high thermal stability. Wang et al [15] also showed that it is possible to understand at the molecular level the interactions and conformation transition of ssDNA upon binding with graphene by MD simulations.
Moreover, recent computational findings suggest that MD simulations can disclose [16][17][18] the favorable orientation and binding location of the analyte in aptasensors, namely sensors which utilize aptamers (single stranded DNA or RNA oligonucleotide) as bioreceptors, as well as the anisotropic distribution of ions [17] which is found to contribute towards the binding affinity of the sensors.
In the present paper, the sensitivity of an acoustic wave (AW) sensor, namely a quartz-crystal-microbalance sensor, is addressed from computational perspectives. To the best of our knowledge this is the first time that the optimization of the analyte binding region of an AW gold-surface biosensor is pursued, by means of atomistic simulations. AW based technology is widely used in biosensors [19,20], where the transduction mechanism is manifested in terms of propagation of AW along surface of the substrate. Any kind of change in physico-chemical environment or any slight variation [21,22] in mass localization near the surface of the sensor can cause differences in nature of AW. Functionalization of transducer surface with proper bio conjugation element should [21,22] maximize coverage of the surface, while minimizing nonspecific bindings for the bio-receptor. Moreover, the detection of the bio-molecule is known to vary as a function [23][24][25][26][27][28] of the thickness of the chain length as well as the flexibility of the linker groups. Thus, optimizing the physical properties of surface conjugated ligand is crucial to enhance output of the sensor. Revealing the precise density of the attached ligands with discrete knowledge of their distribution pattern over the transducer could be tricky for experiments. In this context, atomistic simulations can be considered as an alternative treatment to understand potential role of surface chemistry in biosensors.
In a recent experimental work [19], the detection of streptavidin on a Rayleigh surface AW based resonator biosensor has been achieved involving the immobilization of biotinylated polyethylene glycol (bPEG) on the gold surface of the electrode with a detection limit of 104 picomolar (pM) and a normalized sensitivity of −296 m 2 kg −1 . PEG is a hydrophilic polymer composed of the subunit -(CH2-CH2-O) n -, which has been often covalently attached to surfaces, with a process called PEGylation; its main goal is to passivate the underlying inorganic surface to prevent unspecific binding of biomolecules [1,2].
However, the molecular level understanding of the experimental evidence still remains elusive. Streptadivin [29][30][31], the homotetrameric bacterial protein has strong non covalent binding affinity towards biotin, a B class Vitamin, and thus the protein-ligand system plays a pivotal role in molecular biology and biotechnological applications. The protein has multiple biotin binding sites and can get adsorbed as a top layer on functionalized surface, which can further bind additional biotinylated ligands through their free binding sites [32,33]. In principle, a range of different conformations are possible for the immobilized ligand due to the random distribution of biotin sites towards water. In this work, we disclose the effects of different surface coverage densities on biotin site distribution along with its exposure to water, by tuning the surface coverage densities and linker biotinylation ratio. We provide proof-of-concept that MD can be used as a complementary tool for AW biosensor response prediction. Effects that are conventionally neglected due to inherent atomistic complexiety, such as discrete surface morphology, biotins exposure to water and surface chemistry dynamics, are explicitly modeled.

(PEG)n and biotin−(PEG)n in water set-up
A PEG molecule of 2 KDa molecular weight, i.e. with n=45 repetitive subunits -(CH 2 -CH 2 -O) n -, was built using Schrödinger (Maestro) software [34]. The topology and force field (FF) parameters were generated via Topolbuild 1.2.1 software [35] to be consistent with standard OPLS/AA (all atom) FF [36]. The same procedure was applied for the preparation of the bPEG, after downloading the crystal structure of biotin (PDB: 2BDO) [37] and covalently attaching it to the PEG chain via the peptide bond.
200 ns MD [38,39] were performed for both PEG and bPEG molecules in water with OPLS and explicit SPC/E [40] water model implemented in GROMACS-5.0.4 version [41,42]. The PEG system was simulated in a rectangular box of dimension 4.2×4.2×20.7 nm 3 and it consisted of 318 PEG atoms (figure S1(a), available online at stacks.iop.org/ NANO/32/095702/mmedia) and 12 062 water molecules (figure S1(b)). For bPEG, the set-up comprised of 348 of atoms of bPEG (figure S1(c)) and a total number of 35 608 water molecules (figure S1(d)). In this latter case, the box dimensions were 8.9×5.9×21.1 nm 3 respectively. Periodic boundary condition [38,39] was applied to all simulations. The systems were firstly minimized with the steepest descent algorithm [43] for 50 000 steps. The classical equations of motion was integrated employing leap-frog algorithm [38,39] by using a time step of 2 fs. The temperature was maintained at 300 K in terms of V-rescale modified Berendsen thermostat [44], whereas for pressure coupling (1 bar) Parrinello-Rahman barostat [45] was used. NVT (isothermal isochoric) and NPT (isothermal isobaric) equilibration were performed for 100 pico second (ps). Short range cut off for van der Waals as well as for electrostatic interaction ∼1.2 nm, were applied. Long range electrostatic interaction was taken into account via particle mesh ewald method [46]. Bond lengths were constrained using LINCS algorithm [38,39] and velocities were assigned as per Maxwell Boltzmann distribution at prescribed temperature. Trajectories were saved in every 2 ps. All the analyses were performed on 200 nano second (ns) long isothermal isobaric simulated trajectories.

(bPEG) n -linked gold surfaces set-up
A crystalline slab of atomic gold with the (111) facet exposed, namely Au(111), was constructed, with surface dimensions of 100×100 Å 2 . The slab was four atomic layers thick in the direction normal to the surface and was held rigid during the simulation. Au(111) surface was chosen as the most stable and the most widely targeted gold facet by proteins and peptides [47,48]. Experimental [19] surface density of PEG molecules on top of Au(111) surface was provided, namely 100-200 ng cm −2 corresponding to 0.6 molecule nm −2 . In order to mimic the experimental set-up, two atomistic molecular mechanics models were investigated having a surface density of PEG molecules of 0.4 molecule nm −2 and 0.8 molecule nm −2 , respectively. For each system the effect of having 100% of biotin-PEG (bPEG) as well as, 50% of b-PEG and 50% of PEG, was studied. To summarize, we considered four systems with two different surface functionalization density and different biotinylation ratio: System 49-bPEG. For the lowest density system with 0.4 molecule nm −2 , one bPEG was attached to each gold slab of 15 × 15 Å 2 unit surface area and the unit slab was replicated using gmx genconf command of GROMACS-5.0.4 version [41,42] to construct a surface of 100 × 100 Å 2 .
System 25-bPEG/24-PEG. On top of system 49-bPEG, a surface with the same density of PEG molecules but including a mixture of 50% of bPEG and 50% PEG, was constructed with a 100 × 100 Å 2 surface area.
System 98-bPEG. For the highest density system with 0.8 molecule nm −2 , two bPEG, instead of one, were attached to the gold slab unit of 15 × 15 Å 2 surface and further replicated to obtain the surface of 100 × 100 Å 2 .
System 49-bPEG/49-PEG. On top of system 98-bPEG, a surface with the same density of PEG molecules but including a mixture of 50% of bPEG and 50% PEG, was constructed with a 100 × 100 Å 2 surface area.
In the case of bPEG and PEG mixture, different head groups were homogeneously distributed on top of Au(111) surface.
The description of the Au(111) surface atoms was based on the GolP all-atom classical FF for the protein-Au(111) interaction developed and extensively tested by our group [47][48][49][50][51][52][53]. GolP quantitatively reproduced experimental adsorption energies for small molecules on Au(111) [47] and allowed the prediction of binding affinity scale for amino acids on Au coherent with available experiments [49]. SPC water was used being compatible with OPLS/AA, thus with GolP. In GolP, polarization of the gold surface atoms is taken into account.
The covalent chemical bond between the gold surface atom (AU) and the S atom of the terminal thiol groups of the linkers (either bPEG or PEG) were created using Schrodinger (Maestro) software. The AU-S bond length were restrained, while AU-S-C angle and dihedral AU-S-C-C were free to rotate during the simulations. Bonded and non-bonded LJ FF parameters, were taken according to a previous study [54].
The modeled system was simulated for 500 ns in explicit water model. The box dimensions of the systems were 10.7×10.3×22.0 nm. The analyses were reported for last 400 ns trajectories over 500 ns for each system.

Results and discussion
3.1. MD simulations of (PEG) and (bPEG) in water PEG and bPEG polymers were freely equilibrated in water and the RMSD indicated that the systems were well equilibrated upon the time scale of the simulation (figures S2(a)-(b)). A structure representative of the final conformational ensemble of PEG after simulation is depicted in figure S2(c). The mean value of the radius of gyration for the PEG system resulted equal to∼1.1 nm, (figure S2(d)) in good agreement with experimentally available data [57]. The final conformation of bPEG is shown in figure S2(e). Overall the PEG chain retained similar conformations upon biotin attachment with respect to its native conformation, being the mean value of the radius of gyration of bPEG∼1 nm, as also supported by available NMR data [58]. Additionally, the RMSD of biotin computed with respect to its crystal structure, (figure S3) indicated that the structure of biotin is conserved throughout the entire simulation.  figure 1(a) a top view and in figure 1(B) a side view of the surface structure, is depicted, respectively. RMSD of bPEG generated over the entire simulated trajectory is reported in figure S4(a) of SI.
We remark that the oxygen end-atom of biotin, namely O3 atom, is known [29,31] to form hydrogen bond with the Streptavidin protein, thus representing one of the key element for sensitivity and selectivity of the biosensor. The dynamics of O3 end-atom of bPEG is analyzed through the entire MD to disclose the influence of the PEG chain on biotin exposure to water (and consequently its availability towards the streptavidin analyte protein). In figure 2(a) the distance between O3 end-atom of each biotin and AU of surface gold atoms ( -d O3 AU ) are reported over simulated data and the distance distribution (H( -d O3 AU )), is plotted. This analysis shows that biotin molecules arrange on top of gold surface, assuming different ranges of distances of the end-group towards Au(111) ( figure 2(a)). The distribution of distances are clustered in four groups: (i) set1 in which biotins are largely distant from gold, lying on an average distance of 10-30 Å from Au(111), (ii) set2 in which biotins are very close to the surface, namely with a mean distance of ( -d O3 AU )∼5 Å, (iii) set3 in which biotins lie in an intermediate layer with respect to set1 and set2, namely with an average distance of ( -d O3 AU )∼7 Å and finally (iv) set4 in which biotins exhibit random fluctuations, undergoing back and forth motion between set2 and set3.
In order to correlate the distance fluctuations of the endatom of biotin with the biotin solvent exposure, the SASA was evaluated for each single biotin along the entire trajectory. Using g_sasa command of GROMACS it was possible to collect the SASA of biotins (SASA btn ) and to plot the corresponding histogram (H(SASA btn )). More details on the SASA calculation can be found in the SI. We provided an estimation that the SASA of a fully solvent exposed biotin is equal to∼4 nm 2 , then in figure 2(b) it is shown that biotins  belonging to set1 are maximally exposed to solvent, even if with a broad range of possible values of SASA btn ( figure 2(b)), together with a small fractions of biotins belonging to intermediate layer (set3) which exhibit a similar trend. However, the remaining biotins experienced very small value of SASA btn as depicted by the unimodal nature of H(SASA btn ). Thus, biotin end-atom distance fluctuation resulted to be highly correlated to the relative magnitude of SASA for that biotin. The largest the distance, the largest the magnitude of SASA.
Several studies [24,59] on atomistic and coarse grained simulation of PEG demonstrate that this polymer is prone to form a globular shape in aqueous medium. Thus, in figures 2(c)-(d), the structural compactness of the entire linker was investigated by evaluating the radius of gyration of the single bPEG chains from atomistic simulations, in particular the z component namely R .
g bPEG Initially, both the bPEG/PEG chains are elongated along z axis, the longest dimension of our simulation box; during the course of simulation the linkers achieve crumbled conformation which is significantly reflected in anisotropic nature of radius of gyration along z axis. Thus, we primarily focus on addressing z component of radius of gyration vector in our analysis. The obtained distribution of values were plotted in figure 2(c) (H (R g bPEG )) and they were correlated to the corresponding SASA (SASA bPEG ) in figure 2(d). From figure 2(c) we find that the bPEG with approximately a globular shape but less compactness having a mean value of R g bPEG ∼1.7 nm, experienced a broad distribution of SASA values, as reported for (H(SASA bPEG )) in figure 2(d). On the contrary, the comparatively rigid bPEG linkers with a small mean value of R g bPEG ∼0.7 nm, reported small value of SASA . bPEG This correlation between radius of gyration and SASA highlighted the importance of globular shape of the bPEG linkers. Moreover, the results resemble [60] enhancement of SASA for series of globular proteins in accordance with increased value of their radius of gyration. In most of the cases, once the terminal group of the linker was bound to the Au(111), the biotins attached to those linkers remained buried from the solvent for the rest of the simulation. As we shall discuss in the next sections the collective entanglement as well as the compactness of the linkers are strongly influenced by the variation of the surface density (i.e. grafted linker spacing) of bPEG.

System 25-bPEG/24-PEG.
In order to understand whether the surface chemistry of the biosensor could be affected by the biotinylated ratio of the linkers, a surface with the same surface density of linkers as in section 3.2.1 but with a 1:1 ratio of biotinylated and non-biotinylated linkers was prepared and simulated, and the system was called 25-bPEG/ 24-PEG. In figure S4(b) we report RMSD for the entire trajectory. In figures 1(c)-(d) the final structure of the surface after relaxation are reported from different viewpoints.
For this system it was found that in analogy with System 49-bPEG, there existed an amount of biotins attached to the 25 bPEG linkers, which remained at large distances from the gold surface, and those were also the more exposed to water (figures S5(a)-(b)). Instead, the group of biotins lying close to AU∼5Å, remained buried from the solvent, displaying very low values of SASA. Overall, the amount of non-biotinylated linkers, PEG, did not affect the atomistic behavior of the biotinylated linkers. Further analysis, comparing the corresponding component of radius of gyration of both bPEG and PEG with SASA were conducted. Some of the bPEGs could be classified in two regimes, the first having R g bPEG <1 nm, the latter, having R g bPEG ∼2 nm ( figure S5(c)). The distribution of SASA for both biotinylated and non biotinylated linkers ( figure S5(d)) indicated that the ligands with larger value of radius of gyration are more exposed to water, similarly the more compact ones are buried close to the gold (surface). These distinguishable behaviors associated with the radius of gyration for both bPEG and PEG along with their SASA observed in System 25-bPEG/24-PEG is similar to that for 49-bPEG. However, due to the introduction of more hydrophilic PEG in this mixed system, some of the PEGs also showed a more extended, 'brush like' cylindrical conformation corresponding to smaller value of the radius of gyration ( figure S6(a)). Conversely, they showed larger values of SASA ( figure S6(b)). Therefore, due to the collective effects associated to the coexistence of rigid 'mushroom like' spherical bPEG linkers with more 'brush like' extended PEG, the fraction of linkers exhibiting smaller values of the radius of gyration is increased with respect to 49-bPEG (table S1). Conversely, the overall percentage of linkers with larger SASA is enhanced (table S2).

System 98-bPEG.
We comment here on the effect of increasing the density of bPEG linkers on top of the Au(111) surface. To investigate this aspect we built a system with a surface density which was two times larger with respect to 49-bPEG system, namely system 98-bPEG. The final snapshot in different angles are reported in figures 3(a)-(b) RMSD analysis of the bPEG linkers is reported in figure S4(c).
From the analysis of the plots reporting the distances between O3 end-atom of each biotin and AU of surface gold atoms ( -d O3 AU ) (figure 4(a)) we observed an important change in the range of distances experienced by the biotin end-atoms namely, most of the biotins were lying at larger distance from the gold surface (compare with figure 2(a)). The most distant biotins, set1, stabilized at ranges of distances d O3 AU ∼20-50 Å, whereas other groups of biotins, set3, remained confined between 10 and 20 Å. A third group of biotins instead was experiencing fluctuations back and forth between those two ranges of distances. The computed SASA values were comparable with those obtained for the previous systems ( figure 4(b)). More specifically, biotins belonging to set1 and a most from set2 reported broad values of SASA, H(SASA btn ). For the bPEG, the linkers atoms are distributed at a larger distance over gold surface, but the values of the radius of gyration and SASA (figures 4(c)-(d)) can still be grouped into three categories as remain consistent to system-49bPEG. The key difference is given by an increased number of linkers with a more extended, 'brush like' cylindrical conformation corresponding to smaller value of the radius of gyration (table S1) and figures S6(c)-(d).
The present results are in line with previous lower resolution simulations of PEG chains grafted to a non adsorbing surface [24]. To summarize, at very low linker density, the linker chains behaves like an isolated chain in solution, leading to the so called 'mushroom like' conformation. As the grafted linker density increases, a larger number of the polymer chains become crowded and repel each other, leading to a 'brush like' conformation with a larger distance of the end-groups from the surface. In particular, the scenario is a mixed state of 'mushroom-brush' conformations, however the end group of linkers exhibiting 'mushroom-like' shape are situated far from the gold surface with respect to 49-bPEG system, whereas some bPEGs are in elongated 'brushlike' form. Although, due to the highly crowded linker environment, the collective SASA of the linkers is reduced (table S2), the overall solvent exposure of biotin is increased (table 1) as explained in details in discussion.

Discussion of the results
Here we focus on the discussion of the collective behavior of the linkers' terminal groups, namely biotins. The SASA btn for each single biotins was computed over the simulation trajectories and the values were grouped in different SASA btn domains, corresponding to different degree of solvent exposure. Finally, after counting the number of individual contribution in each domains, the values were normalized by the total number of biotins as well as number of configurations. In table 1, System 49-bPEG and System 25-bPEG/24-PEG, were shown to exhibit more than 60% of biotins buried and unavailable towards solvent. Whereas, the number of partially available biotins is enhanced in System 98-bPEG and the maximum exposure is found for System 49-bPEG/49-PEG. In table 2 the distance between O3 end-atom of each biotin and the Au(111) atoms are grouped with respect to a threshold value∼10 Å for all the systems. From the reported values, a clear trend is outlined as a function of the linkers density, namely at low linkers density the majority of biotins are lying close to the gold surface, namely in System 49-bPEG and System 25-bPEG/24-PEG, whereas increasing the linkers density, the repulsion between linker chains is pushing the biotins end-group far from the gold surface, namely in System 98-bPEG and System 49-bPEG/49-PEG. In the case of mixed linkers 49-bPEG/49-PEG the different hydrophobicity between the bPEG and PEG chains, induced the more hydrophilic PEG linkers to assume a more extended conformations (tables S1, S2) at larger distances between the terminal end-atom (C90) and the gold surface, as reported in table 3. The radial distribution function between oxygen atoms of water molecules and the entire PEGs generated through the simulated trajectory is shown in figure S8. Similar analysis for bPEGs (figure S8) indicate that PEG is more hydrophilic than bPEG. This collective behavior seems to favor  the transition of a larger number of PEG linkers to the so called 'brush like' cylindrical conformation, increasing the amount of biotins of the bPEG linkers lying at larger distances from the surface, which is going from 92.1% to 97.2%, as reported in table 2. Thus, from our simulation we can propose larger number density of linkers and the suitable biotinylation ratio are the two trigger points that can optimize performance of a biosensor.

Conclusions
Fully atomistic MD simulations have been applied to reveal the structure and dynamics of functionalized gold surfaces of a AW biosensor, which cannot be easily captured by experiments. Recent advances in computer power and simulation methods allowed fully atomistic simulations of large PEGylated gold surfaces, functionalized with a different number and types of linkers.
The main goal of the paper was to clarify the role of the grafted linker density of bPEG chains on a gold surface and to disclose the influence of the pegylation ratio on the sensitivity of the biosensor. Results suggested that at lower linker density, namely system 49-bPEG and 25-bPEG/24PEG, the linker chains mostly maintained their high 'mushroom like' conformation typical of the isolated chain in solution.
Increasing the linker density from 0.4 to 0.8 molecule nm −2 , the repulsion between linker chains, cause the extension of the chains leading to an increased number of structural conformations with large distances of the end-groups from the gold surface, namely System 98-bPEG. This latter behavior favors the exposure of biotin and thus the sensitivity of the biosensor. Moreover, including a 1:1 ratio of pegylation in System 49-bPEG/49-PEG, the accessibility of biotin molecules towards solvent improve, due to the increased different hydrophobicity of the linkers chains introduced by the PEG polymer chains, which strongly favor the availability of biotins molecule towards the protein analyte.