Constant-pH Brownian Dynamics Simulations of a Protein near a Charged Surface

We have developed a rigid-body Brownian dynamics algorithm that allows for simulations of a globular protein suspended in an ionic solution confined by a charged planar boundary, with an explicit treatment of pH-dependent protein protonation equilibria and their couplings to the electrostatic potential of the plane. Electrostatic interactions are described within a framework of the continuum Poisson-Boltzmann model, whereas protein-plane hydrodynamic interactions are evaluated based on analytical expressions for the position- and orientation-dependent near-wall friction tensor of a spheroid. The algorithm was applied to simulate near-surface diffusion of lysozyme in solutions having pH in the range 4–10 and ionic strengths of 10 and 150 mM. As a reference, we performed Brownian dynamics simulations in which the protein is assigned a fixed, most probable protonation state, appropriate for given solution conditions and unaffected by the presence of the charged plane, and Brownian dynamics simulations in which the protein probes possible protonation states with the pH-dependent probability, but these variations are not coupled to the electric field generated by the boundary. We show that electrostatic interactions with the negatively charged plane substantially modify probabilities of different protonation states of lysozyme and shift protonation equilibria of both acidic and basic amino acid side chains toward higher pH values. Consequently, equilibrium energy distributions, equilibrium position-orientation distributions, and functions that characterize rotational dynamics, which for a protein with multiple ionization sites, such as lysozyme, in the presence of a charged obstacle are pH-dependent, are significantly affected by the approach taken to incorporate the solution pH into simulations.


INTRODUCTION
Interactions of proteins with charged surfaces and diffusion of proteins in electric fields are important in many applications, such as the chromatographic separation of proteins, 1 protein immobilization, 2 biosensing, 3 design of biocompatible materials for medical applications, 4 and removal of protein fouling from metal surfaces in the food and pharmaceutical industries. 5 Reversible binding of various peripheral proteins to the negatively charged cytoplasmic leaflet of the cell membrane, oftentimes mediated by electrostatic interactions, is an important step in a diverse array of cellular processes, from cell signaling to membrane trafficking. 6 All proteins contain basic and acidic titratable groups capable of exchanging protons with their environment. Solvated protein switches constantly between available protonation states. Consequently, the protein charge distribution and its net charge constantly fluctuate. Magnitudes of these fluctuations, and their equilibria, depend on the solution pH. If a proton uptake/release by a given titratable group in the protein, A, is assumed to adhere to a scheme: , usually used in the form of the negative decimal logarithm, pK a = − log 10 K a . The pK a 's of ionizable groups in proteins and protein complexes are often much different than these measured for isolated groups in solution. 7,8 Shifts in pK a 's between native and denatured states of a protein or between bound and free forms of a protein complex give rise to pH effects on protein activity, stability, binding kinetics, and binding constants. 9−12 One should expect the presence of a charged surface, which is a source of an electrostatic potential, to affect protein protonation equilibria. pK a 's of protein ionizable groups are anticipated to be shifted relative to values observed in the absence of an external electric field due to protein-surface electrostatic interactions.
Variations in the protein charge distribution, coupled to the external electric field through the described above shifting of pK a values, inevitably affect diffusional dynamics of the protein suspended in an ionic solution near a charged boundary. This can be seen, although indirectly, in experiments. For example, single-molecule total internal reflection fluorescence (TIRF) microscopy and dual polarization interferometry experiments showed pH-dependent adsorption kinetics of, respectively, the BSA protein on silica 13 and the Mefp-1 protein on silicon oxynitride. 14 Single-molecule TIRF microscopy studies on the dynamics of α-lactalbumin at the silica surface 15 showed that affinity, immobilized adsorption, and surface diffusion of the protein can be controlled by pH. Solution pH should be thus included in computational modeling and simulations of protein-surface interactions, near surface diffusion, and adsorption kinetics. Dynamical studies may require consideration of changes in the protonation states of ionizable groups as a function of time and position and orientation of the protein in the electric field of the charged surface.
A number of computer simulation studies on different aspects of protein-surface systems, using either atomistic or mesoscopic interaction potentials, have been published. 16−27 Starzyk and Cieplak 16 investigated conformational changes of tryptophan cage induced by a uniform electric field and the electric field near the mica surface, employing all-atom molecular dynamics simulations. Xie and co-authors 17 applied the parallel tempering Monte Carlo approach to simulate preferred orientations of a coarse-grained lysozyme model adsorbed on positively and negatively charged surfaces, taking into account electrostatic and van der Waals interactions between the protein and the surface. Brancolini and coauthors 18 studied the molecular basis of interactions between ubiquitin and gold nanoparticles combining quantum mechanics, molecular dynamics, and Brownian dynamics simulations. Kubiak and co-authors 19,20 performed atomistically detailed molecular dynamics simulations of lysozyme adsorbed on mica and silica surfaces, focusing on identification of residues involved in surface-protein contacts, preferred orientations of the protein, and possible alterations of its secondary and tertiary structure. A somewhat similar study was conducted by Broemstrup and Reuter 21 for the proteinase 3 protein anchored at three different phospholipid bilayers. Yu and co-authors 22 performed coarse-grained molecular dynamics simulations of lysozyme adsorption on different surfaces (hydrophobic, neutral hydrophilic, zwitterionic, negatively charged, and positively charged) and described conformational changes of the protein upon adsorption, its orientations on surfaces, interactions that promote adsorption on different surfaces and key residues playing the role in the anchoring of the protein. The role of lysozyme flexibility upon adsorption to the source 15S cation exchange surface was evaluated by Steudle and Pleiss 23 who used all-atom molecular dynamics. Lumb and Sansom 24 combined rigid-body Brownian dynamics and molecular dynamics to decipher the mechanism of recognition (in particular the role of electrostatic interactions) by pleckstrin homology domains of specific lipids within the cytoplasmic leaflet of the plasma membrane. Ravichandran and co-authors 25 performed Brownian dynamics simulations of initial stages of adsorption of a rigid, atomistically detailed model of lysozyme on a positively charged plane, and described the role of the protein charge distribution on controlling the adsorption behavior. In another work, Ravichandran and Talbot 26 used Brownian dynamics to simulate adsorption of lysozyme molecules, modeled as uniformly charged spheres, on a charged planar surface, and evaluated lysozyme lateral diffusion as a function of the surface coverage. Gorba and co-authors 27 performed Brownian dynamics simulations of protein-surface systems consisting of up to 150 horse heart cytochrome c molecules, modeled as dipolar spheres. Their main objective was to examine the overall behavior of cytochrome molecules in the presence of the charged surface. All of these works neglect variations of protonation states of titratable groups and their coupling to the electrostatic field of the charged surface. Instead, the protonation states of titratable groups of biomolecules are kept fixed during simulations and the solution pH is accounted for implicitly by assigning to acidic and basic residues predefined protonation states, determined for instance using standard empirical algorithms 28 or electrostatic free energy calculations. 29 Different constant-pH molecular dynamics algorithms that treat the pH-dependent protonation equilibria of biomolecules in an explicit fashion are available. 30−37 Even though in principle, some of these methods can be applied to study dynamics of proteins near biological and nonbiological surfaces, either directly or after adding an implementation of the coupling of protein protonation equilibria to the molecular potential of the surface, constant-pH molecular dynamics simulations of diffusion and adsorption kinetics are problematic due to computational complexity, sampling, and convergence issues, 38 significant length of trajectories required to calculate position and orientation correlation functions of a sufficient quality so that they can be used to reasonably estimate diffusion coefficients, 39,40 and finite size effects on diffusion due to periodic boundary conditions. 41−44 As far as we are aware, constant-pH molecular dynamics simulations have not been applied to study diffusional dynamics in proteinsurface systems.
Diffusional dynamics of biomolecules of arbitrary shapes, treated as rigid bodies, can be modeled using Brownian dynamics methods. 45,46 These methods benefit from the continuum representation of the solvent and freezing of internal degrees of freedom of diffusing molecules, which reduce the computational cost allowing for rapid generation of dynamical trajectories over longer time scales than are accessible by molecular dynamics. Rigid-body BD technique has been used in studies on protein adsorption on different surfaces, as already mentioned above. 18,24−27 Algorithms applied in these studies lack the explicit treatment of pHdependent protonation equilibria. Moreover, as they neglect protein-surface hydrodynamic interactions, 47−49 they cannot provide a quantitative description of near-surface diffusion and adsorption kinetics.
The issues that we tackle in this paper are thus the following: how does the solution pH affect electrostatic interactions of proteins with charged surfaces and near-surface diffusional dynamics of proteins, and what are the possible consequences of neglecting pH-dependent protonation equilibria in computational algorithms that are being applied to study protein-surface systems. To answer these questions, we have developed a rigid-body constant-pH Brownian dynamics (cpH BD) algorithm, which allows us to investigate effects of solution pH on diffusive dynamics of globular proteins near charged planar surfaces. In our method, BD sampling of the protein configuration space, which consists of positions and orientations evaluated relative to the charged plane, is augmented by the sampling of discrete protein protonation states through Monte Carlo (MC) moves. MC sampling of protonation states is based on the protein electrostatic free energy in a given protonation state and in a given configuration relative to the surface at a given solution pH and ionic strength.
As a prototypical protein for our studies, we choose hen egg white lysozyme (Figure 1). Lysozyme is an attractive model in a wide spectrum of physical, chemical, and biological problems, including protein adsorption at charged biological 52−55 and nonbiological surfaces. 19,20,23,25,56,57 Lysozyme contains 32 ionizable groups ( Figure 1). Its isoelectric point was reported to be 11.35. 58 It is a very stable and robust protein, whose secondary and tertiary structures have been reported not to undergo significant changes in acidic solutions down to pH 1.0. 59,60 This protein maintains its activity across a broad range of pH conditions, from highly acidic to highly basic. 61 Lysozyme resists significant conformational changes upon adsorbing on solid surfaces. 62,63 In total internal reflection fluorescence studies of lysozyme adsorption on negatively charged silica surfaces, 62 lysozyme behaved as expected for a hard protein, undergoing primarily electrostatic adsorption with minimal conformational changes. There is also no evidence that the conformation of lysozyme changes significantly during initial stages of an adsorption on biological membranes. 54 Thus, our approach to model lysozyme as a rigid conglomerate of atoms seems reasonable. Moreover, the hydrodynamic shape of lysozyme is well described by a prolate ellipsoid, 64 which enables us to evaluate protein-surface hydrodynamic interactions in BD simulations based on analytic expressions for friction and mobility tensors of an axisymmetric body diffusing in the presence of an infinite, planar wall. 49 We simulated Brownian motions of lysozyme suspended in an ionic solution near a homogeneously (negatively) charged infinite plane (Figure 1) at the solution pH in the range 4−10 and at ionic strengths of 10 and 150 mM. We performed three types of Brownian dynamics simulations, which differ in the way the solution pH is accounted for. The first type is cpH BD simulations, in which the protein probes possible protonation states with the pH-dependent probability, and this probability is affected by electrostatic interactions with the charged plane. We will refer to these types of simulations using the term "cpH BD simulations with the biased protonation sampling". The second type are cpH BD simulations, in which the protein probes possible protonation states with the pH-dependent probability, but these variations are not coupled to the electric field of the charged surface; from the protonation standpoint, the protein behaves as it would in an unbounded solution in the absence of any external fields. These are "cpH BD simulations with the unbiased protonation sampling". Finally, the third type are conventional BD simulations, in which the protein is assigned a fixed, most probable protonation state, appropriate for a given solution pH, and unaffected by the presence of the charged plane.
What follows is the analysis and discussion of results of these three approaches to model pH effects in BD simulations of protein-surface systems, preceded by the presentation of the constant-pH BD algorithm.

THEORY AND METHODS
2.1. Rigid-Body Constant-pH Brownian Dynamics. In constant-pH Brownian dynamics simulations, the protein, treated as a rigid body, undergoes translational and rotational diffusion and concurrently probes different protonation states, with the probability that depends on the solution pH and, as protonation equilibria are coupled to the external electrostatic potential, on the position and orientation of the protein relative to the charged plane. Each step of the iterative cpH BD algorithm involves two stages: BD sampling of protein configuration space and MC sampling of discrete protonation states. The central part of the BD sampling is the mobility tensor whose matrix, due to hydrodynamic interactions, depends on the position and orientation of the protein relative to the plane. 65 Brownian dynamics trajectory is propagated with the protein assigned a set of discrete charges, corresponding to the current protonation state. Protein diffuses in an ionic solution confined to a half-space above the charged planar surface, influenced by electrostatic forces and torques resulting from the electrostatic potential of the boundary. This potential is obtained in an analytical form as a solution of the Poisson−Boltzmann equation. 66 Periodically, the BD trajectory is interrupted and a change in the protonation state of the protein is attempted using a Monte Carlo algorithm. 67−69 If the change is accepted, the simulation is restarted with a different set of discrete charges assigned to the protein corresponding to the new protonation state. Otherwise, the simulation continues without modifications in the protein charge distribution. MC sampling of protonation states is based on the protein electrostatic free energy in a given protonation state and at a given solution pH. The free energy of a protonation state is evaluated numerically as described previously in the case of a protein in an unbounded solution 70 using the Poisson−Boltzmann (PB) model, 71 with a correction that we propose for the position-and orientation-dependent protein-surface electrostatic interactions. These two stages are described in detail in the Supporting Information.
2.2. Model of the Lysozyme Molecule. Structure of hen egg white lysozyme (PDB ID 6LYZ, 50 Figure 1) was used as a basis for the model of a protein with multiple ionization sites.  Table S1 for details). The hull of the lysozyme hydrodynamic prolate ellipsoid of revolution is depicted as a black mesh. Center: the excluded volume model of lysozyme. Each heavy atom is assigned a sphere of radii 2 Å. Right: schematics of the protein-plane system considered in this work. Lysozyme is suspended in an ionic solvent confined by a charged, impenetrable boundary (shown in gray, σ is the constant surface charge density). Solvent is treated as a continuum and assigned the viscosity η, the dielectric constant ϵ s , and the Debye length 1/κ. Protein is modeled as a rigid conglomerate of spheres and assigned a distribution of point charges. Hydrodynamic shape of the lysozyme is approximated by a prolate spheroid of semi-axes (a, b). An outer reflective planar boundary limits values of the z coordinate of the protein center. Left and center drawings were done using the UCSF Chimera package. 51 We started by performing a molecular dynamics simulation of lysozyme in the CPT ensemble, with the GROMACS 4.6.5 package, 72 and allowing the crystallographic structure to relax. Details of the simulation protocol are given in the Supporting Information.
Three types of interactions evaluated in Brownian dynamic simulations, namely, hydrodynamic, excluded volume, and electrostatic, require three representations of the lysozyme molecule. Prior to creating these three representations, coordinates of lysozyme atoms, resulting from CPT simulations, were transformed to a molecule fixed frame whose origin coincides with the diffusion center 73 of the molecule, and axes point along principal axes of the lysozyme rotational diffusion tensor. This transformation, which requires calculation of the diffusion tensor of the molecule from its atomic structure, was done using in-house software as described elsewhere. 69 From the hydrodynamic standpoint, the lysozyme molecule is described as a prolate ellipsoid of revolution (Figure 1), in accordance with previous depolarized light scattering measurements. 64 Parameters of this ellipsoid were derived directly from the lysozyme three-dimensional atomic structure using the HullRad package. 74 We obtained values of 20.5 and 13.8 Å for the long and the short semi-axis of the lysozyme hydrodynamic ellipsoid, respectively. The center of the ellipsoid coincides with the diffusion center of lysozyme and its semi-axes are collinear with corresponding principal axes of the lysozyme rotational diffusion tensor. Protein-plane excluded volume interactions (see the Supporting Information) are evaluated with the protein treated as a rigid conglomerate of spheres corresponding to lysozyme heavy atoms ( Figure 1). All spheres are assigned the same radius of 2 Å. Electrostatic interactions are evaluated during BD simulations treating the protein as a set of point charges in an external electrostatic potential (equation S20). These charges correspond to formal charges on lysozyme titratable groups in a given protonation state. Each charge is assigned to a particular atom within a given group (Figure 1). Positions and values of charges are given in Table S1. 2.3. Simulation Setup. We have conducted cpH BD simulations with either the biased (equation S4) or unbiased protonation sampling (equation S2) as well as BD simulations in which the protein is assigned the most probable protonation state (in the absence of the charged plane) at a given pH and the ionic strength; this state is kept fixed during simulations.
All Brownian dynamics simulations were performed using the software developed in-house. During simulations, the motion of the protein in the direction normal to the charged plane is restricted to the interval (0,7a) (Figure 1), where a is the long semi-axis of the lysozyme hydrodynamic ellipsoid. The charged plane, which is impenetrable for the protein (see The same most probable protonation states were obtained at ionic strengths of 10 and 150 mM.

ACS Omega
http://pubs.acs.org/journal/acsodf Article the Supporting Information), is located at z = 0. Additionally, Brownian dynamics steps that result in values of the z coordinate of the protein's center that are greater than 7a are rejected, which effectively creates an outer reflective boundary. 75 No restrictions were imposed on motions in directions parallel, (x, y), to the charged plane. Simulation time step (Δt in equation S5) was set to 20 ps. Probabilities of protein protonation states were evaluated at each cpH BD step, using the MC procedure described in the Supporting Information, with the number of MC steps set to 10 6 . Current protonation state of the protein was saved to a file at each cpH BD step. The solvent dielectric constant was set to 78.54, the solvent viscosity was set to 0.089 cP, the temperature was set to 298.15 K. Surface charge density on the planar boundary was set to −0.015 e/Å 2 . For a comparison, the total number of negative lattice sites on the mica basal plate corresponds to one negative charge per 48 Å 2 resulting in the surface charge density of −0.02 e/Å 2 , 76 whereas area per lipid in the DOPG bilayer at 298.15 K is roughly 71 Å 2 , 77 resulting in the surface charge density of −0.014 e/Å 2 . For each type of simulations and for each value of the solution pH and the ionic strength, we have generated three separate 40.5 μs trajectories. At the beginning of a trajectory, the z coordinate of the protein center was chosen randomly, with an equal probability, from the interval (0,7a), and the orientation of the protein in the laboratory frame was generated using uniformly distributed random rotation matrices. Starting configuration was accepted provided that there was no overlap between the protein and the charged boundary. Position and orientation of the protein and its electrostatic energy in the external potential of the charged plane were saved to a trajectory file at each step of the BD (or cpH BD) simulation. The first 0.5 μs of each trajectory was omitted from the further analysis. The three trajectories generated for a given set of conditions were analyzed separately, and results of this analysis (i.e., curves and distributions) were then averaged; standard deviations are shown in two-dimensional plots as error bars.

RESULTS AND DISCUSSION
3.1. Lysozyme Titration Behavior. Table 1 contains charges of lysozyme ionizable groups in the most probable, at a given pH and an ionic strength, protonation state. These most probable states were found with the Monte Carlo sampling of the free energy levels of protein protonation states. Free energies were computed numerically using the PB approach in the absence of the charged plane, as described in the Supporting Information.
For the whole considered pH range, changing the value of the ionic strength from 10 to 150 mM does not affect the lowest energy protonation state. There are four ionizable We only show charges of groups for which the mean square fluctuation of the charge exceeds 0.01 e 2 . Average values of remaining charges are equal to values given in Table 1. The first number in each cell is a result of cpH BD simulations with the unbiased protonation sampling whereas the second number was obtained from cpH BD simulations with the biased protonation sampling.

ACS Omega
http://pubs.acs.org/journal/acsodf Article groups that change their charge in the lowest energy state upon varying the solution pH between 4 and 10. These are as follows: the N-terminal amino group (NTER), which loses proton between pH 8 and pH 9; the side chain of tyrosine 23 (TYR 23), which loses proton between pH 9 and pH 10; the side chain of lysine 33 (LYS 33), which loses proton between pH 9 and pH 10; and the side chain of glutamic acid 35 (GLU 35), which loses proton between pH 4 and pH 5. We should note that HIS 15 remains neutral (i.e., deprotonated) in the whole considered pH range. According to NMR measurements, 78 the pK a of HIS 15 in hen egg white lysozyme is 5.35 ± 0.07. The reason for this discrepancy is the following. The lysozyme structure that we use in our study to construct the protein model results from molecular dynamics simulations in which amino acids are assigned default protonation states at neutral pH (see the Supporting Information), with histidine residues deprotonated. Consequently, due to the protonation− structure relationship (charge states of protein ionizable groups drive structural changes within the protein 79 ), the deprotonation of HIS 15 is energetically favorable. This holds true in the whole considered pH range, at least within the framework of the applied Poisson−Boltzmann model. We believe this not to be an issue as we are not interested in quantitative predictions of protonation equilibria of particular residues of lysozyme; however, one should be aware of limitations of the described algorithm, resulting from the neglect of the protein flexibility and the coupling between protonation equilibria and conformational changes. Average charges of lysozyme ionizable groups calculated over cpH BD trajectories generated for different values of the solution pH and the ionic strength are shown in Table 2. There are 14 ionizable groups whose protonation state fluctuates significantly in cpH BD simulations.
We first consider results of cpH BD simulations with the unbiased protonation sampling. Unlike in the described above case of lowest energy protonation states, influence of the solution ionic strength is now clearly visible. For a given value of the solution pH, to within the numerical accuracy and roundoff errors, absolute values of average charges of basic ionizable groups (i.e., NTER, LYS) at the 150 mM ionic strength are higher than at the 10 mM ionic strength, whereas absolute values of average charges of acidic groups (i.e., CTER, ASP, GLU, TYR) at the 150 mM ionic strength are lower than at the 10 mM ionic strength. This is a consequence of the increased, at the higher solution ionic strength, screening of electrostatic interactions within the protein. Screening of repulsive interactions between charged basic groups increases the probability of their protonation. Screening of repulsive interactions between charged acidic groups decreases the probability of their protonation. Finally, screening of attractive interactions between charged basic groups and charged acidic groups increases probability of protonation states in which basic and acidic groups are neutral. These three effects, together with the spatial distribution of ionizable groups within the protein and a prevalent number of basic ionizable groups, are responsible for the described dependence of average charges on the solution ionic strength.
A different behavior is observed for results of cpH BD simulations with the biased protonation sampling (Table 2). Free energies of different protonation states of the protein depend now not only on the solution pH and the ionic strength but also on the position and orientation of the protein with respect to the charged plane. The dependence of average charges of lysozyme ionizable groups on the ionic strength is different than in the case of cpH BD simulations with the unbiased protonation sampling. Now, for a given value of the solution pH, absolute values of average charges of basic residues at the 150 mM ionic strength are lower than at the 10 mM ionic strength, whereas absolute values of charges of acidic residues at the 150 mM ionic strength are higher than at the 10 mM ionic strength. This stems from the fact that proteinsurface electrostatic interactions affect values of intrinsic pK a 's of titratable groups. It is visible from Table 2 that the coupling to the external electric field shifts protonation equilibria of both basic and acidic groups toward higher pH values. Taking NTER at the 10 mM ionic strength as an example, we see that simulations with the unbiased protonation sampling lead one to determine that the pH value, which corresponds to 50% probability of deprotonation (or an average charge of 0.5 e), lies between 8 and 9. Results of simulations with the biased protonation sampling give the corresponding pH value between 9 and 10. For GLU 35, at the 10 mM ionic strength, the pH value corresponding to 50% probability of deprotonation (or an average charge of −0.5 e) lies between 4 and 5 based on results of cpH BD simulations with the unbiased protonation sampling and between 5 and 6 based on simulations in which the coupling between protonation equilibria and the external field is accounted for. It follows from the definition of the acidic dissociation constant that the observed shifting of pH, at which populations of protonated and deprotonated forms of a given group are equal, toward higher values, corresponds to an upward-shifting of the group's pK a . Shifts of protonation equilibria toward higher pH values are observed for basic and acidic ionizable groups for both considered values of the solution ionic strength; however, smaller shifts are observed at the higher ionic strength ( Table  2). This can be explained as follows. Any basic group, when protonated, is positively charged and interacts with the negatively charged surface, which results in a lower electrostatic energy in comparison with the situation in which the considered basic group is deprotonated (neutral). It is thus favorable, from the energetic standpoint, for the basic group to be in the protonated state, which translates to an increase in its pK a value. On the other hand, an acidic group is neutral in the protonated state and negatively charged in the deprotonated state. Electrostatic interaction between the deprotonated acidic group and the negatively charged surface results in a higher energy in comparison with the situation in which the group is protonated. So, for the acidic group, it is also favorable to be in the protonated state, which again translates to an increase in its pK a value. As electrostatic interactions with the charged surface are screened by dissolved ions (equation S20), shifts in pK a values decrease with an increasing ionic strength.
Mean square fluctuations of average charges of lysozyme titratable groups (see equation S28 in the Supporting Information) are shown in Table S2. Again, we compare results of the two types of cpH BD simulations, which differ in the way intrinsic pK a 's of lysozyme ionizable groups are evaluated. Generally, cpH BD simulations with the biased protonation sampling result in larger overall fluctuations of charges, which are clearly visible for residues with pK a values falling within the considered pH range, especially at the 10 mM ionic strength. The increase in magnitudes of charge fluctuations, due to the coupling of protonation equilibria to the electric field of the charged surface, is less pronounced at ACS Omega http://pubs.acs.org/journal/acsodf Article the 150 mM ionic strength due to the weakening of lysozymesurface electrostatic interactions. Uptake and release of protons by ionizable groups causes the average net charge of lysozyme to be dependent on the solution pH. In Figure 2, we show such dependencies resulting from the two types of cpH BD simulations. For a comparison, for each pH, we show also the value of the lysozyme total charge in the protonation state having the lowest energy (evaluated in the absence of the charged plane). We will use these total charges as a reference when discussing results of cpH BD simulations with the biased and unbiased protonation sampling. As we have stated above, in the absence of the charged plane, the same most probable protonation states were obtained at 10 and 150 mM ionic strengths and thus pH dependencies of the lysozyme total charge shown in the left and the right panel of Figure 2 are exactly the same. This is a feature of the lysozyme-like model protein considered in this work and cannot be generalized to other proteins.
Let us start with discussing results of cpH BD simulations in which the unbiased protonation sampling is performed. Depending on the solution pH, the protein samples protonation states that result in net charges that can be lower, higher, or equal to the reference total charge. However, as the isoelectric point of lysozyme lies above pH 10, in cpH BD simulations, we observe only positive net charges of the protein. Differences between reference total charges and average net charges resulting from cpH BD simulations for the solution pH values below 6 are much smaller at the ionic strength of 10 mM than at the ionic strength of 150 mM, whereas the opposite is observed for the solution pH above 7 ( Figure 2). An interesting behavior can be observed for solutions having pH 4 and pH 5 and the ionic strength of 10 mM. The average net charge at pH 4 is slightly lower than the reference value, whereas at pH 5, the average net charge is slightly higher than the reference value. At the 150 mM ionic strength, both average net values are substantially larger than corresponding reference values. In the solution having pH 4 and the 10 mM ionic strength, the unbiased protonation sampling gives rise to a population of proteins with net charges lower than the total charge in the most probable state, whereas at pH 5, there exists a population of proteins with net charges higher than the total charge in the most probable state. At the 150 mM ionic strength, sampling of protonation states results in larger populations of proteins with net charges higher than the reference value, both at the solution pH 4 and the solution pH 5. Consequently, at the 150 mM ionic strength, for both pH values, cpH BD simulations with the unbiased sampling result in average net charges of lysozyme that are higher than reference values. At the solution pH 7, the average net charge is almost exactly the same as the reference total charge, regardless the solution ionic strength. This results from the fact that, at this pH, and for both considered values of the ionic strength, the free energy gap between the most probable and the second most probable protonation state is so large that states above the most probable one are rarely accepted in cpH BD simulations.
In cpH BD simulations with the biased protonation sampling electrostatic interactions with the charged plane shift intrinsic pK a 's of lysozyme titratable groups toward higher values than these obtained as a result of the unbiased protonation sampling. Thus, for the whole considered pH range, and regardless the ionic strength, one may expect cpH BD simulations conducted with the biased protonation sampling to result in higher average net charges of lysozyme than cpH BD simulations conducted with the unbiased protonation sampling, with differences observed at the 150 mM ionic strength being smaller than these observed at the 10 mM ionic strength. Such a behavior is indeed visible in Figure  2. Smallest differences between average net charges resulting from the two cpH BD approaches are visible at pH 7 (the reason for that was already given above) and at pH 6 and the ionic strength of 150 mM. In this case, the ionic screening decreases the contribution of protein-surface interactions to values of intrinsic pK a 's (equation S4). Average net charges resulting from cpH BD simulations with the biased sampling conducted at pH 8 and pH 9, and the 150 mM ionic strength, are clearly smaller than corresponding values of total charges in the most probable protonation states. The reason for this is twofold. For these pH values in cpH BD simulations, the most dominant are populations of proteins whose net charges are either equal to charges in most probable protonation states or smaller; additionally, due to a substantial ionic screening, shifts in intrinsic pK a 's are too small to result in an increase in lysozyme average net charge above the reference total value.
In the case of cpH BD simulations with the unbiased protonation sampling, average net charges obtained for the 10 mM ionic strength are lower than those obtained for the 150 mM ionic strength (which, as we have explained above, is a consequence of the increased, at the higher ionic strength, screening of repulsive electrostatic interactions within the protein), whereas the opposite is observed for cpH BD simulations with the biased protonation sampling (Figure 2). , results of cpH BD simulations with the unbiased protonation sampling. Fixed state, sum of charges assigned to titratable groups in the most probable at a given pH (in the absence of the charged plane) protonation state (see data presented in Table 1).

Equilibrium Distributions of Electrostatic Interaction Energies between Lysozyme and the Charged
Plane. Values of energies of protein-surface electrostatic interactions, registered during different types of BD simulations, and for different solution conditions, served us to create distributions presented in Figure 3. Here, we show results of simulations conducted at solutions having the ionic strength of either 10 or 150 mM and pH 5, 7, and 9. Additionally, in Figure S1, we show energy distributions obtained for ionic strengths of 10 and 150 mM and pH 4, 6, 8, and 10. As expected, the introduction of the protonation degree of freedom affects values and distributions of registered energies.
At the 10 mM ionic strength, distributions obtained from BD simulations in which the protein is assigned a fixed protonation state have only a single peak, whose position moves, in a regular fashion, toward less negative energy values as the solution pH increases. Distributions obtained from cpH BD simulations have, depending on the solution pH, one, two, or even three peaks, which are more or less overlapping. Similarly, as in the case of distributions resulting from fixed-state BD simulations, centers of these distributions move toward less negative values with the increasing solution pH. For pH 4, pH 9, and pH 10, distributions obtained from the three types of BD simulations are significantly different from each other. For pH 5, pH 6, and pH 8, there are similarities in pairs of distributions resulting from different simulation methods. Distributions obtained from the three types of simulations for the solution pH 7 are quite similar to each other. When the ionic strength is increased to 150 mM, differences between energy distributions become overall less pronounced.
Below, we discuss effects of protonation equilibria on energy distributions using as examples mostly results of simulations conducted at pH 4 and pH 5.
We begin with the comparison of distributions obtained from BD simulations in which the protonation state of lysozyme is fixed and cpH BD simulations with the unbiased protonation sampling. At the solution pH 5 and the ionic strength of 10 mM (Figure 3), in both energy distributions, there is a main peak located roughly at −15 kcal/mol. The main difference between these two distribution is the much Bottom: pH = 9. pK a, intr mod , distributions obtained from cpH BD simulations with the biased protonation sampling. pK a, intr o , distributions obtained from cpH BD simulations with the unbiased protonation sampling. Fixed state, distributions obtained from BD simulations in which the protein is assigned the most probable at a given pH (in the absence of the charged plane) protonation state (see data presented in Table 1).

ACS Omega
http://pubs.acs.org/journal/acsodf Article smaller peak, at −18 kcal/mol, which is present only in the distribution resulting from cpH BD simulations. The total charge of lysozyme, in its most probable state at pH 5, is +8 e (Figure 2). In line with values of average charges of lysozyme ionizable groups and their fluctuations calculated from cpH BD simulations (Table 2 and Table S2), in these simulations, the dominant population consists of proteins whose total charge is +8 e due to the deprotonation of either one of acidic residues: GLU 7, GLU 35, ASP 18, and ASP 101 (most probably GLU 35 for which at pH 5 we observe the largest charge fluctuations of 0.1 e 2 ). There is also a smaller population of proteins in which these residues are protonated and whose total charge is +9 e. This agrees with the average net charge of lysozyme at pH 5 and the 10 mM ionic strength being only slightly larger than the total charge in the most probable protonation state ( Figure 2). The small peak visible in the distribution obtained from cpH BD simulations with the unbiased sampling corresponds to this smaller population of proteins. At pH 4 and the 10 mM ionic strength, the average net charge of lysozyme is somewhat smaller than the total charge in the most probable state (Figure 2). Again, it is GLU 35 that shows the largest charge fluctuations from all acidic residues (0.48 e 2 , Table S2). In the energy distribution resulting from cpH BD simulations with the unbiased protonation sampling, there are now two peaks ( Figure S1) of nearly equal size. The first one located roughly at −18 kcal/mol coincides with the single peak in the distribution obtained from fixed-state BD simulations; the second one is located at roughly −15 kcal/mol. The first peak corresponds to the population of proteins whose total charge is +9 e (protonated GLU 35), whereas the second peak corresponds to the population of proteins whose total charge is +8 e (deprotonated GLU 35). As the average charge of GLU 35 is −0.48 e, these two populations are comparable in size, hence the similar height of both peaks. One may thus conclude, somewhat naively, as energy distributions are created for ensembles of different protein-plane configurations, that for the 10 mM ionic strength, an addition of one positive charge, shifts protein-surface interactions energy by approximately −3 kcal/mol, whereas removal of one positive charge causes a shift of +3 kcal/mol. Moreover, a single peak in the energy distribution appears to be related to a population of proteins with a particular value of the net charge. While we are able to explain the pH dependence of energy distributions resulting from cpH BD simulations as a result of the protonation behavior of GLU 35, we should consider a minor contribution resulting from, less probable at pH 4 and pH 5, protonations of GLU 7, ASP 18, and ASP 101 ( Table 2 and Table S2). At pH 4, protonations of these residues feeds the population of proteins whose net charge is +10 e, contributing to the tail visible on the left side of the cpH BD energy distribution ( Figure S1). We now move to distributions obtained from cpH BD simulations with the biased protonation sampling. At the solution pH 5 and the ionic strength of 10 mM, the energy distribution obtained from these simulations (Figure 3) has three peaks. The less distinct from the three, at −15 kcal/mol, overlaps with the single peak of the distribution obtained from fixed-state BD simulations (which correspond to the lysozyme net charge of +8 e). The second peak, at −18 kcal/mol, which is also the highest, coincides with the smaller of the two peaks of the distribution resulting from cpH BD simulations with the unbiased protonation sampling (which correspond to the lysozyme net charge of +9 e). The third peak is located slightly below −20 kcal/mol. This third peak corresponds to a population of proteins whose net charge is two units larger (+10 e) than the total charge in the most probable state. From Figure 2, we see that the average net charge of lysozyme resulting from cpH BD simulations with the biased protonation sampling is more than one unit of e higher than the total charge in the most probable state. From Table 2 and  Table S2, we can identify two resides contributing to the observed increase in the lysozyme net charge, GLU 7 (the average charge of −0.68 e and the magnitude of mean square fluctuation of 0.32 e 2 ) and GLU 35 (−0.21 e and 0.79 e 2 ). Protonation of these two residues feeds populations of proteins with net charges of +10 e and +9 e. At the solution pH 4 and the 10 mM ionic strength, the average net charge of lysozyme resulting from cpH BD simulations with the biased protonation sampling is roughly 1.5 e higher that the total charge in the most probable state (Figure 2), which signifies contributions from protonation of at least two residues, which are deprotonated in the most probable state. The energy distribution has two overlapping peaks ( Figure S1): one located at −22 kcal/mol and one located slightly below −20 kcal/mol. The right tail of the distribution overlaps with the distribution resulting from fixed-state BD simulations. Its left tail reaches −26 kcal/mol. Such a distribution can be explained as resulting from contributions of four populations of proteins with total charges of +11 e, +10 e, +9 e, and +8 e. From Table  2 and Table S2, we may identify the four residues involved. Finally, we would like to comment on energy distributions resulting from the three types of BD simulations conducted for the 150 mM ionic strength (Figure 3 and Figure S1). In this case, distributions obtained for a given pH substantially overlap. At pH 5 and the 150 mM ionic strength, the average net charge of lysozyme resulting from cpH BD simulations with the biased protonation sampling is higher than the total charge in the most probable state by roughly 1 e (Figure 2). For cpH BD simulations with the unbiased sampling, the increase relative to the most probable state is roughly 0.5 e. In the case of the cpH BD simulations with the biased sampling, contributions to the increase in the average net charge of lysozyme come from protonations of GLU 35 (−0.57 e and 0.43 e 2 ), GLU 7 (−0.87 e and 0.13 e 2 ), ASP 101 (−0.92 e and 0.08 e 2 ), and ASP 18 (−0.96 e and 0.04 e 2 ) ( Table 2 and Table  S2 Difference between the sampling of protonation states in the two types of cpH BD simulations at the 150 mM ionic strength is rather small. Moreover, populations of proteins with net charges larger than the total charge in the most probable state diminish. Consequently, at a given pH, the three considered energy distributions are not that different from each other. Multimodal character, that is now visible also for distributions obtained from the fixed-state BD simulations, is not due to the existence of two charge states of the protein but reflects some orientational preferences of the protein above the plane. In addition, in Figures S2 and S3, distributions obtained for solutions having pH 4, 6, 8, and 10 are presented. Position is defined as the distance between the diffusion center of the protein (which is also the center of its hydrodynamic ellipsoid) from the charged surface and denoted z. As a partial measure of the orientation, we choose the angle between the normal to the plane and the shorter semi-axis of the lysozyme ellipsoid along one of the axes of the lysozyme-fixed coordinate frame. This angle is denoted α. While such a definition does not contain full information about the orientation of the protein relative to the surface, it is sufficient to illustrate differences in results obtained from different types of BD simulations. Depending on the solution pH, position-orientation distributions have up to three peaks whose maxima are located roughly at (α, z): (0.0,1.0), (1.5,1.0), and (2,7,0.8). At lower pH values, the most prominent peak of the three is the one located at (2.7,0.8). The other two peaks, barely visible at lower pH values, increase upon increasing the solution pH.
We first consider results of BD simulations conducted at the 10 mM ionic strength (Figure 4 and Figure S2). At the solution pH in the range 4−6, shapes of distributions obtained from the three types of BD simulations are qualitatively quite similar to each other, with peaks located at (0.0,1.0) and (1.5,1.0) barely visible. At the solution pH 4, cpH BD simulation with the biased protonation sampling result in a distribution in which the (2.7,0.8) peak is notably higher than in distributions obtained either from BD simulations in which the lysozyme protonation state is fixed or from cpH BD simulations with the unbiased protonation sampling. Moreover, the (2.7,0.8) peak in the distribution obtained from cpH BD simulations with the unbiased protonation sampling is lower than the (2.7,0.8) peak in the distribution obtained from fixed-state simulations. This agrees with energy distributions obtained at pH 4 and the 10 mM ionic strength. As we explained above, at the solution pH 4, the total charge of lysozyme in the most probable protonation state is +9 e. In cpH BD simulations with the unbiased protonation sampling, we observe populations of proteins with net charges of +9 e and +8 e and the average net charge of lysozyme is slightly ACS Omega http://pubs.acs.org/journal/acsodf Article below +9 e. In cpH BD simulations with the biased protonation sampling, we observe populations of proteins with net charges of +11 e, +10 e, +9 e, and +8 e, and the average net charge of lysozyme is nearly +11 e. Thus, the height of the (2.7,0.8) peak in position-orientation distributions is directly related to the net charge of lysozyme and the overall strength of its interactions with the charged surface. At pH 5, we observe that the (2.7,0.8) peak is the highest in the case of cpH BD simulations with the biased protonation sampling and the lowest in the case of fixed-state simulations. Again, this can be related to energy distributions and our analysis of lysozyme charges. As expected, at pH 7, distributions obtained from the three types of BD simulations are practically the same. Largest differences between positionorientation distributions are observed for the solution pH in the range 8−10, where we also observe a significant increase of (0.0,1.0) and (1.5,1.0) peaks, in part as a consequence of the complete deprotonation of almost all acidic groups of lysozyme (with the exception of TYR 20 and TYR 23, see Table 2 and  Table S2). At pH 8, distributions obtained from fixed-state BD simulations and cpH BD simulations with the biased protonation sampling are quite similar, in line with corresponding energy distributions depicted in Figure S1. However, the distribution resulting from cpH BD simulations with the unbiased protonation sampling is different from the other two, with (0.0,1.0) and (1.5,1.0) peaks being more pronounced. Lastly, distributions obtained from the three types of BD simulations at pH 9 and pH 10 are all clearly different from each other. The behavior of position-orientation distributions for the solution pH in the range 8−10 is related to the protonation behavior of basic groups: NTER, LYS 33, and LYS 116, the acidic TYR 23, and to a smaller extent LYS 1, LYS 96, and LYS 97 (see Table 2 and Table S2). When the solution ionic strength is increased to 150 mM ( Figure 5 and Figure S3), differences between distributions obtained from different types of BD simulations become less apparent. It can be noted however that, at pH 9, positionorientation distributions obtained from fixed-state BD simulations and cpH BD simulations with the biased protonation sampling have similar shapes, whereas the distribution obtained from cpH BD simulations without the coupling diverges from the other two as its (0.0,1.0) and (1.5,1.0) peaks are more marked and the (2.7,0.8) peak is lower. This agrees with energy distributions shown in Figure 3. The distribution obtained from cpH BD simulations with the biased protonation sampling, conducted for the solution pH 10, is clearly different from the other two, which again is in agreement with energy distributions obtained for this pH value, presented in Figure S1. 3.4. Lysozyme Rotational Dynamics. Rotational dynamics of lysozyme was investigated in terms of correlation functions C 2, x (t) (equation S25) and orientation distributions G x (θ, t) (equation S26). The ensemble of performed BD simulations was extended as we considered also diffusion of lysozyme near a neutral plane. C 2, x (t) curves resulting from BD simulations performed for solutions having pH 5, 7, and 9 and ionic strengths 10 mM and 150 mM are presented in Figure 6, whereas results of simulations performed for solutions having pH 4, 6, 8, and 10 are shown in Figure S4. In both figures, as a reference, we also show the C 2, x (t) function constructed from BD trajectories of lysozyme near the neutral plane. With the plateau value converging to zero, this curve corresponds to free diffusion. 46 From the comparison of curves obtained from BD simulations in which the plane confining the solution is charged with the curve obtained from simulations in which the plane is neutral, we see a much slower convergence of plateau values to zero in the former case. Electrostatic interactions with the plane keep the diffusing molecule in the vicinity of the boundary and as the shape and charge distribution of the molecule are not isotropic (see Figure 1) its orientational freedom, close to the boundary, becomes restricted as a result of both steric and electrostatic interactions. Changing the orientation may require diffusing away from the charged plane. Consequently, rotational dynamics of the molecule is considerably diminished as evidenced by slower, in comparison with the curve obtained for the neutral surface, decays of autocorrelation functions in Figure 6 and Figure S4. Weakening of protein-plane electrostatic interactions, as a consequence of either increasing the solution pH (which results in a decreasing average net charge of the protein, Figure  2) or increasing the ionic strength, gradually abolishes orientational restrictions, which in turn results in faster dynamics.
At the 10 mM ionic strength, different types of BD simulations in which the plane is charged result in disparate decay curves for solutions having pH 4, 5, 8, 9, and 10. At the solution pH 5, cpH BD simulations with the biased protonation sampling result in the slowest and the most restricted rotational dynamics, whereas fixed-state BD simulations and simulations with the unbiased protonation sampling result in curves that are close to each other and correspond to a faster decay. A similar behavior is observed at pH 9 and pH 10. At pH 4 and pH 8, autocorrelation functions obtained from fixed-state simulations and from cpH BD simulations with the biased protonation sampling are close to each other and correspond to slower decays than curves obtained from cpH BD simulations with the unbiased protonation sampling. C 2, x (t) curves obtained for pH 6 and pH 7 are almost indistinguishable from each other.
At the 150 mM ionic strength, differences between decays are less visible (Figure 6 and Figure S4). As the electrostatic interactions are weaker due to the ionic screening, the protein is able to explore the space further away from the plane more effectively. Its rotational diffusion is less restricted and, to a much smaller extent, affected by electrostatic interactions and the treatment of pH effects adopted in simulations.
In Figure S5, Figure 7, and Figure S6, we present G x (θ, t) functions (obtained for solutions having pH 5, 7, and 9 and ionic strengths of 10 and 150 mM) for time intervals 1, 100, and 1000, respectively.
In the case of BD simulations with the neutral plane we deal with free diffusion and consequently G x (θ, t) distributions are characterized by smooth curves with a single maximum that shifts from 0 to π/2 radians with an increasing t. For t = 100 ns and t = 1000 ns, orientation distributions are symmetric as, after 100 ns, all orientations of the molecule, which was initially in a particular orientation, are equally probable.
Distributions obtained from different types of BD simulations in which the boundary plane is charged, for t = 1 ns ( Figure S5) are, for a given solution pH value and a given ionic strength, almost indistinguishable from each other. As , results of cpH BD simulations with the unbiased protonation sampling. Fixed state, results of simulations in which the protein is assigned the most probable at a given pH (in the absence of the charged plane) protonation state. Neutral, results of BD simulations in which the plane is uncharged. All distributions obtained for the time interval of 100 ns. electrostatic and steric interactions restrict rotational dynamics of the molecule, these distributions are narrower than the one describing diffusion near the neutral plane and their maxima are shifted toward smaller values of θ. Weakening of proteinplane electrostatic interactions, with an increasing solution pH and an increasing ionic strength, results in widening of G x (θ, t) curves. Orientation distributions obtained for t = 100 ns (Figure 7) are all asymmetric, unlike the distribution obtained based on simulations in which the plane is neutral. In agreement with decay curves shown in Figure 6, at pH 5 and pH 9, regardless the ionic strength, rotational dynamics resulting from cpH BD simulations with the biased protonation sampling appears to be more restricted than dynamics observed in BD simulations with the unbiased protonation sampling and in fixed-state BD simulations. Largest deviations between the curve obtained from cpH BD simulations with the biased protonation sampling and the other two is visible for the solution pH 9 and the 10 mM ionic strength; we note that, in this case, all distributions show bimodal structure with peaks located roughly at θ = π/6 and θ = π/2. The first peak corresponds to the restricted diffusion close to the surface whereas the second peak signifies unrestricted diffusion away from the boundary. Such a structure is also visible at the solution pH 9 and the 150 mM ionic strength, although to a much smaller extent. At pH 7, all distributions are close to each other. Overall, similar observations can be made for distributions obtained for t = 1000 ns ( Figure S6), although we should note that the multimodal character of distributions is now visible also at lower values of the solution pH. Moreover, at pH 9 and the 150 mM ionic strength, distributions obtained from simulations with the charged plane nearly converge to the distribution resulting from simulations in which the plane is neutral. However, for considered time intervals, rotational diffusion of lysozyme near the charged plane never appears to be free.

CONCLUSIONS
Our aim here was to answer two questions: the first, how does the solution pH affect electrostatic interactions of proteins with charged surfaces and near-surface diffusive dynamics of proteins, and the second, what are the possible consequences of neglecting pH-dependent protonation equilibria in computational algorithms that are being applied to study protein-surface systems.
Regarding the first question, we have considered that pHdependent variations of protonation states of protein ionizable groups are coupled to positions and orientations of the protein relative to the surface, through electrostatic interactions of protein charges with the charged plane. Employing constant-pH BD simulations of a positively charged protein diffusing in an ionic solution bounded by a negatively charged planar surface, we have shown that protein-plane electrostatic interactions shift values of pK a 's of protein ionizable groups toward higher values. Consequently, the average net charge of the protein is increased in comparison with its average net charge in an unbounded solution having the same pH value. pH-dependent fluctuations of the protein net charge and charge distribution modulate confinement effects of the bounding plane, which elongates the time scale of rotational dynamics of the molecule. Moreover, as the protein switches between different protonation states and charge distributions, its rotational dynamics is influenced by electrostatic torques of fluctuating directions and magnitudes. Equilibrium distributions of the electrostatic energy of protein-plane interactions, as well as distributions of positions and orientations of the protein relative to the plane, are strongly pH-dependent.
Regarding the second question, the comparison of results of cpH BD simulations with results of a conventional approach, in which the protein is assigned a fixed protonation state, chosen to be the most probable one at a given solution pH and an ionic strength for the protein in an unbounded solution, clearly shows that, accounting for the protonation degree of freedom, in a way which is sound from the standpoint of physics, at certain pH values results in a drastically different picture of the studied protein-surface system, which we have illustrated using as examples equilibrium energy distributions, equilibrium position-orientation distributions, rotational correlation functions, and time-dependent orientation distributions.
Differences between results of cpH BD simulations and fixed-state BD simulations presented in our work diminish with the increasing solution ionic strength, and differences observed at the 150 mM ionic strength may, at first glance, appear not to be that significant. One should keep in mind, however, that all analyzed here quantities are evaluated over an ensemble of protein-surface configurations, in which protein-surface distances span the range between zero and a dozen times the Debye−Huckel length at the 150 mM ionic strength (∼8 Å). The Debye−Huckel length, which defines the effective range of electrostatic interactions in an ionic solution, decreases with an increasing ionic strength. In our simulations, at the higher ionic strength, the protein overall diffuses further from the boundary, effects of the protein charge distribution on its dynamics are less pronounced due to the weakening of electrostatic interactions, and differences between results of cpH BD and fixed-state BD simulations appear to be smaller. The closer the protein is to the charged boundary, the more electrostatic interactions and the solution pH are bound to affect protein dynamics. We believe that modeling of proteins diffusing in the presence of charged obstacles, and adsorbing on charged boundaries, requires the treatment of pH effects that extends beyond the fixed-state approach.
Detailed description of the cpH BD algorithm; details of molecular dynamics simulations of lysozyme in the CPT ensemble; details of PB calculations; definitions of correlation functions C 2, x (t) and orientation distributions G x (θ, t); description of calculations of mean square fluctuations of average charges of lysozyme titratable groups; Table S1: location of charges assigned to lysozyme titratable groups; Table S2: mean square fluctuations of charges of lysozyme titratable groups at different solution pH values; Figure S1