### 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 2*h* 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 2*h*
_{
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.