We model the transport of electrically charged solute molecules by a laminar flow within a nanoslit microfluidic channel with electrostatic surface potential. We derive the governing convection–diffusion equation, solve it numerically, and compare it with a Taylor–Aris-like approximation, which gives excellent results for small Péclet numbers. We discuss our results in light of designing an assay that can measure simultaneously the hydrodynamic size and electric charge of single molecules by tracking their motion in such nanoslit channels with electrostatic surface potential.

## I. INTRODUCTION

Structure and electrical charge are fundamental characteristics of biological molecules. Electrostatic effects are ubiquitous and dominate the activities and properties of many biomolecules. The main driving forces behind protein folding are the burying hydrophobic amino acid residues in the interior, the formation of salt bridges, and the formation of favorable energetic contacts of polar groups on the surface with solvent molecules. Ionic or ionizable groups on the surface of biomolecules play crucial roles in fundamental biological processes ranging from ligand recognition, signaling, protein folding, aggregation, enzymatic catalysis to redox reactions.^{1–9} Metamorphic proteins and nucleic acid chains exhibit more than one folded active conformation in native conditions, often with different functions.^{10–12} To investigate the structure–function relation of biomolecules, one requires knowledge of the full energy landscape and the effective charge on multiple pre-existing conformational isomers. At the present time, established techniques such as x-ray crystallography, small-angle x-ray scattering (SAXS), small-angle neutron scattering (SANS), nuclear magnetic resonance spectroscopy (NMR), and cryo-electron microscopy (cryoEM) are capable of delivering detailed atomic information on biomolecules (see, for example, Ref. 13). As a result, theoretical structure-based electrostatic calculations are now routinely performed and used to dissect the structure–function inter-relationships and properties of these biomolecules. However, computational methods that perform structure-based electrostatic calculations face challenges when simulating specific biomolecule–ion interactions or biomolecule–biomolecule interactions in cellular settings where they are immersed in a solution of numerous ions and other small molecules and macromolecules. To handle important effects such as preferential hydration, ion–ion competition, charge regulation, or solvation is computationally expensive, which severely limits the temporal and spatial scales of these calculations, often leading to inaccurate results.^{14} Although recent experimental advances are beginning to explore the energy landscapes of proteins and nucleic acid chains,^{15,16} measuring electric charges of biological molecules in ion-containing buffer solutions at single molecule level is far from trivial.^{17–20}

Solid-state nanopores have been used to characterize the charge and size of biomolecules, and they offer great potential for scalable high-throughput rapid protein and nucleic acid analysis. The core principle here is to measure the dwell times of ionic current blockage while molecules (proteins and nucleic acids) chain through the nanopores. However, this method is limited by the passage time of a single amino or nucleic acid, which is typically on the order of 10 ns to 1 *µ*s, poor signal-to-noise ratio due to thermal noise, interactions with the nanopore walls, and the detrimental effect of the applied high voltages on the protein or oligonucleotide structure itself.^{21,22} Recently, a high-precision method to measure the electrical charge on single nanometer-scale objects using a geometry-induced electrostatic fluidic trap has been introduced.^{23,24} The chief principle of this method is to modulate the spatial interaction potential of a charged object by introducing topographical indentations onto a sub-micrometer channel with like-charged walls and then to stably confine the particles in local interaction potential minima. By measuring the escape rates of fluorescently labeled biomolecules in such traps, which scales exponentially with the potential difference of having a particle inside or outside the trap, their effective electrical charges were back-calculated with $<1$ *e* precision. Although the method has promising potential, there are a few major limitations that impact its accuracy in determining the electrical charges. These are the determination of surface charge properties of the silica walls that remain poorly resolved to date^{25–28} and the required *a priori* knowledge of the molecules’ hydrodynamic radii that necessitates independent measurements. Furthermore, the uncertainty of determining the charge is often defined and limited by the accuracy of the *in silico* numerical simulations of the entire experiment and by the accuracy of the complex fabrication process.

In the present article, we propose a new method for determining *simultaneously* the size and the *effective* electrical charge of fluorescently labeled biomolecules in a nanochannel with thin conductive surfaces. The core idea is to measure the charge-dependent drift velocity and the dispersion of biomolecules under lateral pressure-driven parabolic laminar flow inside a nanochannel (or capillary) having a constant electrostatic surface potential. For sufficiently small Péclet numbers (ratio between characteristic diffusive and convective transport times), the broadening of the concentration peaks can be used for measuring molecular diffusion coefficients. This broadening is perfectly described by the theory developed by Taylor and Aris, who derived an expression for the apparent diffusion constant of the concentration spreading as a function of the actual solute’s diffusion constant, the channel size, and the mean flow velocity.^{29–31} In this work, we consider a similar problem, the convective–diffusive transport of a solute along a planar nanoslit channel but with channel walls having non-zero surface potential and for charged solute molecules. The electrostatic interaction with the channel walls and the Debye layer around them alters the transversal distribution of the charged species across the channel. When a pressure difference is applied along the channel, the mean transport rate of molecules is modified significantly by these interactions in a charge-dependent manner. Although electrokinetic transport has been modeled extensively in the past,^{32–34} very few theoretical works have targeted the convective–diffusive transport of charged solute in a nanochannel and are mostly limited to neutral solute particles.^{35,36} Since Faxén derived the correction factor for the drag coefficients of particles close to an interface in the beginning of the 20th century,^{37} low Reynolds number hydrodynamics and classical Navier–Stokes equation have been used to model the transport of particles within confinements of similar dimensions. Quite recently, this classical framework was successfully applied to find approximate analytical formulas for the effective charge, diffusion coefficients, and mobility of proteins through nanopores.^{38} In Sec. II, we develop the general theory for this problem and describe our approach to numerically solve the corresponding convection–diffusion equation. In Sec. III, we derive a Taylor–Aris-like approximation for this problem. In Sec. IV, we present numerical and approximate results for realistic parameters, highlighting the feasibility of a single-molecule assay for measuring electric charge and diffusion by tracking single molecules in channels with non-zero surface voltage.

## II. THEORETICAL FOUNDATION

*L*and width

*h*(boundaries

*y*= ±

*h*/2) with a stationary laminar flow along the channel’s axis (

*x*-axis). Both surfaces of the slit channel are kept at the same constant surface potential (voltage), and the solute molecules that are transported by convection and diffusion bear the electric charge

*q*. For the sake of simplicity, we consider the case where the ratio of solute size to channel width is small enough to ignore Faxén’s correction factors for the drag coefficient, or of any electroviscous effects.

^{39}This assumption holds true for the majority of biomolecules of interest that are smaller than or equal to 10 nm in diameter inside a slit with a width in the order of 100 nm. In this case, the solute’s flux density is given by

*c*(

**r**,

*t*) is the local concentration of the solute molecules at position $r=x,y$ and time

*t*,

*D*is the diffusion coefficient,

*u*(

*y*) is the

*y*-dependent flow speed along the

*x*-direction,

*ψ*(

*y*) is the

*y*-dependent electric potential, and

*γ*=

*k*

_{B}

*T*/

*D*is the molecular friction coefficient. This leads to the following partial differential equation (PDE) for the solute transport:

^{40}

*η*is the solution’s viscosity and

*∂*

_{x}

*p*< 0 is a constant pressure gradient along

*x*. With the no-slip boundary conditions

*u*(

*y*= ±

*h*/2) = 0, this equation has the analytic solution

*u*

_{0}on the surface, this slip velocity has only to be added to the just found velocity profile. However, slipping should be important only for non-wetting liquids or nanostructured surfaces but less for smooth water/glass or water/metal interfaces, channel widths of 100 nm, and low Reynolds numbers as considered in the present paper.

^{41,42}

*ψ*shall obey the linearized Poisson–Boltzmann equation,

*κ*is the inverse Debye length of the solvent, which is given by

*N*

_{A}being the Avogadro–Loschmidt constant,

*e*being the elementary charge,

*ɛ*and

*I*being the dielectric constant and ionic strength (in mol/l) of the solution, respectively, and

*ɛ*

_{0}being the vacuum dielectric constant. In all numerical simulations below, we use

*ɛ*= 80.1 for water at room temperature (20 °C). With the boundary conditions

*ψ*(

*y*= ±

*h*/2) =

*ζ*, where

*ζ*denotes the surface-potential, this equation has the analytic solution

*x*-direction, i.e.,

*c*(0,

*y*,

*t*) =

*c*(

*L*,

*y*,

*t*), and no-flux conditions at the channel boundaries, i.e.,

*x*= 0 and

*x*=

*L*can essentially be set arbitrarily since their effect is negligible sufficiently away from the boundaries, particularly when

*L*≫

*h*. We will seek this solution in the form

*c*(

**r**,

*t*) is represented by a discrete Fourier series in the coordinate

*x*with discrete Fourier frequencies

*q*

_{k}= 2

*πk*/

*L*, where 2

*K*+ 1 is the number of summed frequency components. Here, we have introduced new functions $w\u0303k(y,t)$ that depend only on the transverse coordinate

*y*and time

*t*. This automatically assures periodic boundary conditions along

*x*. By inserting this

*ansatz*into Eq. (2), one finds that the auxiliary functions $w\u0303k$ obey the elliptic partial differential equations

*c*

_{0}(

**r**) is the initial concentration distribution at time

*t*= 0. When introducing new variables

*k*-dependent complex-valued potentials

## III. TAYLOR–ARIS APPROXIMATION

^{30,31,45}which gives an excellent description of the actual situation for small Péclet numbers. One can extend this approximation to the case of convection–diffusion transport of electrically charged solutes in an electric field that is considered here. For this purpose, we expand the concentration

*c*(

*x*,

*y*,

*t*) into the form $exp(\u2212\mu \psi )a\u0304(x,t)+b(x,y,t)$, where $a\u0304(x,t)$ does not depend on

*y*,

*b*(

*x*,

*y*,

*t*) describes small deviations from $a\u0304$, and the factor exp(−

*μψ*) is proportional to the Boltzmann equilibrium distribution of charge

*q*in an electric potential

*ψ*(

*y*). Inserting this into Eq. (2) and integrating over the transversal coordinate

*y*result in

*δu*=

*u*− ⟨

*u*⟩

_{ψ}. Following the arguments of Taylor and Aris

^{30,31,45}and assuming that the gradients of

*b*along

*y*are much bigger than gradients along

*x*and that the temporal change of

*b*is negligible, the last equation can be approximately written as

*A*and

*B*. Constant

*A*has to be zero to fulfill the boundary conditions, and constant

*B*is determined by requiring that $b\psi =0$ so that we find

*b*from Eq. (6) is proportional to $\u2202xa\u0304$, the resulting equation for $a\u0304$ is a convection–diffusion equation along only one spatial dimension

*x*with drift velocity

*σ*

_{0}along the

*x*-direction and a Boltzmann equilibrium distribution along the

*y*-direction, the approximate concentration distribution at time

*t*is

*u*

_{0}is the fluid velocity at the boundaries.

## IV. NUMERICAL RESULTS AND DISCUSSION

We numerically simulated the transport of a spatially localized concentration of solute along a 1 mm long slit channel with width *h* = 200 nm and with surface potential of *ζ* = 1 V. The initial concentration distribution was assumed to be Gaussian along the *x*-direction with a spatial square-root variance *σ*_{0} = 10 *μ*m. We performed the numerical simulation for two scenarios: one with an electrically neutral solute and one with a solute having an electric charge of *q* = 1*e*. The discretization grid of the numerical calculation had a grid spacing of 250 nm along the *x*-direction and of 1 nm along the *y*-direction. The mean flow velocity was set to $u\u0304=66.67\mu $m/s, and the solute’s diffusion coefficient was assumed to be *D* = 100 *µ*m^{2}/s. The calculation result is shown in Fig. 1, where the panels in the left column show the transport for electrically neutral molecules (i.e., for *q* = 0) and the right panels for solute molecules with a single positive elementary charge. As can be seen, the electric field forces the charged molecules toward the center of the channel, which leads to an on average faster transport along the *x*-direction for the charged solute as compared to that of the neutral solute.

A comparison of this exact numerical calculation with the Taylor–Aris approximation of Eq. (9) shows that the approximation is virtually indistinguishable from the exact result (see Fig. 2), where a cross section averaged broadening of concentrations at various *x*-positions is plotted for comparison. This is due to the extremely small Péclet number of ∼0.018, which is dominated by the small channel width *h*. For all parameter values relevant for the present study (average flow speed, channels size, and viscosity), one remains always in the regime of small Péclet numbers so that the Taylor–Aris approach is a perfect approximation of the exact results but offering a tremendously better computational efficiency of more than five orders of magnitude [for example, the numerical computation time for Fig. 2 was ∼450 s, whereas the computation time for the Taylor–Aris approximation was ∼3 ms using the same personal computer (PC) hardware and Matlab software platform]. The Taylor–Aris approximation also clearly shows that the small Péclet numbers lead to an effective diffusion coefficient, *D*_{eff}, that is virtually equal to *D*. However, the drift velocity *u*_{drift} becomes strongly charge-dependent and differs considerably from $u\u0304$ for charged solute molecules, in stark contrast to the classical Taylor–Aris dispersion, which leads to an enhanced diffusion coefficient but unchanged drift velocity.

Using Eq. (7), we calculated the dependence of the drift velocity as a function of surface potential *ζ* and Debye length *κ*^{−1} (i.e., ionic strength *I*) as shown in Fig. 3. The yellow line in Fig. 3 depicts the line of maximum drift velocity as a function of surface potential. For large surface potential above *ζ* = 0.5 V, this maximum drift velocity saturates and comes close to its maximally possible value of $3u\u0304/2$ (the maximum flow velocity at the channel’s center line). For negative surface potential, one sees a retardation of the drift velocity, which is due to increasing concentration of solute close to the channel walls due to electrostatic attraction. Thus, convection–diffusion of solutes in the slit channel with sufficiently large surface potential will lead to a separation between neutral and charged solutes. However, for large surface potential values as on the right side of Fig. 1, increasing the charge will nearly not change the drift velocity anymore. If one is interested in efficiently separating differently charged solutes, one should work with surface potentials where the slope of the yellow line in Fig. 2 is the largest, i.e., close to *ζ* ∼ 0, but not too small so that one can observe efficient separation of different solutes within reasonable channel length (and time). If one aims for a drift velocity enhancement of ∼10% over $u\u0304$ for a singly charged solute, reasonable surface potential values are around *ζ* = 50 mV for an ionic strength of ∼100 *μ*M (channel width 200 nm). As can be seen from Fig. 3, working at even higher ionic strength is possible when correspondingly working at larger values of surface potential (see stretching of relative drift velocity values along horizontal axis for small Debye length, i.e., high ionic strength). To estimate how the plot of Fig. 3 will change with smaller channel width, we calculated the ionic strength that maximizes the relative drift velocity at *ζ* = 1 V surface potential as a function of channel width. The result is shown in Fig. 4. As can be seen, for technically reasonable channel width values between 40 and 400 nm, this ionic strength varies between 1 *μ*M and 1 mM, and this roughly gives an estimate on how the vertical axes in Fig. 3 re-scale when changing the channel width. For example, when working with a channel having a width of 40 nm instead of 200 nm, one can work at ∼25 times higher ionic strength while observing a similar behavior as shown in Fig. 3.

*single solute molecules*. For this purpose, we calculated, in the Taylor–Aris approximation of Eq. (9), the relative drift velocity as a function of the product of electric charge and surface potential (note that only this product determines the drift velocity, not the separate values of charge or surface potential). The result for a 200 nm wide channel at ∼100

*μ*M ionic strength is presented in Fig. 5. The annotations in Fig. 5 refer to six different relative drift velocities for a solute molecule with a charge between 1 and 6 elementary charges at a surface potential of 50 mV. Let us denote by

*R*the difference of relative drift velocities between two charge states differing by one elementary charge, and let us assume that one tracks the motion of a molecule along a channel of length

*L*; then, the mean difference in traveled distances after run length

*L*is

*RL*, while at the same time, the diffusion-induced spreading of position is $2DL/u\u0304$. If we assume that for distinguishing different charge states based on mean drift velocity, the mean difference in traveled distances should be at least three times bigger than the diffusion-induced spreading, then we find the following estimate for the required minimum run length (channel length

*L*) over which on has to track a molecule:

*µ*m

^{2}/s (hydrodynamic radius of a solute molecule is ∼1 nm) and a mean flow speed $u\u0304$ of 1000

*μ*m/s, one needs a microfluidic channel length of ∼4500

*μ*m for distinguishing between

*q*= 5

*e*and

*q*= 6

*e*(

*R*≈ 0.02, see Fig. 5), with mean travel times (total observation time per molecule) of 3–4.5 s. To prevent photobleaching for such a long time, stroboscopic excitation could be employed. For smaller values of charges (

*q*= 4

*e*and

*q*= 5

*e*,

*q*= 3

*e*and

*q*= 4

*e*, etc.), shorter channel lengths would be sufficient. For example, one would need only a 60

*μ*m long channel for distinguishing between zero and one elementary charge (

*R*≈ 0.17), with a mean travel time of 40–60 ms, which would require, however, imaging at 1000 fps. However, our estimate shows that by tracking molecules along sufficiently long nanoslit channels, one should be able to resolve between 0 and 6 charge states when using a microfluidic chip with a total channel length of 4.5 mm, which could be realized in a compact way by a meandering channel design. By comparing the absolute mean drift velocity

*u*

_{drift}with the mean flow velocity $u\u0304$ for an applied surface potential and given geometry of the slit, the absolute charge of the solute species can be estimated. Moreover, evaluation of a molecule’s diffusive motion around its mean drift position would simultaneously allow for determining its diffusion coefficient and, thus, hydrodynamic size. This can be done using the existing single particle tracking algorithms that employ a stochastic drift-diffusion model for the particle motion along the

*x*-axis.

^{47,48}Note that in the presence of slip boundary conditions, the relative drift velocities

*R*will decrease with increasing fluid velocity at the boundary so that the minimum run length

*L*required for distinguishing different charge states will be longer.

We would like to emphasize that the mean-field based theoretical approach presented here is applicable as long as the channel surfaces can be considered to be in the “far-field” of the charged molecules, i.e., at least one Debye length, *κ*^{−1}, away. The channel width and electrolyte concentration considered in this work satisfy this criterion, as shown in Fig. 4. Furthermore, the considered Debye lengths are substantially longer than the Bjerrum length (0.71 nm for water at room temperature) so that ion–ion correlations are not significant and the Poisson–Boltzmann theory remains generally valid.^{49} The surface potential values that we consider in our paper and which are required for successfully separating charges are so large that the probability to find a (likely) charged molecule close to the walls is negligible as compared to the probability to find it toward the middle of the channel.

Importantly, the interaction of charged molecules with the channel walls is governed by their *effective* and not their formal structural charge. Whereas the latter can be derived from the structure of the molecular species in a straightforward manner, the former is a resultant of charge regulation and renormalization phenomena that take into account the acid–base equilibrium of the molecules’ ionizable groups modified by their local environment and by intramolecular Coulomb interactions, as well as counter-ion condensation on the molecules’ surface.^{50} The free energy of the long-range interactions between a molecule and the channel walls at a given electrolyte concentration is determined by this effective charge.^{51}

Finally, let us mention that a well-defined surface potential on the channel walls can be generated with transparent electrodes on these walls.

## V. CONCLUSION

We have considered the problem of the convection–diffusion of charged solute molecules in a nanoslit microfluidic channel with electric surface potential and transported by a laminar parabolic flow. We have derived the governing equations for this problem, and solved them numerically. We have also found a Taylor–Aris-like approximation of the solution, which yields excellent results for small Péclet numbers. Finally, we have discussed our numerical results in light of designing a single-molecule assay that could simultaneously measure hydrodynamic size and electric charge of biomolecules, providing reasonable estimates for channel width, length, surface potential, and flow speed and for the range of ionic strength values of the solvent for which such an assay could work. With our mean-field Poisson–Boltzmann theory, we currently consider only monovalent electrolytes neglecting ion–ion correlations, ion-specific effects such as ion affinity, and the finite size of ions. A suitable modified Poisson–Boltzmann equation can be incorporated in the future that takes into account most of these effects.

## ACKNOWLEDGMENTS

J.E. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy, Grant No. EXC 2067/1-390729940, and the European Research Council (ERC) via project “smMIET” (Grant Agreement No. 884488) under the European Union’s Horizon 2020 research and innovation program. M.C. would like to acknowledge the ICONIC EPSRC Programme Grant (No. EP/P020720/1) for support. N.K. acknowledges EPSRC funding for RFI.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors declare no conflict of interest.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Soft Condensed Matter Physics in Molecular and Cell Biology*

*q*-canonical Monte Carlo sampling for modeling the linkage between charge regulation and conformational equilibria of peptides

*Textbook of Structural Biology*

*Computational Approaches for the Prediction of pKa Values*

*Theoretical Microfluidics*