Governing equation and boundary conditions
The time-dependent Smoluchowski equation can be written as:
$$ \frac{d p(\mathbf{x},t)}{dt}=\nabla\cdot \mathbf{J} (\mathbf{x},t), ~~~~\mathbf{x}\in\Omega, $$
((2))
where p(x,t) is the concentration distribution of the reactants, and the concentration flux J(x,t) is defined as:
$$ \mathbf{J} (\mathbf{x},t) = D (\mathbf{x})[\nabla p(\mathbf{x},t) + \beta p(\mathbf{x},t) \nabla W(\mathbf{x})], $$
((3))
where D(x) is the diffusion coefficient; for simplicity, it is assumed to be constant. β=1/k
B
T is the inverse Boltzmann energy with the Boltzmann constant k
B
and kinetic temperature T. W(x) is the potential mean force (PMF) for the diffusing particle due to solvent-mediated interactions with the target molecule. The equation is solved in a three-dimensional domain Ω, subject to the following boundary conditions. First,
$$ p(\mathbf{x},t)= p_{\text{bulk}} ~\text{for} ~\mathbf{x}\in\Gamma_{b}, $$
((4))
specifying a Dirichlet BC on the outer boundary Γ
b
where the concentration is equal to a bulk concentration p
bulk. The outer boundary is often a spherical surface with a radius chosen to ensure that the ligand-protein potential is spherically symmetric and/or can be approximated analytically [6]. For the current study with mAChE, this outer boundary has radius R
2≈128 as determined following a procedure similar to Song et al. and Chen et al. [6,23]. Also following Song et al. and Chen et al., p is normalized such that p
bulk=1.
The active site boundary Γ
a
was modeled using either reactive Robin or absolute absorbing Dirichlet BC:
$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t) = \alpha p(\mathbf{x},t) ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}, $$
((5))
or
$$ p(\mathbf{x})= 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}, $$
((6))
respectively. The coefficient α is chosen to model an intrinsic reaction rate for the active site. Finally, a reflective Neumann BC is defined on the non-reactive boundary of molecule
$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t) = 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{m}. $$
((7))
Figure 10 shows the simulation domain along with all boundaries.
Given a solution to Eq. 2, the reaction rate is calculated from the integral of the flux across the reactive boundary [37]:
$$ k_{\text{on}}= p_{\text{bulk}}^{-1} \iint\limits_{\Gamma_{a}} \mathbf{n}(\mathbf{x}_{s}) \cdot \mathbf{J} (\mathbf{x}_{s},t) d\mathbf{x}^{\prime}_{s}. $$
((8))
In order to solve the Smoluchowski equation (Eq. 2) subject to the reactive Robin BC (Eq. 5), the simulation domain is extended to include a sub-domain Ω
a
that is separated from Ω by Γ
a
, and we then reformulate Eq. 2 as:
$$ {\fontsize{8.1pt}{9.6pt}\selectfont{ \begin{aligned} {} \frac{dp^{r}(\mathbf{x},t)}{dt}= &\nabla\cdot \left(D (\mathbf{x})[\nabla p^{r}(\mathbf{x},t) + \beta p^{r}(\mathbf{x},t) \nabla W(\mathbf{x})]\right) \\ - &\alpha p^{r} (\mathbf{x}, t) \iiint\limits_{\Omega_{a}} [\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')]\cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}', ~~~\mathbf{x}\in\Omega, \end{aligned}}} $$
((9))
subject to the reflective Neumann BC:
$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J}^{r} (\mathbf{x},t) = 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}. $$
((10))
The derivation of Eq. 9 is detailed in Appendix A, which demonstrates
$$ {\lim}_{h_{r} \to 0} p^{r}(\mathbf{x},t) = p(\mathbf{x},t). $$
((11))
In Eq. 9, the normalized kernel function, w(x), is a positive bell-shaped function with at least first continuous derivative and compact support κ
h
r
such that w(|r|>κ
h
r
)=0. The value of κ depends on the specific functional form of w(x), which is specified in Section ‘SPH discretization of equations and boundary conditions’. In particular, w(x) satisfies the following conditions:
$$ \iiint\limits_{\Omega\cup\Omega_{a}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' =1 $$
((12))
and
$$ {\lim}_{h_{r}\to 0}w(\mathbf{x}-\mathbf{x}', h_{r})=\delta(\mathbf{x}-\mathbf{x}'). $$
((13))
The normal unit vector n in Eq. 9 can be found in terms of a smoothed color function \(\tilde {\phi }\) as defined in Appendix A:
$$ \mathbf{n}(\mathbf{x}) = \frac{\nabla \tilde{\phi}(\mathbf{x})}{|\nabla \tilde{\phi}(\mathbf{x})|}, ~~~ \mathbf{x} \in \Omega\cup\Omega_{a}. $$
((14))
SPH discretization of equations and boundary conditions
In this section, we present SPH discretization of the Smoluchowski equation, using Eq. 2 if the Dirichlet BC is used and Eq. 9 if the Robin BC is assumed. To simplify notation, we omit superscript r for the variables in Eq. 9 in the subsequent derivations.
The domain Ω, and the boundaries Γ
a
and Γ
b
(extended as domains Ω
a
and Ω
b
respectively), are discretized with a set of N points with positions denoted by a vector r
i
(i=1,...,N). The points (which are commonly referred to as particles in SPH) are used to discretize and solve the governing equation. Initially, the particles are distributed uniformly (e.g., placed on a regular cubic lattice) with d
i
as the prescribed number density at r
i
. The discretization is based on a meshless interpolation scheme:
$$ A_{i} \approx \sum\limits_{j} \frac{A_{j}}{d_{j}}w(\mathbf{r}_{ij},h), $$
((15))
where, A
i
=A(r
i
) is a function defined at particle i, A
j
=A(r
j
) is the function defined at particle i’s neighboring particles j with distances r
ij
=r
i
−r
j
, and w(r
ij
,h) is the weighting kernel function. The interpolation scheme assumes a summation over all neighboring SPH particles but, due to the compact support of w, only particles within distance κ
h from r
i
have a non-zero contribution to the summation. Spatial derivatives of A can be calculated as
$$ \nabla_{i} A_{i} \approx \sum\limits_{j} \frac{A_{j}}{d_{j}} \nabla_{i} w(\mathbf{r}_{ij},h). $$
((16))
In the present work, we use a cubic spline kernel as the weighting function
$$ w(\mathbf{r},h) = \frac{1}{\pi h^{3}} \begin{cases} 1-\frac{3}{2}q^{2}+\frac{3}{4}q^{3} & ~~~0 \leq q \leq 1 \\ \frac{1}{4}(2-q)^{3} & ~~~1 < q \leq 2 \\ 0 & ~ q>2, \end{cases} $$
((17))
where q=|r|/h. With this form of weighting function, only particles within 2h distance from particle i contribute to the summations in the SPH equations. We have chosen h=1.3Δ
x where Δ
x is the size of the cubic lattice.
The SPH approximation of functions and their spatial derivatives allows the Smoluchowski equation subject to the Dirichlet BC (Eq. 2) to be written as a ODE governing the evolution of concentration on particle i as:
$$\begin{array}{@{}rcl@{}} \frac{d p_{i}}{dt} &=& \sum_{j\in \Omega\cup\Omega_{b}\cup\Omega_{a}} \frac{D_{i}+D_{j}}{d_{j}} (p_{i}-p_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&+ \beta \sum_{j\in \Omega} \frac{D_{i} p_{i}\,+\,D_{j} p_{j}}{d_{j}} (W_{i}\,-\,W_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}. \end{array} $$
((18))
The derivations of the first and second terms on the right-hand side of Equation 18 can be found in Monaghan et al. [27] where r
ij
is the magnitude of the vector r
ij
. If the reactive Robin BC on reactive boundary is imposed, Equation 9 is then solved instead and its corresponding SPH discretization form is:
$$\begin{array}{@{}rcl@{}} \frac{d p_{i}}{dt} &=& \sum_{j\in \Omega\cup\Omega_{b}} \frac{D_{i}+D_{j}}{d_{j}} (p_{i}-p_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&+\beta \sum_{j\in \Omega} \frac{D_{i} p_{i}+D_{j} p_{j}}{d_{j}} (W_{i}-W_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&- \alpha p_{i} \sum_{k\in \Omega_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h_{r})}{d r_{ik}}. \end{array} $$
((19))
To integrate the SPH Eqs. 18 and 19, an explicit Verlet scheme [38] is employed. The last term in Eq. 19 is obtained by discretizing the integral in Eq. 9 as a Riemann sum:
$$ \begin{aligned} {} \alpha p(\mathbf{x}, t) &\iiint\limits_{\Omega_{a}} [\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')]\cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' \\ &= \alpha p(\mathbf{x}, t) \sum_{k \in \Omega_{a}} \Delta V_{k} [\mathbf{n}(\mathbf{x}) + \mathbf{n}_{k}] \cdot \nabla_{\mathbf{x}} w(\mathbf{x} - \mathbf{r}_{k}, h_{r}), \end{aligned} $$
((20))
where \(\Delta V_{k} = \frac {1}{d_{k}}\) is the volume of particle k and \(\sum _{k \in \Omega _{a}}\) is the summation over the reactive boundary particles within 2h
r
distance from particle i. The SPH expression for calculating the normal unit vector is obtained as:
$$ \mathbf{n}_{i} = \frac{\sum\limits_{j \in \Omega\cup\Omega_{a}}\frac{1}{d_{j}}(\phi_{j}-\phi_{i})\nabla_{i} w(\mathbf{r}_{ij}, h_{r})}{\left |\sum\limits_{j \in \Omega\cup\Omega_{a}}\frac{1}{d_{j}}(\phi_{j}-\phi_{i})\nabla_{i} w(\mathbf{r}_{ij}, h_{r}) \right |}. $$
((21))
In the simulations presented below, we set h
r
=h but it could be set differently in future applications.
Note that the reflective Neumann BC (Equation 7 or 10) can be simply enforced in SPH by excluding the contribution from the boundary particles in the summation. The Dirichlet BC (Equation 4 or 6) is enforced by assigning the fixed boundary value of concentration on the boundary particles. If the Robin BC is imposed, the reaction rate k
on(t) can be calculated by Equation 42 in Appendix B and its corresponding SPH discretization form is:
$$ k_{\text{on}}= \sum_{i \in \Omega} \alpha p_{i} \left[\sum_{k\in \Gamma_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h)}{d r_{ik}}\right]. $$
((22))
Otherwise, when the Dirichlet BC is enforced on the reactive boundary, the discretization of Eq. 41 in Appendix B is:
$$ k_{\text{on}}= \sum_{i \in \Omega} (\mathbf{n}_{i} \cdot \mathbf{J}_{i}) \left[\sum_{k\in \Gamma_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h)}{d r_{ik}}\right], $$
((23))
where
$$ \begin{aligned} \mathbf{J}_{i} &= D_{i} \sum_{j\in \Omega\cup\Omega_{b}\cup\Omega_{a}} \frac{p_{j}-p_{i}}{d_{j}}\frac{\mathbf{r}_{ij}}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}\\ &\quad+ \beta D_{i} p_{i} \sum_{j\in \Omega} \frac{W_{j}-W_{i}}{d_{j}} \frac{\mathbf{r}_{ij}}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}. \end{aligned} $$
((24))
Calculation of potentials of mean force
We calculated the potential of mean force W(x) using the recently published 2.1 Å resolution structure of mAChE [39]. To prepare this structure for the calculation, we assigned titration states of ionizable residues using PROPKA [40] at pH 7, and we used PDB2PQR v1.8 [41,42] to assign atomic radii and charges. APBS v1.4 was used to perform a nonlinear Poisson-Boltzmann multi-grid calculation of the electrostatic potential over the entire simulation domain [43]. The small and large domains were set to 600 Å and 400 Å, respectively, with a fine grid spacing of 0.600 Å. For APBS calculations, we used the single Debye-Hückel boundary condition, a smoothed molecular surface, and protein and solvent dielectrics of 2 and 78.54, respectively. Atomic charges were mapped onto the grids using cubic B-spline discretization. The calculated potential was mapped onto the SPH discretization points of protein and solvent via trilinear interpolation.