# Brownian Dynamics Simulations of Proteins in the Presence of Surfaces: Long-
range Electrostatics and Mean-field Hydrodynamics

Martin Reinhardt  [  Neil J. Bruce  [  Daria B. Kokh  [  Rebecca C. Wade  [
<EMAIL_ADDRESS>
###### Abstract

Simulations of macromolecular diffusion and adsorption in confined
environments can offer valuable mechanistic insights into numerous biophysical
processes. In order to model solutes at atomic detail on relevant time scales,
Brownian Dynamics simulations can be carried out with the approximation of
rigid body solutes moving through a continuum solvent. This allows the
precomputation of interaction potential grids for the solutes, thereby
allowing the computationally efficient calculation of forces. However,
hydrodynamic and long-range electrostatic interactions cannot be fully treated
with grid-based approaches alone. Here, we develop a treatment of both
hydrodynamic and electrostatic interactions to include the presence of
surfaces by modeling grid-based and long-range interactions. We describe its
application to simulate the self-association and many-molecule adsorption of
the well-characterized protein Hen Egg-White Lysozyme to mica-like and silica-
like surfaces. We find that the computational model can recover a number of
experimental observables of the adsorption process and provide insights into
their determinants. The computational model is implemented in the Simulation
of Diffusional Association (SDA) software package.

###### keywords:

Macromolecular Diffusion, Hydrodynamics, Electrostatics, Confinement,
Macromolecular Crowding, Biomolecular Association

HITS] Molecular and Cellular Modeling Group, Heidelberg Institute for
Theoretical Studies (HITS), Schloß-Wolfsbrunnenweg 35, 69118 Heidelberg,
Germany [MPI] Theoretical Biophysics Group, Max Planck Institute for
Biophysical Chemistry, Am Fassberg 11, 37077 Göttingen, Germany HITS]Molecular
and Cellular Modeling Group, Heidelberg Institute for Theoretical Studies
(HITS), Schloß-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany HITS]Molecular
and Cellular Modeling Group, Heidelberg Institute for Theoretical Studies
(HITS), Schloß-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany HITS]Molecular
and Cellular Modeling Group, Heidelberg Institute for Theoretical Studies
(HITS), Schloß-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany Center for
Molecular Biology (ZMBH), University of Heidelberg, Im Neuenheimer Feld 282,
69120 Heidelberg, Germany Interdisciplinary Center for Scientific Computing
(IWR), Heidelberg University, Im Neuenheimer Feld 368, 69120 Heidelberg,
Germany +49-6221-533 247 +49-6221-533 298 IR,NMR,UV

##  1 Introduction

The diffusion of macromolecules, such as proteins, is essential for cellular
processes, including signal transmission and the transport of matter 1, 2, 3.
Diffusion-related properties are often highly dependent on the environment of
the molecules4, 5, 6, 7, 8, 9, with numerous biological consequences 10, 11.
Furthermore, interactions with surfaces influence the diffusive behavior of
these molecules. Firstly, in the case of attractive interactions, adsorption
processes occur, including protein aggregation at confining12 or artificial
surfaces such as implants 13, and these are of interest for their implications
for drug delivery and biosensing 14. Secondly, inside cells, surfaces restrict
diffusive motion, e.g., due to the plasma membrane surrounding the cellular
volume as a whole, or because of internal cytoskeletal elements. These
restrictions have been identified as a major determinant of macromolecular
kinetics and reactivity, and have been analyzed theoretically 15, 16, 4, as
well as experimentally 17, 18, 19. Predicting the effects of such surfaces on
protein diffusion is, however, a highly complex task 20, 21, 22.

A number of factors contribute to the diffusive behavior of molecules in the
presence of surfaces. Aside from direct steric obstruction, hydrodynamic
interactions (HI) with other molecules 23 and with surfaces have been shown to
lead to a significant reduction in macromolecular self-diffusion 24, 25, 26,
27. In addition, occurrences of non-Gaussian mean-squared displacements near
surfaces have been observed for colloidal solutions28, 29, 27. The adsorption
processes of molecules to surfaces are affected by long-range electrostatic
interactions and, therefore, depend on the charge density of the given surface
and the surrounding salt concentration30. Furthermore, short-range
interactions, such as van-der-Waals forces, also affect the orientation and
structure of molecules on surfaces 31, 32, 33.

Computational approaches can provide mechanistic insights for systems that are
often hard to access experimentally with techniques such as nuclear magnetic
resonance 34 or fluorescence spectroscopy 35. However, while molecular
dynamics (MD) simulations of dense protein solutions in atomic detail with
explicit solvent models have been performed in recent years 36, 37, 38, they
require very large computational resources, and the short lengths of the
simulated timescales often make it hard to obtain statistically significant
results on diffusion-related processes. Therefore, it is reasonable to employ
different levels of detail in computer simulations, depending on the length
and time scales of the system of interest. While MD simulations are suitable
for capturing short-range macromolecule-surface interactions 39, 40, Brownian
Dynamics (BD) simulations can be used to study the properties defining
processes occurring on longer time scales, such as the the kinetic and
structural properties of macromolecular adsorption41, 42. BD simulations have
been performed of spherical particle models of the macromolecules with a
charged surface43, 44, 45. So far, however, BD simulations with structurally
detailed molecular models have mostly been restricted to systems consisting of
one solute and a surface, although they have been used to simulate the
diffusion of many hydrophobin molecules to an uncharged surface46 and a few
diffusing hen egg white lysozyme (HEWL) molecules to a charged surface 47.

The Simulation of Diffusional Association (SDA) BD software package48, 49, 50
can be used to simulate the motion of pairs of macromolecules, of (dense)
solutions of macromolecules, and of macromolecules in the presence of
surfaces. SDA uses a rigid-body model, that permits intermolecular interaction
forces to be computed efficiently by precomputation of their interaction
potentials on three-dimensional discretized grids while at the same time
calculating the interactions on an atomic basis rather than approximating the
solute by, e.g., a sphere or an ellipsoid. HI between solutes is modelled by a
mean-field approximation described by Mereghetti and Wade 51. Here, we first
extend this approach to introduce a treatment of solute-surface HI by a first-
order approximation that can be used to simulate the adsorption of either a
single solute or of multiple solutes to a surface. Then, for long-range
electrostatic interactions that extend beyond the dimensions of the
electrostatic potential grids, we introduce a Debye-Hückel (DH) approximation
for surfaces. This treatment complements the DH approximation that was
previously introduced for long-range electrostatic interactions between
solutes 52, for which we here describe an improved treatment of the transition
between grid and DH formulations at the grid edges. We describe the validation
and application of these new computational models to test systems containing
the experimentally well-characterized protein, HEWL. We first compute and
analyze the rate of approach of two HEWL proteins, and then we simulate the
adsorption process of over a hundred HEWL molecules to attractive mica- and
silica-like surfaces for different bulk protein concentrations. These
simulations allow us to investigate the contributions of HI and long-range
electrostatic interactions to these protein-surface adsorption processes.

##  2 Theory and Methods

###  2.1 Brownian Dynamics Simulations

In BD simulations, a mesoscopic model is employed to simulate the motion of
solutes over length and time scales relevant for Brownian motion, that is the
random motion of solutes in fluids where the solutes move much more slowly
than the solvent molecules. An implicit solvent model is used: besides the
systematic forces, the effect of the solvent is included through stochastic
sampling of the collisions with solvent molecules. In SDA, solute trajectories
are propagated according to the algorithm described by Ermak and McCammon 53.
The trajectories are composed of successive displacements of the solutes, each
taken over a short time step $\Delta t$. The translational displacement of a
Brownian particle $i$ is calculated as

$\displaystyle\begin{aligned}
\mathrm{\mathbf{r}}_{i}^{n+1}=\mathrm{\mathbf{r}}_{i}^{n}&+\sum_{j}\frac{\partial\mathrm{\mathbf{D}}_{ij}}{\partial\mathrm{\mathbf{r}}_{j}^{n}}\>\Delta
t\\\
&+\sum_{j}\frac{\mathrm{\mathbf{D}}_{ij}}{k_{\textnormal{B}}T}\>\mathrm{\mathbf{F}}_{i}^{n}\>\Delta
t+\mathrm{\mathbf{R}}_{i}\;\end{aligned}$ (1)

where $\mathrm{\mathbf{r}}_{i}^{n}$ denotes the position of the center of
geometry of particle $i$. The superscript $n$ indicates that the variable is
to be evaluated at the beginning of the time step, $n+1$ is the result after
the time step. $\mathrm{\mathbf{F}}_{i}$ is the total systematic force acting
on particle $i$ and $\mathrm{\mathbf{D}}_{ij}$ is a $3\times 3$ subtensor of
the hydrodynamically-coupled diffusion tensor $\mathrm{\mathbf{D}}$ of the
system of Brownian particles, where the diagonal subtensor
$\mathrm{\mathbf{D}}_{ii}$ is the infinite dilution diffusion tensor of
particle $i$ and the off-diagonal subtensors account for the configuration-
dependent hydrodynamic coupling between particles $i$ and $j$.
$\mathrm{\mathbf{R}}_{i}$ is a stochastic displacement vector that takes into
account the collisions with the solvent molecules. It is drawn from a Gaussian
distribution with mean $<\mathrm{\mathbf{R}}_{i}>=0$ and covariance
$<\mathrm{\mathbf{R}}_{i}\mathrm{\mathbf{R}}_{j}>=2\mathrm{\mathbf{D}}_{ij}\Delta
t$ for all $i$ and $j$.

Propagating a system of Brownian solutes that each consist of a large number
of particles using Eq. 1 is computationally expensive, as the calculation of
the hydrodynamically-coupled stochastic term of the BD propagation step
requires Cholesky factorization of the tensor $\mathrm{\mathbf{D}}$ at every
time step, which scales as $\mathcal{O}(N^{3})$ for $N$ Brownian solutes,
although with approximations this can be reduced to $\mathcal{O}(N^{2.25})$54,
55 or $\mathcal{O}(N^{2})$56. This compares with $\mathcal{O}(N^{2})$ for the
calculation of the interparticle systematic forces, which can be reduced to
$\mathcal{O}(N)$ through the use of distance cutoffs.

As SDA was initially developed to simulate the association of a pair of
solutes,48 approximations were made to simplify Eq. 1. It was assumed that
each solute diffuses as a single spherical hydrodynamic bead and that
hydrodynamic interactions are negligible, due to the low solute concentration.
This meant that all off-diagonal terms in Eq. 1 could be ignored and the
diagonal subtensors could be could be replaced with scalar isotropic diffusion
translational coefficients $D^{\textnormal{t}}_{i}$, resulting in :

$\mathrm{\mathbf{r}}_{i}^{n+1}=\mathrm{\mathbf{r}}_{i}^{n}+\frac{D^{\textnormal{t}}_{i}}{k_{\textnormal{B}}T}\>\mathrm{\mathbf{F}}_{i}^{n}\>\Delta
t+\mathrm{\mathbf{R}}_{i}\;$ (2)

Importantly, the stochastic vector $\mathrm{\mathbf{R}}_{i}$ is no longer
configuration-dependent and can be drawn from a precomputed distribution. The
propagation of the rotation is calculated through an analogous equation :

$\displaystyle\begin{aligned}
\mathrm{\mathbf{w}}_{i}^{n+1}=\mathrm{\mathbf{w}}_{i}^{n}+\frac{D_{i}^{\textnormal{r}}}{k_{\textnormal{B}}T}\>\mathrm{\mathbf{T}}_{i}^{n}\>\Delta
t+\mathrm{\mathbf{W}}_{i}\end{aligned}$ (3)

where $\mathrm{\mathbf{w}}_{i}$ and $\mathrm{\mathbf{W}}_{i}$ describe the
orientation and the stochastic rotation vector, respectively, of solute $i$.
$D_{i}^{\textnormal{r}}$ denotes the rotational diffusion coefficient, and
$\mathrm{\mathbf{T}}_{i}^{n}$ describes the sum of torques acting on solute
$i$.

###  2.2 Mean-Field Intersolute Hydrodynamic Interactions

As SDA was extended to allow the simulation of protein solutions57, it became
apparent that the assumption that intersolute hydrodynamic interactions could
be ignored was no longer valid for increasing concentration of solutes.
Therefore a mean-field hydrodynamic model51, 58, 59 was developed in which
$D_{i}^{\textnormal{t}}$ is replaced by a local occupied volume fraction-
dependent diffusion coefficient $D_{i}(V_{i}^{\textnormal{frac}})$ and the
translational displacement equation becomes 51:

$\mathrm{\mathbf{r}}_{i}^{n+1}=\mathrm{\mathbf{r}}_{i}^{n}+\frac{D_{i}(V_{i}^{\textnormal{frac}})}{k_{\textnormal{B}}T}\>\mathrm{\mathbf{F}}_{i}^{n}\>\Delta
t+\mathrm{\mathbf{R}}_{i}$ (4)

where $V_{i}^{\textnormal{frac}}=\sum_{j}\nu_{j}/V_{i}$ denotes the
dimensionless local occupied volume fraction around solute $i$. It is obtained
by summing over the volumes $\nu_{j}$ of the surrounding solutes $j$,
calculated by approximating these as spheres of radius $a_{j}$. The sum only
includes the solutes within a sphere with radius $R_{\textnormal{cut}}$ with
volume $V_{i}=(4\pi/3)\>(R_{\textnormal{cut}})^{3}$ centered on the solute
$i$. The volume fraction dependent short-time translational diffusion
coefficient $D_{i}(V_{i}^{\textnormal{frac}})$ is then obtained using the
Tokuyama model 58, 59, derived for a concentrated hard-sphere suspension of
solutes interacting with both direct and hydrodynamic interactions. An
equation analogous to Eq. 4 is used for the rotational motion with the volume
fraction dependent short-time rotational diffusion coefficient obtained by
using the model derived by Cichocki et al. 60, which includes the lubrication
forces as well as two- and three-body expansions of the mobility functions.
For a larger number of solutes, the approach correctly reproduces the short
and long-time diffusion coefficients51 while avoiding the computationally
expensive Cholesky factorization.

###  2.3 Hydrodynamic Interactions in the Presence of a Surface

While the intersolute mean-field approach allows concentrated solute solutions
to be simulated, it is not able to model the effects that obstructing surfaces
have on the diffusion of nearby solutes. In this work, we extend this model to
account for these effects. In doing so, we distinguish between two different
HI effects. Firstly, a single solute near a surface experiences a large
reduction in its diffusion due to HI with the surface, i.e., the interaction
with its own reflected flow field. Secondly, a solute in a crowded environment
also interacts with the reflected flow fields from the surrounding solutes.
Naturally, these two effects are correlated. However, computing HI correctly
up to a high order becomes difficult both analytically and computationally for
as few as two24 or three spheres in the presence of a surface, and, therefore,
is often not accounted for in BD simulations at all. In order to estimate the
magnitude of the effect of HI on many solutes close to a surface, we combine
the above two effects in a linear first-order approach, thereby approximating
them as uncoupled.

The effect of HI on the diffusive behavior of a single sphere in the presence
of a planar surface is analytically accessible and, due to direct coupling of
its movement with the change of its image flow field, is known to lead to a
large reduction in the effective diffusion coefficient close to the surface.
Due to the symmetries of the problem, the translational propagation step can
be easily split into a step in the plane parallel to the surface and a step in
the direction perpendicular to the surface. The relative reduction,
$\zeta_{\perp}$, in the short-time diffusion coefficients of a sphere
perpendicular to the surface has been determined by Lorentz 61 and Faxen 62 :

$\displaystyle\begin{split}\zeta_{\perp}(z)=&\;\frac{9}{16}\frac{a}{z}+\frac{1}{8}\left(\frac{a}{z}\right)^{3}\\\
&-\frac{45}{256}\left(\frac{a}{z}\right)^{4}-\frac{1}{16}\left(\frac{a}{z}\right)^{5}\end{split}$
(5)

where $a$ is the sphere’s radius and $z$ the height of the sphere’s center
above the surface. They also calculated a result for the relative reduction in
the diffusion coefficient parallel to the surface
($\zeta_{\parallel}(z)=1-\frac{9}{8}\frac{a}{z}$) which is, however, only
valid at $a/z<1/4$, i.e., when the sphere’s center is farther away from the
surface than four times the sphere’s radius. More refined calculations were
conducted by Brenner et al. 63, 64, leading to the analytical result:

$\begin{split}&\zeta_{\parallel}(z)=\left[\frac{4\sinh\alpha}{3}\sum_{n=0}^{\infty}\frac{n(n+1)}{(2n-1)(2n+3)}\right.\\\
&\left.\left(\frac{2\sinh[(2n+1)\alpha]+(2n+1)\sinh[2\alpha]}{(2\sinh[(n+1/2)\alpha])^{2}-((2n+1)\sinh\alpha)^{2}}-1\right)\right]^{-1}\end{split}$

(6)

where $\alpha=\cosh^{-1}((a/z)^{-1})$. A numerical integration at every time
step in the simulation would be too costly computationally. Therefore, we
numerically integrate the expression for $\zeta_{\parallel}(z)$ in Eq. 6 once
in $a/z$ and determine a third order polynomial fit. As can be seen from
Figure 1, the first and second order fits still deviate visibly from the
analytical result obtained through full integration, but the third order fit
matches the result with high accuracy.

Figure 1: Relative height-dependent reduction in a solute’s diffusion parallel
to a surface. Polynomial fits up to third order to the analytical solution of
Eq. 6 describing the normalized self-diffusion coefficient of a sphere
parallel to a plane are shown. As the calculation of the analytical solution
is computationally demanding, and the often employed first order approximation
leads to unphysical negative values when a sphere is in close proximity to the
plane, a third order fit was implemented in this work.

The solute is now first assumed to diffuse with a short-time diffusion
coefficient obtained by the mean-field approach. In the presence of a surface,
this coefficient is further lowered by the relative reduction of Eq. 5
perpendicular to the surface, and the third order approximation of Eq. 6
parallel to the surface. The resulting short-time diffusion coefficient is
then used in the propagation step for BD simulations of Eq. 4. The resulting
diffusion coefficient for motion in the plane parallel to the surface is given
by

$D_{i}^{\parallel}(V_{i}^{\textnormal{frac}})=D_{i}(V_{i}^{\textnormal{frac}})\zeta_{\parallel}(z)\;$
(7)

and, equivalently, the diffusion coefficient for motion in the direction
perpendicular to the surface is obtained through multiplying by
$\zeta_{\perp}(z)$. The reduction of the rotational short-time diffusion
coefficient due to HI between a single solute and the surface is not included,
as it is much smaller than that for the translational short-time diffusion65
and is only apparent at very small surface-solute separations. The larger
reduction of the rotational diffusion of a solute as a result of its crowded
environment is, however, accounted for as described above.

The method of images is employed in hydrodynamic theory66 to satisfy the
boundary condition of zero flow field perpendicular to a surface by mirroring
the flow field of the solutes at the surface and, due to the linearity of the
Stokes equation, the superposition of the flow field yields the result. For a
non-zero flow field parallel to the surface, higher order terms arise. For a
discussion of the accuracies of the mobility matrices for a single sphere see,
for example, reference 60. However, for cases without collective motion
parallel to the surface, we restrict ourselves to the force monopole. Using
this assumption, we can extend the mean-field approach beyond the surface by
assuming an image flow field created by a mirrored configuration of solutes on
the other side of the surface. Figure 2 shows a snapshot of a simulation with
spheres diffusing in the presence of a surface. For each sphere, the local
occupied volume fraction is calculated within a cutoff sphere with
$R_{\textnormal{cut}}$, now also including the image solutes and thereby
accounting for the reflected flow field of the surrounding solutes. In
Mereghetti and Wade 51, $R_{\textnormal{cut}}=4a$, i.e., four times the radius
of the solute, was determined to be a good choice for crowded environments.
Note that the height of the atomically detailed surface at which the boundary
condition of zero flow in the direction perpendicular to the surface is set
has to be taken into account, and can be set as an input parameter in SDA.

Figure 2: Image solute method for computing HI for solutes diffusing in the
presence of a surface. The color indicates the local occupied volume fraction
$V_{\textnormal{frac}}$ within the cutoff sphere with
$R_{\textnormal{cut}}=4a$ with the solute with radius $a$ of interest placed
at the center. To resemble HEWL in the adsorption simulation, $a$ was set to
15 Å in this figure. The simulated volume is viewed as an orthographic
projection viewed along the plane of the surface. The HI of a solute with a
surface can be included by considering the interactions with image solutes
with mirrored positions on the other side of the surface. The occupied volume
fraction is then used with the hydrodynamic mean-field approach.

###  2.4 Grid-based Interactions

In SDA, the forces between a pair of macromolecules $1$ and $2$ are calculated
as finite-difference derivatives of the pairwise interaction, $\Delta G$,
defined as:

$\displaystyle\begin{aligned} \Delta
G&=\frac{1}{2}\sum_{i_{2}}\Phi_{\textnormal{el}_{1}}(\mathrm{\mathbf{r}}_{i_{2}})\cdot
q_{i_{2}}+\frac{1}{2}\sum_{i_{1}}\Phi_{\textnormal{el}_{2}}(\mathrm{\mathbf{r}}_{i_{1}})\cdot
q_{i_{1}}\\\
&+\sum_{i_{2}}\Phi_{\textnormal{ed}_{1}}(\mathrm{\mathbf{r}}_{i_{2}})\cdot
q_{i_{2}}^{2}\\\
&+\sum_{i_{1}}\Phi_{\textnormal{ed}_{2}}(\mathrm{\mathbf{r}}_{i_{1}})\cdot
q_{i_{1}}^{2}\\\
&+\sum_{n_{2}}\Phi_{\textnormal{np}_{1}}(\mathrm{\mathbf{r}}_{n_{2}})\cdot
SASA_{n_{2}}\\\
&+\sum_{n_{1}}\Phi_{\textnormal{np}_{2}}(\mathrm{\mathbf{r}}_{n_{1}})\cdot
SASA_{n_{1}}\\\
&+\sum_{n_{2}}\Phi_{\textnormal{rep}_{1}}(\mathrm{\mathbf{r}}_{n_{2}})+\sum_{n_{1}}\Phi_{\textnormal{rep}_{2}}(\mathrm{\mathbf{r}}_{n_{1}})\end{aligned}$
(8)

where $\mathrm{\mathbf{r}}$ refers to the atomic coordinates of either atoms
($n_{1,2}$) or charges ($i_{1,2}$). A detailed description and
parameterization can be found in references 67, 57, 50 and 68. The first two
terms in Eq. 8 represent the interaction energies of the charges ($q_{i_{2}}$
or $q_{j_{1}}$) of one solute with the electrostatic potential
($\Phi_{\textnormal{el}_{1}}$ or $\Phi_{\textnormal{el}_{2}}$) of another
solute 68.

To calculate the electrostatic interactions, the effective charge model 68
(ECM) is employed. These charges are calculated such that in a uniform
dielectric they reproduce the electrostatic potential in a shell around the
solute that was previously computed by solving the Poisson-Boltzmann (PB)
equation. Thereby, the effective charges implicitly account for
inhomogeneities in the dielectric that do not have to be considered further
during the BD simulation. Importantly, the required number of effective
charges necessary to maintain a high accuracy is commonly much smaller than
the number of partial atomic charges.

The third and fourth terms describe the electrostatic desolvation energy and
account for the effects of the low dielectric cavity of one solute on the
interaction energies of another 68. These terms are computed as the
interaction of the charges and the electrostatic desolvation potential of one
another ($\Phi_{\textnormal{ed}_{1}}$ or $\Phi_{\textnormal{ed}_{2}}$) 69
using the parameterization of reference 67. The fifth and sixth terms
correspond to the non-polar interactions due to the burial of the solvent
accessible surface areas (SASAs) of the surface atoms 67. The last two terms
of Eq. 8 describe the soft-core repulsive potential with an inverse power
function to avoid overlaps of solutes while at the same time preventing abrupt
force changes at close contact.

In general, the potentials could be recomputed at every step of a simulation
for the newly propagated configuration of atoms. In practise, we treat the
macromolecules as rigid bodies (each having a single conformation or an
ensemble of rigid conformations), allowing us to map all interaction
potentials onto grids centered on each of the macromolecules that are only
calculated once before the start of the BD simulation.

###  2.5 Debye-Hückel Model for Long-range Solute-Solute Electrostatic
Interactions

For short-range interactions, such as van der Waals and desolvation forces,
the interaction potentials decay quickly with distance, requiring only small
grids to obtain a negligible truncation error. In addition, efficient storage
algorithms, such as DT-Grid 70, 71, exist that further lower the memory
requirements for a given grid size. In contrast, for long-range electrostatic
interactions, algorithms such as DT-grid are not applicable and using
electrostatic potential grids sufficiently large to ensure the truncation
error is small can result in both very high memory requirements and slow
memory access times.

To alleviate the problem of the truncation of the electrostatic interaction
grid at the boundaries, the long-range interactions can be captured by the
Debye-Hückel (DH) approximation that was implemented in SDA by Mereghetti et
al. 52 and Martinez et al. 50. The method is well suited in combination with
the effective charge approximation, since beyond the region where
electrostatic interactions are treated on a grid-point charge basis, the
energy term, $\Delta G_{DH}$, can simply be added to the energies given in Eq.
8. According to the DH theory of dilute electrolyte solutions, all ions in the
solvent are treated as point charges while any pair of solutes $1$ and $2$ is
treated as spheres with radii $a_{1}$ and $a_{2}$ and net formal charges
$Q_{1}$ and $Q_{2}$. At a center-center separation of
$r=|\mathrm{\mathbf{r}}_{1}-\mathrm{\mathbf{r}}_{2}|$, the additional DH
energy term can be calculated as:

$\displaystyle\Delta G_{\textnormal{DH}}(r)=\begin{cases}\infty,&\text{if
}r<a\\\
\frac{Q_{1}Q_{2}}{4\pi\epsilon_{0}\epsilon_{\textnormal{r}}}\frac{\exp(-\kappa(r-a))}{r(1+\kappa
a)}&\text{if }r\geq a\end{cases}$ (9)

where ${a=a_{1}+a_{2}}$ is the sum of the radii of the two interaction solutes
and ${\kappa=(\lambda_{D})^{-1}}$ is defined as the inverse Debye length and
is calculated as:

$\kappa^{2}=\frac{2e^{2}}{k_{\textnormal{B}}T\epsilon_{0}\epsilon_{\textnormal{r}}}I$
(10)

where $\epsilon_{0}$ denotes the vacuum permittivity,
$\epsilon_{\textnormal{r}}$ the relative permittivity of the solvent, $e$ the
elementary charge, $k_{\textnormal{B}}$ the Boltzmann constant, $T$ the
temperature, and $I$ the salt ionic strength.

The full electrostatic potential grid is isotropically cut off to fulfill the
assumption of centrosymmetry at the grid edges and to switch to the analytical
DH potential beyond the cutoff distance. However, in the model described by
Mereghetti et al., when a solute is on the border of the electrostatic
potential of another solute, the grid-point charge interactions are calculated
using only the fraction of its effective charges that is inside the grid,
while the other charges are not considered. Furthermore, discontinuities in
the forces may arise once the solute moves fully outside the potential grid,
as the electrostatic calculation is switched to the DH treatment in a single
step. To improve the treatment of the transition between the full grid-based
treatment and the region with only monopole interactions, a hybrid approach
was introduced by Martinez et al. 50. This affected the region where only a
fraction of the effective charges of the first solute is inside the
electrostatic potential grid of the second solute, which we refer to as the
’partial’ region. This region is illustrated for HEWL in Figure 3. The
intersolute forces inside the partial region are computed from two
contributions 50: For the effective charges located inside the grid, the
partial interaction force $F_{\textnormal{partial}}(r)$ is calculated based on
the effective charges interacting with the potential grid. Outside the grid,
the DH treatment is applied but with the term computed for a solute charge
equal to the sum of the effective charges lying outside the spherical grid
boundary.

Here, we further modify the interaction model to improve the treatment of the
partial region. This improved model is aimed at ensuring force and torque
continuity, as well as consistency with the effective charge model. As before,
for the effective charges located inside the grid, the partial interaction
force $F_{\textnormal{partial}}(r)$ is calculated based on the effective
charges interacting with the potential grid. However, outside the grid, the
use of a single charge (defined as a fraction of the effective rather than the
formal net charge of the solute) for the DH calculation can lead to
inconsistencies, e.g., in the calculation of the torque. Therefore, here, the
effective charges are still employed for the calculation of the DH force
$F_{\textnormal{DH}}(r)$, as they represent the geometry of the solute, but a
scaling factor of total net formal charge / total effective charge is applied
to each charge. As such, the scaled value of the effective charges leads to
the correct formal net charge with the DH approximation accounting for the
dielectric medium.

Figure 3: Illustration of the’partial’ region of a solute electrostatically
interacting with another solute or surface. By way of example, HEWL is shown
in cartoon represntation with its molecular surface and with red and blue
spheres representing the effective charges computed at pH 7. The effective
charges in the ’partial’ region that overlap with the potential grid of
another solute or surface positioned to the left of HEWL are used to compute
the grid-point charge interaction forces.

This procedure gives a smooth transition from the ’inside’ region, in which
one solute is fully inside the grid of another, to the ’partial’ region, where
it is partially outside the grid. The transition to the outside region (when
the solute is fully outside the grid of the other solute) will, however, lead
to a discontinuity in the forces and potentials, as the interaction is being
reduced from a representation of many effective charges to a single charge in
one step. Furthermore, on approach of two solutes, the abrupt appearance of
effective charges may lead to artificial repulsive forces as, in the absence
of intersolute torques, the solutes will not have been able to adjust to a
favorable orientation with respect to each other. Therefore, we here introduce
a linear scaling from the partial force $F_{\textnormal{partial}}(r)$ to the
monopole interaction $F_{\textnormal{DH}}(r)$, within part of the partial
region,

$\displaystyle
F(r)=\lambda_{\textnormal{F}}\,F_{\textnormal{DH}}(r)+(1-\lambda_{\textnormal{F}})\,F_{\textnormal{partial}}(r)\;$
(11)

where the prefactor $\lambda_{\textnormal{F}}$ scales from zero at
${r=s_{1}-a_{2}}$, to one at $r=s_{1}$, i.e.
${\lambda_{\textnormal{F}}=(r-s_{1}+a_{2})/a_{2}}$ for ${s_{1}-a_{2}\leq r\leq
s_{1}}$, where $s_{1}$ denotes the half-length of one side of the cubic
electrostatic grid of solute 1. The split of the partial region into two parts
ensures that the electrostatic force present upon a solute entering or leaving
the grid of another is maintained while simultaneously avoiding
discontinuities in the forces and numerical instabilities.

###  2.6 Debye-Hückel Model for Long-Range Solute-Surface Electrostatic
Interactions

Here, we describe our extension of the DH model for use with homogeneously
charged surfaces. We assume an infinitely extended surface, which is locally
well fulfilled under experimental conditions and a given in simulations if
periodic boundary conditions are employed.

Placing the surface without loss of generality at $z=0$, the corresponding
potential at height $z$ is known from continuum electrostatics by solving the
Debye-Hückel equation:

$\frac{d^{2}}{dz^{2}}\Phi(z)=\kappa^{2}\Phi(z)\;$ (12)

As the system is assumed to be uniform in the x and y directions, the
potential only depends on the distance from the surface, $z$. Through an
exponential ansatz under the boundary condition that the potential
$\Phi(z)\rightarrow 0$ for $z\rightarrow\infty$, the solution

$\displaystyle\Phi(z)=\Phi_{0}\>e^{-\kappa z}$ (13)

is obtained, where $\Phi_{0}=\Phi(z=0)$. To determine $\Phi_{0}$, two
different cases can be considered, depending on the system to be modelled.

Firstly, we consider an infinitely thick surface with its charges located at
the solute-solution interface and aqueous solution on one side and an infinite
low dielectric region on the other. In this case, there is a surplus of
counterions in solution to balance the net charge of the surface, i.e., the
charge of the surface is equal and opposite to the charge of the solution. For
this case, Grahame 72 calculated :

$\displaystyle\sigma=-\int_{0}^{\infty}\rho_{\textnormal{e}}\>\mathrm{d}z$
(14)

where $\sigma$ denotes the surface charge density and $\rho_{\textnormal{e}}$
is the volume charge density of the solution, with the condition of

$\displaystyle\frac{\mathrm{d}\Phi}{\mathrm{d}z}(z\rightarrow\infty)=0\;$ (15)

He then used the differential form of Gauss’ law to calculate $\rho_{e}$ and
arrived at:

$\displaystyle\sigma$
$\displaystyle=\epsilon_{\textnormal{r}}\epsilon_{0}\kappa\Phi_{0}\>e^{-\kappa
z}|_{z=0}$ (16) $\displaystyle\Rightarrow\Phi_{0}$
$\displaystyle=\frac{\sigma}{\epsilon_{\textnormal{r}}\epsilon_{0}\kappa}$
(17)

This result is referred to as the Gouy-Chapman model.

Now, we again consider the solute to be an extended sphere of radius $a$ with
a low-dielectric region inside which the electrostatic potential does not
decay exponentially. Therefore, accounting for the sphere radius and
substituting Eq. 17 into Eq. 13, we obtain:

$\Phi=-\frac{\sigma}{\epsilon_{\textnormal{r}}\epsilon_{0}\kappa}\exp(-\kappa(z-a))\;$
(18)

Secondly, we consider a thin surface with solution on both sides, i.e., with
negligible thickness for long-range electrostatics. Here, it is assumed that
the excess counterions balancing the charge of the surface in solution are
equally distributed on both sides of the thin surface. Integrating Eq. 14 from
$-\infty$ to $\infty$, and again considering the solute as an extended sphere
gives:

$\Phi=-\frac{\sigma}{2\epsilon_{\textnormal{r}}\epsilon_{0}\kappa}\exp(-\kappa(z-a))\;.$
(19)

Here, $\Phi_{0}$ differs from Eq. 18 by a factor of one half. Whether a thick
or a thin surface is employed can be set in the input files of the SDA
simulation package.

##  3 Simulation Details

###  3.1 Approach Rate Constant for Two Proteins

First, we assessed the effects of truncating the precomputed electrostatic
potential grids used to model intersolute electrostatic interactions, and the
performance of the intersolute DH model for correcting for this truncation. We
simulated the diffusion of a pair of HEWL molecules and monitored the rate at
which they first approach each other to within separations $d$ sampled at 1 Å
intervals in the range from 30 to 80 Å. We compared the convergence of the
approach rate constant computed from these three sets of simulations, with
differing treatments of the long-range electrostatic interactions between
solutes that lie partially or fully outside the potential grids of their
interacting solutes, as the extent of the electrostatic potential grid was
increased. Namely, we first used a simple truncation of the potential at the
grid boundary with no DH correction at longer distances. Secondly, we used a
truncation of the grids and the DH correction that considers the full
effective charge distribution in the transition treatment at the grid
boundary. Thirdly, we used a truncation of the grids with the full transition
treatment that includes the scaling of Eq. 11.

For each grid extent and treatment of long-range electrostatics, 5,000,000 BD
trajectories were performed in which a pair of HEWL molecules were initially
separated by a center-to-center distance $b=600$ Å, with both molecules
rotated in different random orientations at the beginning of each trajectory.
The trajectories were simulated until the proteins reached a separation of
$c=900$ Å. The fraction of trajectories $\beta(d)$ in which the two molecules
diffused to a separation $d$ was monitored, allowing a rate constant for the
diffusional approach of the two molecules to a separation $d$, which we name
the approach rate constant $k_{\textnormal{app}}(d)$, to be calculated from:
73

$k_{\textnormal{app}}(d)=\frac{k(b)\beta(d)}{1-(1-\beta(d))\frac{k(b)}{k(c)}}$
(20)

where $k(b)$ and $k(c)$ are the rate constants for diffusion to the relative
separations $b$ and $c$ under the assumption that the interaction force
between the two solutes is centrosymmetric at separations of $b$ and larger.
If the force acting between two solutes at a distance $x$ can be assumed to be
negligible, $k(x)$ can be calculated using the Smoluchowski approximation of
$k(x)\approx 4\pi Dx$ where $D$ is the sum of the diffusion coefficients of
the two solutes. To check that the intersolute interaction energy is indeed
negligible at the assigned $b$ and $c$ distances, we also calculated the
results by determining $k(x)$ by numerical integration of

$k(x)=\left(\int_{x}^{\infty}\mathrm{d}r\left[\frac{\exp(U(r)/k_{\textnormal{B}}T)}{4\pi
r^{2}D}\right]\right)^{-1}\;$ (21)

We found that the interaction energy at ${b=600}$ Å and ${c=900}$ Å was
negligible.

The crystal structure of HEWL74 was taken from the RCSB Protein Data Bank (PDB
code: 1HEL). The protonation states of ionizable amino acid residues were
assigned with PDB2PQR75 at pH 7.0 using the atomic partial charges and radii
from the Amber force field, 76 giving a resulting net charge of +8 e on HEWL.
Cubic electrostatic potential grids of length 96, 128, 160, 224, 352 and 480 Å
and grid spacing 1 Å were created by solving the linearized Poisson-Boltzmann
(PB) equation using the finite difference multigrid method implemented in APBS
1.477 with an ionic strength of 5 mM, a solvent dielectric constant of 78.4,
and a protein interior dielectric of 4.0. The dielectric interface was
generated using a smoothed molecular surface and the potential at the boundary
grid points was assigned using a multiple DH sphere model. Effective charges
were calculated using the ECM module of SDA 768, 50 by placing charge sites on
the side chains of charged amino acid residues and the protein chain terminii.
Potential grids describing electrostatic desolvation interactions were created
using the make_edhdlj_grid module in SDA 750 following the parameterization
described by Gabdoulline et al.67

The BD simulations were performed with a modified version of SDA 7 using a
HEWL infinite dilution translational diffusion coefficient of 0.01232
Å${}^{2}/$ps and a rotational diffusion coefficient of $2.3057\times 10^{-5}$
rad${}^{2}/$ps. Both diffusion coefficients were calculated with HYDROPRO 78.
During the BD simulations, the intersolute repulsive interactions were
modelled using an excluded volume grid of spacing 1 Å generated for the HEWL
structure using a probe radius of 1.77 Å. Any attempted BD steps that resulted
in an atom of the other HEWL molecule entering this volume were rejected. A
HEWL radius of 15 Å, the radius of gyration calculated with HYDROPRO,78 was
used to describe the low dielectric cavity of the protein in the DH model.

###  3.2 Adsorption of HEWL to a Mica Surface

To assess the performance of the new solute-surface DH electrostatic model and
the solute-surface HI models, we simulated the adsorption of HEWL to a mica
surface at a range of HEWL concentrations in aqueous solutions of 5 and 50 mM
ionic strength. The mica surface was approximated using a homogeneously
charged graphite lattice surface with a charge density of $\sigma=-0.0215$
$e/$Å${}^{2}=-0.331$ C/m2, corresponding to that of a mica surface at pH 779.
This resulted in a charge of $-0.053$ _e_ on each graphite atom. The
electrostatic potentials of this surface when surrounded by aqueous solvent of
ionic strengths of 5 and 50 mM were calculated by solving the linearized PB
equation with APBS 1.4.77 To approximate the low dielectric interior of a
macroscopically-sized surface, 80 additional neutral graphite layers were
stacked at separations of 1 Å below the charged surface layer. All layers had
planar dimensions of 502.68 Å $\times$ 501.26 Å. For each ionic strength, the
linearized PB equation was solved using a coarse cubic grid of length 1056 Å
with 353 grid points in each dimension. The other calculation parameters were
as described for the HEWL calculation above. Following this, two finely-spaced
electrostatic potential grids of different sizes were computed using the
coarse grid to provide a focusing boundary condition. Both grids had spacings
of 1 Å and lengths of 352 Å in the dimensions in the plane of the surface. The
larger grid had a total length of 352 Å in the dimension perpendicular to the
surface and extended 336 Å above the surface. The smaller grid had a total
length of 128 Å and extended 114 Å above the surface. For simulations at 50
mM, additional HEWL electrostatic potential grids and effective charges
corresponding to this ionic strength were created following the procedure
described above. In all the simulations performed in the presence of a
surface, we used the HEWL electrostatic potential grid with sides of length 96
Å. We show below that, when combined with the intersolute DH correction, this
grid size is sufficient to accurately represent HEWL electrostatic
interactions.

In SDA, the electrostatic interactions between a pair of interacting species
are approximated using the effective charge model68 as the mean of the
interactions due to the set of effective charges on each interacting partner
with the PB-derived electrostatic potential on the other (the first two terms
in Eq. 8). These two terms should be approximately equal. To account for the
desolvation of solute-facing charges on each interacting partner as the two
species approach each other closely, additional correction terms are applied69
(the third and fourth terms in Eq. 8). This effective charge model is able to
approximate the PB-calculated electrostatic interaction forces between the
partners at each time step of the simulation, without requiring repeated
solution of the PB equation. However, it is more difficult to define
appropriate effective charges on an infinite surface by fitting the
electrostatic potential within a given shell around the surface. For a planar
surface, Romanowska et al.33 found that the effective charges calculated on
the surface had to be scaled so that the first two terms in Eq. 8 were similar
in magnitude for docked HEWL configurations in which the protein was close to
the surface. Here, to describe the adsorption process, we need to compute the
electrostatic interactions between HEWL molecules and the surface to good
accuracy at all heights above the surface. As the scaling used in Romanowska
et al. is not able to do this, we instead calculated the surface – HEWL
electrostatic interaction in one direction only, i.e. using the electrostatic
potential grids on the surface and the effective charges on the HEWL
molecules, effectively replacing the first two terms of Eq. 8 with a single
term with no $\frac{1}{2}$ prefactor. Therefore, effective charges were not
computed for the mica or silica surfaces modelled in this work.

The interaction grids for defining the electrostatic desolvation potential of
the surface were calculated with make_edhdlj_grid from SDA 750 following the
parameterization described by Gabdoulline et al.67. The electrostatic
desolvation grid of HEWL described in the previous section was used. The
additional grids describing the non-polar desolvation potentials of the
surface and HEWL were calculated following the parameterization used in the
graphite – hydrophobin simulations described by Mereghetti and Wade.46 Unlike
the two solute simulations described in the previous section, repulsive
interactions were modelled using a soft-core repulsive model, again following
the parameterization used for the graphite – hydrophobin simulations46. The
parameterization of all interaction terms used in these simulations were
chosen to be consistent with those used in previous simulations of HEWL
solutions.52

BD simulations, each of a duration of 12 $\mu$s, were performed with a time
step of 0.4 ps at HEWL concentrations of 5, 10, 20 and 30 mg/ml. Note that, as
the experiments were performed on timescales from minutes to hours and at
concentrations in the range of ${\mu\mathrm{g}/\mathrm{ml}}$80, 81, it is not
feasible to simulate these systems in atomic detail at the experimental time
and length scales. Therefore, we used higher concentrations than in the
experiments, that were chosen with the aim of reproducing the qualitative
trends observed in experiments and providing insights into the mechanisms of
the adsorption processes and the contributions of the different components of
the interaction forces.

For each HEWL concentration, an initial configuration of 150 HEWL molecules
was placed in a box with periodic boundaries and length 322.8 Å in both
dimensions parallel to the surface. These dimensions were chosen to be
significantly smaller than the surface used for generating the PB
electrostatic grid in order to minimize the curvature of the isoelectric field
contours, so that the electrostatic potential calculated in the non-periodic
PB calculation provided a good approximation of that of an infinite charged
plane in the simulated volume. The extents of the simulated volumes above the
surface were chosen to give the desired HEWL concentrations (Table
LABEL:tab:adsorption_box_sizes). The region within 120 Å, i.e., four times the
diameter of HEWL, above the surface was left empty to avoid an immediate
collapse of a number of proteins onto the surface that would not occur at
experimental conditions at lower concentrations.

A protein was considered to be adsorbed if its geometric center was located
less than 50 Å above the surface. This criterion was chosen because,
considering the size and the ellipsoidal shape of HEWL, it captures two
adsorption layers (see Figure 7 and corresponding discussion). We tested that
the results were robust with respect to the choice of this criterion, and that
the proteins positioned below this threshold were stably associated with the
surface and not diffusing freely.

At each HEWL concentration, four simulations were performed for a solution of
5 mM salt concentration. The first three, in which HI were ignored, varied in
the treatment of long-ranged electrostatic interactions: (1) the larger
electrostatic potential grid, which extended 336 Å above the surface, with the
charged surface DH correction above this and intersolute DH correction for
solute – solute interactions that extend beyond grid boundaries; (2) the
smaller electrostatic potential grid that extended 114 Å above the surface and
both DH corrections; and (3) the smaller potential grid with no corrections.
In the fourth simulation, HI were accounted for and the smaller electrostatic
potential grid was used with both DH corrections.

Two further simulations were performed at each concentration assuming a
solution of 50 mM salt concentration, one neglecting and one accounting for
HI. In both of these simulations, the smaller electrostatic potential grids
and both DH corrections were used. As we wish to model a macroscopically-sized
charged surface, we represented it as an infinitely thick solvent-excluding
surface, i.e. using Eq. 17. The HEWL radius used for calculating HI was 15 Å
which we also used for the solvent excluded region in the DH models. All other
simulation parameters were as described for the two-HEWL simulations.

Lastly, as a comparison, another set of simulations with a silica surface was
conducted at the HEWL concentrations mentioned above. Again, silica was
approximated using a homogeneously charged graphite lattice surface, but with
a charge density of ${\sigma=-0.0013}$ e$/$
Å2=${-0.02~{}\mathrm{C}/\mathrm{m}^{2}}$, corresponding to that of a silica
surface at pH 780. This assignment resulted in a charge of $-$0.0032 e on each
graphite atom. Here, the simulations were only conducted at 5 mM salt
concentration using an electrostatic potential grid that extended 114 Å above
the surface and was solved as described above for the mica surface. The DH
correction and HI were included. All other interactions were calculated as for
the mica surface.

##  4 Results

###  4.1 Convergence of the approach rate constant for two HEWL proteins with
increasing electrostatic potential grid size

The rate constants for the approach of two HEWL proteins to separations in the
range from 30 to 80 Å were calculated using three different models for
electrostatic interactions, and the convergence of these models with
increasing electrostatic potential grid size was determined (Figure 4). Using
the first model, a simple truncation of the electrostatic interactions when
the effective effective charges of one protein lie outside the electrostatic
potential grid of the other, large truncation errors were seen when using the
smaller grids (Figure 4a). Only when using the second largest grid, which
included the potential within a radius 176 Å from the protein center, did the
rate constant curve converge to that of the larger grid with a 240 Å radius.
When using smaller grids, considerably higher rate constants were computed due
to the lack of electrostatic repulsion between the HEWL proteins, which have
identical net charges, at longer distances.

Figure 4: Approach rate constant for two HEWL protein solutes. Simulations
were performed with six sizes of the electrostatic potential grids centered on
each of the two proteins as indicated by the colors and the corresponding
radii in Å. (a) No long-range electrostatic interactions beyond the grid
boundaries were considered. (b) Using the DH approximation beyond the grid
radius with the transition treatment as in Martinez et al. 50, but modified so
that the charge $Q$ for each solute in Eq. 9 is defined by the formal charge
of the protein, rather than the sum of its effective charges. (c) Using the
full improved transition treatment with a smooth scaling switch (Eq. (11))
between the grid-based and Debye-Hückel regimes.

While the simulations performed with a simple truncation of the electrostatic
interactions showed a clear convergence of the approach rate constant for all
separations as the grid size increased (Figure 4a), this was not the case when
the transition treatment described in Martinez et al. 50, with the corrected
assignment of solute net charges $Q$ but not the scaling of Eq. (11), was
used, although the errors in these simulations were generally of smaller
magnitude than obtained with simple truncation (Figure 4b). The errors in
these simulations had the opposite effect to those obtained with the
truncation model, leading to a decrease in the approach rate constant at all
separations. We again observed that only the rate constants from simulations
using the second largest grid (with a 176 Å radius) agreed with those obtained
using the largest grid at all separations. The smallest grid produced the
largest error in the predicted rate constants at smaller separations (Figure
4b), whereas the third largest grid (with a 112 Åradius) produced the largest
error in the predicted rate constants at larger separations. The reason for
these errors is apparent when examining the rate constants predicted for the
simulations performed using the two smallest grid sizes. At larger
separations, the rate constants match well with those of the simulations
performed with the 240 Å radius grids, showing that, when corrected to use the
formal charge of the proteins in place of the sum of their effective charges,
the Debye-Hückel model accurately models electrostatic interactions at these
longer distances. It is only when the effective charges on the two HEWL
proteins begin to enter into the electrostatic potential grids of the other
protein, at separations slightly larger than the grid radius, that the
predicted rate constants diverged. The abrupt entrance of the effective
charges into the electrostatic grid of the interacting protein means that the
proteins did not have the opportunity to rearrange into more favorable
orientations, leading to unphysical repulsions in this region, and resulting
in a reduction in the calculated rate constants.

The simulations run with the smoothed transition between the grid-based and
Debye-Hückel regimes showed much reduced errors for all grid sizes at all
separations (Figure 4c). Even with the smallest 48 Å radius grids, the
divergence in the predicted rate constants from those obtained from the
simulations with the largest grid size was minimal. For this reason, the 48 Å
radius grid was used in all subsequent simulations, as it can accurately
describe the electrostatic interactions of HEWL when coupled with the new
Debye-Hückel model.

###  4.2 Adsorption of Multiple HEWL Proteins to a Mica Surface

Figure 5: Simulation of HEWL adsorption to a mica surface. Protein adsorption
in simulations with different initial bulk protein concentrations (5 - 30
mg$/$ml) are shown by color and are compared for two different salt
concentrations (5 (a-d) and 50 (e-f) mM) and different treatments of
electrostatic interactions and HI. (a) - (c) show the results from different
treatments of electrostatic interactions, comparing the use of a large
electrostatic potential grid (336 Å height above the surface)(a) to a smaller
one (114 Å) with (b) and without an additional long-range DH treatment (c).
The effect of including solute-solute and solute-surface HI is shown in (d)
and (f).

In all the simulations of systems consisting of multiple HEWL proteins and a
mica surface, which were performed with several simulation conditions and
interaction treatments, we observed that both the number of adsorbed proteins
and the corresponding surface coverage in ${\mathrm{mg}/\mathrm{m^{2}}}$
increased with simulation time, with the adsorption rate being faster
initially and then gradually levelling off, see Figure 5. However, the number
of adsorbed proteins was highly dependent on the initial bulk protein
concentration.

The simulations at 5 mM salt concentration with the electrostatic interactions
calculated using a grid with a height of 114 Å above the surface and the
interactions beyond the grid captured with the DH model (Figure 5b) yielded
similar surface coverage values at all HEWL concentrations to those obtained
when using a larger electrostatic potential grid extending 336 Å above the
surface together with the DH approximation beyond the grid (Figure 5a). When
the DH treatment to correct for the truncation error was not used, the
resulting surface coverage was greatly reduced (Figure 5c).

At the higher salt concentration of 50 mM (Figure 5e) and otherwise the same
simulation conditions as in Figure 5b, the results remain similar for low HEWL
concentrations but differ for higher concentrations. The number of adsorbed
proteins after 12 $\mathrm{\mu}$s of the simulations at 30 mg$/$ml protein
concentration is almost half that at 5 mM. In contrast, at the lowest protein
concentration of 5 mg$/$ml, the number of proteins adsorbed at the two salt
concentrations is very similar.

Figure 6: Short-time diffusion coefficients of the HEWL proteins computed from
BD simulations with different HI contributions included. Each dot represents a
protein. For the simulation conditions of Figure 5d, the short-time
coefficients of the individual proteins are plotted against the height of
their center above the surface. The red dots indicate the results when HI
between the proteins only were considered, green when, in addition, the image
flow field of the surrounding proteins was considered, and blue when also
accounting for the HI of the protein with its own reflected flow field.

Considering HI between the proteins, as well as between the proteins and the
surface (Figure 5d and f, at 5 mM and 50 mM salt concentration, respectively),
the results are similar to those from corresponding simulations with HI
neglected (Figure 5b and e, respectively).

####  4.2.1 Adsorption Patterns

Examination of the distribution of the proteins revealed that two protein
layers can form on the mica surface during the simulations. This distribution
is shown in Figure 7 at the end of the simulation with a 30 mg$/$ml protein
concentration and the simulation conditions of Figure 5d.

We next evaluated the short-time diffusion coefficients of the proteins during
the assembly of the first adsorption layer. For the simulations at a protein
concentration of 30 mg$/$ml (Figure 5d), scatter plots of the diffusion
coefficients parallel and perpendicular to the surface are shown in Figure 6a
and b, respectively, against the height of the protein’s center above the
surface, with one point representing one solute. When using the mean-field
approach to account for the HI between the solutes only, the short-time
diffusion coefficient increases as the proteins approach close to the surface
(red dots). The reason is that this method considers the part of the space
occupied by the surface as empty. Extending the mean-field approach by taking
into account image solutes beyond the surface, also accounts for the flow
field reflection of the surrounding solutes, and resulted in the diffusion
coefficients of the proteins being roughly constant at all distances from the
surface, as shown by the green dots. However, when the HI computation included
the direct coupling effects of the solutes with their own reflected flow
field, there was a considerable decrease of the short-time diffusion
coefficents close to the surface, as shown by the blue dots.

Figure 7: Adsorption patterns of HEWL on the mica surface at the end of a BD
simulation (i.e. after 12 $\rm{\mu}s$), at 30 mg$/$ml HEWL concentration. The
circles indicate the positions of the centers of geometry of the proteins and
are viewed (a) from the side and (b) from the top. The circles are smaller
than the size of HEWL to faciliate visualization. In (a), the dashed line
indicates the projection of a sphere with the radius of gyration of HEWL. It
is distorted, as the distribution of the proteins in the $x$ and $z$
directions is shown on different scales.

###  4.3 Adsorption to a Silica Surface

We next conducted simulations with a silica surface at 5 mM salt concentration
using an electrostatic potential grid extending 114 Å above the surface and
the DH approximation beyond. The full HI model was used for the simulations.
In comparison to mica, the number of adsorbed proteins at all four HEWL
concentrations simulated was lower for the silica surface (Figure 8), due to
the lower surface charge of silica compared to mica. Furthermore, in
comparison to mica, the adsorption curves levelled off much earlier,
especially at the higher protein concentrations. The number of proteins in the
plateau region of the simulations at 20 and 30 mg$/$ml bulk protein
concentrations is more similar, although it was slightly higher at the higher
concentration.

Figure 8: Simulation of HEWL adsorption to a silica surface at four bulk
concentrations of HEWL (5, 10, 20 and 30 mg/ml). Apart from replacing the mica
surface with a silica surface, the simulation conditions and the treatment of
interactions, as well as the color scheme, are the same as for Figure 5d.

##  5 Discussion

We first consider our results with respect to the newly introduced methods,
and then consider the insights into the systems studied that can be obtained
from the simulations.

We computed the approach rate constants characterizing the rate at which two
HEWL molecules reach a certain distance from one another by diffusion. The
introduction of the new treatment of the transition of the electrostatic
interactions between the grid based full charge representation and a one-
dimensional DH approximation leads to accurate results for three dimensional
interaction grids that are much smaller — by almost a factor five in the case
of HEWL — than previously required for the same accuracy.

The reproduction of the adsorption of HEWL proteins to a mica surface required
the simulation of a high number of molecules — here, 150 HEWL proteins — over
several microseconds to cover the full adsorption process. The DH
approximation for homogeneously charged surfaces, together with the improved
treatment of the solute electrostatic grid transition, was used to capture the
long-range electrostatic interactions. The agreement between the results
obtained by using grids of 336 and 114 Å height (Figure 5a and 5b,
respectively), indicate that also in this case, a full charge representation
is only necessary at close distances when using the improved DH approximation
for long-range interactions. Furthermore, even though the electrostatic
interaction decays exponentially, the much smaller numbers of adsorbed
proteins observed when the interactions beyond the grid are neglected (Figure
5c) shows that the long-range electrostatic forces still critically influence
the long-time diffusion behavior.

From a computational perspective, the reduction in memory and storage
requirements achieved by using the improved long-range DH tratement becomes
important when simulating the diffusion of large macromolecules or when
carrying out simulations with large numbers of different solute molecules,
each of which has its own precomputed electrostatic potential grid.

We furthermore extended the computationally efficient mean-field treatment of
HI51 for many molecules to include confining surfaces. Both the relative
anisotropic reduction in diffusion due to solute-surface HI, as well as the
reduction based on the HI with the reflected flow-field of the surrounding
solutes are included, in addition to the direct effect of HI between solutes.
It is important to note that, due to the complexity of HI, these effects are
combined in a first-order approximation, thereby neglecting higher order
terms.

For HEWL adsorption, although they largely decreased the mobility of the
solutes in the vicinity of the surface (Figure 6), the HI showed little to no
effect on the adsorption curves at 5 mM (Figure 5d) and 50 mM (Figure 5f) salt
concentration. Similarly to other studies, the results show that while HI
predictably change the local properties of a system, it is generally hard to
predict the effect of HI on averaged global observables a priori. Antosiewicz
et al. 83 analyzed the kinetics of diffusional encounters of, among others,
Barnase and Barstar, and found that, even though the magnitudes of the torques
resulting from the hydrodynamic coupling of the associating molecules were
comparable with the magnitudes of the torques from electrostatic interactions,
the overall effects of the hydrodynamic torques on the association kinetics
were rather small. For studies of solutes in the presence of a surface, there
is agreement that the diffusivity of the solutes is reduced, although there is
a strong dependence on the properties and assumptions of the individual
systems studied as regards whether this reduction also influences the
adsorption properties. For irreversible adsorption processes of spheres to an
attractive surface, Pagonabarraga and Rubí 84 found that although HI
influences local properties such as the pair correlation, the macroscopic
quantities were largely unaffected. In contrast, studies on similar models
revealed that the time dependence of the surface coverage near saturation85
and the saturation level86 were highly influenced by the HI between the free
solutes and the adsorbing surface.

Importantly, these models often only consider adsorption as a purely diffusive
process, e.g., as in the studies based on the Random Sequential Adsorption
models 87. Therefore, changes in diffusivity influence both the approach to
the surface and the task of finding a free spot on the surface. While some
models include an attractive force to the surface 84, 85, they generally do
not consider interactions between the solutes beyond exclusion and HI. As
these neglected interactions are often strong for protein interactions, the
results for proteins can deviate highly from the predictions of these models.

The simulation results can be compared to a number of experimental studies 80,
81, 88, 89. However, it should be recalled that, firstly, we simulate bulk
concentrations of several ${\mathrm{mg}/\mathrm{ml}}$ compared to
${\mathrm{\mu g}/\mathrm{ml}}$ in the experimental studies. Secondly, the
results often vary between different experimental techniques88, 90, and even
between different setups for the same technique 91, 89, 90. Therefore, the
simulations described here are only intended to reproduce qualitative trends,
but can thereby provide insights into the adsorption processes.

We observe that the simulated adsorption kinetics strongly depend on the
protein concentration, which is in agreement with experimental studies 80, 81,
88, 89. Furthermore, the adsorption kinetics are dependent on the salt
concentration, as, due to the change in the decay of the electrostatic
potentials, the repulsive forces between already adsorbed proteins and newly
approaching ones becomes more dominant in comparison to the attractive forces
with the surface. This observation is again in agreement with experiments80.

The conditions of the adsorption simulation to a mica surface match
experimental studies conducted at 5 mM ionic strength by Daly et al. 80 using,
among other techniques, streaming current measurements to investigate the
saturated surface coverage, based on the adsorption mechanisms proposed by
Robeson and Tilton 81 under the same conditions. The surface coverage obtained
of 2.3 ${\mathrm{mg}/\mathrm{m^{2}}}$ is slightly higher than in the
experimental studies (2.0 ${\mathrm{mg}/\mathrm{m^{2}}}$). However, the
reported amount was observed in the first layer on the surface. If we only
consider the first layer (adsorption curves not shown), we obtain around the
same value of 2.0 ${\mathrm{mg}/\mathrm{m^{2}}}$. Especially for the mica
surface, however, it is not possible to identify whether the adsorption curves
for different concentrations will converge to the same level, as the levelling
off is very slow at the lower concentrations. Simulating even higher
concentrations would, however, further increase the artificial effect of
repulsion between the proteins in the bulk that pushes them towards the
surface, which does not occur to such a degree in experiments at more dilute
conditions.

The qualitative features of the adsorption kinetics to the silica surface
match those for the mica surface. The adsorption curves for the two highest
concentrations appear converged at the end of the simulations. However, it is
unclear if the small remaining difference in the number of adsorbed proteins
between the two plateaus is due to additional repulsion between the proteins
in bulk, or if it would vanish after a much longer simulation time. Two
experimental studies using neutron reflection88 and dual polarization
interferometry89 reported surface coverages of 3.29, 2.10 and 0.74
${\mathrm{mg}/\mathrm{m^{2}}}$ for protein concentrations of 4, 1 and 0.03
${\mathrm{mg}/\mathrm{ml}}$, respectively, so our result is slightly higher
than that for the lower protein concentration. As the final surface coverage
is much lower than that obtained by just considering the first layer at the
mica surface, this shows that the adsorption to the silica surface is not
limited by available space on the surface, but rather regulated by the
electrostatic interactions between the proteins and the surface. In contrast,
for the higher protein concentrations, the neutron reflection study reports
several adsorption layers, which we, however, only observe for the more highly
charged mica surface.

Two additional approximations underlying the BD simulations presented here
should be emphasized again. Firstly, HEWL is treated as a rigid body. While
HEWL is not known to show large conformational transitions upon adsorbing to a
surface92, flexible tail and side-chain movements may facilitate binding to
the surface93, and hence, may result in a denser packing of the proteins on
the surface. For a small number of flexible proteins, the BD simulations could
be combined with a MD simulation technique to account for induced fit in the
last step of the adsorption, as shown in Ref. 32. Secondly, in the BD
simulations of systems of many solutes and a surface, the forces acting on one
solute are obtained by summing the contributions from the interactions the
surrounding solutes and the surface in an additive manner. It is thereby
assumed that the linearized PB equation holds, which may not be fully
justified for highly charged species, e.g., for the mica surface without any
oppositely charged proteins adsorbed to it. With these assumptions, the
combined methodological advances described here will enable the
computationally efficient study on other adsorption processes or of systems
with confined crowded environments that are hard to access in vivo
experimentally. A BD simulation of 12 $\mu$s with over 150 HEWL molecules
treated in atomic detail for one bulk protein concentration took less than a
week on a compute node with 16 CPU cores.

##  6 Conclusions

We have here presented new computational models to efficiently capture
hydrodynamic and long-range electrostatic interactions for BD simulations of
proteins in the presence of a surface. These models have been implemented in
the SDA software package. Proteins and surfaces are modelled at full atomic
detail while approximating them by one or more rigid bodies. Considering the
self-association of HEWL and the adsorption process of many HEWL molecules to
mica and silica surfaces, long-range electrostatic interactions were shown to
critically influence the diffusional kinetics at low ionic strength, but, with
the improved methodology described here, only had to be calculated in full
detail at close distances. In contrast, while HI lowered the mobility of the
proteins close to the surface, they had almost no effect on the observed
protein adsorption kinetics. The simulations were able to reproduce trends in
protein-surface adsorption properties observed in different experimental
conditions. In conclusion, the methodology presented here enables the
simulation of over a hundred proteins in the presence of a surface modelled in
atomic detail at very low computational cost compared to, for example,
atomistic explicit solvent MD simulations.

##  7 Code and Data Availability

The SDA software package is available at http://mcm.h-its.org/sda7 and the new
methodology described here is available in SDA version 7.3, which also
includes example files for the systems simulated here.

This research has received funding from the European Union Seventh Framework
Programme (FP7/2007-2013) under grant agreement number 604102 (HBP Ramp-Up
Phase), the European Union Horizon 2020 Framework Programme for Research and
Innovation under grant agreement numbers 720270, 785907 and 945539 (SGA1, SGA2
and SGA3), and the Klaus Tschira Foundation. The authors thank Dr. Stefan
Richter for software support, Abraham Muniz-Chicharro for testing of the
software, and Dr. Julia Romanowska for initial studies on the protein-surface
adsorption systems.

## References

  * Mika and Poolman 2011 Mika, J. T.; Poolman, B. Macromolecule diffusion and confinement in prokaryotic cells. _Current Opinion in Biotechnology_ 2011, _22_ , 117–126
  * Mourão et al. 2014 Mourão, M. A.; Hakim, J. B.; Schnell, S. Connecting the dots: The effects of macromolecular crowding on cell physiology. _Biophysical Journal_ 2014, _107_ , 2761–2766
  * Luby-Phelps 2013 Luby-Phelps, K. The physical chemistry of cytoplasm and its influence on cell function: An update. _Molecular Biology of the Cell_ 2013, _24_ , 2593–2596
  * Minton 2001 Minton, A. P. The Influence of Macromolecular Crowding and Macromolecular Confinement on Biochemical Reactions in Physiological Media. _Journal of Biological Chemistry_ 2001, _276_ , 10577–10580
  * Rivas et al. 2004 Rivas, G.; Ferrone, F.; Herzfeld, J. Life in a crowded world. _EMBO reports_ 2004, _5_ , 23–27
  * Minton 2006 Minton, A. P. How can biochemical reactions within cells differ from those in test tubes? _Journal of Cell Science_ 2006, _119_ , 2863–2869
  * Klumpp et al. 2013 Klumpp, S.; Scott, M.; Pedersen, S.; Hwa, T. Molecular crowding limits translation and cell growth. _Proceedings of the National Academy of Sciences of the United States of America_ 2013, _110_ , 16754–16759
  * Kuznetsova et al. 2015 Kuznetsova, I. M.; Zaslavsky, B. Y.; Breydo, L.; Turoverov, K. K.; Uversky, V. N. Beyond the excluded volume effects: Mechanistic complexity of the crowded milieu. _Molecules_ 2015, _20_ , 1377–1409
  * Rivas and Minton 2016 Rivas, G.; Minton, A. P. Macromolecular Crowding In Vitro, In Vivo, and In Between. _Trends in Biochemical Sciences_ 2016, _41_ , 970–981
  * Zhou 2013 Zhou, H.-X. Influence of crowded cellular environments on protein folding, binding, and oligomerization: Biological consequences and potentials of atomistic modeling. _FEBS Letters_ 2013, _587_ , 1053–1061
  * Guin and Gruebele 2019 Guin, D.; Gruebele, M. Weak Chemical Interactions That Drive Protein Evolution: Crowding, Sticking, and Quinary Structure in Folding and Function. _Chemical Reviews_ 2019, _119_ , 10691–10717
  * Vieira et al. 2009 Vieira, E. P.; Rocha, S.; Carmo Pereira, M.; Möhwald, H.; Coelho, M. A. N. Adsorption and Diffusion of Plasma Proteins on Hydrophilic and Hydrophobic Surfaces: Effect of Trifluoroethanol on Protein Structure. _Langmuir_ 2009, _25_ , 9879–9886
  * Wang et al. 2002 Wang, L.; Schönhoff, M.; Möhwald, H. Lipids Coupled to Polyelectrolyte Multilayers: Ultraslow Diffusion and the Dynamics of Electrostatic Interactions. _The Journal of Physical Chemistry B_ 2002, _106_ , 9135–9142
  * Gray 2004 Gray, J. J. The interaction of proteins with solid surfaces. _Current Opinion in Structural Biology_ 2004, _14_ , 110–115
  * Minton 1992 Minton, A. Confinement as a determinant of macromolecular structure and reactivity. _Biophysical Journal_ 1992, _63_ , 1090–1100
  * Minton 1995 Minton, A. Confinement as a determinant of macromolecular structure and reactivity. II. Effects of weakly attractive interactions between confined macrosolutes and confining structures. _Biophysical Journal_ 1995, _68_ , 1311–1322
  * Zhou et al. 2008 Zhou, H.-X.; Rivas, G.; Minton, A. P. Macromolecular crowding and confinemment: biochemical, biophysical, and potential physiological consequences. _Annu Rev Biophys_ 2008, _37_ , 375–397
  * Hoffmann et al. 2009 Hoffmann, E. K.; Lambert, I. H.; Pedersen, S. F. Physiology of Cell Volume Regulation in Vertebrates. _Physiological Reviews_ 2009, _89_ , 193–277
  * Konorty et al. 2009 Konorty, M.; Brumfeld, V.; Vermeglio, A.; Kahana, N.; Medalia, O.; Minsky, A. Photosynthetic System in Blastochloris viridis Revisited. _Biochemistry_ 2009, _48_ , 4753–4761
  * Dix and Verkman 2008 Dix, J. A.; Verkman, A. Crowding Effects on Diffusion in Solutions and Cells. _Annual Review of Biophysics_ 2008, _37_ , 247–263
  * Cohavi et al. 2010 Cohavi, O.; Corni, S.; De Rienzo, F.; Di Felice, R.; Gottschalk, K. E.; Hoefling, M.; Kokh, D.; Molinari, E.; Schreiber, G.; Vaskevich, A.; Wade, R. C. Protein-surface interactions: challenging experiments and computations. _Journal of Molecular Recognition_ 2010, _23_ , 259–262
  * Ozboyaci et al. 2016 Ozboyaci, M.; Kokh, D. B.; Corni, S.; Wade, R. C. Modeling and simulation of protein-surface interactions: achievements and challenges. _Quarterly Reviews of Biophysics_ 2016, _49_
  * Ando and Skolnick 2013 Ando, T.; Skolnick, J. Importance of Excluded Volume and Hydrodynamic Interactions on Macromolecular Diffusion In Vivo. _Quantum Bio-Informatics V_ 2013, _30_ , 375–387
  * Dufresne et al. 2000 Dufresne, E. R.; Squires, T. M.; Brenner, M. P.; Grier, D. G. Hydrodynamic Coupling of Two Brownian Spheres to a Planar Surface. _Physical Review Letters_ 2000, _85_ , 3317–3320
  * Doster and Longeville 2007 Doster, W.; Longeville, S. Microscopic diffusion and hydrodynamic interactions of hemoglobin in red blood cells. _Biophysical Journal_ 2007, _93_ , 1360–1368
  * Lisicki and Nägele 2016 Lisicki, M.; Nägele, G. In _Lecture Notes in Physics_ ; Lang, P. R., Liu, Y., Eds.; Lecture Notes in Physics; Springer International Publishing: Cham, 2016; Vol. 917; pp 313–386
  * Czajka et al. 2019 Czajka, P.; Antosiewicz, J. M.; Długosz, M. Effects of Hydrodynamic Interactions on the Near-Surface Diffusion of Spheroidal Molecules. _ACS Omega_ 2019, _4_ , 17016–17030
  * Skaug et al. 2013 Skaug, M. J.; Mabry, J.; Schwartz, D. K. Intermittent Molecular Hopping at the Solid-Liquid Interface. _Physical Review Letters_ 2013, _110_ , 256101
  * Matse et al. 2017 Matse, M.; Chubynsky, M. V.; Bechhoefer, J. Test of the diffusing-diffusivity mechanism using near-wall colloidal dynamics. _Physical Review E_ 2017, _96_ , 1–12
  * Antosiewicz and Długosz 2020 Antosiewicz, J. M.; Długosz, M. Constant-pH Brownian Dynamics Simulations of a Protein near a Charged Surface. _ACS Omega_ 2020, _5_ , 30282–30298
  * Buijs and Hlady 1997 Buijs, J.; Hlady, V. Adsorption Kinetics, Conformation, and Mobility of the Growth Hormone and Lysozyme on Solid Surfaces, Studied with TIRF. _Journal of Colloid and Interface Science_ 1997, _190_ , 171–181
  * Ozboyaci et al. 2016 Ozboyaci, M.; Kokh, D. B.; Wade, R. C. Three steps to gold: mechanism of protein adsorption revealed by Brownian and molecular dynamics simulations. _Physical Chemistry Chemical Physics_ 2016, _18_ , 10191–10200
  * Romanowska et al. 2015 Romanowska, J.; Kokh, D. B.; Wade, R. C. When the Label Matters: Adsorption of Labeled and Unlabeled Proteins on Charged Surfaces. _Nano Letters_ 2015, _15_ , 7508–7513
  * Li and Liu 2013 Li, C.; Liu, M. Protein dynamics in living cells studied by in-cell NMR spectroscopy. _FEBS Letters_ 2013, _587_ , 1008–1011
  * Ignatova and Gierasch 2004 Ignatova, Z.; Gierasch, L. M. Monitoring protein stability and aggregation in vivo by real-time fluorescent labeling. _Proceedings of the National Academy of Sciences_ 2004, _101_ , 523–528
  * Yu et al. 2016 Yu, I.; Mori, T.; Ando, T.; Harada, R.; Jung, J.; Sugita, Y.; Feig, M. Biomolecular interactions modulate macromolecular structure and dynamics in atomistic model of a bacterial cytoplasm. _eLife_ 2016, _5_ , 1–22
  * Feig et al. 2017 Feig, M.; Yu, I.; Wang, P. H.; Nawrocki, G.; Sugita, Y. Crowding in Cellular Environments at an Atomistic Level from Computer Simulations. _Journal of Physical Chemistry B_ 2017, _121_ , 8009–8025
  * von Bülow et al. 2019 von Bülow, S.; Siggel, M.; Linke, M.; Hummer, G. Dynamic cluster formation determines viscosity and diffusion in dense protein solutions. _Proceedings of the National Academy of Sciences of the United States of America_ 2019, _116_ , 9843–9852
  * Wei et al. 2012 Wei, T.; Carignano, M. A.; Szleifer, I. Molecular dynamics simulation of lysozyme adsorption/desorption on hydrophobic surfaces. _Journal of Physical Chemistry B_ 2012, _116_ , 10189–10194
  * Kubiak et al. 2016 Kubiak, K.; Mulheran, P. A.; Kubiak, K.; Mulheran, P. A. Molecular Dynamics Simulations of Hen Egg White Lysozyme Adsorption at a Charged Solid. 2016, 12189–12200
  * Brancolini et al. 2012 Brancolini, G.; Kokh, D. B.; Calzolai, L.; Wade, R. C.; Corni, S. Docking of Ubiquitin to Gold Nanoparticles. _ACS Nano_ 2012, _6_ , 9863–9878
  * Brancolini et al. 2015 Brancolini, G.; Corazza, A.; Vuano, M.; Fogolari, F.; Mimmi, M. C.; Bellotti, V.; Stoppini, M.; Corni, S.; Esposito, G. Probing the Influence of Citrate-Capped Gold Nanoparticles on an Amyloidogenic Protein. _ACS Nano_ 2015, _9_ , 2600–2613
  * Ravichandran and Talbot 2000 Ravichandran, S.; Talbot, J. Mobility of adsorbed proteins: a Brownian dynamics study. _Biophysical Journal_ 2000, _78_ , 110–20
  * Carlsson et al. 2004 Carlsson, F.; Hyltner, E.; Arnebrant, T.; Malmsten, M.; Linse, P. Lysozyme Adsorption to Charged Surfaces. A Monte Carlo Study. _J. Phys. Chem. B_ 2004, _108_ , 9871–9881
  * Gorba et al. 2004 Gorba, C.; Geyer, T.; Helms, V. Brownian dynamics simulations of simplified cytochrome c molecules in the presence of a charged surface. _The Journal of Chemical Physics_ 2004, _121_ , 457–464
  * Mereghetti and Wade 2011 Mereghetti, P.; Wade, R. C. Diffusion of hydrophobin proteins in solution and interactions with a graphite surface. _BMC Biophysics_ 2011, _4_ , 9
  * Ravichandran et al. 2001 Ravichandran, S.; Madura, J. D.; Talbot, J. A Brownian Dynamics Study of the Initial Stages of Hen Egg-White Lysozyme Adsorption at a Solid Interface. _The Journal of Physical Chemistry B_ 2001, _105_ , 3610–3613
  * Gabdoulline and Wade 1997 Gabdoulline, R.; Wade, R. Simulation of the diffusional association of barnase and barstar. _Biophysical Journal_ 1997, _72_ , 1917–1929
  * Gabdoulline and Wade 1998 Gabdoulline, R. R.; Wade, R. C. Brownian dynamics simulation of protein-protein diffusional encounter. _Methods (San Diego, Calif.)_ 1998, _14_ , 329–341
  * Martinez et al. 2015 Martinez, M.; Bruce, N. J.; Romanowska, J.; Kokh, D. B.; Ozboyaci, M.; Yu, X.; Öztürk, M. A.; Richter, S.; Wade, R. C. SDA 7: A modular and parallel implementation of the simulation of diffusional association software. _Journal of Computational Chemistry_ 2015, _36_ , 1631–1645
  * Mereghetti and Wade 2012 Mereghetti, P.; Wade, R. C. Atomic Detail Brownian Dynamics Simulations of Concentrated Protein Solutions with a Mean Field Treatment of Hydrodynamic Interactions. _The Journal of Physical Chemistry B_ 2012, _116_ , 8523–8533
  * Mereghetti et al. 2014 Mereghetti, P.; Martinez, M.; Wade, R. C. Long range Debye-Hückel correction for computation of grid-based electrostatic forces between biomacromolecules. _BMC Biophysics_ 2014, _7_ , 4
  * Ermak and McCammon 1978 Ermak, D. L.; McCammon, J. A. Brownian dynamics with hydrodynamic interactions. _The Journal of Chemical Physics_ 1978, _69_ , 1352–1360
  * Schmidt et al. 2011 Schmidt, R. R.; Cifre, J. G. H.; de la Torre, J. G. Comparison of Brownian dynamics algorithms with hydrodynamic interaction. _The Journal of Chemical Physics_ 2011, _135_ , 084116
  * Saadat and Khomami 2014 Saadat, A.; Khomami, B. Computationally efficient algorithms for incorporation of hydrodynamic and excluded volume interactions in Brownian dynamics simulations: A comparative study of the Krylov subspace and Chebyshev based techniques. _The Journal of Chemical Physics_ 2014, _140_ , 184903
  * Geyer and Winter 2009 Geyer, T.; Winter, U. An O(N2) approximation for hydrodynamic interactions in Brownian dynamics simulations. _The Journal of Chemical Physics_ 2009, _130_ , 114905
  * Mereghetti et al. 2010 Mereghetti, P.; Gabdoulline, R. R.; Wade, R. C. Brownian dynamics simulation of protein solutions: Structural and dynamical properties. _Biophysical Journal_ 2010, _99_ , 3782–3791
  * Tokuyama and Oppenheim 1994 Tokuyama, M.; Oppenheim, I. Dynamics of hard-sphere suspensions. _Physical Review E_ 1994, _50_ , 6–9
  * Tokuyama et al. 2011 Tokuyama, M.; Moriki, T.; Kimura, Y. Self-diffusion of biomolecules in solution. _Physical Review E_ 2011, _83_ , 051402
  * Cichocki et al. 2000 Cichocki, B.; Jones, R. B.; Kutteh, R.; Wajnryb, E. Friction and mobility for colloidal spheres in Stokes flow near a boundary: The multipole method and applications. _The Journal of Chemical Physics_ 2000, _112_ , 2548
  * Lorentz 1907 Lorentz, H. Abhandlungen über theoretische Physik. _Monatshefte für Mathematik und Physik_ 1907, _18_ , A32–A32
  * Faxén 1922 Faxén, H. Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist. _Annalen der Physik_ 1922, _373_ , 89–119
  * Brenner 1961 Brenner, H. The slow motion of a sphere through a viscous fluid towards a plane surface. _Chemical Engineering Science_ 1961, _16_ , 242–251
  * Goldman et al. 1967 Goldman, A.; Cox, R.; Brenner, H. Slow viscous motion of a sphere parallel to a plane wall-I Motion through a quiescent fluid. _Chemical Engineering Science_ 1967, _22_ , 637–651
  * Cichocki and Jones 1998 Cichocki, B.; Jones, R. Image representation of a spherical particle near a hard wall. _Physica A: Statistical Mechanics and its Applications_ 1998, _258_ , 273–302
  * Ainley et al. 2008 Ainley, J.; Durkin, S.; Embid, R.; Boindala, P.; Cortez, R. The method of images for regularized Stokeslets. _Journal of Computational Physics_ 2008, _227_ , 4600–4616
  * Gabdoulline and Wade 2009 Gabdoulline, R. R.; Wade, R. C. On the contributions of diffusion and thermal activation to electron transfer between Phormidium laminosum plastocyanin and cytochrome f: Brownian dynamics simulations with explicit modeling of nonpolar desolvation interactions and electron transfer event. _Journal of the American Chemical Society_ 2009, _131_ , 9230–9238
  * Gabdoulline and Wade 1996 Gabdoulline, R. R.; Wade, R. C. Effective Charges for Macromolecules in Solvent. _The Journal of Physical Chemistry_ 1996, _100_ , 3868–3878
  * Elcock et al. 1999 Elcock, A. H.; Gabdoulline, R. R.; Wade, R. C.; McCammon, J. Computer simulation of protein-protein association kinetics: acetylcholinesterase-fasciculin. _Journal of Molecular Biology_ 1999, _291_ , 149–162
  * Nielsen and Museth 2006 Nielsen, M. B.; Museth, K. Dynamic Tubular Grid: An Efficient Data Structure and Algorithms for High Resolution Level Sets. _Journal of Scientific Computing_ 2006, _26_ , 261–299
  * Ozboyaci et al. 2016 Ozboyaci, M.; Martinez, M.; Wade, R. C. An Efficient Low Storage and Memory Treatment of Gridded Interaction Fields for Simulations of Macromolecular Association. _Journal of Chemical Theory and Computation_ 2016, _12_ , 4563–4577
  * Grahame 1947 Grahame, D. C. The Electrical Double Layer and the Theory of Electrocapillarity. _Chemical Reviews_ 1947, _41_ , 441–501
  * Northrup et al. 1984 Northrup, S. H.; Allison, S. a.; McCammon, J. A. Brownian dynamics simulation of diffusion-influenced bimolecular reactions. _The Journal of Chemical Physics_ 1984, _80_ , 1517
  * Wilson et al. 1992 Wilson, K.; Malcolm, B.; Matthews, B. Structural and thermodynamic analysis of compensating mutations within the core of chicken egg white lysozyme. _J Biol Chem._ 1992, _267_ , 10842–10849
  * Dolinsky et al. 2004 Dolinsky, T. J.; Nielsen, J. E.; McCammon, J. A.; Baker, N. A. PDB2PQR: An automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. _Nucleic Acids Research_ 2004, _32_ , 665–667
  * Cornell et al. 1996 Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. _Journal of the American Chemical Society_ 1996, _118_ , 2309–2309
  * Baker et al. 2001 Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: application to microtubules and the ribosome. _Proceedings of the National Academy of Sciences of the United States of America_ 2001, _98_ , 10037–41
  * Ortega et al. 2011 Ortega, A.; Amorós, D.; García de la Torre, J. Prediction of Hydrodynamic and Other Solution Properties of Rigid Proteins from Atomic- and Residue-Level Models. _Biophysical Journal_ 2011, _101_ , 892–898
  * Shapiro 2015 Shapiro, Y. E. _Encyclopedia of Surface and Colloid Science, Third Edition_ ; Taylor & Francis, 2015; pp 5701–5717
  * Daly et al. 2003 Daly, S. M.; Przybycien, T. M.; Tilton, R. D. Coverage-Dependent Orientation of Lysozyme Adsorbed on Silica. _Langmuir_ 2003, _19_ , 3848–3857
  * Robeson and Tilton 1996 Robeson, J. L.; Tilton, R. D. Spontaneous Reconfiguration of Adsorbed Lysozyme Layers Observed by Total Internal Reflection Fluorescence with a pH-Sensitive Fluorophore. _Langmuir_ 1996, _12_ , 6104–6113
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. _Journal of Molecular Graphics_ 1996, _14_ , 33–38
  * Antosiewicz et al. 2017 Antosiewicz, J. M.; Kamiński, K.; Długosz, M. Hydrodynamic Steering in Protein Association Revisited: Surprisingly Minuscule Effects of Considerable Torques. _The Journal of Physical Chemistry B_ 2017, _121_ , 8475–8491
  * Pagonabarraga and Rubí 1994 Pagonabarraga, I.; Rubí, J. M. Influence of hydrodynamic interactions on the adsorption process of large particles. _Physical Review Letters_ 1994, _73_ , 114–117
  * Wojtaszczyk and Avalos 1998 Wojtaszczyk, P.; Avalos, J. B. Influence of hydrodynamic interactions on the kinetics of colloidal particles’ adsorption. _Physical Review Letters_ 1998, _80_ , 754–757
  * Pérez et al. 2015 Pérez, C. A.; Moncho-Jordá, A.; Hidalgo-Álvarez, R.; Casanova, H. A comparative study on the effect of hydrodynamic interactions in the non-sequential deposition of concentrated colloidal dispersions: Stochastic rotation dynamics and Brownian dynamics simulations. _Molecular Physics_ 2015, _113_ , 3587–3597
  * Schaaf and Talbot 1989 Schaaf, P.; Talbot, J. Kinetics of Random Sequential Adsorption. _Physical Review Letters_ 1989, _62_ , 175–178
  * Su et al. 1998 Su, T. J.; Lu, J. R.; Thomas, R. K.; Cui, Z. F.; Penfold, J. The Adsorption of Lysozyme at the Silica-Water Interface: A Neutron Reflection Study. _Journal of Colloid and Interface Science_ 1998, _203_ , 419–429
  * Lu et al. 2004 Lu, J. R.; Swann, M. J.; Peel, L. L.; Freeman, N. J. Lysozyme Adsorption Studies at the Silica/Water Interface Using Dual Polarization Interferometry. _Langmuir_ 2004, _20_ , 1827–1832
  * Lundin et al. 2010 Lundin, M.; Elofsson, U. M.; Blomberg, E.; Rutland, M. W. Adsorption of lysozyme, $\beta$-casein and their layer-by-layer formation on hydrophilic surfaces: Effect of ionic strength. _Colloids and Surfaces B: Biointerfaces_ 2010, _77_ , 1–11
  * Kim et al. 2002 Kim, D. T.; Blanch, H. W.; Radke, C. J. Direct Imaging of Lysozyme Adsorption onto Mica by Atomic Force Microscopy. _Langmuir_ 2002, _18_ , 5841–5850
  * Daly et al. 2007 Daly, S. M.; Przybycien, T. M.; Tilton, R. D. Aggregation of lysozyme and of poly(ethylene glycol)-modified lysozyme after adsorption to silica. _Colloids and Surfaces B: Biointerfaces_ 2007, _57_ , 81–88
  * Steudle and Pleiss 2011 Steudle, A.; Pleiss, J. Modelling of lysozyme binding to a cation exchange surface at atomic detail: the role of flexibility. _Biophysical journal_ 2011, _100_ , 3016–3024