# WeatherBench Probability: A benchmark dataset for probabilistic medium-range
weather forecasting along with deep learning baseline models

Sagar Garg
Technical University of Munich, Germany Stephan Rasp
ClimateAi, Inc
<EMAIL_ADDRESS>
Nils Thuerey
Technical University of Munich, Germany

###### Abstract

WeatherBench is a benchmark dataset for medium-range weather forecasting of
geopotential, temperature and precipitation, consisting of preprocessed data,
predefined evaluation metrics and a number of baseline models. WeatherBench
Probability extends this to probabilistic forecasting by adding a set of
established probabilistic verification metrics (continuous ranked probability
score, spread-skill ratio and rank histograms) and a state-of-the-art
operational baseline using the ECWMF IFS ensemble forecast. In addition, we
test three different probabilistic machine learning methods—Monte Carlo
dropout, parametric prediction and categorical prediction, in which the
probability distribution is discretized. We find that plain Monte Carlo
dropout severely underestimates uncertainty. The parametric and categorical
models both produce fairly reliable forecasts of similar quality. The
parametric models have fewer degrees of freedom while the categorical model is
more flexible when it comes to predicting non-Gaussian distributions. None of
the models are able to match the skill of the operational IFS model. We hope
that this benchmark will enable other researchers to evaluate their
probabilistic approaches.

##  1 Introduction

WeatherBench [Rasp et al., 2020] introduced a benchmark dataset for medium-
range, data-driven weather forecasting, focused on 3 and 5 day predictions of
global geopotential, temperature and precipitation fields. The three core
components of WeatherBench are 1) a preprocessed training and validation
dataset hosted in a convenient format on a public server, 2) a set of
precisely defined evaluation metrics and 3) a set of baselines, ranging from
simple ones such as climatologies and persistence to dynamical models run at
different resolutions and simple machine learning models.

The idea behind WeatherBench is to provide a framework in which to evaluate
progress in data-driven weather forecasting. Since its publication,
WeatherBench has been used in several papers. Weyn et al. [2020] built a
convolutional neural network (CNN) trained using a cubed-sphere coordinate
transformation. They then extended this approach for subseasonal prediction
[Weyn et al., 2021]. Rasp and Thuerey [2021] used a deep Resnet architecture
pretrained on climate model simulations that achieves comparable skill to
dynamical baselines at similar resolutions. Gneiting and Walz [2019] used the
WeatherBench dataset to compare different evaluation metrics. Keisler [2022]
trained a graph neural network and achieved impressive skill on metrics very
similar to those suggested in WeatherBench.

One short-coming of the original WeatherBench benchmark is that it exclusively
focuses on deterministic prediction. However, in medium-range weather
forecasting is was long ago recognized that the inherent unpredictability of
the atmosphere necessitates a probabilistic approach [Bauer et al., 2015]. For
this reason, virtually all operational weather forecasting centers run
ensemble forecasts, where each ensemble member represents one realistic
potential future outcome. Here we aim to extend WeatherBench by adding a
probabilistic dimension to the verification metrics and baselines. Building on
this we then train a range of probabilistic neural networks.

There is some previous work in building probabilistic data-driven methods.
Scher and Messori [2020] tested several methods to create a data-driven
ensemble forecast: random initial condition perturbations, singular vector
initial condition perturbations and different random seeds for the neural
networks. The latter results in the best probabilistic scores. Clare et al.
[2021] explored dropout ensembles, similar to this work. We will contrast our
results to theirs below.

We begin by briefly summarizing the original WeatherBench data and describing
the new evaluation metrics and baselines in Section 2. Then, in Section 3 we
review relevant literature in probabilistic machine learning before describing
the three methods we will use. This is followed by the results of our methods
in Section 4, which we discuss further in Section 5, before summarizing our
findings in Section 6.

##  2 Data, Evaluation and Baselines

###  2.1 Data

Here we will provide a brief summary of the WeatherBench dataset. For a
detailed description refer to Rasp et al. [2020]. The latest version of the
data is available at https://github.com/pangeo-data/WeatherBench. WeatherBench
contains regridded ERA5 [Hersbach et al., 2020] data from 1979 to 2018 at
hourly resolution with 2017 and 2018 set aside for evaluation. Here we use the
coarsest horizontal resolution, 5.625∘ (32 $\times$ 64 grid points), at 7
vertical levels (50, 250, 500, 600, 700, 850, and 925 hPa).

As a probabilistic baseline we use the 50-member operational ECMWF IFS
ensemble, archived in the THORPEX Interactive Grand Global Ensemble (TIGGE)
archive [Bougeault et al., 2010]. The target variables 500 hPa geopotential
(Z500), 850 hPa temperature (T850), 2 m temperature (T2M) and 6h accumulated
precipitation (PR) for 2017 and 2018 regridded to match WeatherBench are now
also available at Rasp [2022].

###  2.2 Evaluation

For evaluating the probabilistic forecasts, we chose several well-established
metrics. First, we evaluate the deterministic skill of the ensemble mean using
the root mean squared error (RMSE) defined as

$\textrm{RMSE}=\frac{1}{N_{\textrm{forecasts}}}\sum_{i}^{{N_{\textrm{forecasts}}}}\sqrt{\frac{1}{N_{\textrm{lat}}N_{\textrm{lon}}}\sum_{j}^{N_{\textrm{lat}}}\sum_{k}^{N_{\textrm{lon}}}L(j)(\bar{f}_{i,j,k}-t_{i,j,k})^{2}}$
(1)

where $\bar{f}$ is the ensemble mean forecast and $t$ is the ERA5 truth.
$L(j)$ is the latitude weighting factor for the latitude at the $j$th latitude
index:

$L(j)=\frac{\cos(\textrm{lat}(j))}{\frac{1}{N_{\textrm{lat}}}\sum_{j}^{N_{\textrm{lat}}}\cos(\textrm{lat}(j))}$
(2)

As a first-order measure of the reliability of the ensemble, we look at the
spread-skill ratio ($\textrm{Spread}/\textrm{RMSE}$), where Spread defined as

$\textrm{Spread}=\frac{1}{N_{\textrm{forecasts}}}\sum_{i}^{{N_{\textrm{forecasts}}}}\sqrt{\frac{1}{N_{\textrm{lat}}N_{\textrm{lon}}}\sum_{j}^{N_{\textrm{lat}}}\sum_{k}^{N_{\textrm{lon}}}L(j)\textrm{var}(f_{j,k})}$
(3)

where $\textrm{var}(f_{j,k})$ is the variance in the ensemble dimension. A
perfectly reliable ensemble should have a spread-skill ratio of 1. Smaller
values indicate underdispersion (i.e. the probabilistic forecast is
overconfident in its forecasts), larger values overdispersion [Fortin et al.,
2014].

As a more detailed visual evaluation of reliability we use rank histograms,
described in detail in Wilks [2006](p. 316). A perfect rank histogram is
uniform, while a U-shape indicates underdispersion and a reverse U-shape
overdispersion.

Finally, calibration and sharpness of the ensemble together is evaluated using
the continuous ranked probability score (CRPS) Wilks [2006](p. 302), defined
as:

$\textrm{CRPS}=\int_{-\infty}^{\infty}[F(f)-1(t\leq z)]^{2}dz,$ (4)

where $F$ denotes the CDF of the forecast distribution and $1(t\leq z)$ is an
indicator function that is 1 if $t\leq z$ and 0 otherwise. In the case of a
deterministic forecast the CRPS reduces to the mean absolute error. For
computing the CRPS for ensemble forecasts, we use the implementation in the
xskillscore111https://xskillscore.readthedocs.io Python package.

##  3 Probabilistic neural networks

We explore three different methods to create probabilistic predictions using
neural networks. The basic network architecture of all networks is a deep
Resnet, which mirrors that of Rasp and Thuerey [2021] with only the final
layer changed. The input variables are geopotential, temperature, zonal and
meridional wind and specific humidity at seven vertical levels (50, 250, 500,
600, 700, 850 and 925hPa), 2-meter temperature, 6-hourly accumulated
precipitation, the top-of-atmosphere incoming solar radiation, all at the
current time step $t$, $t-6h$ and $t-12h$, and, finally three constant fields:
the land-sea mask, orography and the latitude at each grid point. All
variables, levels and time steps are stacked to create a 114 channel input.
All features are normalized
$(x\leftarrow\frac{x-\mathrm{mean}}{\mathrm{std}})$, except precipitation
where $\mathrm{mean}=0$ to keep values positive. Precipitation is also log-
transformed
$(\mathrm{PR}\leftarrow\log(\mathrm{PR}+\epsilon)-\log(\epsilon))$, where
$\epsilon=0.001$. For evaluation the transformation is reversed again. All
models were trained using ERA5 data from 1979 to 2015, with 2016 used for
validation and 2017/18 for testing. Predictions are made for a lead-time of 3
days.222Ideally we would have also trained models for 5 days lead time and
continuous models. However, because the first two authors moved on to
different positions shortly after the initial work was finished, this was not
possible. We would expect qualitatively similar results for other lead times.
The training procedure is the same as in Rasp and Thuerey [2021]. No
pretraining with climate models was used for any of the experiments in this
paper.

###  3.1 Monte-Carlo Dropout (MC Dropout)

Dropout [Srivastava et al., 2014] was introduced as a simple way of combating
overfitting in neural networks. Typically, dropout is only applied during
training but then turned off during inference. In fact, our base architecture
described above uses dropout with a rate of 0.1. However, it is also possible
to use dropout during inference to obtain a stochastic prediction [Gal and
Ghahramani, 2015]. Here we train our model with different dropout rates, 0,
0.1, 0.2 and 0.5, and create 50 random realizations for the test dataset to
match the 50-member IFS ensemble.

###  3.2 Parametric prediction

Another way to create a probabilistic forecast using neural networks is to
directly predict the parameters of a distribution. This is a common technique
in ensemble post-processing. Here we follow the methodology of Rasp and Lerch
[2018]. The first step is to prescribe a distribution for the target
variables. For Z500, T850 and T2M we use a Gaussian distribution parameterized
by its mean $\mu$ and variance $\sigma$. This means that instead of a single
output channel, the network now has two channels for $\mu$ and $\sigma$. To
optimize for the distributional parameters we chose the CRPS as a
probabilistic loss function. For a Gaussian distribution, there exists a
closed-form, differentiable solution of the CRPS:

$\mathrm{}{CRPS}\left(F_{\mu,\sigma},y\right)=\sigma\left\\{\frac{y-\mu}{\sigma}\left[2\Phi\left(\frac{y-\mu}{\sigma}\right)-1\right]+2\varphi\left(\frac{y-\mu}{\sigma}\right)-\frac{1}{\sqrt{\pi}}\right\\}$
(5)

For precipitation, a Gaussian distribution would be a poor fit because
precipitation tends to be very skewed and is censored at zero. Several
distributions have been proposed in literature such as a left-censored
generalized extreme value distribution [Scheuerer, 2014a]. We struggled to
implement these for neural network training. The closed-forms of these
distributions are quite complicated with limited valid ranges for some
parameters. Unfortunately, we did not end up figuring out how to achieve a
stable training using these distributions. We will discuss this in Section 5.

For evaluating the parametric forecasts we computed the CRPS and the spread-
skill ratio directly from the predicted parametric distributions but created a
randomly sampled 50-member ensemble for computing the rank histogram to match
the other methods.

###  3.3 Categorical prediction

Finally, we test a non-parametric, probabilistic method based on discretized
bins. This has recently been used by Sønderby et al. [2020] and Espeholt et
al. [2021] for precipitation nowcasting. The idea is to divide the value range
to be predicted into discrete bins and predict the probability for each of
these bins. The problem is then simply a multi-class classification problem.
The number of channels in the last layer are changed to match the number of
bins, followed by a softmax layer. The targets are one-hot encoded, with a one
indicating the correct bin and zeros otherwise. As a loss function we use the
standard categorical cross-entropy, also called log-loss

There are some subtleties in how to construct the bins. Generally, there is a
trade-off between the bin width, which can be thought of as the probabilistic
resolution, and the number of bins. Higher resolution will lead to finer
probability distributions but will also make training harder. Further, the
residual blocks of the core architecture have 128 channels which sets an upper
limit to the number of channels. However, taking temperature prediction as an
example, temperature values can fluctuate by around 100C from the poles to the
deserts, which means that the bins would be rather large. For this reason, for
geopotential and temperature, we predict the change in time
($X_{t=3d}-X_{t=0}$; after normalization) rather than the absolute values.
These differences will have a much smaller spread. Before evaluation, we then
simply add the difference.

Note that we train our model separately for each variable, unlike in the other
experiments where Z, T, T2M are trained together and precipitation trained
separately. Empirically, this has given better results compared to training
all variables together. For Z, T and T2M, we use 50 bins ranging from -1.5 to
1.5. For precipitation, we use 50 bins from 0 to 8 (after the log-transform).

##  4 Results

Table 1: RMSE, Spread-skill ratio and CRPS for different methods for 3 day
forecasts.

| RMSE of ensemble mean | Spread-Skill Ratio | CRPS
---|---|---|---
| Z500 (m2s-2) | T850 (K) | T2M (K) | TP (mm) | Z500 | T850 | T2M | TP | Z500 (m2s-2) | T850 (K) | T2M (K) | TP (mm)
| MC Dropout
---
(Dr=0.1)
312.96 | 1.80 | 1.52 | 2.2 | 0.40 | 0.34 | 0.39 | 0.13 | 155.70 | 1.03 | 0.77 | 0.57
Parametric | 315.30 | 1.82 | 1.55 | - | 0.87 | 0.92 | 0.95 | - | 142.67 | 0.90 | 0.70 | -
Categorical | 327.48 | 1.80 | 1.49 | 2.2 | 0.91 | 0.94 | 0.98 | 1.24 | 142.59 | 0.87 | 0.65 | 0.47
Tigge (3/5 days) | 145/297 | 1.20/1.73 | 1.26/1.57 | 2.02/2.15 | 1.05/1.00 | 0.93/0.96 | 0.69/0.80 | 0.84/0.85 | 65.6/127 | 0.60/0.83 | 0.58/0.70 | 0.41/0.47
Determinstic | 313.70 | 1.79 | 1.53 | 2.2 | - | - | - | - | 194.90 | 1.24 | 0.96 | 0.65
Clare et al. [2021] | 375/627 | 2.11/2.91 |  |  |  |  |  |  | 211/1500 | 1.22 / 1.69 |  |

Figure 1: (top) Ensemble mean RMSE for TIGGE baseline and different data-
driven methods. (bottom) Same for CRPS.

The results for the different verification metrics are shown in Table 1 and
Fig. 1. For MC dropout we also tested the sensitivity to the dropout rate,
shown in Fig. 2 (only Z500 is shown; other variables behave qualitatively
similarly). The ensemble mean RMSE and the CRPS is lowest for a dropout rate
of 0.1. The spread-skill ratio shows that the dropout ensemble is severely
underdispersive with the spread being less than half of what it should be.
These results are consistent with Clare et al. [2021]. Increasing the dropout
rate beyond 0.1 has a slight positive effect on the spread-skill ratio but
this comes at the cost of severely decreasing skill measured by the ensemble
mean RMSE and CRPS. This will be further discussed in Section 5. For now, we
use a dropout rate of 0.1 for all further evaluation.

The ensemble mean RMSE is reasonably similar between all deep learning methods
with slight fluctuations between the variables. We also added the RMSE of the
deterministic deep learning model without pretraining described in Rasp and
Thuerey [2021] as well as the scores reported in Clare et al. [2021]. The
operational TIGGE ensemble outperforms the data-driven methods quite
significantly. The only exception is precipitation, but, as already discussed
in Rasp and Thuerey [2021], the RMSE is not a particularly fitting metric to
evaluate a highly intermittent and skewed field like precipitation.

The spread-skill ratio shows that MC dropout is severely underdispersive for
all variables. The parametric and categorical networks have similar ratios
just under one for geopotential and temperature, while the ratio is above one
for categorical precipitation. TIGGE has ratios very close to one for Z500 and
T850 but is somewhat underdispersive for T2M. It is important to remember here
though that the TIGGE results are not post-processed, which would likely
improve calibration.

Further insight into calibration and reliability can be gained from looking at
the rank histograms in Fig. 3. Here, we see a very strong U-shape for MC
dropout. TIGGE only has a slight U-shape for Z500 and T850 but strong peaks on
the left-hand side for T2M and PR. For T2M, this hints at a bias issue but
would require further investigation. The peak for PR could be a manifestation
of the well-known drizzle bias. The parametric predictions are skewed which is
a sign of a high-bias. We did not further investigate the cause for this. The
categorical predictions do not have a bias and are only moderately U-shaped,
suggesting good calibration.

The CRPS evaluates both calibration and sharpness. Again we see similar scores
for the parametric and categorical networks, with slightly better results for
categorical. MC dropout is worse on CRPS, a consequence of its
underdispersion. Further, we see that the deterministic forecast (for which
the CRPS is simply the MAE) naturally scores worse since it does not provide
any probabilistic information. Again, TIGGE has significantly better CRPS
scores compared to the data driven methods.

Figure 2: RMSE, CRPS and spread-skill ratio for different dropout rates for
Z500. Figure 3: Rank histograms for different variables and methods. All rank
histograms have 51 bins.

##  5 Discussion

###  5.1 Verification caveats

There are several aspects to consider in the above verification. First, the
verification is done at a very coarse grid, 5.625 degrees. The operational IFS
ensemble forecast runs at a resolution of around 13 km. The required
regridding could potentially add a source of error, which might partly account
for the initial error at $t=0h$ of the operational ensemble. However, another
cause of initial error is the fact that the operational forecasts are started
from different initial conditions, the operational analysis as opposed to the
ERA5 reanalysis. For a fair comparison of pure forecasting skill, each method
should be compared against its own analysis, starting from zero error at
$t=0h$.

For precipitation, the verification metrics used, RMSE and CRPS, are poor
choices as they suffer heavily from the double penalty problem. In the
original WeatherBench paper and here, we still decided to include the metrics
but urge everyone to interpret these results with the appropriate caution.

###  5.2 Parametric versus categorical prediction

Both, parametric and categorical forecasts performed well in this study, with
similar verification scores. They also both have been used successfully in
previous research, even if parametric methods have a much longer history,
especially in post-processing. Parametric approaches have the advantage of
simplicity. Only a few (usually 2-3) parameters are to be estimated to produce
a full distribution. For parameters where the underlying distribution is
thought to be well known, e.g. temperature and geopotential, a parametric
approach is a great choice. However, for more unusually distributed variables,
e.g. precipitation, problems can arise. First, often the "true" distribution
is not know which requires assumptions to be made. If those assumptions are
not true, an artificial error is introduced. Second, more complex distribution
functions and their derivatives required for gradient descent optimization can
be difficult to implement. In the case of precipitation, we attempted to
implement a censored extreme value distribution described in Scheuerer
[2014b]. However, we encountered exploding gradients during optimization,
likely a result of the limited valid ranges of some of the distributional
parameters. One could likely fix this by adding constraints on the parameter
values and tinkering with the implementation of the loss. However, it shows
some of the difficulty with the parametric approach.

The categorical approach, on the other hand, requires no assumptions on the
distribution. In theory, any distribution can be learned. This is potentially
a big advantage of this method, especially for more complex distributions. The
trade-off is that instead of learning 2-3 values, now $>$100 values have to be
estimated, which can be problematic for small sample sizes. In fact, we have
observed that initially the estimated "distributions" are not smooth and only
become so after training. Further, one introduces a probabilistic resolution,
i.e. the bin width. As described in the Methods, this required a reframing of
the problem for geopotential and temperature to predict the change rather than
the absolute value. Nevertheless, the categorical approach, given enough
samples, turns out to work very well empirically and is easy to implement.

###  5.3 Spatially and temporally correlated forecasts

A final topic to discuss is spatial correlation in the forecast fields. The
parametric and categorical methods predict separate distributions at each grid
point and at each point in time. While this is fine for some applications,
others also require knowledge about spatial and temporal patterns. For flood
forecasting, for example, one is not simply interested in the distribution of
rain at a single location in 5 days time but rather in the distribution of
cumulative rain in the entire catchment area. For rain, which can have very
different characteristics, e.g. convective versus stratiform, adding the
distributions in space and time will not lead to the correct answer. Further,
visual inspection of weather maps is still a very common use case for decision
makers. With point-wise distribution methods, it is impossible to create
"realistic" weather maps. Monte-Carlo dropout produces an ensemble of
spatially coherent outputs but produces less reliable results compared to the
other methods.

To solve this issue, generative methods could be an option. Generative
adversarial networks (GANs) [Goodfellow et al., 2014] are one way of producing
realistic weather output. Recently this has been used for precipitation
nowcasting [Ravuri et al., 2021] and downscaling [Price and Rasp, 2022, Harris
et al., 2022]. While GANs are, in principal, a well-suited method to produce
coherent output, in practice they can be hard to train and empirically don’t
always produce the correct distribution.

##  6 Conclusion

Here we introduced WeatherBench Probability, a probabilistic extension to the
WeatherBench benchmark dataset. We added commonly used probabilistic
verification metrics along with an operational numerical weather prediction
baseline. Further, we trained three different sets of deep learning models to
produce probabilistic predictions for 3 days ahead, Monte-Carlo dropout,
parametric and categorical networks. Our version of Monte-Carlo dropout
produced severely underdispersive results. The other two methods produced
similarly good results across all scores, with trade-offs between the two
discussed in Section 5.

We hope that this paper complements the WeatherBench dataset and enables fast
and easy comparison between different data-driven forecasting methods. In the
last year or two, a number of promising approaches have been proposed and
others are in development which makes for an exciting future for this research
field.

## References

  * Rasp et al. [2020] Stephan Rasp, Peter D. Dueben, Sebastian Scher, Jonathan A. Weyn, Soukayna Mouatadid, and Nils Thuerey.  WeatherBench: A Benchmark Data Set for Data-Driven Weather Forecasting.  _Journal of Advances in Modeling Earth Systems_ , 12(11), 11 2020.  ISSN 1942-2466.  doi: 10.1029/2020MS002203.  URL https://onlinelibrary.wiley.com/doi/10.1029/2020MS002203.
  * Weyn et al. [2020] Jonathan A. Weyn, Dale R. Durran, and Rich Caruana.  Improving Data-Driven Global Weather Prediction Using Deep Convolutional Neural Networks on a Cubed Sphere.  _Journal of Advances in Modeling Earth Systems_ , 12(9), 9 2020.  ISSN 1942-2466.  doi: 10.1029/2020MS002109.  URL https://onlinelibrary.wiley.com/doi/10.1029/2020MS002109.
  * Weyn et al. [2021] Jonathan A. Weyn, Dale R. Durran, Rich Caruana, and Nathaniel Cresswell-Clay.  Sub-Seasonal Forecasting With a Large Ensemble of Deep-Learning Weather Prediction Models.  _Journal of Advances in Modeling Earth Systems_ , 13(7):e2021MS002502, 7 2021.  ISSN 1942-2466.  doi: 10.1029/2021MS002502.  URL https://onlinelibrary.wiley.com/doi/full/10.1029/2021MS002502https://onlinelibrary.wiley.com/doi/abs/10.1029/2021MS002502https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2021MS002502.
  * Rasp and Thuerey [2021] Stephan Rasp and Nils Thuerey.  Data-Driven Medium-Range Weather Prediction With a Resnet Pretrained on Climate Simulations: A New Model for WeatherBench.  _Journal of Advances in Modeling Earth Systems_ , 13(2):e2020MS002405, 2 2021.  ISSN 19422466.  doi: 10.1029/2020MS002405.  URL https://github.com/pangeo-data/WeatherBench.
  * Gneiting and Walz [2019] Tilmann Gneiting and Eva-Maria Walz.  Receiver operating characteristic (ROC) movies, universal ROC (UROC) curves, and coefficient of predictive ability (CPA).  _arXiv_ , 11 2019.  URL http://arxiv.org/abs/1912.01956.
  * Keisler [2022] Ryan Keisler.  Forecasting Global Weather with Graph Neural Networks.  2 2022.  doi: 10.48550/arxiv.2202.07575.  URL https://arxiv.org/abs/2202.07575v1.
  * Bauer et al. [2015] Peter Bauer, Alan Thorpe, and Gilbert Brunet.  The quiet revolution of numerical weather prediction.  _Nature_ , 525(7567):47–55, 9 2015.  ISSN 0028-0836.  doi: 10.1038/nature14956.  URL http://www.nature.com.emedien.ub.uni-muenchen.de/nature/journal/v525/n7567/full/nature14956.html.
  * Scher and Messori [2020] Sebastian Scher and Gabriele Messori.  Ensemble neural network forecasts with singular value decomposition.  2 2020.  URL https://arxiv.org/abs/2002.05398.
  * Clare et al. [2021] Mariana Clare, Omar Jamil, and Cyril Morcrette.  A computationally efficient neural network for predicting weather forecast probabilities.  3 2021.  URL http://arxiv.org/abs/2103.14430.
  * Hersbach et al. [2020] Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, Adrian Simmons, Cornel Soci, Saleh Abdalla, Xavier Abellan, Gianpaolo Balsamo, Peter Bechtold, Gionata Biavati, Jean Bidlot, Massimo Bonavita, Giovanna Chiara, Per Dahlgren, Dick Dee, Michail Diamantakis, Rossana Dragani, Johannes Flemming, Richard Forbes, Manuel Fuentes, Alan Geer, Leo Haimberger, Sean Healy, Robin J. Hogan, Elías Hólm, Marta Janisková, Sarah Keeley, Patrick Laloyaux, Philippe Lopez, Cristina Lupu, Gabor Radnoti, Patricia Rosnay, Iryna Rozum, Freja Vamborg, Sebastien Villaume, and Jean-Noël Thépaut.  The ERA5 Global Reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , page qj.3803, 5 2020.  ISSN 0035-9009.  doi: 10.1002/qj.3803.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qj.3803.
  * Bougeault et al. [2010] Philippe Bougeault, Zoltan Toth, Craig Bishop, Barbara Brown, David Burridge, De Hui Chen, Beth Ebert, Manuel Fuentes, Thomas M. Hamill, Ken Mylne, Jean Nicolau, Tiziana Paccagnella, Young-Youn Park, David Parsons, Baudouin Raoult, Doug Schuster, Pedro Silva Dias, Richard Swinbank, Yoshiaki Takeuchi, Warren Tennant, Laurence Wilson, and Steve Worley.  The THORPEX Interactive Grand Global Ensemble.  _Bulletin of the American Meteorological Society_ , 91(8):1059–1072, 8 2010.  doi: 10.1175/2010BAMS2853.1.  URL http://journals.ametsoc.org/doi/10.1175/2010BAMS2853.1.
  * Rasp [2022] Stephan Rasp.  TIGGE data for WeatherBench Probability.  4 2022.  doi: 10.5281/ZENODO.6480846.  URL https://doi.org/10.5281/zenodo.6480846#.YmVG6z5Nf1o.mendeley.
  * Fortin et al. [2014] V. Fortin, M. Abaza, F. Anctil, and R. Turcotte.  Why Should Ensemble Spread Match the RMSE of the Ensemble Mean?  _Journal of Hydrometeorology_ , 15(4):1708–1713, 8 2014.  ISSN 1525-755X.  doi: 10.1175/JHM-D-14-0008.1.  URL http://journals.ametsoc.org/doi/abs/10.1175/JHM-D-14-0008.1.
  * Wilks [2006] Daniel S Wilks.  _Statistical Methods in the Atmospheric Sciences_.  Elsevier, 2006.  ISBN 0127519661.  URL http://cds.cern.ch/record/992087.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A Simple Way to Prevent Neural Networks from Overfitting.  _Journal of Machine Learning Research_ , 15:1929–1958, 2014\.  URL http://www.jmlr.org/papers/volume15/srivastava14a/srivastava14a.pdf?utm_content=buffer79b43&utm_medium=social&utm_source=twitter.com&utm_campaign=buffer.
  * Gal and Ghahramani [2015] Yarin Gal and Zoubin Ghahramani.  Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.  _33rd International Conference on Machine Learning, ICML 2016_ , 3:1651–1660, 6 2015.  URL http://arxiv.org/abs/1506.02142.
  * Rasp and Lerch [2018] Stephan Rasp and Sebastian Lerch.  Neural Networks for Postprocessing Ensemble Weather Forecasts.  _Monthly Weather Review_ , 146(11):3885–3900, 11 2018.  doi: 10.1175/MWR-D-18-0187.1.  URL http://journals.ametsoc.org/doi/10.1175/MWR-D-18-0187.1.
  * Scheuerer [2014a] M. Scheuerer.  Probabilistic quantitative precipitation forecasting using Ensemble Model Output Statistics.  _Quarterly Journal of the Royal Meteorological Society_ , 140(680):1086–1096, 4 2014a.  ISSN 00359009.  doi: 10.1002/qj.2183.  URL http://doi.wiley.com/10.1002/qj.2183.
  * Sønderby et al. [2020] Casper Kaae Sønderby, Lasse Espeholt, Jonathan Heek, Mostafa Dehghani, Avital Oliver, Tim Salimans, Shreya Agrawal, Jason Hickey, and Nal Kalchbrenner.  MetNet: A Neural Weather Model for Precipitation Forecasting.  3 2020.  URL http://arxiv.org/abs/2003.12140.
  * Espeholt et al. [2021] Lasse Espeholt, Shreya Agrawal, Casper Sønderby, Manoj Kumar, Jonathan Heek, Carla Bromberg, Cenk Gazen, Jason Hickey, Aaron Bell, and Nal Kalchbrenner.  Skillful Twelve Hour Precipitation Forecasts using Large Context Neural Networks.  11 2021.  doi: 10.48550/arxiv.2111.07470.  URL https://arxiv.org/abs/2111.07470v1.
  * Scheuerer [2014b] M. Scheuerer.  Probabilistic quantitative precipitation forecasting using Ensemble Model Output Statistics.  _Quarterly Journal of the Royal Meteorological Society_ , 140(680):1086–1096, 4 2014b.  ISSN 00359009.  doi: 10.1002/qj.2183.  URL http://doi.wiley.com/10.1002/qj.2183.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Advances in neural information processing systems_ , pages 2672–2680, 2014.
  * Ravuri et al. [2021] Suman Ravuri, Karel Lenc, Matthew Willson, Dmitry Kangin, Remi Lam, Piotr Mirowski, Megan Fitzsimons, Maria Athanassiadou, Sheleem Kashem, Sam Madge, Rachel Prudden, Amol Mandhane, Aidan Clark, Andrew Brock, Karen Simonyan, Raia Hadsell, Niall Robinson, Ellen Clancy, Alberto Arribas, and Shakir Mohamed.  Skillful Precipitation Nowcasting using Deep Generative Models of Radar.  4 2021.  URL https://arxiv.org/abs/2104.00954v1.
  * Price and Rasp [2022] Ilan Price and Stephan Rasp.  Increasing the accuracy and resolution of precipitation forecasts using deep generative models.  page 2022, 3 2022.  doi: 10.48550/arxiv.2203.12297.  URL https://arxiv.org/abs/2203.12297v1.
  * Harris et al. [2022] Lucy Harris, Andrew T. T. McRae, Matthew Chantry, Peter D. Dueben, and Tim N. Palmer.  A Generative Deep Learning Approach to Stochastic Downscaling of Precipitation Forecasts.  4 2022.  doi: 10.48550/arxiv.2204.02028.  URL https://arxiv.org/abs/2204.02028v1.


# Surface-Facilitated Trapping by Active Sites: From Catalysts to Viruses

Mikita M. Misiura  Department of Chemistry and Center for Theoretical
Biological Physics, Rice University, Houston, TX 77005, USA  Alexander M.
Berezhkovskii  Mathematical and Statistical Computing Laboratory, Office of
Intramural Research, Center for Information Technology, National Institutes of
Health, Bethesda, MD 20892  Sergey M. Bezrukov  Section on Molecular
Transport, Eunice Kennedy Shriver National Institute of Child Health and Human
Development, National Institutes of Health, Bethesda, MD 20892  Anatoly B.
Kolomeisky  Department of Chemistry and Center for Theoretical Biological
Physics, Rice University, Houston, TX 77005, USA<EMAIL_ADDRESS>
###### Abstract

Trapping by active sites on surfaces plays important roles in various chemical
and biological processes, including catalysis, enzymatic reactions, and viral
entry into host cells. However, the mechanisms of these processes remain not
well understood, mostly because the existing theoretical descriptions are not
fully accounting for the role of the surfaces. Here we present a theoretical
investigation on the dynamics of surface-assisted trapping by specific active
sites. In our model, a diffusing particle can be occasionally reversibly bound
to the surface and diffuse on it before reaching the final target site. An
approximate theoretical framework is developed, and its predictions are tested
by Brownian Dynamics computer simulations. It is found that the surface
diffusion can be crucial in mediating the association to active sites. Our
theoretical predictions work reasonably well as long as the size of the active
site is much smaller than the overall surface area. Potential applications of
our method are discussed.

[Rice University]Department of Chemical and Biomolecular Engineering, Rice
University, Houston, TX 77005, USA [Rice University]Department of Physics and
Astronomy, Rice University, Houston, TX 77005, USA

## Introduction

There are multiple chemical and biological processes that involve association
of particles to specific reactive sites that are located on surfaces. Examples
include catalytic reactions on heterogeneous catalysts1 and enzymatic
processes on active sites located on the protein surfaces.2 Analysis of these
systems typically considers the molecular surfaces as inactive and not
participating in these processes. However, the surfaces of solid-state
catalysts and enzymatically active proteins are clearly not inert, and the
reacting species interact with them affecting the overall outcome of the
process. These interactions might involve hydrogen bonds, electrostatic forces
and hydrophobic interactions, and they might achieve significant strengths
that could strongly modify the overall dynamics of the process. For instance,
recent all-atom molecular dynamics simulations of association of small ligands
with dihydrofolate reductase proteins explicitly demonstrated the effect of
surface facilitation on the kinetics.3

A relatively new example of the reaction facilitated by surface diffusion is
channel-catalyzed transport of membrane-bound peripheral proteins.4, 5 It was
recently shown that $\alpha$-synuclein – the cytosolic protein implicated in
Parkinson disease – enters mitochondrial periplasmic space through a process
involving reversible binding to the outer membrane of mitochondria as the
first crucial step. Protein then diffuses on the membrane surface and either
escapes into the bulk or is captured by the $\beta$-barrel nanopore of the
voltage-dependent anion channel, VDAC. The capture of a membrane-bound
$\alpha$-synuclein molecule happens via its highly negatively charged
C-terminal tail entering the nanopore with a net-positive charge. The on-rate
of the process is defined by the surface concentration of $\alpha$-synuclein
molecules and their conformation, both of which are affected by the membrane
lipid composition.6, 7 Since the lipid effects are strong it was suggested
that changes in the lipid composition of the host membrane could constitute a
potent regulatory mechanism of the $\alpha$-synuclein interaction with VDAC,
including its VDAC-catalyzed translocation into the periplasmic space where
$\alpha$-synuclein obstructs the normal functioning of the machinery of the
inner membrane integral proteins of mitochondria.

The crucial involvement of the membrane surface in the $\alpha$-synuclein
capture was clearly demonstrated in experiments with another $\beta$-barrel
nanopore,8 formed by $\alpha$-hemolysin which is known to have a pronounced
structural asymmetry.9 On one side of the membrane, it protrudes above the
membrane surface for about 5 nanometers (the cap side of the channel
structure), while on the other side it is flush with the surface (the stem
side). It was found that the protruding, cap-side opening of the channel
displays several orders of magnitude lower capture rates, presumably capturing
$\alpha$-synuclein only from the bulk solution, compared to the rates of
capture from the membrane surface by the other, stem-side opening. This huge
difference in the capture on-rates at otherwise symmetric conditions, as well
as the highly expressed sensitivity of the on-rate to the membrane lipid
composition, support the picture according to which the $\alpha$-synuclein
molecule first binds to the membrane and then diffuses over its surface before
being captured by the nanopore or released back to the bulk solution.

Another example of the complex process where interactions with the surfaces
are important is the entry of viruses into host cells.10 Recent experimental
studies indicated that the viral penetration is a very complex process that
involves multiple pathways. It was shown that viruses often bind first to a
random spot on the cell’s surface and then diffuse along the surface before
binding to the appropriate receptor that allow them to enter into the cell.
11, 12 High-resolution fluorescence microscopy and single particle tracking
techniques have illustrated how $\lambda$ phage is able to first bind to a
random location on the E. coli surface and then use complicated target-
searching process to reach one of the cell’s poles.11 It was also shown that
sometimes multiple receptors are needed for a virus to enter the cell, which
might lead to even longer and more complex diffusion processes.13

The problem of association of small ligands with specific active sites on the
surfaces has been investigated before using various theoretical tools.14, 15,
16, 17 A diffusion limited association, when the reaction of binding to the
surface (both specific and non-specific) is instantaneous, has been considered
by Berg.14 In the elegant approach by Zhou and Szabo, the effect of the
surface was modelled as a short-range attraction potential, leading to
analytical results for different systems.15, 16 The surface-assisted
association has been analyzed as a dynamic search process by one of us,17 but
the analytical results have been obtained only in certain limiting cases. It
should be noted here that the enhancement of the association rates due to the
presence of the surface is an example of more general dynamic target search
facilitation which is accomplished by switching the participating particles
between different dynamic regimes with different dimensionalities. In this
case, particles are alternating between 3D bulk motion and 2D surface sliding.
It is worth mentioning that significant progress has been achieved earlier in
understanding the mechanisms of facilitated diffusion,18, 19, 20, 21, 22, 23,
24, 25, 26 where the dynamics switches between 1D and 3D regimes. This 3D/1D
facilitation mechanism has been extensively discussed due its importance for
protein target search on DNA that initiates all major biological processes.
These ideas were initially stimulated by pioneering work of Adam and Delbruck
on dimensionality reduction and diffusion in biological systems.27

In this work, we propose a new theoretical framework to describe the surface-
assisted molecular/particle associations to specific sites. It can be viewed
as a generalization of the Smoluchovskii-Collins-Kimball method that accounts
for nonspecific interactions in the system. We model one of the binding
partners as an immobile partially absorbing sphere with a small perfectly
absorbing patch that corresponds to the active site, and the second binding
partner is viewed as a diffusing point particle that might bind to the
surface. It can diffuse then along the surface until it finds the active site,
or it can dissociate back into the bulk solution. If the size of the active
site is small, the overall binding process is considered as trapping by the
uniform partially absorbing sphere. This allows us to derive explicit formulas
for dynamic properties of the process. Extensive Brownian Dynamics (BD)
simulations are utilized then to test our theoretical predictions and the
range of parameters where our approximate method works.

## Theoretical Method

To obtain a theoretical description of the surface-assisted trapping of
particles to active sites, let us consider a system schematically shown in
Fig. 1. We have an immobile sphere of radius $R$ surrounded by the medium with
diffusing point particles. When the sphere is perfectly absorbing, a particle
touching its surface is instantly trapped. In this case, the steady-state flux
of trapped particles is given by the Smoluchovski formula,28

$J_{Sm}=k_{Sm}c_{\infty},\quad k_{Sm}=4\pi DR,$ (1)

where $k_{Sm}$ is the Smoluchowski rate constant, $D$ is the particle
diffusivity and $c_{\infty}$ is the concentration of particles far away from
the sphere.

If the sphere is only partially absorbing, the steady-state flux of trapped
particles is lower than $J_{Sm}$ since a particle touching the sphere has a
chance to escape back to the bulk solution. This flux is given by the Collins-
Kimball theory,29

$J_{CK}=k_{CK}c_{\infty}=J_{Sm}P_{tr}^{CK},$ (2)

where $k_{CK}$ is the Collins-Kimball rate constant, and $P_{tr}^{CK}$ is the
trapping probability for a particle that starts from the surface of the
sphere. The Collins-Kimball rate constant is given by

$k_{CK}=\frac{k_{Sm}k_{tr}}{k_{Sm}+k_{tr}}=k_{Sm}P_{tr}^{CK},$ (3)

where $k_{tr}=4\pi R^{2}\kappa_{tr}$ is the rate constant in the trapping-
controlled regime (when $D\rightarrow\infty$), with $\kappa_{tr}$ denoting the
surface trapping rate ($\kappa_{tr}=0$ and $\infty$ for perfectly reflecting
and absorbing surfaces, respectively). As follows from Eqs. (1), (2) and (3),
the trapping probability, defined as
$P_{tr}^{CK}=J_{CK}/J_{Sm}=k_{CK}/k_{Sm}$, is equal to

$P_{tr}^{CK}=\frac{k_{tr}}{k_{Sm}+k_{tr}}=\frac{1}{1+D/(R\kappa_{tr})}.$ (4)

As $D/(R\kappa_{tr})\rightarrow 0$, the trapping probability tends to unity
and $k_{CK}$ approaches $k_{Sm}$. In the opposite limit when
$D/(R\kappa_{tr})\rightarrow\infty$, the probability $P_{tr}^{CK}$ tends to
zero as $R\kappa_{tr}/D$ and $k_{CK}$ approaches $k_{tr}$.

Figure 1: A schematic view of dynamic processes during the particle
association to an active site located on a partially absorbing spherical
surface. Details are explained in the text.

However, in a more realistic situation when the trapping ability is not
uniform over the surface of the sphere (Fig. 1), we can introduce an effective
trapping probability $P_{tr}$ and an effective rate constant $k$, such that
the steady-state flux of trapped particles can be written as

$J_{ss}=J_{Sm}P_{tr}=kc_{\infty},\quad k=k_{Sm}P_{tr}.$ (5)

By analogy with Eq. (4), we can use $P_{tr}$ to introduce an effective uniform
trapping rate of the surface, $\kappa_{eff}$, by the following relation,

$P_{tr}=\frac{J_{ss}}{J_{Sm}}=\frac{k}{k_{Sm}}=\frac{1}{1+D/(R\kappa_{eff})},$
(6)

which leads to

$\kappa_{eff}=\frac{DP_{tr}}{R(1-P_{tr})}.$ (7)

This allows us to map the problem of nonuniform trapping by a spherical
surface to an effective Collins-Kimball problem of trapping by the uniform
partially absorbing sphere characterized by the effective trapping rate
$\kappa_{eff}$. To determine $\kappa_{eff}$, we have to find the steady-state
flux, from which one can calculate the effective rate constant and trapping
probability using the relations given in Eq. (5) and eventually the expression
for $\kappa_{eff}$ by Eq. (7).

Applying the "boundary homogenization" approach outlined above to trapping of
diffusing particles by a small absorbing disk of radius $a$ located on the
surface of otherwise reflecting sphere of radius $R$ ($a\ll R$: see Fig. 1),
we can write the steady-state flux as

$J_{ss}=k_{HBP}c_{\infty}=J_{Sm}P_{tr}^{HBP}.$ (8)

Here $k_{HBP}$ is the Hill-Berg-Purcell rate constant that describes trapping
of diffusing particles by an absorbing disk of radius $a$ located on the
otherwise reflecting flat wall,30, 31

$k_{HBP}=4Da=k_{Sm}P_{tr}^{HBP},$ (9)

and $P_{tr}^{HBP}=J_{ss}/J_{Sm}=k_{HBP}/k_{Sm}$ is the probability to be
trapped by the disk for a particle whose starting position is uniformly
distributed over the surface of the sphere, including the disk,

$P_{tr}^{HBP}=\frac{a}{\pi R}.$ (10)

Using this, we can introduce the effective trapping rate $\kappa_{eff}^{HBP}$
which, according to Eq. (7), is given by

$\kappa_{eff}^{HBP}=\frac{DP_{tr}^{HBP}}{R(1-P_{tr}^{HBP})}=\frac{Da}{R(\pi
R-a)}\simeq\frac{Da}{\pi R^{2}}$ (11)

as $a\ll R$.

Now we generalize the above considerations to a more realistic case where
diffusing particles can reversibly bind to the surface of the sphere outside
the disk and diffuse on this surface with a diffusivity $D_{s}$ (Fig. 1). This
allows particles to reach the disk not only via direct association from the
bulk solution but also by the surface diffusion, leading to the overall
increase in the steady-state flux and hence the effective rate constant. The
steady-state flux of particles trapped by the disk is the sum of the fluxes
coming to the disk from the bulk and from the surface, denoted by
$J_{ss}^{(b)}$ and $J_{ss}^{(s)}$, respectively. Since the disk is small,
$a\ll R$, its presence affects the particle concentration only in close
vicinity of the disk. With this in mind,, we assume that the flux
$J_{ss}^{(b)}$ can be well approximated as

$J_{ss}^{(b)}=k_{HBP}c(R),$ (12)

where $c(R)$ is the steady-state concentration of unbound particles near the
surface of the sphere.

In addition, we assume that the flux $J_{ss}^{(s)}$ is the product of the mean
number $N_{s}$ of bound particles diffusing on the surface at the stationary
state and the rate constant $k_{s}$ that describes the trapping of these
particles by the disk,

$J_{ss}^{(s)}=k_{s}N_{s}.$ (13)

Since the disk is small, $k_{s}$ is the inverse mean fist-passage time
$\tau_{s}$ of a particle diffusing on the surface to the disk boundary,
conditional on that the particle starting point is uniformly distributed over
the surface. This mean first-passage time is given by32

$\tau_{s}=\frac{R^{2}}{D_{s}}f(a/R),\quad f(z)=2\ln(2/z)-1.$ (14)

Thus, the rate constant $k_{s}$ is

$k_{s}=\frac{1}{\tau_{s}}=\frac{D_{s}}{R^{2}f(a/R)}.$ (15)

The steady-state number of particles diffusing on the surface is proportional
to the stationary concentration $c(R)$ of unbound particles near the surface.
Let $k_{b}$ and $k_{d}$ be the rate constants for the particle binding to and
dissociation from the surface, respectively, as illustrated in Fig. 1. We use
these rate constants and the trapping rate constant $k_{s}$ to write the
balance equation for the steady-state number of particles diffusing on the
surface,

$(k_{d}+k_{s})N_{s}=k_{b}c(R).$ (16)

Thus, we have

$N_{s}=\frac{k_{b}}{k_{d}+k_{s}}c(R).$ (17)

This allows us to write the flux $J_{ss}^{(s)}$, Eq. (13), in terms of the
concentration $c(R)$,

$J_{ss}^{(s)}=k_{s}N_{s}=\frac{k_{s}k_{b}}{k_{d}+k_{s}}c(R).$ (18)

Summing up the two contributions to the steady-state flux of the particles
trapped by the disk, we obtain

$J_{ss}=J_{ss}^{(b)}+J_{ss}^{(s)}=\left(k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}\right)c(R).$
(19)

This steady-state flux is maintained by the bulk diffusion of particles to the
sphere from infinity. Assuming that the angular anisotropy of the particle
concentration in the bulk is small and, therefore, can be neglected, and
denoting this concentration at distance $r$ from the center of the sphere
($r>R$) by $c(r)$, one can write the steady-state flux in the bulk as

$J_{ss}=4\pi r^{2}\frac{dc(r)}{dr},\quad r>R.$ (20)

Solving this equation subject to the boundary condition
$c(r)|_{r\rightarrow\infty}=c_{\infty}$, we find that

$c(r)=c_{\infty}-\frac{J_{ss}}{4\pi Dr}.$ (21)

This expression is now used to find the concentration $c(R)$ of unbound
particles near the surface of the sphere:

$c(R)=c_{\infty}-\frac{J_{ss}}{4\pi DR}=c_{\infty}-\frac{J_{ss}}{k_{Sm}}.$
(22)

Substituting the above expression for $c(R)$ into Eq. (19), we arrive at

$J_{ss}=\left(k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}\right)\left(c_{\infty}-\frac{J_{ss}}{k_{Sm}}\right).$
(23)

Eventually, we determine the steady-state flux by solving this equation. This
leads to

$J_{ss}=\frac{k_{Sm}\left(k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}\right)}{k_{Sm}+k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}}c_{\infty}=J_{Sm}\frac{k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}}{k_{Sm}+k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}}.$
(24)

Now we take advantage of the relations in Eq. (5) to find the effective rate
constant and trapping probability:

$k=\frac{J_{ss}}{c_{\infty}}=\frac{k_{Sm}\left(k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}\right)}{k_{Sm}+k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}},$
(25)

and

$P_{tr}=\frac{J_{ss}}{J_{Sm}}=\frac{k}{k_{Sm}}=\frac{k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}}{k_{Sm}+k_{HBP}+\frac{k_{s}k_{b}}{k_{d}+k_{s}}}.$
(26)

Finally, we find the effective trapping rate from Eq. (7)

$\kappa_{eff}=\kappa_{eff}^{HBP}+\frac{1}{4\pi
R^{2}}\frac{k_{s}k_{b}}{k_{d}+k_{s}}.$ (27)

A particle diffusing on the surface of the sphere either dissociates and
becomes unbound or is trapped by the disk. The probability of the second
outcome, denoted by $P_{tr}^{(s)}$, is given by

$P_{tr}^{(s)}=\frac{k_{s}}{k_{s}+k_{d}}.$ (28)

Introducing the notation $\kappa_{s}$ for the binding rate of unbound
particles to the surface of the sphere outside the absorbing disk,

$\kappa_{s}=\frac{k_{b}}{4\pi R^{2}},$ (29)

we can rewrite the effective trapping rate, Eq. (27), as

$\kappa_{eff}=\kappa_{eff}^{HBP}+\kappa_{s}P_{tr}^{(s)}.$ (30)

We use the trapping probability $P_{tr}^{(s)}$ to write Eqs. (25) and (26) in
the form convenient for further analysis:

$k=\frac{k_{Sm}\left(k_{HBP}+k_{b}P_{tr}^{(s)}\right)}{k_{Sm}+k_{HBP}+k_{b}P_{tr}^{(s)}},$
(31)

and

$P_{tr}=\frac{k_{HBP}+k_{b}P_{tr}^{(s)}}{k_{Sm}+k_{HBP}+k_{b}P_{tr}^{(s)}}.$
(32)

To summarize, main results of our analysis are the expressions for the
effective rate constant, Eq. (31), the trapping probability, Eq. (32), and the
trapping rate, Eq. (30).

The effective rate constant $k$ is a function of the geometric parameters $a$
and $R$ as well as the parameters that characterize the dynamics of the
system: $D$, $D_{s}$, $k_{b}$ (or $\kappa_{s}$), and $k_{d}$. From Eqs. (31)
and (32), one can see that the effect of particle binding to the surface on
the effective rate constant and trapping probability (at fixed values of $a$,
$R$ and $D$) is determined by the term
$k_{b}P_{tr}^{(s)}=k_{b}k_{s}/(k_{d}+k_{s})$. The rate constant $k$ is a
monotonically increasing function of this parameter. As $k_{b}P_{tr}^{(s)}$
increase from zero (no effect of binding) to infinity (very strong binding),
the effective rate constant grows from its minimum value
$k_{min}=k_{Sm}k_{HBP}/(k_{Sm}+k_{HBP})\simeq k_{HBP}$ to its maximum value
$k_{max}=k_{Sm}$. Correspondingly, the effective trapping probability
increases from $k_{HBP}/k_{Sm}$ to unity.

Note that the effective rate constant in Eq. (31) approaches zero, as
$a\rightarrow 0$, much more slowly than $k_{HBP}$ given in Eq. (9). This is a
consequence of the fact that the surface diffusion of the particles to the
disk plays a dominant role in this limiting case. According to Eq. (31), we
have

$k|_{a\rightarrow
0}\simeq\frac{k_{b}k_{s}}{k_{d}}=\frac{k_{b}D_{s}}{k_{d}R^{2}\left[2\ln(2R/a)-1\right]},$
(33)

whereas $k_{HBP}=4Da$ is proportional to the disk radius.

## Numerical Tests of the Theory

To test our analytical theory for a wide range of parameters, we employ
Brownian Dynamics (BD) simulations. More specifically, our goal is to compute
the mean lifetime of a Brownian particle diffusing in a spherical layer
between two concentric spheres. The outer sphere of radius $R_{out}$ is a
reflecting boundary for the particle. The inner sphere of radius $R$
($R<R_{out}$) contains a small absorbing disk of radius $a$ ($a\ll R$). The
particle can reversibly bind to the surface of the inner sphere outside the
disk, and we assume that it diffuses on the surface with the same diffusivity
$D$ as in the bulk. Thus the particle can reach the disk and be trapped either
coming to the disk from the bulk or by surface diffusion. At $t=0$, the
particle is not bound to the surface and its starting point is uniformly
distributed over the surface of the inner sphere, including the disk area.

Figure 2: Comparison between numerical simulations and theoretical predictions
for mean lifetimes of Brownian particles for $R_{out}=2,3,4$ and $6$.

The system is characterized by 6 dimensional parameters such as $a$, $R$,
$R_{out}$, $D$, $k_{d}$, and $\kappa_{s}$. It seems convenient then to use $R$
as the unit of length and $R^{2}/D$ as the unit of time. So in our computer
simulations we choose $R=D=1$. The mean particle lifetime $\tau_{sim}$
obtained from the computer simulations can be compared with its counterpart
predicted by our analytical theory, $\tau_{theory}$. Then the measure of the
success of our theoretical analysis is how close are these two times, i.e., we
need to calculate the ratio $\tau_{theory}/\tau_{sim}$.

The main idea of our approach is to replace the non-uniform surface of the
inner sphere by homogeneous partially absorbing surface with the effective
trapping rate $\kappa_{eff}$. The mean lifetime of a particle diffusing
between partially absorbing and fully reflecting surfaces of radii $R$ and
$R_{out}$, respectively, that starts from the partially absorbing surface, is
given by33

$\tau_{theory}=\frac{1}{3\kappa_{eff}R}(R_{out}^{3}-R^{3}).$ (34)

This time can be rewritten in the dimensionless form as

$\widetilde{\tau}_{theory}=\tau\frac{D}{R^{2}}=\frac{1}{3\widetilde{\kappa}_{eff}}\left[\left(\frac{R_{out}}{R}\right)^{3}-1\right],$
(35)

where $\widetilde{\kappa}_{eff}$ is the dimensionless effective trapping rate
of the surface given by

${\widetilde{\kappa}}_{eff}=\kappa_{eff}\frac{R}{D}.$ (36)

For numerical tests, it is convenient to present the dimensionless effective
trapping rate as

$\widetilde{\kappa}_{eff}=\frac{a}{\pi
R-a}+\frac{\varepsilon_{2}}{1+\varepsilon_{3}},$ (37)

where new dimensionless parameters are given by

$\varepsilon_{2}=\frac{k_{b}}{k_{Sm}}=\frac{\kappa_{s}R}{D},\quad\varepsilon_{3}=\frac{k_{d}}{k_{s}}=k_{d}\frac{R^{2}}{D}\left[2\ln(2R/a)-1\right].$
(38)

Our idea is to test analytical predictions in computer simulations by varying
$\varepsilon_{2}$ and $\varepsilon_{3}$ since they contain all relevant
parameters of the process.

The comparison between theoretical predictions and the results of computer
simulations for the surface-assisted binding process are presented in Fig. 2.
One can see that pushing the outer boundary away from the surface (increasing
$R_{out}$) significantly improves the agreement between the theory and
simulations. This is an expected result because the theoretical arguments for
substituting the originally heterogeneous surface with an effective
homogeneous one are applicable for $R_{out}/R\gg 1$. One should also notice
that theory works quite well for all ranges of parameters already for small
$\varepsilon_{2}$ and large $\varepsilon_{3}$ that corresponds to weak
bindings to the surface and frequent dissociations from the surface. In this
case, the surface plays a relatively small role in the overall dynamics. The
opposite limit of large $\varepsilon_{2}$ and small $\varepsilon_{3}$ presents
the situation where the presence of the surface is crucial. In this case, the
agreement between theory and simulations is not as good, although increasing
$R_{out}$ clearly shows the right tendency.

To emphasize the important role of the surface in our system, we also consider
an effective acceleration parameter $\gamma$ which is defined as a ratio of
the effective trapping rate constant $k$, Eq. (31), to the Hill-Berg-Purcell
rate constant $k_{HBP}$, Eq. (9), that corresponds to the fully inert sphere.
This ratio provides a quantitative measure of how the presence of the surface
influences the overall trapping rate by the specific site. It can be shown
that

$\gamma=\frac{k}{k_{HBP}}=\frac{1+\frac{\pi
R}{a}\frac{\epsilon_{2}}{1+\epsilon_{3}}}{1+\frac{a}{\pi
R}+\frac{\epsilon_{2}}{1+\epsilon_{3}}}.$ (39)

Fig. 3 presents this acceleration as a function of two dimensionless
parameters $\epsilon_{2}$ and $\epsilon_{3}$ for the system with $a=0.05$. As
given in Eqs. (38), the parameter $\epsilon_{2}$ is proportional to the
absorption constant ($\kappa_{s}$) and the parameter $\epsilon_{3}$ is
proportional to the desorption rate constant ($k_{d}$). From the Fig. 3, one
can see that the highest values of the acceleration are achieved for high
values of $\epsilon_{2}$ and low values of $\epsilon_{3}$. The maximum
theoretically possible acceleration can be simply estimated as

$\gamma_{max}\simeq\frac{k_{Sm}}{k_{H}}=\frac{\pi R}{a}.$ (40)

This means that the surface effect is most pronounced when the active site is
very small and the surface binding is strong.

Figure 3: Acceleration of binding due to presence of nonspecific interactions
for $R=1$, $a=0.05$ as a function of dimensionless parameters $\epsilon_{2}$
and $\epsilon_{3}$. Maximum theoretical acceleration for this system is
$\gamma\sim 62.8$.

## Summary and Conclusions

We developed a new theoretical framework to describe trapping of particles by
an active site located on the surface. It generalizes a Smoluchovskii-Collins-
Kimball approach by taking into account nonspecific interactions between the
particles and the surface where the active site is located. Our main idea is
to replace the inhomogeneous surface by an effective homogeneous one. This
allowed us to obtain a full description of the processes in the system.
Importantly, we present an explicit analytical expression for the effective
trapping rate. Our approximate theoretical method gives correct predictions in
the limiting cases of strong binding and no binding to the surface. It was
also tested using extensive BD computer simulations. The analysis suggests
that nonspecific interactions can significantly accelerate the trapping
process in the systems with small active sites.

## Acknowledgments

MM and ABK acknowledge the support from the Welch Foundation (C-1559), from
the NSF (CHE-1953453 and MCB-1941106), and from the Center for Theoretical
Biological Physics sponsored by the NSF (PHY-2019745). AMB and SMB were
supported by the Intramural Research Program of the National Institutes of
Health, Eunice Kennedy Shriver National Institute of Child Health and Human
Development, and the Center for Information Technology.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Somorjai and Li 2010 Somorjai, G. A.; Li, Y. _Introduction to surface chemistry and catalysis_ ; John Wiley & Sons, 2010
  * Copeland 2000 Copeland, R. A. _Enzymes: a practical introduction to structure, mechanism, and data analysis_ ; John Wiley & Sons, 2000
  * Nerukh et al. 2012 Nerukh, D.; Okimoto, N.; Suenaga, A.; Taiji, M. Ligand diffusion on protein surface observed in molecular dynamics simulation. _The journal of physical chemistry letters_ 2012, _3_ , 3476–3479
  * Rostovtseva et al. 2015 Rostovtseva, T. K.; Gurnev, P. A.; Protchenko, O.; Hoogerheide, D. P.; Yap, T. L.; Philpott, C. C.; Lee, J. C.; Bezrukov, S. M. $\alpha$-Synuclein shows high affinity interaction with voltage-dependent anion channel, suggesting mechanisms of mitochondrial regulation and toxicity in Parkinson disease. _Journal of Biological Chemistry_ 2015, _290_ , 18467–18477
  * Hoogerheide et al. 2018 Hoogerheide, D. P.; Gurnev, P. A.; Rostovtseva, T. K.; Bezrukov, S. M. Real-time nanopore-based recognition of protein translocation success. _Biophysical journal_ 2018, _114_ , 772–776
  * Jacobs et al. 2019 Jacobs, D.; Hoogerheide, D. P.; Rovini, A.; Jiang, Z.; Lee, J. C.; Rostovtseva, T. K.; Bezrukov, S. M. Probing membrane association of $\alpha$-synuclein domains with VDAC nanopore reveals unexpected binding pattern. _Scientific reports_ 2019, _9_ , 1–14
  * Hoogerheide et al. 2020 Hoogerheide, D. P.; Rostovtseva, T. K.; Jacobs, D.; Gurnev, P. A.; Bezrukov, S. M. Tunable Electromechanical Nanopore Trap Reveals Populations of Peripheral Membrane Protein Binding Conformations. _ACS nano_ 2020,
  * Gurnev et al. 2014 Gurnev, P. A.; Yap, T. L.; Pfefferkorn, C. M.; Rostovtseva, T. K.; Berezhkovskii, A. M.; Lee, J. C.; Parsegian, V. A.; Bezrukov, S. M. Alpha-synuclein lipid-dependent membrane binding and translocation through the $\alpha$-hemolysin channel. _Biophysical journal_ 2014, _106_ , 556–565
  * Song et al. 1996 Song, L.; Hobaugh, M. R.; Shustak, C.; Cheley, S.; Bayley, H.; Gouaux, J. E. Structure of staphylococcal $\alpha$-hemolysin, a heptameric transmembrane pore. _Science_ 1996, _274_ , 1859–1865
  * Sobhy 2017 Sobhy, H. A comparative review of viral entry and attachment during large and giant dsDNA virus infections. _Archives of virology_ 2017, _162_ , 3567–3585
  * Rothenberg et al. 2011 Rothenberg, E.; Sepúlveda, L. A.; Skinner, S. O.; Zeng, L.; Selvin, P. R.; Golding, I. Single-virus tracking reveals a spatial receptor-dependent search mechanism. _Biophysical journal_ 2011, _100_ , 2875–2882
  * Marsh and Helenius 2006 Marsh, M.; Helenius, A. Virus entry: open sesame. _Cell_ 2006, _124_ , 729–740
  * Gibbons et al. 2010 Gibbons, M. M.; Chou, T.; D’orsogna, M. R. Diffusion-dependent mechanisms of receptor engagement and viral entry. _The Journal of Physical Chemistry B_ 2010, _114_ , 15403–15412
  * Berg 1985 Berg, O. Orientation constraints in diffusion-limited macromolecular association. The role of surface diffusion as a rate-enhancing mechanism. _Biophysical journal_ 1985, _47_ , 1–14
  * Zhou and Szabo 2004 Zhou, H.-X.; Szabo, A. Enhancement of association rates by nonspecific binding to DNA and cell membranes. _Physical review letters_ 2004, _93_ , 178101
  * Zhou 2010 Zhou, H.-X. Rate theories for biologists. _Quarterly reviews of biophysics_ 2010, _43_ , 219
  * Shin and Kolomeisky 2018 Shin, J.; Kolomeisky, A. B. Surface-assisted dynamic search processes. _The Journal of Physical Chemistry B_ 2018, _122_ , 2243–2250
  * Mirny et al. 2009 Mirny, L.; Slutsky, M.; Wunderlich, Z.; Tafvizi, A.; Leith, J.; Kosmrlj, A. How a protein searches for its site on DNA: the mechanism of facilitated diffusion. _Journal of Physics A: Mathematical and Theoretical_ 2009, _42_ , 434013
  * Kolomeisky 2011 Kolomeisky, A. B. Physics of protein–DNA interactions: mechanisms of facilitated target search. _Physical Chemistry Chemical Physics_ 2011, _13_ , 2088–2095
  * Sheinman et al. 2012 Sheinman, M.; Bénichou, O.; Kafri, Y.; Voituriez, R. Classes of fast and specific search mechanisms for proteins on DNA. _Reports on Progress in Physics_ 2012, _75_ , 026601
  * Koslover et al. 2011 Koslover, E. F.; de la Rosa, M. A. D.; Spakowitz, A. J. Theoretical and computational modeling of target-site search kinetics in vitro and in vivo. _Biophysical journal_ 2011, _101_ , 856–865
  * Bauer and Metzler 2013 Bauer, M.; Metzler, R. In vivo facilitated diffusion model. _PloS one_ 2013, _8_ , e53956
  * Shvets et al. 2018 Shvets, A. A.; Kochugaeva, M. P.; Kolomeisky, A. B. Mechanisms of protein search for targets on DNA: theoretical insights. _Molecules_ 2018, _23_ , 2106
  * Kochugaeva et al. 2016 Kochugaeva, M. P.; Shvets, A. A.; Kolomeisky, A. B. How conformational dynamics influences the protein search for targets on DNA. _Journal of Physics A: Mathematical and Theoretical_ 2016, _49_ , 444004
  * Felipe et al. 2021 Felipe, C.; Shin, J.; Kolomeisky, A. B. DNA Looping and DNA Conformational Fluctuations Can Accelerate Protein Target Search. _The Journal of Physical Chemistry B_ 2021, _125_ , 1727–1734
  * Cencini and Pigolotti 2018 Cencini, M.; Pigolotti, S. Energetic funnel facilitates facilitated diffusion. _Nucleic acids research_ 2018, _46_ , 558–567
  * Adam et al. 1968 Adam, G.; Delbrück, M.; Rich, A.; Davidson, N. Structural chemistry and molecular biology. _San Francisco, USA_ 1968,
  * Smoluchowski 1916 Smoluchowski, M. v. Drei vortrage uber diffusion, brownsche bewegung und koagulation von kolloidteilchen. _Zeitschrift fur Physik_ 1916, _17_ , 557–585
  * Collins and Kimball 1949 Collins, F. C.; Kimball, G. E. Diffusion-Controlled Reaction Rates. _J. Colloid Sci._ 1949, _4_ , 425–437
  * Hill 1975 Hill, T. L. Effect of Rotation on the Diffusion-Controlled Rate of Ligand-Protein Association. _Proc. Natl. Acad. Sci. U. S. A._ 1975, _72_ , 4918–4922
  * Berg and Purcell 1977 Berg, H. C.; Purcell, E. M. Physics of Chemoreception. _Biophys. J._ 1977, _20_ , 193–219
  * Bloomfield and Prager 1979 Bloomfield, V.; Prager, S. Diffusion-controlled reactions on spherical surfaces. Application to bacteriophage tail fiber attachment. _Biophysical journal_ 1979, _27_ , 447–453
  * Dagdug et al. 2016 Dagdug, L.; Vázquez, M.-V.; Berezhkovskii, A. M.; Zitserman, V. Y. Boundary homogenization for a sphere with an absorbing cap of arbitrary size. _The Journal of chemical physics_ 2016, _145_ , 214101