Protein-Ligand Interaction Detection with a Novel Method of Transient Induced Molecular Electronic Spectroscopy (TIMES): Experimental and Theoretical Studies.

Protein-ligand interaction detection without disturbances (e.g., surface immobilization, fluorescent labeling, and crystallization) presents a key question in protein chemistry and drug discovery. The emergent technology of transient induced molecular electronic spectroscopy (TIMES), which incorporates a unique design of microfluidic platform and integrated sensing electrodes, is designed to operate in a label-free and immobilization-free manner to provide crucial information for protein-ligand interactions in relevant physiological conditions. Through experiments and theoretical simulations, we demonstrate that the TIMES technique actually detects protein-ligand binding through signals generated by surface electric polarization. The accuracy and sensitivity of experiments were demonstrated by precise measurements of dissociation constant of lysozyme and N-acetyl-d-glucosamine (NAG) ligand and its trimer, NAG3. Computational fluid dynamics (CFD) computation is performed to demonstrate that the surface's electric polarization signal originates from the induced image charges during the transition state of surface mass transport, which is governed by the overall effects of protein concentration, hydraulic forces, and surface fouling due to protein adsorption. Hybrid atomistic molecular dynamics (MD) simulations and free energy computation show that ligand binding affects lysozyme structure and stability, producing different adsorption orientation and surface polarization to give the characteristic TIMES signals. Although the current work is focused on protein-ligand interactions, the TIMES method is a general technique that can be applied to study signals from reactions between many kinds of molecules.

Introduction

Protein–ligand interaction is a subject of great interest in the biochemical field due to its scientific significance and practical applications in drug discovery.[1,2] Functional proteins, especially those surrounded by the liquid environment, have extraordinary complexities and degrees of freedom to form 3D structures, while their biological functions are sensitive to and can be modified substantially by their binding with molecules (i.e., ligands) that are much smaller than themselves. The abilities to quantitatively and precisely characterize protein–ligand interactions are essential to understanding and controlling protein’s properties. Existing techniques, such as isothermal calorimetry (ITC),[3] surface plasmon resonance (SPR),[4] biologically modified field effect transistors (BioFET),[5] fluorescence resonance energy transfer (FRET),[6] differential optical scattering,[7] electrophoretic mobility shift assays (EMSA),[8] and small molecule microarray,[9] require the formation of protein crystal or aggregate, fluorescent labeling, or surface immobilization of molecules. Given the small size of ligand molecules and the importance of protein folding in 3D space for the reactions, fluorescent labeling and molecular immobilization can introduce significant disturbances to the reactions, producing potentially incorrect or misleading results in key parameters such as reaction coefficients (e.g., dissociation constant, Kd). On the other hand, these label-free and immobilization-free methods that are currently available, such as ITC and differential optical scattering techniques, render low throughput and limited temporal resolution, and often work only under special conditions (e.g., protein crystallization or exothermic reactions).

Transient induced molecular electronic spectroscopy (TIMES) is a new technique to characterize protein–ligand interactions without the above-mentioned constraints. TIMES is a label-free, immobilization-free technique, and produces accurate and repeatable results with high temporal resolution. In TIMES, the readout is related to molecular interactions with the electrode surface, whereas the reaction itself is performed in the bulk space. As a method of signal readout, the TIMES signal shows the electric response of the reaction products approaching the electrode surface connected to a low-noise electric amplifier. In this paper, we make original contributions in four areas through experiments and physical computations: (a) We demonstrate the accuracy of the TIMES technique by measuring the dissociation constant of lysozyme protein[10,11] with N-acetyl-d-glucosamine (NAG)[11−14] and its trimer, N,N′,N″-triacetylchitotriose (NAG3)[15,16] ligands, and showing that the dissociation coefficient of protein–ligand complex made of the same type of molecule can differ by 3 orders of magnitude. (b) Aided by an analytical model and detailed computational fluid dynamics (CFD) calculations, we show that the measured TIMES signal is directly proportional to the induced charge of a protein molecule (or protein–ligand complex) approaching the electrode. (c) We relate the macroscopic level molecular transfer in a microfluidic channel to the microscopic molecular interfacial mass transfer by incorporating the effect of hydraulic forces and surface’s biofouling (i.e., protein adsorption and desorption from the electrode surface subject to the flow induced shear stress). (d) We perform full-atom molecular dynamics (MD) simulation combined with binding free energy computation to elucidate the fundamentals of the electric signal, which is related to the adsorbed protein’s charge distribution (such as net charge, dipole moment, etc.) and surface polarization at the microscopic level. Through these efforts, we have demonstrated the feasibility and established the physical foundation of the TIMES technique as a method to investigate protein–ligand interactions without labeling or immobilization.

Results and Discussion

TIMES Experimental Measurements of Protein–Ligand Binding Dissociation Constant (Kd)

Figure shows the measured results of the dissociation constants between lysozyme and two ligands: NAG and NAG3. According to the analytical model we developed previously,[17] the measured signal of TIMES system can be represented aswith the electrode area A, volume concentration of molecules (protein, ligand, or protein–ligand complex) no,, rate of molecular adsorption to the electrode surface K+,, the coefficient with ζ being the zeta potential, molecular induced charge as a function of time q(t – to), Q(t – to) = γq(t – to), and diffusion time for the molecule to transport transversely toward the electrode.[17]

Figure 1

Binding between lysozyme and NAG, tested in 1× PBS buffer pH = 7.4 by TIMES: (a) current vs time, (b) calculation of Kd, and (c) histogram of Kd. Binding between lysozyme and NAG3, tested in 1× PBS buffer pH = 7.4 by TIMES: (d) current vs time, (e) calculation of Kd, and (f) histogram of Kd.

Utilizing electric signals for ligand, protein, and protein–ligand complex, protein–ligand dissociation constant can be obtained to estimate protein–ligand dissociation coefficient (Kd),where nL, nP, and nC represent the bulk concentration of ligand, protein, and complex, respectively, in the unit of mol/L. Figure a shows TIMES signals for different ratios of lysozyme and NAG. We have a total of 4 unknowns to be found: the time-dependent induced charge response by the protein, ligand, and protein–ligand complex, namely, QP(t), QL(t), and QP–L(t) in eq , and the dissociation constant, Kd. To find all these values, we perform an experiment by flowing four samples through the device: for example, samples containing protein only, ligand only, and 1:2 and 2:1 protein to ligand concentration ratios before reaction. From the measured TIMES signals of the 4 samples and using eqs and 2 one can find unique solutions for QP(t), QL(t), QP–L(t), and Kd at each time point. Since out of the 4 unknowns only Kd is time independent, we will obtain a histogram for Kd found at each time point. This histogram can produce not only the value of Kd but also the quality of the measurement, since a reliable measurement should yield a tight distribution of the Kd value. In other words, in one single set of measurements, we essentially measure Kd 1000 times over a duration of 1 s at a sampling rate of 1000 s–1. Figure b shows the dissociation constant, Kd, and Figure c shows the histogram of Kd obtained from the method described above. We have found that the most likely value of dissociation constant for lysozyme and NAG is 12.59 mM. The experiment was repeated three times, and the averaged dissociation constant estimated from TIMES is summarized in Table , which also includes previously published value(s)[12,18,19] for comparison. The large dissociation constant (i.e., on the order of 10 mM) suggests that the binding between lysozyme and NAG is very weak. The binding can be strengthened significantly by using the trimer of NAG (NAG3), as shown in Figure d–f. The histogram shows that the value of the dissociation constant between lysozyme and NAG3 is 39.81 μM (Figure f), which is nearly 3 orders of magnitude lower than the value between lysozyme and NAG. The measured value is found to be close to the reported value by other detection methods[16,20] (see Table ). It is noted that, besides the complex of lysozyme–NAG or NAG3, we have also applied the TIMES technique to measure the binding constant of different complexes, e.g., trypsin and p-aminobenzamidine, and thermolysin and phosphoramidon, and obtained very precise values.

Table 1

Comparison of Lysozyme–Ligand Dissociation Constant, Kd, between the TIMES Method and the Previously Published Results Achieved by Other Methods

	Kd
ligand	TIMES	lit. reps
NAG	12.59 mM	16 mM,[12] 24 mM,[18] 7 mM[19]
NAG3	39.81 μM	39 μM,[16] 38.3 μM,[20] 19.6 μM[16]

By observing the TIMES signal waveform (see Figure a and Figure d), one realizes that the signals produced by protein and protein–ligand appear to be obviously different even though the size, molecular weight, and dipole moment of ligand molecule are orders of magnitude smaller than those of protein. A reasonable explanation for why a protein bound with a small ligand can produce significant change in the TIMES signal is that ligand binding can alter the folding and/or orientation of protein. We will present atomistic simulations to elucidate this point later in this paper.

CFD Computations

The TIMES electric signal arises from molecular interactions with the electrode surface and is also affected by external hydraulic forces and surface fouling due to protein adsorption. To further investigate the mechanism, we use both experiments and CFD computations. Figure a presents the measured induced currents under different flow rates. By changing the flow rate of lysozyme in the microfluidic channel, the shear stress is changed and so is the driving force to pull the protein away from the electrode. Therefore, we anticipate that the average dwelling time for a protein molecule on the electrode surface is reduced with increased flow rate, which is consistent with the data in Figure a. Different lysozyme concentrations are also examined as shown in Figure b. By flowing different protein concentrations from 200 μM to 800 μM into the microfluidic channel, the signal intensity increases linearly with the concentration as shown in eq and then saturates as the concentration becomes very high. It is believed that the signal intensity saturation is caused by Coulomb repelling and steric hindrance of molecules near the electrode surface. In other words, the surface adsorption rate K+ in Equation is no longer constant, but decreases with increasing molecular concentration. Finally, it is noted that all the temporal profiles of the TIMES signal (see Figure a,b) exhibit a similar characteristic waveform: starting with a fast increase in the signal intensity and followed by a decrease that often displays a negative overshoot before returning to the zero value. Such a general pattern of the waveforms suggests that the waveform follows the protein flux toward the electrode surface. Each time a protein reaches and leaves the electrode, the TIMES signal is produced.

Figure 2

Experimental current measured by TIMES with different lysozyme flow rates (30−300 μL/min) at constant lysozyme concentration of 200 μM (a); different lysozyme concentrations (200−800 μM) at constant flow rate of 300 μL/min (b); the computation of temporal profiles of scaled flux (F) in CFD simulation as a function of dimensionless groups, Da, ε and Pe (corresponding to (c), (d) and (e) labeling on Figures, respectively). The profile (c) is computed at (ε = 40 and Pe = 5000); the profile (d) is calculated at (ε = 40 and Da = 10); the profile (e) is at Pe = 5000. The dimensionless time (τ) is defined as .

We investigate this inference using fluidic dynamic simulation of the mass transfer process inside a microfluidic channel. CFD computation is performed with a fluid dynamics model, which consists of diffusion, fully developed laminar flow convection, and surface reactions. Protein–surface interactions are generally significantly stronger than protein–protein interactions, e.g., the lysozyme–Au(111) surface binding free energy (∼59 kT) measured from the potential of mean field profile by the umbrella sampling method,[21] which will be reported in our future publication, compared to lysozyme–lysozyme interaction energy (∼0.93 kT) incorporating the hydration and ion effects implicitly through Debye–Hückel theory.[22] Due to the strong surface–protein interactions, it is conceived that the substrate gold surface is covered with a layer of tightly adsorbed proteins and then floppy multilayer adsorption is built up. To simplify the analysis, we adopt a Langmuir adsorption model, in which the effect of surface jamming limit packing is incorporated and only the first-layer adsorption is considered (see eq in Experimental Section). Most of the previous studies[23−25] focused on the steady-state adsorption behavior inside a microfluidic channel involve surface adsorption or reactions; whereas in this work, particular emphasis is placed on the transition state to interpret the result from TIMES experiments. The equations as well as the initial and boundary conditions are scaled in order to reveal the dimensionless parameters governing the system and to explain the general mechanism. A detailed description of the simulation model and its scaled form are provided in the Experimental Section and Supporting Information.

To analyze the effect of the surface reactions, convection, and diffusion on the scaled surface flux (F) in the transition state, dimensionless groups, Damköhler (Da), relative concentration between the bulk and the fully saturated surface (ε), and Péclet (Pe) are introduced,with the initial concentration CA0 of protein solution before entering the microchannel, the bulk protein concentration CA inside the microchannel, protein diffusion coefficient DAB, microchannel height h, adsorption rate Kads, maximum surface adsorption amount CA(max)*, and the average bulk velocity U. The reduced number Da represents the ratio of the surface reaction rate and the diffusion rate, and Pe stands for mass transfer rate ratio of convection and molecular diffusion.

Figure c shows that when there is no protein adsorption (i.e., Kads = 0 or Da = 0), no negative overshooting occurs. As Da increases and the value of Pe(=5000) is fixed, both the intensities of flux peak and negative overshooting increase. The results suggest that the surface adsorption can enhance the mass flux toward the surface. When the surface adsorption amount goes beyond the threshold amount, the extra accumulated amount joins the bulk phase, which leads to the signal’s negative overshooting. Figure e presents the concentration effects on the scaled surface flux (F) temporal profiles at different Da. The literature[23] shows that Da can vary from the order of ∼0.1 to ∼10. If the surface fouling reaction is strong (e.g., Da = 20), at a constant Pe, the surface flux can be scaled by the concentration. However, for the weak surface fouling reactions (e.g., Da = 0.1), the surface flux is not proportional to the concentration increase and the flux reaches its maximum value and declines since the mass transport is limited by surface reactions. Figure d shows the effect of Pe on the scaled flux at a constant ε and Da. Here Pe varies in the range of 102–104, which is the order of magnitude suggested by previous study.[23] One can observe that when Pe increases (i.e., higher flow speed), the period of transition state shortens, whereas flux intensity and overshooting phenomena become magnified. It indicates that strong convection can enhance mass transfer flux to the surface, due to the large concentration gradient normal to the surface. Our results demonstrate that the reduced-unit flux (Figure d,e) from CFD is qualitatively consistent with the experimental results of electric current (Figure a,b). The similarity of their trends illustrates that TIMES electric signal is introduced in the transition state. It is also noteworthy that the CFD simulation offers good explanations for the negative overshooting of the TIMES signal by attributing this phenomenon to protein adsorption kinetics. The overshooting on protein adsorption kinetics can also be visible in the previous experimental measurements[26] of lysozyme adsorption on the C16 hydrophobic self-assembled surface (SAMs) by using total internal reflectance fluorescence (TIRF).

Atomistic Simulations

To interpret the surface’s electric response at the molecular level, hybrid atomistic simulations are employed in this study. MD simulations are first carried out to simulate the solvation structures of lysozyme and lysozyme–NAG complex. Next, with the solvation structures, the adsorbed lysozyme and lysozyme–NAG complex are predicted by using hybrid MM/PBSA[27−31] and full-atom MD simulation according to our previously established protocol.[32] Previous studies[32−34] show that it is computationally expensive to perform full-atom MD simulations to predict protein adsorption in explicit water, and the simulation results are highly dependent on the protein initial orientation, due to the large molecular size and slow rotational motion. Therefore, we perform MM/PBSA to predict the initial orientation of lysozyme on Au(111) surface at a fixed protein–surface distance by treating protein as a rigid body in an implicit water environment. Then full-atom MD simulations are carried out to relax protein conformation on polarizable Au(111) surfaces.

MM/PBSA computations show that, before binding with ligand NAG, the lysozyme molecule most likely “lies down” with residues (Ile78–Asn93) contacting the gold surface, whereas the lysozyme–NAG complex most possibly “stands up” with residues (Asn65–Asn74) close to the Au(111) substrate surface. After 25 ns further relaxation with MD simulations, both protein and protein–ligand complex undergo slight conformational changes. Figures a and 3b compare the final different conformations of pure protein and the complex. Consequently, their dipole moment directions are very different leading to distinct gold surface polarization. For a pure lysozyme, the angle between the dipole moment vector and positive Z-axis normal to the surface is (90.36 ± 4.98)°, while for NAG–lysozyme, the angle is (27.86 ± 8.10)°. To characterize the gold surface polarization, the X–Y plane above the gold substrate above the top of the gold surface (8.075 × 7.992 nm2) is divided into 17 × 17 grids. At each grid point of the position , the electric potential ϕ (kJ mol–1 e–1) of the gold surface’s image charges is calculated by summing up the Coulombic interactions over all surface atoms N,with electric conversion factor f = 138.935 kJ mol–1 nm e–2 and the surface atom i position . The electric potential contour shows the overall effect of protein’s charge distribution as well as all other contributions from solvent environment, i.e., hydration water and ions. As shown in Figure c,d, different contours of induced electric potential resulting from gold substrate surface image charges are detected for the adsorbed pure lysozyme and lysozyme–NAG complex. To further investigate the reason for different lysozyme adsorption orientations before and after ligand-binding, root-mean-square displacement (RMSD) of protein’s heavy backbone excluding hydrogen atoms (see Figure e) is introduced to quantify protein structural stability. By using MD simulation trajectories of both pure lysozyme and lysozyme-NAG complex, RMSD is computed after an optimal overlap,[35] where each instantaneous structure is translated and rotated to superimpose the reference structure. From RMSD profile of pure lysozyme and NAG–lysozyme complex, one can observe that, compared to the lysozyme solvation structure, lysozyme bound with the NAG ligand displays smaller structural fluctuations, particularly around the area (atom #514–587) corresponding to the adsorption region (Asn65–Asn74) for the complex. For further verification, we also examine experimental structure data from the Protein Data Bank. RMSD is computed to compare structures of both pure protein and protein of a complex. The same large variance is also detected at that particular area, i.e., Asn65–Asn74 (see Figure S2). Previous studies[32] showed that a gold surface has a large surface energy. Protein’s slight structural rearrangement can affect protein adsorption orientation due to the strong dehydration free energy of the gold surface.[35] It should be noted, that for a lysozyme molecule, sulfur atoms are not exposed to the solvent environment. Therefore, the possibility of thio (Au–S) interaction between protein and gold surface can be excluded, particularly at the earlier stage when a protein molecule reaches the gold surface. Given the aforementioned analysis, it can be concluded that the binding of NAG results in the variance in lysozyme structural stability and hence affects protein dipole moment direction, which induces different electric responses as the detected TIMES signal.

Figure 3

Snapshots of lysozyme before (a) and after (b) binding with NAG on Au(111) surface from atomistic simulations for 20 ns, their corresponding contours (c, d) of induced gold surface electric potential ϕ on the X–Y plane above the surface of 0.3 nm by taking average with the configurations of the last 2 ns, and the comparison of RMSD of simulated solution structures of pure lysozyme and lysozyme of the complex binding with NAG ligand (e). RMSD profiles for pure lysozyme and lysozyme of the lysozyme–NAG complex solution structures. In RMSD computation, the first conformation of the sampled trajectories was used as a reference. The area around atom #514–587 is indicated with red arrows on the snapshots (b) and RMSD curves (e) respectively. Note: atom number is only for heavy backbone atoms. For the purpose of clarity, water and ions are excluded from the snapshots (a, b). The induced surface electric potential was computed by fixing lysozyme (or lysozyme–NAG complex) configuration to relax the system for 2 ns. The last 1 ns data including total 200 configurations was used to take the average of ϕ.

Conclusion

In summary, the current paper presents the method of transient induced molecular electronic spectroscopy (TIMES), which is capable of detecting protein–ligand interactions in aqueous phase without the need for surface immobilization and fluorescent labeling. The method has been characterized by experiment and physical computations. We have shown experimentally that the TIMES method can accurately measure the dissociation constant for lysozyme–NAG and lysozyme–NAG3 interactions and demonstrated that, compared to monomer NAG, trimer NAG3 can enhance the binding with lysozyme by nearly 1000 times. The theory presented here was further examined experimentally with different concentrations and flow rates of protein in a microfluidic channel as well as theoretically by performing CFD and atomistic simulations. The CFD simulations suggest that attractive protein–electrode force and repulsive shear force in a microfluidic channel determine the surface flux of protein, which gives the general waveform of the TIMES signal characterized by a positive peak followed by a negative overshoot before returning to the baseline. With the induced electrical signal from proteins dwelling on the electrode surface, we systematically investigated the effects of protein and ligand concentration and hydraulic shear stress on protein–ligand binding and proteins’ kinetic transport at the water–surface interface inside a microfluidic channel. The efficient hybrid MM/PBSA calculations and MD simulations predict the most probable adsorption orientations of protein and protein–ligand complex and the subsequent surface polarization. The results support our theory that protein configuration change due to ligand binding contributes to the TIMES signal and enables the method to detect protein–ligand interactions and find the reaction dissociation constant. Although we have focused our studies to protein–ligand interactions here, the technique and general principle of TIMES can be easily extended to study interactions of different kinds of molecules as an effective tool to characterize biomolecular reactions in conditions closest to their native environments. In this paper we demonstrate the TIMES method for characterization of protein ligand interactions with a single binding site. More complicated systems involving multiple-ligand binding will be investigated in the future.

Experimental Section

TIMES Experiment Setup

The setup of the TIMES system is shown in Figure . It consists of a microfluidic channel to allow the biomaterials to flow through, a pair of gold electrodes on the floor of the microfluidic channel as the sensing electrode and ground electrode, two inlets to inject biomolecules of interest and buffer respectively, and a transimpedance amplifier with its input connected to the gold sensing electrode. The microfluidic device was fabricated on a 1 mm thick glass slide (VWR). Before fabrication, the glass slide was cleaned in acetone, methanol, and isopropyl alcohol (IPA) sequentially with sonication. A thin layer of NR-1500PY (Futurrex, USA) photoresist was spun on the glass slide. After lithographic patterning of the photoresist, 100 nm thick titanium and 200 nm gold layers were deposited by sputtering (Denton Discovery 18, Denton Vacuum, LLC), followed by the lift-off process in acetone. The fabricated Ti/Au pattern has a sensing area of 1 × 1 mm2, and the electrode pattern is extended outside of the channel for wire connection by soldering. The microfluidic channel was fabricated using soft lithography process. The mold was fabricated on a 4 inches Si wafer. After the standard wafer cleaning process, a layer of 30 μm thick SU8-2050 (Microchem) photoresist was spun on the silicon wafer and patterned photolithographically to form the SU8 mold. Uncured polydimethylsiloxane (PDMS) was poured onto the SU8 mold and cured in a 65 °C oven. After curing, holes were punched on the PDMS to form inlets and outlets. Then the PDMS and the glass slide with patterned electrodes were both treated with UV ozone before they were aligned and bonded together. The sensing electrode was connected to the input of a low-noise transimpedance amplifier (SR570, Stanford Research System, Inc.), and the ground electrode was connected to the instrument ground. The transimpedance of the amplifier was set at 100 MΩ, and the voltage output of the amplifier was digitized by a DAQ board (National Instrument) and recorded by Labview at 1 kHz sampling rate. The primary sources of noise are interference from the environment and thermal noise of transimpedance amplifier. Using an electromagnetic shielded chamber, an improved amplifier with lower thermal noise, and digital filters (Labview), the signal-to-noise ratio can be improved to allow us to measure reactions with very low (e.g., pM) dissociation constant.

Figure 4

Schematic of transient induced molecular electronic spectroscopy. The setup consists of a microfluidic channel with electrodes and fluid inlets and outlet. Biomolecule solution and pure buffer are introduced into the microchannel through two separate inlets, and the transient current signal from the gold electrode is amplified and recorded by the external circuit.

Biomolecule Test on TIMES

Before the test, the microfluidic channel is first filled with buffer from inlet. Chosen amounts of protein and ligand are dissolved in buffer solution before the test. For samples that contain both protein and ligand molecules, the samples are set aside for 3 h before the test to ensure that the reaction has reached the equilibrium state. All the measurements are conducted at room temperature, and each test is repeated three times to confirm repeatability. After each test, the device is washed with buffer to remove any biomolecule residues in the microfluidic channel or on the electrode surface. The protein and ligand binding experiments are performed in 1× PBS buffer at pH = 7.4. The data obtained after the amplifier and ADC are low pass filtered digitally in Matlab to remove noise, and the dissociation constant Kd is calculated and plotted with Matlab. The detailed analysis is shown in the Supporting Information. The protein experiment under different flow rates and concentrations is performed in 50 mM Tris-HCl buffer at pH = 7.4.

CFD Computation

The transport of lysozyme solution inside a microchannel is modeled with the mass transfer equation (eq ), which includes convection and diffusion, coupled with surface adsorption (eq ),with lysozyme bulk concentration CA, bulk solution concentration near the surface CAS, surface adsorption concentration CA*, and maximum surface adsorption amount CA(max)*. Due to large protein–surface binding free energy, Kdes is ignored in our computation to simplify the analysis. A fully developed laminar velocity profile is adopted. The microfluidic channel is initially filled with pure water before being flushed with protein solution of concentration CA0. The governing equations, initial and boundary conditions in dimensional and scaled forms are shown in the Supporting Information. COMSOL multiphysics software (Version 5.1, COMSOL Inc. USA) is used to solve the partial differential equations.

MD simulations are performed with Gromacs software package (version 4.6.5)[36] in NVT ensemble by using Charmm36 force field[37−39] for protein and NAG molecules, and tip3p water model. Lysozyme crystal structure (pdb code: 3TXJ) is obtained from the Protein Data Bank and is with a net charge of +8e at pH 7 (see the Supporting Information). The system is neutralized by adding Cl– ions. In addition, 40 pairs of Na+ and Cl– are added into the system to keep the ion concentration equal to 120 mM. To obtain solvation structures, pure lysozyme and lysozyme–NAG complex are first equilibrated in water environment for 50 ns. Detailed discussion about MD simulation is shown in the Supporting Information.

A two-step procedure is adopted to predict protein adsorption and consequent surface polarization. First, hybrid molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) computations[27−31] are performed to predict the protein initial orientation on Au (111) surface based on protein binding free energy, which consists of protein–surface interactions and hydration or dehydration free energy, according to our previously established protocol[32] to serve as an initial value for the following MD simulations. In MM/PBSA, the solvated protein and protein–NAG complex are treated as rigid bodies respectively, and are rotated around their center of mass on Au (111) surface (8.075 × 7.992 nm2) while fixing protein–surface minimum distances (i.e., 0.3 nm) to search for the most energetically favorable orientations. To simplify the computation in MM/PBSA, the nonpolarizable Au surface parameters[40,41] are used for the Au–protein and Au–water interactions. We also change the protein–surface distance to 0.26 nm, which is closer to the surface, and the same most top-ranking orientations for both pure protein and complex were identified from MM/PBSA. The surface tension of Au (111) (γ = 1.41 J/m2) is adopted from the literature report[41] in MM/PBSA computations. Second, a full relaxation of the initial adsorbed lysozyme and lysozyme–NAG configurations is performed with full-atom MD simulations for 20 ns with polarizable force field parameters[42] of Au (111) surfaces, which was developed by Walsh et al. and accounts for the interactions between peptides or protein and the induced surface image charges by introducing dummy atoms to form rigid-rod dipoles free to rotate around atomic sites. Surface atoms are aligned with periodic boundary image atoms in accordance with the gold crystal lattice to mimic a large surface without boundary effects. Two repulsive walls are built on the top and bottom layers of the z-direction to confine solvent molecules. At the top of the gold surface, a water box of 7.6 nm height is built.