key: cord-0901622-0tqp4ar5
authors: Chen, Jixin
title: Simulating stochastic adsorption of diluted solute molecules at interfaces
date: 2022-01-11
journal: AIP Adv
DOI: 10.1063/5.0064140
sha: 1cf6721d40bf6850736bb7436baeddfd3de423d0
doc_id: 901622
cord_uid: 0tqp4ar5

This report uses Monte Carlo simulations to connect stochastic single-molecule and ensemble surface adsorption of molecules from dilute solutions. Monte Carlo simulations often use a fundamental time resolution to simulate each discrete step for each molecule. The adsorption rate obtained from such a simulation surprisingly contains an error compared to the results obtained from the traditional method. Simulating adsorption kinetics is interesting in many processes, such as mass transportation within cells, the kinetics of drug–receptor interactions, membrane filtration, and other general reaction kinetics in diluted solutions. Thus, it is important to understand the origin of the disagreement and find a way to correct the results. This report reviews the traditional model, explains the single-molecule simulations, and introduces a method to correct the results of adsorption rate. For example, one can bin finer time steps into time steps of interest to simulate the fractal diffusion or simply introduce a correction factor for the simulations. Then two model systems, self-assembled monolayer (SAM) and biosensing on the patterned surface, are simulated to check the accuracy of the equations. It is found that the adsorption rate of SAM is highly dependent on the conditions and the uncertainty is large. However, the biosensing system is relatively accurate. This is because the concentration gradient near the interface varies significantly with reaction conditions for SAMs while relatively stable for the biosensing system.

Molecules in diluted solution are moving constantly and randomly (Fig. 1) . When still, a major contributor of this motion, diffusion, is fundamentally relevant in many fields such as cell biology, biosensing, separation, fluidic dynamics, reaction kinetics, catalysis, and batteries. [1] [2] [3] [4] [5] [6] [7] [8] [9] The ensemble kinetics of diffusion has been summarized by Fick's laws of diffusion. 2, 3 For example, the diffusion of materials in a diluted solution from a high-concentration reservoir into a tubing space forms a time-and space-dependent concentration gradient function C(x, t) 2,3 derived from the heat equation formulated by Joseph Fourier in 1822, 3

where t is time (unit s), x is the distance away from the interface (m), C(x, t) is the concentration gradient (molecules m −3 ), and D(x, t) is the time-and-space-dependent diffusion coefficient (m 2 s −1 ), a constant under many conditions. Note that these units are for the 1D simplification of a 3D diffusion when the diffusions of the other two dimensions do not matter. For a real 1D diffusion, C will have the unit molecules m −1 . The special case with constant D has a formal analytical solution, 2

where C(x, t) (unit mol/m 3 ) is the concentration of the molecules at space x (m) from origin x = 0 and time t from initial time t = 0 when all molecules are at the origin. Equation ( 2) is a 1D Gaussian distribution function with the standard deviation σ = √ 2Dt and is normalized to the unit in all x spaces, i.e., at any giving snapshot of t, A ∫ this concentration profile represents the probability density function (PDF) of a single molecule diffusing from the origin into space over time. It is a very important achievement in history that Stokes and Einstein come up with a single-particle random walking model to predict the diffusion constant of particles doing Brownian motion. 10 The diffusion constant D for molecules, colloids, proteins, or much bigger particles in a solution can be estimated by the Stokes-Einstein equation, 11, 12 

where k B is the Boltzmann constant, T is temperature, η is the viscosity of the solution, and r is the radius of the particle. For a molecule approximated to a small ball, r (m) can be estimated from the molecular weight Mw = 4 3 πr 3 ρ (kg), where ρ (kg/m 3 ) is the density of the neat molecule in the solid or liquid state (all in SI units).

An interesting application of this diffusion theory is to predict the adsorption rate of molecules in a diluted solution to a solid surface. In 1937, Langmuir and Schaefer came up with an equation to predict the adsorption rate at the short-time limit (a continuous model). 13 Langmuir and Schaefer obtained adsorption kinetics by directly integrating Fick's second law equation at a surface, assuming to absorb any molecules that have "crossed" it. 13, 14 The time-dependent concentration gradient at the surface is

Integrating this equation over time, Langmuir and Schaefer gave the following equation: 13

where Γ(t) (unit mol) is the number of molecules adsorbed on an area of surface A (unit m 2 ) at time t (s), c b (mol m −3 ) is the concentration of the adsorbate in the bulk solution, and D is the diffusion constant (m 2 s −1 ). In 1946, Ward and Tordai added a back-diffusion term in the equation to account for the adsorption during the longer period ( Fig. 1 ), 14

where c(τ) is the sub-surface molar concentration of the adsorbate near the surface (mol m −3 ) and τ is a dummy variable with the unit of time (s). Because our current measuring techniques, especially various chemical imaging methods, have discrete integration times, 4, 6, 9, [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] finding the correlation between the continuous and discrete solutions of adsorption is important. The difficulty of understanding lies in the self-similar fractal nature of diffusion and its probability density function. From a finite-difference point of view, the broadening of the ensemble diffusion profile is a combination of the broadening of a smaller fraction of the earlier profiles (Fig. 2) . This creates a question on directly using the combined ensemble curves [Eq. (2) ] to calculate the adsorption probability of a molecule from the bulk to the surface. Non-Gaussian diffusion is often observed in a confined space or near a surface further complicated the problem. 7, 16 The adsorption of the molecules is greatly affected by the reabsorption of the reflected molecules. 8 It is becoming important for us to re-evaluate the theory of adsorption under ideal conditions to 

scitation.org/journal/adv reduce the difficulty in building or understanding more complicated models.

In this report, we are comparing this theory with the results of Monte Carlo simulations to revisit the correlation between the continuous theoretical models with the discrete simulation results. 25 An additional task of the comparison is to gain insights into a gap between the theory and experiments. 24, [26] [27] [28] [29] [30] [31] [32] The Ward-Tordai (WT) equation is widely used to measure the effective diffusion constant from the experimental data. 14, 32 However, the effective diffusion constant is often several orders of magnitude different from the value predicted by the Stoke-Einstein equation. 14,32 These complexities have limited the applications of the Ward-Tordai equation [Eq. (5)] in many fields, such as chemistry, biochemistry, biophysics, biotechnologies, and chemical engineering. 29, 30 For example, it is critical to predict how long it takes for a drug molecule to be adsorbed on the surface of a cell and how long it takes to diffuse to the target binding site inside a cell. It is also important to predict how long we should wait for a typical biosensing platform, such as a glucose sensor and a surface plasmon resonance sensor. Predicting the diffusion of a molecule to the surface is also essential to calculate the corrosion rate of a pipe and the reaction rates in a heterogeneous catalytic system.

In a typical measurement, we often monitor the accumulated molecule adsorbed on the surface with a time resolution τ as the integration time of each measurement step, e.g., a frame of a movie. Assuming the locations of the molecules are known at the beginning of each frame, whose diffusion probabilities will create broadening profiles described by Eq. (2) by replacing the t with ∆t = τ (Fig. 3 ). If we integrate the error functions of all these profiles in the bulk solution ( Fig. 3 and supplementary material), this will represent the adsorption simulated using the same discrete-time resolution,

Comparing Eq. (5) with Eq. (7), we can see that the discrete solution has a two-fold reduction in predicting the adsorption by ignoring the fractal diffusion that happens within ∆t. When we zoom in the time of the frame at when a molecule hits the surface, the probability function collapses into a binary state: it either adsorbs, which truncates the average PDF, or it reflects, which increase the probability of the adsorption in the next moment by the same amount as the truncated error function. We can call this doubling a "mirror effect," which analogies to Zeno's paradoxes. As such, the overall probability of adsorption is twice the sum of the error functions shown in Fig. 3 ,

Equation (5) assumes that the adsorption continues from frame to frame and the concentration gradient near the surface continuously evolves (the concentration decreases over time since the molecules are absorbed by the surface), while Eq. (8) assumes that the concentration gradient evolves during the single frame and then recovers to the original at the beginning of the next frame. Equation (8) can be used to predict the number of collisions of the molecules to a nonstick surface with care being taken for choosing the length of ∆t. The average rate of solute colliding the wall measured in ∆t is

Equation (9) is confusing that the rate of collision depends on the observation time. As such, we run Monte Carlo simulations to reproduce the result of Eq. (8) on a reflective surface and Eq. (5) on an adsorptive surface. Figure 4 shows Monte Carlo simulations of molecules moving in a 1D space that bounced/reflected from both ends. Equation (7) is reproduced from the Monte Carlo simulation using a single PDF. Equation (8) is reproduced by binning >1000 simulation steps to mimic the finer fractal diffusion by combining the PDFs with subframe resolution. The time-dependent collision rate is confirmed to cause by the repetitive collision of the same molecule in a measuring cycle if observed at a finer time resolution. A real mirror effect, doubling crossing probability, is also observed at an imaginary interface in the bulk. See the supplementary material for a discussion of Fig. 4 results. Figure 5 shows Monte Carlo simulations of molecules moving in a 1D space that bounced/reflected from one end and adsorbed on the other. Because one adsorption per simulated step has been satisfied in such a simulation, no binning is required to reproduce the adsorption kinetics predicted by Eq. (5). We can also visualize the evolution of the concentration gradient near the adsorptive surface and the bulk over time. This gradient is originally rationalized in developing the Langmuir-Schaefer and Ward-Tordai equations that the adsorption rate drops over time.

Comparing the two simulations in Figs We can do a mind experiment with molecules aligned perfectly in space shown in Fig. 6 . To maintain the same molecular distribution within the time Δt, the average location of the molecule should be the same as the distance between two molecules, i.e., the net effect is just switch positions during this time,

where d 0 and z are shown in Fig. 6 . Thus, the characteristic integration time ∆tc to calculate the average adsorption rate with no sub-surface concentration gradient is

Thus, the average initial adsorption rate can be calculated by Eq. (9), This equation has the correct unit s −1 for the dimensional analysis and is consistent with another calculation, assuming the molecular exchanging time ∆tc is the average adsorption time for the characteristic surface area d 0 2 . That is,

In short, ∆tc is the characteristic diffusion-limit time for the adsorption. Shorter than this time, the high average hitting rate and the low total number of adsorptions predicted by Eqs. (8) and (9) are due to the isolated but repetitive collision of the same molecules on the surface. Longer than this time, the same molecule might have diffused away beyond the first nearby layer and diffused back to the surface, which should have been counted as "different" molecules if the sub-surface concentration has been maintained, i.e., Eq. (9) has lower estimated the collision frequency longer than this time. Thus, at the very beginning of the adsorption, ∆tc is the time between the molecules that have just collided with the surface and diffused away to the first nearest neighbor layer.

Equation (13) derived from the single-molecule approach is consistent with the results from the ensemble method. Equation (13) predicts values that are consistent with the diffusion-limited steady-state flux of Pt nanoparticles on an ultramicroelectrode (UME). 33 When diameter of UME is r = 10 μm, c b = 25 pM, and ARTICLE scitation.org/journal/adv FIG. 6. Scheme of finding characteristic Δt to calculate the adsorption frequency when the sub-surface concentration is the same as the bulk concentration, i.e., the short-time limit right after a fresh surface is exposed to the solution. D = 1 × 10 −12 m 2 s −1 , the critical adsorption rate is calculated as 0.37 s −1 . This rate is consistent with the calculated value as 0.4 s −1 using the semiempirical steady-state flux equation. 33, 34 The 4/3-order dependence of rate on concentration is very weird, and I initially also think it is wrongly derived. For the adsorption, Eq. (9) already shows that it should be the first-order dependent on the concentrations; thus, the 4/3-order dependence in Eq. (13) is only for the initial rate at a short period when the symmetry is just broken by introducing the plane in the solution. It is unlikely to observe this dependence in a real adsorption experiment when a concentration gradient will develop. It is a useful initial assessment of the order of rate that is independent of the observation integration time, especially for experiments on a system with a small absorption area surrounded by a much larger inactive area and/or with a fast flow rate. 4, 6, 7, 9, 23, 25, [35] [36] [37] [38] However in the bulk solution, the 4/3-order dependence of the diluted concentration on the collision order is surprising. In the physical chemistry textbooks, first-order dependence for high concentration gas reaction is supported by the collision theory, where collision frequency is a first-order function of concentration. This fraction order is consistent with previous simulations on the fractal nature of reaction kinetics that obtains abnormal reaction orders for simple reactions. 39 In the literature, the Ward-Tordai equation often gives several orders of magnitude differences between the measured effective diffusion constant and that predicted by the Stokes-Einstein equation. Thus, it is interesting to simulate the adsorption of molecules on a surface using the typical Langmuir surface adsorption model. We are simulating two systems, self-assembled monolayers (SAMs) and binding of molecules to a patterned surface. SAMs are used in many fields for surface functionalization and the latter is often used in biosensing systems, where both need a simple equation to predict the adsorption kinetics.

The Langmuir adsorption of SAM without energy barriers is simulated using the Gillespie algorithm, a computer-based Monte Carlo simulation method. 40, 41 Specifically, the bulk solute molecules 

where A is the total surface area simulated, a is the average size of a binding site, t is time, and k is the effective rate constant on the surface. The effective rate constant can be pulled out from fitting and labeled as k 1 , and the fitted curve is labeled Γ 1 [ Fig. 7(b) , orange curve].

If an ideal stirring has been introduced that maintains the bulk concentration with no gradient at the front edge of the surface, an exponential decay curve Γ 2 with the same model as Γ 1 but with much faster rate constant k 2 than k 1 would have been observed [ Fig. 7(b) , red curve] with the rate constant,

Surprisingly, Eq. (5) is often considered short-time limited in the literature and overlaps with the simulated data most of the time until ∼80% of saturation [Figs. 7(b), Γ 3 ]. Equation (5) has not considered the shrinking active binding area of the SAMs (due to increasing adsorption coverage) and the slowing-down evolution speed of the concentration gradient (due to the back-diffusion) near the surface over time. Coincidentally, these two effects cancel with each other during most of the adsorption time.

The difference larger than three orders of magnitude between k 1 and k 2 [ Fig. 7(c) ] explains the over six orders of magnitude variations of the effective diffusion constant calculated from the Ward-Tordai equation on experimental data, which is still confusing nowadays in the literature. 14, 32 We can speculate an answer from the simulation results that because convection, flow, and stirring naturally occur under typical experimental conditions, greatly ARTICLE scitation.org/journal/adv changing the formation of the concentration gradient near the surface. We can simulate the initial slope predicted by Eq. (9) with a much higher consistency with experimental results than SAMs for a special case. We simulate a typical biosensing system on a patterned surface and/or under significant flow (Fig. 8) when the sub-surface concentration of the solution does not change much throughout adsorption. This is simulated with a 3D diffusion model (see the supplementary material for detailed derivation and simulation conditions). Figure 8(a) shows the simulation scheme with a 10-nm-radius adsorption spot (with 150 binding sites, each occupies ∼2 nm 2 ) on the wall of a 1000-μm 3 cubic volume. The rate constants of Fig. 8(b) are fitted with Eq. (14) and shown in Fig. 8(c) .

By replacing the ensemble area with the single-molecule binding area, Eq. (9) predicts rate constants, Twice the simulated binding constant k 1 [Eq. (14) ] is compared to the theoretical values k 3 [Eq. (16) ] [ Fig. 8(d) ]. Twice is chosen to correct the fractal diffusion being ignored during the simulation when the majority of the surface is reflective. We can conclude from the consistency between the two values that Eq. (9) holds for such simulations [ Fig. 8(d) ], suggesting that this equation is applicable in biosensing and drug delivery fields. These simulations are consistent with experimental results reported in the literature where single-molecule adsorption of dye molecules on immobilized DNA is measured 43 or protein-surface adsorption within an order of magnitude. 24 These simulations also suggest that the reverse argument is true, that is, for binding of diffusive solute molecules to isolated targets on a surface, Eq. (5) can be used to measure the diffusion constant of the solute, overcoming the several orders of magnitude variation observed in the measurements of SAMs.

In summary, the surface adsorption kinetics represented by the Langmuir-Schaefer equation and the Ward-Tordai equation has been reproduced using Monte Carlo simulations. The fractal nature of diffusion is examined and simulated. A unique conclusion suggested by the single-molecule approach is that there is a characteristic integration time ∆tc for the equations, which has never been a consideration for the ensemble approaches, a missing piece of the Ward-Tordai equation. That is, the characteristic time distinguishes the overestimated multiple collisions from the same molecule shorter than this time and the lower estimation of the collisions longer than this time. We can draw some interesting specific conclusions from the Monte Carlo simulations. For example, when we measure the collision events of probe molecules in the bulk solution to a small target area on a flat surface, the frequency of seeing such events is dependent on the measuring time resolution and how the surface reacts with the probes. With the results obtained from the simulations, we may start to use the simple continuous and discrete equations carefully in various fields, such as calculating the collision frequency of molecules in a diluted gas or liquid solution, membrane penetration, self-assembly, and biosensing. Wide use of these equations in these fields has not yet been seen in the literature. It may also find applications in finer simulations such as molecular dynamics simulations to skip non-interested mass transportation steps among the solvent. More careful simulations and experiments should be carried out in the future to further test the limit of these equations. Nevertheless, I hope you agree that the single-molecule approach pictured in this report is easier to understand than the ensemble equations of diffusion.

The Monte Carlo simulation is carried out on a laptop equipped with an Intel i7 CPU (2.2 GHz) and 16 Gb of memory. A basic version of MATLAB 2014b is used for all simulations. A single central processing unit (CPU) is used for all simulations. A previously coded exponential regression fitting algorithm fitNguess 44 (Github) is used to fit the curves. The two major random functions to generate the step motion of each molecule are from MATLAB, rand() ARTICLE scitation.org/journal/adv creating evenly distributed random numbers, and randn() generating Gaussian distributed random numbers. Detail parameters and settings for the simulations are listed in the supporting information.

The single-molecule diffusion of the solute molecules in confined semi-infinite volumes and their collisions to the walls are simulated.

Inter-molecular collisions are not simulated. Typical computer time used for the simulations is from a few seconds to a few hours. See the last sections of the supplementary material for two example source codes in MATLAB.

See the supplementary material for simulation details, derivations, and example MATLAB source codes for 1D and 3D simulations. Other source codes and data are available upon request.

The author thanks his family and institute for their support during the COVID-19 pandemic. The research reported in this publication was supported by the National Human Genome Research Institute of the National Institutes of Health under Award No. R15HG009972. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Durability of Concrete in Saline Environment

Brownian Motion

Electrochemical Methods Fundamentals and Applications

Diffusion in Nanoporous Materials

The author has no conflicts to disclose.

MATLAB source codes have been included in the supplementary material. Additional codes and data that support the findings of this study are available from the corresponding author upon reasonable request.