Dynamic behaviour of the silica-water-bio electrical double layer in the presence of a divalent electrolyte

Abstract

Electronic devices are becoming increasingly used in chemical- and bio-sensing applications and therefore understanding the silica-electrolyte interface at the atomic scale is becoming increasingly important. For example, field-effect biosensors (BioFETs) operate by measuring perturbations in the electric field produced by the electrical double layer due to biomolecules binding on the surface. In this paper, explicit-solvent atomistic calculations of this electric field are presented and the structure and dynamics of the interface are investigated in different ionic strengths using molecular dynamics simulations. Novel results from simulation of the addition of DNA molecules and divalent ions are also presented, the latter of particular importance in both physiological solutions and biosensing experiments. The simulations demonstrated evidence of charge inversion, which is known to occur experimentally for divalent electrolyte systems. A strong interaction between ions and DNA phosphate groups was demonstrated in mixed electrolyte solutions, which are relevant to experimental observations of device sensitivity in the literature. The bound DNA resulted in local changes to the electric field at the surface; however, the spatial- and temporal-mean electric field showed no significant change. This result is explained by strong screening resulting from a combination of strongly polarised water and a compact layer of counterions around the DNA and silica surface. This work suggests that the saturation of the Stern layer is an important factor in determining BioFET response to increased salt concentration and provides novel insight into the interplay between ions and the EDL.

---

1 Introduction

Silica and water form some of the most abundant chemical systems and understanding the interface between the two is important for a large range of applications such as biosensing,^1–3 drug-delivery,⁴ prebiotic chemistry,⁵ improved fundamental understanding of geochemical processes (e.g. dissolution reactions⁶) and chemical engineering (e.g. water-desalination⁷). The precise structure and dynamics of this interfacial region, including the Electrical Double Layer (EDL) remains elusive, despite over a century of extensive study.^8–12

Addition of charged macromolecules to oxide surfaces results in a perturbation of the EDL which cannot be accurately described by conventional mean-field models. In this work, electrolyte and biomolecule dynamics were studied with atomistic resolution, providing a detailed description of the electric field generated at the interface. DNA was chosen as an example of a highly-charged macromolecular polyelectrolyte which is both well-characterised and has relevance to a range of biotechnology applications. In addition, divalent ions were included, which are known to have a strong influence on the structure of the EDL, important to silica dissolution processes¹³ and prominent in physiological solutions¹⁴ but despite this, have received surprisingly little attention in the atomistic simulation literature.

One application of this work is in improving understanding of the mechanism-of-action of a promising class of biosensors, termed biologically-sensitive field-effect transistors (BioFETs). These sensors operate by detecting changes in the electric field within the EDL as a result of biomolecule binding, as shown schematically in Fig. 1. Reliable and quantitative prediction of changes in the electric field due to biomolecule adsorption, and hence BioFET response, is currently difficult primarily due to the complexity of the EDL.

Fig. 1 Schematic of BioFET operation. Biomolecules can alter the electric field at the interface, resulting in a measurable change in conductance of the channel. Many factors: surface charge; biomolecule charge; biomolecule orientation; surface dipole; ionic strength; and pH, can affect the interfacial electric field.

Not only has the presence of divalent ions recently been shown to increase BioFET sensitivity,^15–17 but the electric field and ion-dynamics at the interface are thought to be crucial in determining field effect biosensor response.

1.1 Importance of EDL structure and ion dynamics in the interfacial region

The most commonly discussed hypothesis for BioFET response is via detecting changes in the electric field due to changes in the “surface” concentration of ionised groups forming the EDL or charges around biomolecules.§¹⁸ The interfacial region is thought to be significantly affected by biomolecules; for example the orientation of biomolecules is thought to be important in determining sensor response.^3,21 The EDL structure in turn can be affected by dense biomolecule layers through ion-exclusion; mathematical models incorporating ion-exclusion effects have been shown to describe experimental signal measurements of DNA hybridisation better than more conventional EDL models.¹⁵

Recent experimental work recognises the importance of ion dynamics in BioFET engineering, with deterministic information extracted from BioFET signals in the frequency domain of the response^15,22,23 which has been explained as a result of adsorption–desorption noise of biomolecules²⁴ and perturbed charge fluctuations in the EDL.²⁵ Experiments have shown a decrease in low frequency noise with increased ionic strength due to increased screening competition between the EDL and the semiconductor device.²⁶ Heitzinger et al. suggested a different trend for DNA-sensing, in which they calculated an increase in the standard deviation of the FET channel current with ionic concentration due to an increasingly variable orientation of the DNA.²¹ These studies show that addition of a biomolecule, such as DNA, can affect EDL dynamics to the extent that a response can be observed in the frequency domain that is not apparent in the time-domain, even under high ionic strength conditions.

Experiments are not able to unambiguously decouple the signal noise originating from the semiconductor device and the EDL region in the electrolyte solution. Most current EDL theories used in the BioFET engineering field are based on equilibrium, mean-field solutions of the Poisson–Boltzmann equation. These models offer the advantage of low computational cost and can be accurate for low ionic strength and low surface-charge systems. However, for BioFET systems, this is rarely the case and finite-size steric effects render the Poisson–Boltzmann equation inaccurate without modification. Modern advances in computational power have enabled the exploration of more detailed atomistic models of the structure and dynamics of the silica-water(-bio) interface²⁷via both classical^28–36 and ab initio Molecular Dynamics (MD).^37–41

This work presents MD simulations of EDL structure and ion dynamics in the interfacial region to investigate how (a) increased ionic strength and (b) addition of DNA perturbs the electric field and charge density at the silica-water interface.

1.2 Divalent ions

Physiological samples often contain divalent cations such as Mg²⁺ and Ca²⁺ which serve important biological functions. For example, diffusely associated divalent cations are thought to have a significant effect on reducing the internal stress in DNA/RNA due to screening of the negative charges on the phosphate backbone, as evidenced by experiment^14,42 and simulation.^43,44

Divalent ions are also known to be important in the phenomenon of charge inversion, in which the first diffuse layer in the EDL contains more counterions than needed to compensate for surface charge, which is then balanced by a second co-ion layer. This phenomenon has been attributed to two (non mutually-exclusive) mechanisms. One mechanism is via ‘specific adsorption’ of ions, caused by forces such as chemical bonding or water-mediated interactions.⁴⁵ The other mechanism is via many-body ion–ion correlations, in which the chemical potential near the surface is reduced due to spatial correlations between discrete ions, with the electrostatic interactions outweighing the entropic cost of forming such a highly correlated system.^45,46 Atomic Force Microscopy (AFM) has measured the effects of charge inversion at the silica-water interface using trivalent and quadrivalent ions at low concentrations (≤1 μM).⁴⁷ These experiments did not demonstrate charge inversion for Mg²⁺; however, it would be expected that divalent ions would require a higher concentration than trivalent ions. Edel and de Mello have measured the streaming current for silica nanochannels and Mg²⁺ counterions, observing charge inversion for concentrations exceeding approximately 400 mM.⁴⁸

From the perspective of BioFET biosensing, the importance of adding divalent and multivalent salts has also been shown. Jayant et al. have recently shown a significant enhancement in DNA-hybridisation sensitivity upon addition of trace amounts of multivalent salt (Mg²⁺ or Co³⁺) to a low-concentration NaCl background buffer.^15,49,50 In their work, a background of 1 mM NaCl was used for hybridisation of ssDNA, and it was found that addition of the complementary strand in 1 mM NaCl with 100 μM Mg²⁺ produced a 350 ± 150% increase in potential shift relative to the control. They found that this effect was only significant when the initial concentration of monovalent salt was low, which supports an ion-competition mechanism. Modeling the DNA as a membrane, they putatively assigned the signal to a combination of (a) increased ion-exclusion of multivalent ions and (b) increased DNA-condensation in the presence of multivalent ions.¹⁵

Other authors support the notion that multivalent salts can have a significant effect on BioFET signal, for example, Rica et al. have observed an increase in DNA hybridisation on addition of multivalent salt (spermidine) and observed a corresponding FET signal indicative of charge inversion at 10 μM spermidine. The signal changed polarity at higher concentrations of spermidine, which they attribute to increased screening of the charge inverted DNA molecule.¹⁶ Shul'ga et al. have reported an almost 100% increase in glucose-sensitive enzyme-FET signal on addition of 0.1 M MgCl₂ which they attribute to divalent cations affecting the rate of charge transfer of the enzyme substrate oxidation.¹⁷

Despite the importance of divalent ions, a monovalent electrolyte is typically assumed in MD simulations of hydrated surface-biomolecule systems or mathematical modelling of BioFET signals. Sakata et al.^28–31 have recently used MD simulations to investigate the EDL structure and dynamics for hydrated silica-water^28,29 and hydrated silica-DNA systems.³¹ This paper extends this work to investigate the effect of divalent Mg²⁺ ions upon the structure and dynamics of the EDL at this technologically important interface.

2 Computational methods

Fig. 2(a) shows a schematic diagram of the simulation domain used in this work. The solid base was modelled as the (100) surface of alpha-quartz (SiO₂) with dimensions of 49.130 Å × 54.050 Å, and a depth of 16.5 Å. At open-circuit potential and biosensing conditions (usually 5 ≲ pH ≲ 9), silica-water interfaces are negatively charged and therefore the upper surface was defined with a ratio of one fifth SiO⁻/SiOH (0.2 C m⁻²), as discussed further in the ESI 1.† The surface charge was then neutralised with Na⁺ to produce an electroneutral unit-cell and a solvent box was put into contact with this surface, similar to the method of Zhang et al.³⁶

Fig. 2 System summary. (a) Side-on schematic of the simulation cell. (b–d) Initial configuration of the 0 mM, 200 mM and 1 M ionic strength systems with DNA. Mg²⁺ = purple, Na⁺ = green, Cl⁻ = yellow, O = red, Si = gold. Mg-coordinated water drawn as transparent. DNA drawn as stick representation. Videos of each trajectory can be found in the ESI.†

Solvent boxes were prepared with an initial density of 1 g cm⁻³ and a vertical height of approximately 73 Å. Three different solvent boxes were considered: 0 M electrolyte (salt free, corresponding to deionised water); approximately 0.2 M ionic strength electrolyte; and approximately 1 M ionic strength electrolyte, each containing a 1 : 1 molar ratio of NaCl to MgCl₂. The system was geometry optimised for 5000 steps and then NVT dynamics⁵¹ were performed at 300 K for 3 ns. Dynamics were performed with the Nosé–Hoover thermostat using a Q ratio of 0.01.⁵¹

For comparison, three further systems were prepared at each ionic strength, incorporating DNA neutralised with Na⁺. DNA was constructed and chemically bonded to the surface following the method of Maekawa et al.³¹ and consisted of a d(AAAAAAAAAA) decamer with a complementary base-T strand and a net charge of −19e. The DNA was superimposed onto the solvent box, the DNA was kept fixed and a 5000 step geometry optimisation was performed for the three cases. Then, water molecules within the DNA were removed and electrolyte ions from inside to outside of the DNA. A further 5000 step geometry optimisation and 100 ps of NVT dynamics were performed in order to further relax the system. Similarly to Luan et al.,⁵² ionic strengths were calculated using the number of electrolyte ions counted beyond those required to neutralise the silica and DNA. The volume of the liquid system (without DNA and after geometry optimisation) was used for this calculation with the result that stated ionic strengths (0.2 M/1 M) are only estimates. Systems referred to as 0 M contain no added electrolyte in the solvent box, but are electroneutral due to Na⁺ associated with the surface layer and DNA. A summary of each the composition of each model is given in Table 1 and images of the initial configurations of the DNA systems at varying ionic strength are shown in Fig. 2.

Table 1 Total number of atoms of each element in each simulation. For the ‘Na’ row, the first number shown is the number of Na⁺ ions in the solvent box, the number of round brackets is the number of ions initialised as ion pairs with silanolate ions at the surface, and the number in square brackets is the number of ions initialised to neutralise the DNA phosphate groups

Element	0 mM	200 mM	1000 mM	0 mM & DNA	200 mM & DNA	1000 mM & DNA
O	8902	8850	8654	8724	8653	8494
H	13 171	13 067	12 675	12 830	12 688	12 370
Si	1200	1200	1200	1200	1200	8494
Cl	0	18	87	0	6	87
Mg	0	6	29	0	6	29
Na	0 (33)	6 (33)	29 (33)	0 (33)[19]	6 (33)[19]	29 (33)[19]

---

In the simulations, the COMPASS II 1.2 forcefield was used. This forcefield has been parameterised predominately using ab initio data for a wide range of condensed systems; both organic and inorganic systems and on a range of ionic liquids.^53,54 In this forcefield, the charges are the same as in the COMPASS I forcefield (see for reference the previous work by Maekawa et al.^28,30,31), with the exception of the phosphate group on DNA molecules. For this functional group, the COMPASS II defaults were used; −0.3 for the sugar-linking oxygens, −0.822 for the non-linking oxygens (=O and –O⁻) and +0.9246 for phosphorus atoms resulting in a total charge of −1 for each [RCH₂PO₄C(H)R]⁻ substructure. Ewald summation was used for the electrostatic interactions with a 4.184 J mol threshold and an atom-based summation with a 12.5 Å cutoff for the van der Waals interactions. Unless otherwise specified, all analysis was performed over the mean of the last 1 ns with 1 ps windows.

The diffusion coefficient, D, was calculated using the Einstein relation (1), from the Mean Squared Displacement (MSD) as a function of time, t:

The residence times of molecules coordinated with ions was calculated as per the definition of Impey et al.,⁵⁵ based on the rate of decay of the time-correlation function with the parameter t* = 0.

The electric field was calculated in this work from the electrostatic forces on a test charge evaluated using two methods: coulomb summation and Ewald summation. In the Ewald sum, the calculation is periodic in x, y and z, whereas in the coulomb summation method, a finite supercell is used. Details of these methodologies can be found in the ESI 2,† each of which have advantages and disadvantages. Briefly, the coulombic sum method has error from truncation of the full periodicity of the system to a finite sum, and the Ewald sum has error from the contribution of periodic interactions in the z-direction and from the introduction of a non-uniform compensating background charge for orthorhombic unit-cells.⁵⁶

3 Forcefield verification

An accurate description of ion dynamics and structure is fundamental to simulations of charge dynamics in the EDL. The COMPASS II forcefield⁵³ was utilised in this study as it has been parameterised using a wide range of experimental data, including both organic compounds and biochemically relevant ions in the condensed phase. One important requirement of the forcefield is that it describes the solvation of ions accurately. In order to validate the forcefield, initial simulations of the interaction between the Mg²⁺ and Na⁺ cations, and water were performed.

Fig. 3 shows a plot of the coordination number for each ion with respect to pure water (oxygen atoms), as well as snapshots of the equilibrium solvation sphere around the (a) Mg²⁺ and (b) Na⁺ cations. The calculated coordination number of Mg²⁺ was approximately 6, in good agreement with neutron scattering data,⁵⁷ whereas for Na⁺ the coordination number was approximately 5 which is in agreement with neutron scattering of 4.9 ± 1⁵⁸ and ab initio MD of 4.6.⁵⁹ The Mg²⁺ cation showed a more ordered solvation sphere than Na⁺, demonstrated by the initial steep rise and flat region of the curve for the first compared to the second. The secondary peak in the curve shows evidence of a structured secondary solvation sphere. The difference in structure between the ions is due to the smaller ionic radii and stronger coulombic attraction of Mg²⁺ resulting in a much tighter, more ordered solvation sphere.

Fig. 3 Comparison of the water structure around Mg²⁺ and Na⁺ in pure water. Ion–O coordination number (cumulative number of oxygen atoms) for a single Mg²⁺ (blue) and Na⁺ (green) in pure water. There is an initially steep rise for Mg²⁺ followed by an extended flat region, indicating that there is a highly ordered first solvation shell resulting in a radial region where oxygen is absent. Na⁺ by comparison has a more labile first solvation shell, with less distinct regions occupied by the oxygen atoms. The coordination number for Mg²⁺ is 6.2 in agreement with the expected octahedral coordination. For Na⁺ the coordination number is in the range 4.8–5.5, the lower coordination to Na⁺ is expected due to the less efficient packing around the molecule. The curves also show evidence of a second ordered region for Mg²⁺ compared to Na⁺. Shown insets are representative snapshots of the first hydration shell of Mg²⁺ (a) and Na⁺ (b) also demonstrating the difference in ordering of the water molecules. Mg²⁺ shows a tighter binding octahedral structure, as opposed to the more diffuse Na⁺ ion. Snapshots were taken from the system after 2 ns of NVT MD. All molecules shown with atoms within 3.5 Å of the ion, bond lengths (Å) shown in black.

Another test parameter for the forcefield is the diffusion coefficient; the values for Cl⁻, Na⁺, Mg²⁺ and H₂O are presented in Table 2. The calculated diffusion coefficient for water and ions overestimates the experimental value by a factor of approximately 2, which is consistent with other common water forcefields such as TIP3P and is a consequence of the well-known difficulty of accurately capturing water dynamics in empirical forcefields. Agreement within 2- to 3-fold is considered reasonable for diffusion coefficients.⁶⁰ It is promising that the relative diffusivity of the ions is in good agreement with experiment, suggesting qualitatively correct dynamics.

Table 2 Simulated COMPASS II forcefield diffusion coefficients compared to literature data. The Na⁺/Mg²⁺/Cl⁻ diffusion coefficients were calculated using the COMPASS II forcefield at 300 K (with at least 500 ps NPT equilibration) extracting the MSD gradient with respect to time over 8 ps and substituting this into eqn (1). The simulation cell was approximately 46.5 ų, and contained 20 Mg²⁺ and 40 Cl⁻. The COMPASS II forcefield, like the widely used TIP3P model of water, overestimates the diffusion coefficient of water by approximately two fold from experiment, likely resulting in the overestimation of D for the ions. This disagreement reflects the difference in concentration between simulation results and experiment, and the difficulty of accurately parameterising the dynamics of water in a simple empirical forcefield; however, the relative ion diffusivities are in good agreement with experiment

System	Simulated D 1 × 10^−5 cm^2 s^−1	Literature D 1 × 10^−5 cm^2 s^−1
a Extrapolated to infinite dilution. b Simulated using TIP4P and OPLS forcefield at very low ionic strength. c Simulated using TIP3P forcefield.
Na^+ 0.2 M in water	4.1	1.3 (expt.^a)^61
Mg^2+ 1 M MgCl2	1.3	0.71 (expt.^a)^61 0.50 (expt.)^62 0.60–0.79^b ^63
Cl^− 1 M MgCl2	3.6	2.0 (expt.^a)^61 1.4 (expt.)^62 2.4–2.6^b ^63
H2O bulk water	5.9^30	2.3 (expt.)^64 5.2–7.0^c ^65

---

Lastly, the calculated radial distribution function of Mg²⁺/Na⁺ ions to H/O atoms in bulk water is shown in Fig. 4. The coordination number previously shown in Fig. 3 was calculated as the integral of this RDF, assuming a fixed box-volume based on the last frame of analysis. Excellent agreement is seen between the simulated RDF peak position for Mg²⁺–O and Mg²⁺–H versus neutron scattering experiments; however, exact agreement is not expected due to the higher concentration of the experimental data.

Fig. 4 Simulated (red/blue) Mg²⁺–O (left) and Mg²⁺–H (right) Radial Distribution Function (RDF) in bulk water, compared to experimental (green) data for Mg²⁺. For the simulations, 268 water molecules and 1 ion (Mg²⁺/Na⁺ respectively) were prepared in an approximately 20 Å box, NPT molecular dynamics were performed with 2 ns equilibration and a 8 ns production period over which the RDF was calculated. The neutron scattering data was obtained from Bruni et al.⁵⁷ and is measured at 1 : 83 MgCl₂ : H₂O concentration. Excellent agreement is seen between the simulated RDF peak position for Mg²⁺ RDFs versus experiment; however, exact agreement is not expected due to the higher concentration of the experimental data. Na⁺ shows a smaller first peak due to its lower water coordination number and shows a more diffuse RDF for its second peak due to its weaker interaction with water. Na⁺ is also in good agreement with experiments; neutron diffraction data for Na⁺ puts the first RDF peak at 2.35 Å for oxygen and the first peak at 2.91 Å for H.⁵⁸ These results suggest the COMPASS II forcefield is accurately describing the equilibrium ion–water structural properties.

4 Results and discussion

The main aim of the work presented in this paper was to investigate the interfacial EDL structure in a range of ionic strengths, with and without the presence of DNA as an example biomolecule. Simulations were performed using 0 mM, 200 mM & 1000 mM ionic strength 1 : 1 MgCl₂ to NaCl electrolyte, with and without DNA. Videos of each of these molecular dynamics simulations can be found in the ESI.† The simulation results were analysed in sections to examine: ion dynamics in the interfacial region and comparison with accepted double layer models; variation of EDL structure with ionic strength; and the effect of the inclusion of the DNA molecule.

4.1 Ion dynamics at the silica-water interface

Understanding ion dynamics at the silica-water interface is vital not only for improved biosensor design but also in other fields such as geochemistry, where Na⁺ and Mg²⁺ may be important in dissolution reactions,⁶ and chemical engineering for the improved design of water-desalination processes.⁷

4.1.1 Si–O⁻[Na⁺(H₂O)_n] interfacial structure. The simplest interfacial system, representing deionised bulk water or an ionic strength of 0 mM, contains the silica surface with sodium ions neutralising the negatively charged silanols. In order to investigate the structure of the interface at this and higher ionic strengths, as in the previous section, the radial distribution functions (ESI 3†) and coordination numbers (Fig. 5) were calculated for both silanolate–Na⁺ and silanolate–water(H).

Fig. 5 Silanolate-coordination number as a function of 0 mM, 200 mM and 1000 mM ionic strength for silica-water systems. (left) Silanolate–sodium coordination number 0.95, 1.03 and 1.10 for 0 mM, 200 mM and 1000 mM respectively (right) silanolate–water (hydrogen) coordination number 2.83, 2.67 and 2.5 for 0 mM, 200 mM and 1000 mM ionic strength respectively. Increasing ionic strength resulted in an increase in sodium ion coordination to silanolate ions and corresponding decrease in water coordination.

For the 0 mM case, the RDF demonstrated that for the silanolate groups, the mean O⋯H_water hydrogen bond length was 1.25–2.0 Å and each silanolate was coordinated to 2–3 water molecules. The silanolate–Na⁺ coordination number increased with increasing ionic strength, resulting in coordination numbers of 0.95, 1.03 and 1.10 for 0 mM, 200 mM and 1000 mM respectively, as shown in Fig. 5(a). The increased concentration of cations did not affect the ionic bond length or result in multiple sodium ions per silanolate, but did result in a slight reduction in silanolate–water (hydrogen) coordination resulting in coordination numbers of 2.83, 2.67 and 2.5 for 0 mM, 200 mM and 1000 mM ionic strength respectively, as shown in Fig. 5(b). One explanation for this is that Na⁺ accumulation near the interface results in increased displacement of interfacial water.

The residence time of water molecules to silanolate groups was calculated to be 80 ps (based on a silanolate to hydrogen distance cutoff of 3.5 Å) and showed no dependence on ionic strength. This indicates that with increasing interfacial ionic concentration, the water dissociation rate is not significantly affected (ESI 4†).

In conclusion, increasing ionic strength reduced the equilibrium water coordination to silanolate groups, but the water–silanolate kinetics were not significantly affected.

4.1.2 Silica–Na⁺ dissociation. The dissociation of Na⁺ ions, initialised in contact with the silanolate ions, was examined by observing individual ions and calculating their MSD over time. In all three simulations of the silica-electrolyte interface (9 ns simulation time total), only 2 desorption events were observed, one for the 0 mM system after 600 ps and one for the 1 M system after 1700 ps. In both cases, the dissociation mechanism was the same. Examining the desorption event in the 0 mM system in more detail, the calculated diffusion coefficient for the dissociating ion and the roughly linear increase in MSD with time were typical of a unbound stochastic ion (shown in ESI 5†). The mechanism of dissociation is shown in Fig. 6 and was a result of Na⁺ displacement by a fourth water molecule. This resulted in a highly solvated silanolate (coordination number of 4) compared to the average silanolate–water coordination number of 2.5. Sodium desorption kinetics may therefore require a two-step mechanism involving hypercoordination of the silanolate followed by desorption into the bulk.

Fig. 6 Mechanism of dissociation of Na⁺ from the silica surface. Snapshots taken from the 1 M system. At 646 ps water molecules (labeled 1–3) and the Na⁺ (green) were bound to the silanolate ion (large red sphere). Water 4 was loosely associated with a SiO⁻⋯H₂O distance of 2.39 Å. Over the next picosecond, water 4 bound to the silanolate (SiO⁻⋯H distance of 1.49 Å) resulting in dissociation of the Na⁺. The Na⁺ then remained within 5 Å of the silanolate (second solvation sphere) for 40 ps before diffusing into the bulk.

4.1.3 Residence time of the first hydration shell of sodium. The ionic strength of the solution can affect the solvation dynamics of the ions and therefore the structure of the EDL. In order to quantify the characteristic timescale that a water molecule remained coordinated with Na⁺, two systems were considered: the condensed surface layer (≲5 Å from the surface) and the bulk (≳5 Å from the surface). The calculated residence time was approximately 20 ps and 12 ps for the surface and bulk respectively, with a slight dependence on ionic strength (ESI 6†). The demonstrated increase in residence time at the surface suggests a more kinetically stable solvation sphere for surface coordinated ions, and is likely a result of the structuring of water and ions found at the interface. Residence times on the order of picoseconds are consistent with other studies of Na⁺ hydration.⁶⁶

4.1.4 Magnesium ion dynamics. When free magnesium ions approach DNA, it is currently unknown whether the phosphate groups of the DNA displace any of the six Mg²⁺-coordinated water molecules. Several experiments suggest that Mg²⁺ retains its solvation shell,⁶⁷ for example, if Mg²⁺–DNA direct binding was strongly favourable, then Mg²⁺ would be expected to be resolved in X-ray crystal structures of DNA. Furthermore, fluorescence and thermal melting experiments show no indication of Mg²⁺ directly bound to the DNA.⁶⁸

For the duration of all simulated systems presented in this paper, the Mg_(aq)²⁺ ions retained their octahedral water coordination sphere, essentially remaining a single hexahydrated cluster, [Mg(H₂O)₆]²⁺, therefore supporting the notion that Mg²⁺ does not directly coordinate to the DNA. However, NMR experiments have measured a mean water–Mg²⁺ residence time of 1.5 μs⁶⁹ and therefore, microsecond timescale molecular dynamics simulations would be required to sample the full configuration space of the solvent shell.

As discussed, magnesium ions are particularly important in nature and biological systems. Mg²⁺ is known to specifically adsorb to some oxides surface with the extent being highly surface and pH dependent.^70,71 In this work, only on a few occasions did the magnesium ions remain near (<10 Å) the surface. In the 1 M simulation, a Mg²⁺ that was initialised near the surface remained near the surface for the first 1.5 ns and then adsorbed to a specific site, as shown in Fig. 7. This was the only example observed of stable (>500 ps residence time) adsorption of Mg²⁺ to the silica surface, and is the result of the formation of a hydrogen bonded network between two silanolate groups. This result contrasts with the CP/MAS NMR experiments of d'Espinose de la Caillerie et al. which suggested that Mg²⁺ forms direct Si–O–Mg bonds to the surface;⁷² this disagreement may be because their system does not contain sodium ions at the surface. These results suggest that magnesium ions do not readily displace silanolate-bound sodium ions, and that Mg²⁺ interacts with the surface via its solvation shell.

Fig. 7 Adsorption trajectory of an Mg²⁺ cation at the silica-electrolyte (1 M ionic strength) interface. Silica and Mg²⁺ coordinating waters are shown from a snapshot taken at 3 ns. Silanolate groups are shown as black spheres. Hydrogen bonds to the surface are shown with black dotted lines. The trajectory of the Mg²⁺ is shown as a time-coloured line (red-white-blue) from approximately 750 ps (red) to 3 ns (blue). The blue cluster represents the Mg²⁺ reaching a stable adsorption site, the ion reached the site after 1.7 ns and remained there until the end of the simulation (1.3 ns duration). The Mg²⁺ did not bond directly with the negatively charged silanolate groups, instead forming a hydrogen bonded complex via its hexahydrate solvation sphere.

4.2 Poisson–Nernst–Planck (PNP) double layer model

Direct experimental measurement of the distribution of interfacial charge is not available,³⁶ however simulation of the EDL using a continuum model provides a theoretical comparison for expected ion distribution. Due to the high-surface potentials expected at oxide surfaces at the silica-water interface, linearisation of the Poisson–Boltzmann equation is invalid and so the full Poisson–Nernst–Planck equation⁷³ is solved here with a Stern layer, in a simple Gouy–Chapman–Stern (GCS) model⁷⁴ for a mixture of Mg²⁺, Na⁺ and Cl⁻. The details of this model are described in ESI 7.†

The surface potential boundary condition was set to be equivalent to the surface charge density in the MD simulation of 0.2 C m⁻² (see ESI 1 and 7† for further discussion); for example, for the same ionic concentrations as the 1 M mixed-valency MD system, a surface potential of 190 mV was used. This surface potential is in quantitative agreement with experimental measurements at the silica-water interface.⁷⁵ At this surface potential, the ions have reached their bulk concentration within ∼1 nm from the surface, independent of ionic strength (ESI 7†).

The GCS has many limitations in high concentration systems. The dielectric constant of water will not be 80 near the surface due to the strong local interactions. Also, correlated motion between ions is not incorporated, which can be particularly important for systems containing divalent ions.⁴⁵ A further, well-known limitation of the GCS model is that it does not describe finite-size effects and cannot describe the adsorption of ions to specific sites on the surface or ion–water interactions.

In the GCS, high surface potentials can result in extremely high concentrations <1 Å from the surface, which is physically unrealistic¶ due to steric constraints. By increasing the Stern layer thickness, the maximum concentration in the system is reduced, however it is possible to consider extensions to this models to better treat finite-size effects, for example, Kilic and Bazant⁷³ have presented a model which replaces the Stern layer with a layer of c_max cations, or Brown et al.⁷⁶ have presented a model incorporating hydration repulsion interactions which produced a Stern-like layer. These approaches share some findings in common, namely there is expected to be a high cationic concentration within ∼1 nm of the surface under these surface conditions, followed by the smooth decay into the bulk.³⁶

In Fig. 8, the 1 M MD simulation results are compared to the GCS model. As expected, both give the result that the bulk concentration is reached within approximately 1 nm and there is a high cationic concentration within a few angstroms of the surface in a Stern-like layer. The simulations showed that, as discussed previously, strongly favourable solvation of the Mg²⁺ ions resulted in them being distributed roughly evenly through the solvent and not displacing the sodium ions at the surface, with only a small accumulation of ions at the surface (Fig. 8). This result is in good agreement with experimental observations for Mg²⁺ ions around DNA molecules in which the bound state is characteristic by almost complete hydration and free translation and rotational mobility.⁷⁷ These results provide a description that cannot be obtained from the GCS formalism.

Fig. 8 Concentration as a function of z-distance from the surface for the 1 M system (solid lines) compared to the Gouy–Chapman–Stern (GCS) model with a 0.5 Å Stern-layer at 190 mV surface potential (dotted lines). The horizontal lines show the uniform bin size for which the concentration was calculated, and the vertical error bars show the Standard Error of the Mean (SEM). The maximum Na⁺ concentration was 9.5 M and 1.6 M for the MD and GCS results respectively. The atomistic simulation corresponds well with the theory for the bulk concentration, and predicts a high concentration Stern-layer within a several angstroms of the surface. The Mg²⁺ concentration however is much more diffuse than predicted by the GCS because Mg²⁺ did not displace sodium ions from the Stern-layer; such atomistic-chemical effects cannot be described within the GCS formalism.

These results might help to interpret the experimental results of Jayant et al. demonstrating an increase in DNA hybridisation sensitivity upon addition of trace amounts of divalent salt to a monovalent electrolyte system.¹⁵ In their paper, this effect was modelled using a Poisson–Boltzmann model modified to include the effect of variable ion-permittivity due to a biomolecule layer. The results presented here suggest that the double-layer structure for Mg²⁺ containing electrolyte may not be adequately described by the Poisson–Boltzmann model. This discussion will be extended to the effect of DNA later.

4.3 Effect of ionic strength on electrical double layer

The effect of increasing the ionic strength on the interfacial charge distribution is shown in Fig. 9(a). The cumulative charge is plotted as a function of the distance from the unit-cell origin in the z-direction, this figure was obtained as an integral of the average charge distribution (shown in ESI 8†). By plotting the cumulative charge, the zero-charge value corresponds to the position at which the double layer has fully compensated the surface charge. The charge due to silanolate ions is shown at z = 17–18 Å. There is strong similarity between the charge profiles for all ionic strengths and shows positive peaks at 18.5 Å and 21.5 Å, and negative peaks at 20.7 Å and 23.0 Å. These peaks are caused by oriented water molecules; with positive/negative pairs corresponding to layers of oriented water. The similarity between the different ionic strengths indicates that the water is orienting primarily as a result of the surface charge Si–O⁻/Na⁺ layer, as opposed to as a result of the diffuse layer of ions.

Fig. 9 (a) Cumulative charge as a function of distance normal to the silica surface for 0 M, 0.2 M and 1 M systems, this figure is calculated based on cumulative sum of charges in 0.01 Å thick slabs parallel to the xy plane. The negative cumulative charge 16–19 Å was due to silanolate groups at the surface. The subsequent peaks were due to oriented water around the Stern-like layer of electrolyte at the surface. The positive peaks were dominated by hydrogen atoms from water, and the negative peaks by oxygen atoms from water. The inset shows the same graph but with an extended x-axis. (b) The same calculation with water charges are excluded. At ∼3 Å from the surface (z = ∼19 Å), zero cumulative charge was reached due to sodium and magnesium cations neutralising the surface charge. At high salt ionic strength there was a net positive cumulative charge 4–15 Å from the surface (z = 20–30 Å), sometimes called ‘charge inversion’ or ‘concentration polarisation’. The inset shows the same graph but with an extended x-axis. For the 1 M case, the increased positive accumulated charge near the surface induced a negative layer at 30–40 Å due to chloride ions. (c) The mean orientation of water dipoles (cos(θ)) relative to the silica-surface normal, as a function of the z-distance from the surface, using 1 Å bins. A negative value indicates the water hydrogens are pointing towards the silica surface, positive that they are pointing away, and 0 indicates either parallel to the surface or isotropic orientation. Within a few angstroms of the surface (z = 15–19 Å) the water molecules are oriented H-down towards the silanolate groups. Further from the surface, with increasing with ionic strength, the water increasingly orientates H-up towards the chloride ions in the double layer and O-down towards the cations in the Stern-like layer. For high ionic strengths, at ∼18–23 Å from the surface (z = 35–40 Å), the accumulation of chloride anions was sufficient to orientate water H-down, as shown in the inset, which shows the same graph but with an extended x-axis. As expected, the water became isotropic as the bulk was reached.

The distribution of the diffuse-layer charge (≳3 Å from the surface) due to the ions can be seen more clearly by plotting the same function but excluding charges from water atoms, as seen in Fig. 9(b). With increasing ionic strength, an increase in cationic charge at the interface was observed which overcompensates the surface charge up to around ∼1.5 nm from the surface (z = 30 Å). This was due to a significant increase of Na⁺ and Mg²⁺ ions within a few angstroms of the surface, outweighing the cumulative charge from increased chloride ion density in the bulk. This effect is sometimes referred to as charge inversion and is known to occur for multivalent electrolytes near highly charged surfaces and highly charged molecules such as DNA.⁴⁴ In this work, the charge inversion 1.5 nm from the surface was observed to be roughly proportional in magnitude to the ionic strength change (five-fold increase in ionic strength showed a five-fold increase in charge).

The orientation of the water with respect to the normal of the surface is shown in Fig. 9(c). With increasing ionic strength, orientational water polarisation increased. The water orientated H-down towards the negatively charged silica surface near the surface; this oxide-surface induced water polarisation is a well-known phenomenon.^78,79 As expected, the water becomes isotropic as the bulk is reached.

In the high ionic strength simulations (0.2 M & 1 M), the water reoriented at 19–30 Å so as to be H-up towards chloride ions with oxygen towards the interfacial Mg²⁺ and Na⁺. MD simulations of the wet-charged interface in the literature have demonstrated water polarisation in monovalent electrolytes.^28,80 However, the simulations presented here showed that for high ionic strengths, the accumulation of negative charge in the 35 Å to 40 Å region resulted in a secondary layer of H-down orientated water.

In conclusion, with increasing ionic strength, the equilibrium interfacial charge distribution was found to be primarily determined by the water structure around the silanolate and sodium ions at the surface, rather than as a result of the diffuse layer of ions. Charge accumulation and inversion were observed, however, water polarisation lowered the electrostatic energy of the system so as to produce a charge distribution (and therefore potential profile and electric field) which was independent of ionic strength.

4.3.1 Local electric field in the electrical double layer. In the literature, it has been shown in simulations that increased salt concentration led to positive charge accumulation near to the surface due to rearrangement of ions, partially compensated by oriented water at the Stern-like layer directly above the surface.²⁸ By taking the average charge density of atoms and solving the Poisson equation, a decrease in the calculated potential in the interfacial layer was observed. A uniform compensating background charge maintained electroneutrality.

In contrast, in the work presented here, the surface was initialised with a layer of compensating sodium ions; which represents an electroneutral system on the length-scale of the simulation domain. At high ionic strengths this is supported by both theory (Section 4.2) and experiment.⁸¹ For the low ionic strength (“0 mM”') system, this assumption may no longer hold, however the system provides a control against which to contrast the effect of increasing ionic strength. No increase in the interfacial charge density was observed with increasing the ionic strength of the bulk (Fig. 9(a) and ESI 8†); this can be explained by the fact that the surface layer was initialised saturated with a 1 : 1 ratio of counterions. The presence of counterions in the Stern-like layer meant that cation accumulation (Fig. 9(b)) was necessarily weaker than the systems of Maekawa et al., and was completely compensated by water polarisation. This observation has significance for interpreting the response of field-effect transistor-sensors which demonstrate a change in signal with changing ionic strength;³⁰ if many silanolate groups are ion-paired with cations or sterically obstructed from cation binding e.g. proteins, an amorphous surface, then an increase in ionic strength is predicted to correspond to a weaker change in surface-charge accumulation and therefore device response.

A competing factor that could influence device response that is not considered in this work is the effect of variable surface charge. Titration experiments suggest that increasing ionic strength results in an larger apparent negative surface charge, which may be due to altering the chemical equilibrium of the silanol groups.⁸²

4.4 DNA & the double layer

A key aim of this work was to investigate the effect of including DNA molecules on the structure of the interfacial EDL, the electric field and therefore BioFET response. It should be noted that this work does not attempt to provide a detailed analysis of DNA–ion pairing and DNA conformational dynamics, since this topic has received much attention within the literature to date.^43,67,83,84

4.4.1 Effect of DNA on electrolyte structure. The counterion atmosphere around DNA is an area of extensive research, in which it has been proven that there is a closely associated layer of counterions bound to the DNA regardless of their bulk concentration. The ions in this layer are referred to as “condensed ions” in the theory of Manning⁷⁷ (Onsager–Manning–Oosawa condensation). In this theory, for concentrations less than approximately 1 M excess NaCl, 76% of the phosphate groups of β-DNA are calculated to be compensated Na⁺ within ∼10 Å of the DNA surface⁸⁵ which has been confirmed to within 10% by NMR experiments.⁸⁶ For β-DNA in excess MgCl₂ at low concentrations, 44% of phosphate groups are calculated to be compensated by Mg²⁺ (88% charge neutralisation),⁷⁷ this percentage charge neutralisation is supported experimentally by dialysis-monitor titration experiments in which addition of Mg²⁺ to 1 : 1 salt resulted in 85–85% neutralisation,^87,88 and ion condensation has been observed via NMR for divalent ions such as cobalt– and manganese–polyphosphate systems.⁸⁹ In mixed electrolyte, NMR experiments have shown Mg²⁺ can displace DNA-bound Na⁺.⁹⁰

For the MD simulations presented in this paper, there was significant structuring of the water surrounding the DNA when compared to the simulations in absence of DNA (ESI 9†). The RDF of the 1 M DNA–electrolyte system is shown in Fig. 10, and shows that the DNA phosphate backbone attracted a structured counterion cloud in which sodium ions were associated closely with the phosphate groups at approximately 3 Å distance, and formed a secondary more diffuse layer at approximately 6 Å. Consistent with the literature,⁷⁷ the magnesium ions were bound to phosphate groups via hydrogen bonding through the solvation shell; the exception was a single magnesium ion that was initialised in contact with a phosphate group.

Fig. 10 RDF for the phosphate groups of DNA to Mg²⁺ (blue) and Na⁺ (green) respectively, taken from the 1 M silica-DNA system. Na⁺ bound to the phosphate group directly, as seen by the peak at 3 Å. The coordination number is shown as an inset, Mg²⁺ interacted primarily with the DNA through hydrogen bonding of their hexahydrate solvation shell (∼5 Å distance).

The percentage of ions per phosphate group was calculated by inspection of the phosphate–ion coordination number shown in the inset of Fig. 10. Values of 76% and 42% for Na⁺ and Mg²⁺ respectively, were calculated following the methodology of Young et al.,^91,92 in which the second inflection point in the RDF curve was used. These values are in excellent agreement with both the Manning condensation theory predictions and experiment, for non-mixed electrolytes NaCl and MgCl₂ solutions at low concentrations. It can be noted that, given the divalence of magnesium ions, each phosphate group has +1.69e counterion charge within approximately 6 Å, which is evidence of charge inversion around the DNA. This is a phenomenon that is expected for multivalent systems at high concentrations⁴⁶ and is not described by Manning condensation theory, although more recent revisions of the theory have attempted to incorporate these effects.⁹³

Continuing from before, in the work of Jayant et al.¹⁵ the enhanced FET-signal was attributed to a combination of (a) increased ion-exclusion of multivalent ions and (b) increased DNA-condensation in the presence of multivalent ions. The simulations presented here show no evidence of Mg²⁺ ion-exclusion in the DNA layer, suggesting that the previous observed increase in signal in multivalent salt was due to other effects, such as Mg²⁺ induced DNA-condensation onto the surface as supported by other MD studies, which discussed Mg²⁺ induced DNA aggregation.⁴⁴

4.4.2 Effect of DNA on the surface potential and electric field. The experimental response of BioFET devices is still poorly understood, due to a lack of understanding of the interfacial electric field. BioFETs are capable of detecting single molecules suggesting that even fluctuations in the electric field over nanoscale dimensions can be detected;^94,95 despite the atomistic length-scale, there have been few atomistic studies which investigated this behaviour.

The full set of data from the simulations of the charge, potential and electric field in the interfacial region for all six cases, with and without DNA, is shown in ESI 8† for ease of comparison. In the simulations, DNA did not produce a strong effect on the time-average charge distribution of the systems. This might be seen as a counter-intuitive result given that the DNA has a negative charge, however the mechanical flexibility of the DNA–Na⁺ system means that the time-averaged charge of DNA with respect to the surface normal, is expected to be small at any bulk ionic strength.

The electrostatic potential (relative to the silica substrate at 0 V) calculated from this charge distribution, demonstrates that the potential change at the surface (Δ_DNAψ = ψ_DNA(z_surf) − ψ_noDNA(z_surf)) on the addition of DNA was on the scale of millivolts, for example, at the position of the silanolate groups (z_surf = 17 Å), Δψ = −18 mV, −0.5 mV and −37 mV for 0 mM, 200 mM and 1 M systems respectively. This is of the same order of magnitude as surface potential change measurements for biomolecule-oxide systems,^96,97 however these changes were highly sensitive to the choice of surface coordinate and therefore cannot be taken as accurate prediction of surface potential change due to DNA.

Due to natural thermal fluctuations, this mean potential is variable. In order to explore this variability over time due to DNA, the long-range electric field in to the EDL was calculated by measuring the electric field on a test charge 1 Å below the base of the silica. Fig. 11 shows the z-component of the electric field as a function of time for all six cases. No significant difference was found between the electric field for the 0 mM (M = −0.00604, SD = 0.00986) and the 0 mM DNA system (M = −0.00596, SD = 0.0102) based on an independent samples t-test (t(2999) = −0.306, p = 0.759). For the higher ionic strength systems, a small but statistically significant change (200 mM system: p = 5 × 10⁻⁵¹; 1 M system: p = 7 × 10⁻²⁷) in electric field was observed upon addition of DNA ([E(DNA) − E(noDNA)] ≈ 0.002 V Å⁻¹). If this change was a result of the DNA itself and not a result of noise, it would be expected that this response would be strongest in the 0 M system due to lowest ionic screening. Interestingly, the control 0 mM systems showed a 30–40% greater standard deviation in the electric field than the higher ionic strength systems; this indicates that bulk electrolyte may play a role in dampening transient fluctuations in the surface potential.

Fig. 11 500 ps moving average of the electric field (z-component) at a test charge 1 Å below the silica substrate as a function of time, calculated using Ewald summation. For each system, the mean electric field, Ē_z, is shown in the legend and drawn as a colored line, and the standard deviation, σ, is shown in the legend. The addition of DNA did not produce a significant change in Ē_z field for the 0 mM systems but demonstrated a small, statistically significant change for higher ionic strength.

As discussed previously (Section 4.3.1), increasing ionic strength is expected to increase the surface potential.²⁸ The incorporation of an unsaturated Stern-like layer (silanolate groups at the surface without ion paired cations), produced an electric field which was screened as a result of bulk electrolyte and a compensating background charge introduced by the Ewald summation.²⁸ In order to compare these simulations with the work presented in this paper, the mean electric field was calculated (shown in the figure in ESI 10†).

In addition, for the simulations presented in this paper, there was no compensating background charge and the system was neutralised by a Stern-like layer and the bulk electrolyte, resulting in a more compact double layer and therefore a weaker electric field by comparison. No trend in electric field change upon ionic strength increase was observed, in contrast to the strong changes in electric field for the unsaturated Stern layer systems of Maekawa et al. (ΔE_[1M]–[0M] ≈ 0.02 V Å⁻¹).²⁸ This comparison suggests that the electrolyte structure within several angstroms of highly-charged interfaces has a far more significant effect on the electric field, and therefore BioFET response, than biomolecule net-charge/orientation. This also emphasises the importance of developing atomistic models with a realistic description of the Stern layer in order to obtain quantitative atomistic prediction of surface potential.

Finally, for nanowire BioFETs, with small cross sections in the semiconducting region, and a correspondingly high sensitivity, the spatial variation of the field is crucial in understanding the response rather than a smeared out average approximation of the behaviour of an artificial one-dimensional system.^3,98,99 The spatial variation in the electric field at the surface of the silica was investigated for the control system (“0 mM”) and is shown in Fig. 12. This demonstrated that the DNA is having an effect on the local electric field at the surface; however screening from the condensed sodium ions and polarised water reduced the field such that the field perturbation due to DNA was weak compared to thermal noise.

Fig. 12 Weighted electric field (z-component) as a function of space in the xy plane at z = 17 Å for the “0 mM” system without DNA (left) and with DNA (right). The weighted electric field is the electric field divided by its standard deviation, which acts to filter out highly variable thermal noise. The electric field is calculated using the coulomb summation method, taking the mean of the z-component of the field (1 ps frames over the last 100 ps) on a grid of test charges and displayed as a linear interpolated heat map. In the each figure, regular patterns of negative and positive field are due to the regular arrangement of silanolate and ion-paired sodium ions respectively. In the DNA system (right), DNA atoms close to the surface are shown as a transparent overlay, with negative phosphate atoms shown as orange spheres. The positive region at y ≈ 35 Å is due to the positive hydrogen atom on the carbon linker and the negative regions at y ≈ 18 Å and y ≈ 25 Å are due to the negative charges of the DNA phosphate atoms, showing the DNA is having a local effect on the electric field. The magnitude of the electric field due to the DNA is low due to screening from condensed sodium ions and polarised water.

These conclusions suggest that, with a compact neutralising Stern-like layer on the silica surface and around the DNA phosphate groups, that the electric field due to addition of DNA and ionic strength changes is negligible compared to thermal noise. The most likely explanation for BioFET response is that the double layer is more diffuse than determined by these simulations. Over timescales not currently reachable within this MD model, a more diffuse layer might be formed via sodium ions in the Stern-like layer dissociating due to water- or DNA-induced displacement.

5 Conclusions

To the best of the knowledge of the authors, this work presents the first classical molecular dynamics investigation of the bare silica-water interface incorporating magnesium ions, and provides a novel atomistic analysis of the effect of ionic strength and DNA on the electric field and EDL structure at these technologically important interfaces. In this work, molecular dynamics simulations were performed using 0 mM, 200 mM & 1000 mM ionic strength 1 : 1 MgCl₂ to NaCl electrolyte, with and without DNA.

As discussed, understanding ion-dynamics at the silica-water interface is important for a range of systems. The simulations presented here demonstrate that increased ionic strength reduces the equilibrium water coordination to silanolate groups without significantly affecting the water–silanolate kinetics. Sodium ion surface-desorption kinetics required a two-step mechanism involving hypercoordination of the silanolate followed by desorption into the bulk. Sodium ions demonstrated a more kinetically stable solvation sphere at the silica-surface relative to the bulk due to the structuring of water at the interface.

Direct experimental measurement of the distribution of interfacial charge is not available³⁶ and therefore the MD simulation results were compared to a continuum Poisson–Boltzmann model and revealed good agreement with regard to double-layer thickness and sodium ion accumulation. Mg²⁺ did not accumulate significantly at the interface, instead distributing more diffusely than predicted by the Poisson–Boltzmann model. This was as a result of the strong solvation of Mg²⁺ meaning it could not readily displace Na⁺ bound to the surface not described by the Poisson–Boltzmann formalism.

Fundamental understanding of the interfacial charge distribution and electric field is vital to understanding the mechanism of action of field-effect transistor (FET)-sensors. Increasing ionic strength was shown to result in charge inversion due to cation accumulation, an effect which is observed experimentally for divalent ions. The charge inversion 1.5 nm from the surface was observed to be roughly proportional in magnitude to the ionic strength change (five-fold increase in ionic strength showed a five-fold increase in charge). This suggests that charge inversion can begin at lower ionic strengths than those measured by Edel and de Mello.⁴⁸

The results demonstrate that the equilibrium interfacial charge distribution was primarily determined by the water structure around the silanolate and sodium ions at the surface, rather than as a result of the diffuse layer of ions. Furthermore, the electrolyte structure within several angstroms of highly-charged interfaces had a far more significant effect on the electric field, and therefore FET-sensor response, than biomolecule net-charge/orientation. This supports the theory that the mechanism of action of FET-sensor is via modification of the surface chemistry (e.g. altering silanol/silanolate chemical equilibria) rather than the traditional picture based on directly sensing the electric field of the biomolecule.²

As discussed, modelling of the FET-sensor response due to biomolecules is inhibited by a lack of understanding of the interfacial electric field and ion distribution in the presence of biomolecules. The results showed that with DNA present, there was minimal effect due to water polarisation and strong screening by the condensed layer of electrolyte. The first calculation of the time-varying electric field for these systems was presented and, by comparison to a low ionic strength control, showed that bulk electrolyte plays a role in dampening transient fluctuations in the electric field and therefore device response.

These results have relevance to interpreting the experimental results of Jayant et al.¹⁵ which demonstrated an increase in DNA hybridisation sensitivity upon addition of trace amounts of divalent salt to a monovalent-electrolyte system, attributed to ion-exclusion from the DNA region and/or DNA aggregation. The simulations presented here showed no evidence of ion-exclusion from the DNA region and suggests that both (a) the Poisson–Boltzmann model may not be capable of accurately describing the EDL in the presence of mixed electrolyte, and (b) a mechanism other than ion-exclusion, such as DNA aggregation, explains the observed increase in response.

The results also emphasise the role of the Stern-like layer in understanding the response of FET sensors. Changes in the surface charge density (density of silanolate and condensed ions) would be expected to alter the electric field significantly, but in this work, surface charge is effectively fixed due to the long timescale of cation desorption and no chemical reactions. Future work will address this limitation by varying the surface charge on the surface of the model, which can be compared directly with existing experimental titration data.⁷⁵ This will allow quantification of the extent to which the Stern-like layer modulates the electric field and therefore improve fundamental understanding of sensor response.

Acknowledgements

The authors would like to thank Yasushi Shibuta for providing High Performance Computing (HPC) Facilities at the University of Tokyo. The authors acknowledge the use of the IRIDIS HPC Facility, and associated HPC and linux support services at the University of Southampton, in the completion of this work. This work was supported by an EPSRC Doctoral Training Centre grant (EP/G03690X/1) and the JSPS 2015 Summer Fellowship. Data supporting this study are available from the University of Southampton repository at http://dx.doi.org/10.5258/SOTON/401018.

References

1. M.-Y. Shen, B.-R. Li and Y.-K. Li, Biosens. Bioelectron., 2014, 60, 101–111.
2. G. Shalev, Y. Rosenwaks and I. Levy, Biosens. Bioelectron., 2012, 31, 510–515.
3. L. De Vico, L. Iversen, M. H. Sørensen, M. Brandbyge, J. Nygård, K. L. Martinez and J. H. Jensen, Nanoscale, 2011, 3, 3635.
4. C. Barbé, J. Bartlett, L. Kong, K. Finnie, H. Q. Lin, M. Larkin, S. Calleja, A. Bush and G. Calleja, Adv. Mater., 2004, 16, 1959–1966.
5. J. V. Smith, Proc. Natl. Acad. Sci. U. S. A., 1998, 95, 3370–3375.
6. J. D. Kubicki, J. O. Sofo, A. A. Skelton and A. V. Bandura, J. Phys. Chem. C, 2012, 116, 17479–17491.
7. T. Metke, A. S. Westover, R. Carter, L. Oakes, A. Douglas and C. L. Pint, Sci. Rep., 2016, 6, 24680.
8. O. Stern, Z. Elektrochem. Angew. Phys. Chem., 1924, 30, 508–516.
9. H. Helmholtz, Ann. Phys., 1853, 165, 211–233.
10. G. Gouy, J. Phys., 1910, 9, 457–467.
11. D. L. Chapman, Philos. Mag. Series 6, 1913, 25, 475–481.
12. D. Grahame, Chem. Rev., 1947, 41, 441–501.
13. P. M. Dove and C. J. Nix, Geochim. Cosmochim. Acta, 1997, 61, 3329–3340.
14. D. E. Draper, D. Grilley and A. M. Soto, Annu. Rev. Biophys. Biomol. Struct., 2005, 34, 221–243.
15. K. Jayant, K. Auluck, S. Rodriguez, Y. Cao and E. C. Kan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 052817.
16. R. de la Rica, C. Fernández-Sánchez, C. Jiménez-Jorquera and A. Baldi, J. Phys. Chem. B, 2008, 112, 7614–7617.
17. A. A. Shul'ga, M. Koudelka-Hep and N. F. de Rooij, Sens. Actuators, B, 1995, 27, 432–435.
18. L. De Vico, M. H. Sørensen, L. Iversen, D. M. Rogers, B. S. Sørensen, M. Brandbyge, J. Nygård, K. L. Martinez and J. H. Jensen, Nanoscale, 2011, 3, 706.
19. D. Cahen, R. Naaman and Z. Vager, Adv. Funct. Mater., 2005, 15, 1571–1578.
20. C. Heitzinger, N. J. Mauser and C. Ringhofer, SIAM J. Appl. Math., 2010, 70, 1634–1654.
21. C. Heitzinger, Y. Liu, N. J. Mauser, C. Ringhofer and R. W. Dutton, J. Comput. Theor. Nanosci., 2010, 7, 2574–2580.
22. G. Zheng, X. P. A. Gao and C. M. Lieber, Nano Lett., 2010, 10, 3179–3183.
23. S. Rumyantsev, G. Liu, M. S. Shur, R. A. Potyrailo and A. A. Balandin, Nano Lett., 2012, 12, 2294–2298.
24. I. Jokić, M. Frantlović, Z. Djurić, K. Radulović and Z. Jokić, Microelectron. Eng., 2015, 144, 32–36.
25. K. Georgakopoulou, A. Birbas and C. Spathis, J. Appl. Phys., 2015, 117, 104505.
26. M.-P. Lu, E. Vire and L. Montès, Nanotechnology, 2015, 26, 495501.
27. A. Rimola, D. Costa, M. Sodupe, J.-F. Lambert and P. Ugliengo, Chem. Rev., 2013, 113, 4216–4313.
28. Y. Maekawa, Y. Shibuta and T. Sakata, ChemElectroChem, 2014, 1, 1516–1524.
29. Y. Maekawa, Y. Shibuta and T. Sakata, J. Surf. Finish. Soc. Jpn., 2014, 65, 251–256.
30. Y. Maekawa, Y. Shibuta and T. Sakata, Jpn. J. Appl. Phys., 2013, 52, 127001.
31. Y. Maekawa, Y. Shibuta and T. Sakata, Chem. Phys. Lett., 2015, 619, 152–157.
32. C. D. Lorenz, P. S. Crozier, J. A. Anderson and A. Travesset, J. Phys. Chem. C, 2008, 112, 10222–10232.
33. S. V. Patwardhan, F. S. Emami, R. J. Berry, S. E. Jones, R. R. Naik, O. Deschaume, H. Heinz and C. C. Perry, J. Am. Chem. Soc., 2012, 134, 6244–6256.
34. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Chem. Mater., 2014, 26, 2647–2658.
35. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Chem. Mater., 2014, 26, 5725–5734.
36. H. Zhang, A. A. Hassanali, Y. K. Shin, C. Knight and S. J. Singer, J. Chem. Phys., 2011, 134, 024705.
37. B. M. Lowe, C.-K. Skylaris and N. G. Green, J. Colloid Interface Sci., 2015, 451, 231–244.
38. Á. Cimas, F. Tielens, M. Sulpizi, M.-P. Gaigeot and D. Costa, J. Phys.: Condens. Matter, 2014, 26, 244106.
39. K. Leung, I. M. B. Nielsen and L. J. Criscenti, J. Am. Chem. Soc., 2009, 131, 18358–18365.
40. F. Musso, P. Mignon, P. Ugliengo and M. Sodupe, Phys. Chem. Chem. Phys., 2012, 14, 10507.
41. A. Tilocca and A. N. Cormack, ACS Appl. Mater. Interfaces, 2009, 1, 1324–1333.
42. A. M. Soto, V. Misra and D. E. Draper, Biochemistry, 2007, 46, 2973–2983.
43. W. Li, L. Nordenskiöld and Y. Mu, J. Phys. Chem. B, 2011, 115, 14713–14720.
44. B. Luan and A. Aksimentiev, J. Am. Chem. Soc., 2008, 130, 15754–15755.
45. E. Wernersson, R. Kjellander and J. Lyklema, J. Phys. Chem. C, 2010, 114, 1849–1866.
46. A. Y. Grosberg, T. T. Nguyen and B. I. Shklovskii, Rev. Mod. Phys., 2002, 74, 329–345.
47. K. Besteman, M. A. G. Zevenbergen, H. A. Heering and S. G. Lemay, Phys. Rev. Lett., 2004, 93, 170802.
48. J. B. Edel and A. de Mello, Nanofluidics: Nanoscience and Nanotechnology, Royal Society of Chemistry, 2009, p. 7.
49. K. Jayant, K. Auluck, M. Funke, S. Anwar, J. B. Phelps, P. H. Gordon, S. R. Rajwade and E. C. Kan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 012801.
50. K. Jayant, K. Auluck, M. Funke, S. Anwar, J. B. Phelps, P. H. Gordon, S. R. Rajwade and E. C. Kan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 012802.
51. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2004.
52. B. Luan and A. Aksimentiev, Soft Matter, 2010, 6, 243–246.
53. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
54. Accelrys Software Inc., Accelrys Material Studio 8.0, 2007.
55. R. W. Impey, P. A. Madden and I. R. McDonald, J. Phys. Chem., 1983, 87, 5071–5083.
56. E. Wasserman, J. R. Rustad and A. R. Felmy, Surf. Sci., 1999, 424, 19–27.
57. F. Bruni, S. Imberti, R. Mancinelli and M. A. Ricci, J. Chem. Phys., 2012, 136, 064520.
58. N. T. Skipper and G. W. Neilson, J. Phys.: Condens. Matter, 1989, 1, 4141.
59. S. B. Rempe and L. R. Pratt, Fluid Phase Equilib., 2001, 183–184, 121–132.
60. N. Korolev, A. P. Lyubartsev, A. Laaksonen and L. Nordenskiöld, Nucleic Acids Res., 2003, 31, 5971–5981.
61. E. L. Cussler, Diffusion: Mass Transfer in Fluid Systems, Cambridge University Press, 2009, p. 162.
62. S. E. Ingebritsen and W. E. Sanford, Groundwater in Geologic Processes, Cambridge University Press, 1999.
63. A. Chatterjee, M. K. Dixit and B. L. Tembe, J. Phys. Chem. A, 2013, 117, 8703–8709.
64. K. Krynicki, C. D. Green and D. W. Sawyer, Faraday Discuss. Chem. Soc., 1978, 66, 199.
65. D. van der Spoel, P. J. van Maaren and H. J. C. Berendsen, J. Chem. Phys., 1998, 108, 10220–10230.
66. I. Waluyo, C. Huang, D. Nordlund, U. Bergmann, T. M. Weiss, L. G. M. Pettersson and A. Nilsson, J. Chem. Phys., 2011, 134, 064513.
67. C. Maffeo, J. Yoo, J. Comer, D. B. Wells, B. Luan and A. Aksimentiev, J. Phys.: Condens. Matter, 2014, 26, 413101.
68. M. Sotomayor, V. Vásquez, E. Perozo and K. Schulten, Biophys. J., 2007, 92, 886–902.
69. A. Bleuzen, P.-A. Pittet, L. Helm and A. E. Merbach, Magn. Reson. Chem., 1997, 35, 765–773.
70. D. J. Wesolowski, M. L. Machesky, M. K. Ridley, D. A. Palmer, Z. Zhang, P. A. Fenter, M. Predota and P. T. Cummings, ECS Trans., 2008, 11, 167–180.
71. R. Jolsterå, L. Gunneriusson and A. Holmgren, J. Colloid Interface Sci., 2012, 386, 260–267.
72. J.-B. d'Espinose de la Caillerie, M. Kermarec and O. Clause, J. Phys. Chem., 1995, 99, 17273–17281.
73. M. S. Kilic, M. Z. Bazant and A. Ajdari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 021502.
74. O. Stern, Z. Elektrochem. Angew. Phys. Chem., 1924, 30, 508–516.
75. M. A. Brown, Z. Abbas, A. Kleibert, R. G. Green, A. Goel, S. May and T. M. Squires, Phys. Rev. X, 2016, 6, 011007.
76. M. A. Brown, G. V. Bossa and S. May, Langmuir, 2015, 31, 11477–11483.
77. G. S. Manning, Q. Rev. Biophys., 1978, 11, 179.
78. S. Ong, X. Zhao and K. B. Eisenthal, Chem. Phys. Lett., 1992, 191, 327–335.
79. V. Ostroverkhov, G. A. Waychunas and Y. R. Shen, Phys. Rev. Lett., 2005, 94, 046102.
80. J. N. Glosli and M. R. Philpott, Electrochim. Acta, 1996, 41, 2145–2158.
81. T. Baimpos, B. R. Shrestha, S. Raman and M. Valtiner, Langmuir, 2014, 30, 4322–4332.
82. M. Kosmulski, Encyclopedia of Surface and Colloid Science, CRC Press, 2nd edn, 2006, vol. 3, p. 1863.
83. S. Y. Ponomarev, K. M. Thayer and D. L. Beveridge, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 14771–14775.
84. M. Rueda, E. Cubero, C. A. Laughton and M. Orozco, Biophys. J., 2004, 87, 800–811.
85. G. S. Manning, Acc. Chem. Res., 1979, 12, 443–449.
86. C. F. Anderson, M. T. Record Jr. and P. A. Hart, Biophys. Chem., 1978, 7, 301–316.
87. J. Skerjanc and U. P. Strauss, J. Am. Chem. Soc., 1968, 90, 3081–3085.
88. K. R. Bhat, PhD thesis, Rutgers University, 1974.
89. P. Spegt and G. Weill, Biophys. Chem., 1976, 4, 143–149.
90. W. H. Braunlin, C. F. Anderson and M. T. Record, Biopolymers, 1986, 25, 205–214.
91. G. S. Manning and J. Ray, J. Biomol. Struct. Dyn., 1998, 16, 461–476.
92. M. A. Young, B. Jayaram and D. L. Beveridge, J. Am. Chem. Soc., 1997, 119, 59–69.
93. B. I. Shklovskii, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 5802–5811.
94. F. Patolsky, G. Zheng, O. Hayden, M. Lakadamyali, X. Zhuang and C. M. Lieber, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 14017–14022.
95. X. Duan, Y. Li, N. K. Rajan, D. A. Routenberg, Y. Modis and M. A. Reed, Nat. Nanotechnol., 2012, 7, 401–407.
96. E. Souteyrand, J. P. Cloarec, J. R. Martin, C. Wilson, I. Lawrence, S. Mikkelsen and M. F. Lawrence, J. Phys. Chem. B, 1997, 101, 2980–2985.
97. F. N. Ishikawa, M. Curreli, H.-K. Chang, P.-C. Chen, R. Zhang, R. J. Cote, M. E. Thompson and C. Zhou, ACS Nano, 2009, 3, 3969–3976.
98. M. O. Noor and U. J. Krull, Anal. Chim. Acta, 2014, 825, 1–25.
99. D. Passeri, A. Morozzi, K. Kanxheri and A. Scorzoni, Biomed. Eng. Online, 2015, 14, S3.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6cp04101a

‡ These authors contributed equally to this work.

§ However, another more recently hypothesised mechanism of detection is via detecting changes in EDL dipole moment in absence of changes to the surface concentration of ionised groups (e.g. surface charge) or free charges (e.g. electrolyte ions or charged biomolecules like DNA) at the surface. This notion is supported by the experiments of Cahen et al. which showed FET response on addition of neutral organic molecules, or simply oxygen–water vapour.¹⁹ The importance of the dipole moment of the surface has been supported by the simulations of Heitzinger et al. based on the mean-field solution of the Poisson equation.²⁰ This effect could be caused by electronic polarisation of neutral molecules at the silica-water interface.

¶ The maximum concentration possible (c_max = a⁻³ where a is the ionic radius) given steric constraints is 25 M to 207 M for 2–4 Å cations respectively.

---