A Fast Stable Discretization of the
Constant–Convection–Diffusion–Reaction Equations
of Kinetic Capillary Electrophoresis (KCE)

The aim of this paper is to develop a numerical algorithm for solving a general system of equations – the Constant–Convection–Diffusion–Reaction (CCDR) equations (Sec. 1.2) – applicable to the experimental models of both Kinetic Capillary Electrophoresis (KCE; Sec. 1.2) and KSEC (Sec. 4.2). Both are reversible binding reactions in a long and narrow tube: a capillary or a separation column. The CCDR equations model chemical reactions between any number of substances, with an arbitrary reaction term, convected in a fluid and propagated by a constant electric field, describing the spacetime evolution of reactant concentrations. Though the convection is constant in our model, this does not alone eliminate the instability of standard discretization schemes which convection–diffusion equations are notorious for, but it does allow a change of variables which induces stability (Sec. 2.1).

The presented mathematical endeavour is fundamentally important for two directions of our ongoing scientific research in kinetic separation: (1) the accurate computer simulation of experiments under appropriate conditions; (2) the resolution of an experimental–computational inverse problem for determining kinetic rate constants – both within the KCE framework. For these applications, a practical numerical solution method must be both stable and accurate, though for the inverse problem, a low runtime of this direct solver is also critical.

Previously, our simulations of KCE were generated with COMSOL, which employs a streamline upwind Petrov–Galerkin method to solve the general Nernst–Planck Equations (Secs. 1.2 and 2.1.1). This method is neither globally stable however, as evidenced by its output, nor fast enough for resolving the aforementioned inverse problem, described below. Furthermore, the use of a readily available black box simulation software, like COMSOL, has obvious limitations, such as the lack of flexibility for modification or extension to a broader framework, like an inverse problem. By their very generality, such software are likely to be suboptimal for a specific purpose. Thus the introduced numerical algorithm is not only an accurate simulation tool, but also an efficient direct solver serving a larger vision within the KCE framework, as follows.

The KCE system of convection–diffusion–reaction equations includes kinetic rate constants in the reaction term, which can be viewed as parameters that the solution functions of this system of partial differential equations depend on. The efficient and accurate generation of such parametric solutions to the direct problem is necessary for the evaluation of the error target function at each iteration of some optimization algorithm, which minimizes this target function to resolve the inverse problem of determining the parameters that induce a given target signal [30].

Finding the accuracy of a numerical solution is difficult in general, but for the KCE model, it can be compared with an exact solution for a simplified case, derived in Sec. 2.2. The error is analyzed in Sec. 3.3, and the runtimes of the direct solver in Sec. 3.1 are compared in Sec. 3.2.

The equations of our physical model are deduced from the Nernst–Planck Equations, which express the conservation of mass of ions in a fluid medium under the influence of an electric field along a single spatial direction, while accounting for convection, diffusion, and reaction between the ions.

Define $c=(c_1,\dots ,c_N):R_{+}^2\to R_{+}^N$ (where $R_{+}:=[0,+\infty )$) as the spatiotemporal concentration of $N\in N$ convected ions over $(t,x)$ points, $V:R^2\to R^N$ as the velocities of their convection, $D:R_{+}^2\to R_{+}^N$ as their diffusion, $R:R_{+}^N\to R^N$ as a reaction term between the ions, $K\in R^N$ as coefficients arising from various electromagnetic and thermodynamic constants [33], and $E:R^2\to R^N$ as the electric fields influencing the motion of the ions, while denote with $\cdot$ the Hadamard product of two vectors. To arrive at an intermediate form towards our equations, the Nernst–Planck equations [23, 33] reduce to the following vector PDE

with appropriate initial and boundary conditions that ensure the existence and uniqueness of a concentration vector solution.

For the experimental model of Kinetic Capillary Electrophoresis (KCE) [5, 21, 19, 15], the above equation reduces further due to certain features of the experimental setup. Specifically the convection velocities are constant $V\in R^N$, the diffusion coefficients are constant $D\in R_{+}^N$, and the electric field strength $E$ is constant as well. Thus a constant vector $v:=V-K\cdot E\in R_{+}^N$ is introduced, and its components can be thought of as a case of constant convection in a convection–diffusion equation, suggesting the name Constant–Convection–Diffusion–Reaction (CCDR) equations. Therefore, they can be stated as

again with with appropriate initial and boundary conditions [25, 20].

In the formulation of the KCE model, there are three reactants and their concentration vector is denoted either as $c=(L,T,C):R_{+}^2\to R_{+}^3$ representing the experimental ligand, target, and complex, or as $c=(A,B,C)$ in the simplified case of MASKE [17]. The velocities and diffusion coefficients are denoted as $v=(v_L,v_T,v_C)$ and $D=(D_L,D_T,D_C)$ respectively. In the general KCE model, the reaction term is

Here $k_{on},k_{off}\in R_{+}$ are the rate constants of complex formation and dissociation respectively. In MASKE, the concentration of the second ion $T$ is assumed to be constant, so the second equation is omitted and the reaction term becomes $R\left(c\right)=(-k_{on}LT+k_{off}C,k_{on}LT-k_{off}C)$. The point of this simplification is to enable us to find an exact solution [17]. Another simplification for the same reason, is to keep all three concentration components and simplify the reaction term only as $R\left(c\right)=(k_{off}C,k_{off}C,-k_{off}C)$, as in Sec. 2.2 [18].

Since our goal is to give both a numerical and an exact solution to the above equations, the latter must be conveniently derivable from a Green function via function convolutions. So the initial and boundary conditions (IBC) of the CCDR equations are chosen to be idealizations of the initial concentration profiles measured in KCE experiments.

The initial conditions for the KCE equations specify the concentration profiles of the injected reactant plugs. A plug can be represented using various density functions $\varrho :R\to R_{+},{\int }_{-\infty }^{+\infty }\varrho =1$. An identical substance amount is assumed within the plug – regardless of the density function used – to ensure that the amount of reacting molecules remains consistent for various injected concentration profiles. The total amount of molecules at concentration $c\left(x\right)$ in a capillary of constant radius $r$, at locations $x$ within cylinders of infinitesimal height $dx$, is given by the integral

So for two concentration functions to give the same amount of molecules, the areas under them must be equal. To be consistent with our earlier work [15], this amount is set at $\pi r^{2}l\overline{c}$, which arises from the special case of rectangular concentration profiles of heights $\overline{c}\in R^N$ and plug length $l>0$. To rescale a density function $\varrho$ as another density, the proper transformation is $\left(1/l\right)\varrho \left(x/l\right)$, which retains a unit area under the curve. To ensure the standard amount of molecules, the area under the initial concentration profile $IC:=c(0,\cdot ):R_{+}\to R_{+}$ must equal $\overline{c}l$, implying the transformation $\overline{c}\varrho \left(x/l\right)$.

This way our earlier table of initial and boundary conditions [15] can be generalized to arbitrary plug densities as follows.

Note that $x_{det}$ denotes the location of the experimental detector at the end of the capillary, and $\vartheta$ denotes the Heaviside function

We plan to explicitly solve a simplified case of the CCDR vector PDE equation

with reaction function $R\left(c\right)=(k_{off}C,k_{off}C,-k_{off}C):R_{+}^2\to R^3$ where $c=(L,T,C)$, via the method of fundamental solutions as demonstrated earlier for the case of a rectangular plug [18]. Hereby the exact solutions are derived for a plug represented by a Gaussian initial condition, defined via the following density function

The NECEEM vector initial condition according to Sec. 1.3 is

while the boundary conditions remain unspecified for now, for the sake of simplicity.

Introducing some fundamental solutions will aid us in our derivation of an explicit solution. It can be shown by substitution that the equation

with the initial condition $F(0,x)=\delta \left(x\right)$, is solved by

On the other hand, the fundamental solution satisfying

with a spatiotemporal Dirac delta function, is coincidentally the function $F\left[0\right]$.

Thus to find the third concentration component $C$ for some initial condition function $C(0,x)={IC}_C\left(x\right)$, we must convolve it with the first fundamental solution above, resulting in

where $F_C$ is defined similarly as above, with parameters $v_C$ and $D_C$.

In order to get the other two concentration components $L$ and $T$, we break them up into an equilibrium term (the solution of the $k=0,F(0,x)=IC\left(x\right)$ case above) and a dissipation term (the solution with right-hand side $k_{off}C$), resulting in

Since the above formulas include several function convolutions (single and double integrals), it is desirable to simplify our derivation with the following well-known identity

Next, notice that for some ${\mu }_0,{\sigma }_0>0$ and $\mu :={\mu }_{0}l,\sigma :={\sigma }_{0}l$ we have

Employing the above general solution formulas for this special initial condition, we get

With the above, we may now derive the other concentration components as follows

and $T_{eq},T_{dis},T$ are defined similarly.

Despite the extensive simplification of the exact solutions for the Gaussian case above – otherwise containing even more computationally expensive convolution integrals – we will show in Sec. 3.2 that somewhat surprisingly, the numerical method of the previous section is much faster according to our computational tests. It is not only faster, but is capable of solving the CCDR equations for any reaction function and initial-boundary conditions, with a reasonable error according to Sec. 3.3.

The numerical solver algorithm of Sec. 2.1.2 has been implemented in MATLAB for the KCE equations, though can be easily generalized to more reactants and any reaction mapping. The package is available on GitHub [29]. Currently the package contains only the direct solver for KCE, but the inverse solver built on it is also under development, serving our research described in the introduction.

The direct solver sub-package, contains an implementation of the numerical algorithm that is functional for all initial and boundary conditions, and plug types described in Sec. 1.3. Exact solutions have only been derived and implemented for the simplified NECEEM case in Sec. 2.2 with Gaussian and rectangular plug types [18], and for MASKE [17] with any plug type.

Therefore, in order to compare the computational runtimes and errors between the numerical and exact solutions, we are restricted to either of the above two simplified methods. We chose simplified NECEEM with a Gaussian plug for the analysis below, since its equations are closer to the full KCE equations than those of MASKE, consisting of only two lines. Furthermore, the solution formulas can be minimized computationally as much as possible, by eliminating expensive convolution integrals which arise for other density functions. This way we achieve a lower bound on the computational complexity of the exact case, which can be compared with that of the numerical one. Specifically, three single convolution integrals are eliminated for $C,L_{eq},T_{eq}$ and two double convolution integrals are reduced to single ones for $L_{dis},T_{dis}$, as shown in Sec. 2.2.

The inverse problem of our research requires fast and accurate simulation of the solutions to the KCE equations for varying $k_{on},k_{off}$ values, modeling experimental electropherogram signals [30]. Thus we analyze our direct solver in terms of both runtime and error below.

In order to compare the runtime of our numerical solver (Sec. 2.1.2) to that of computing the exact solutions of simplified NECEEM (Sec. 2.2), we ran the solver for each temporal mesh size $I$ five times, and plotted the average of runtimes in Fig. 1. We are only interested in varying $I$ and not $J$ (the spatial mesh size), because the experimental signals are available only at the detector location $x_{det}$, restricting the inverse problem to this single location in space. Test runs were done with an AMD A8-5550M processor, 8.00 GB of RAM, 64-bit Windows 10, and MATLAB R2016a, with the parameters

The runtimes seem to closely follow quadratic and linear relationships for the numerical and exact solvers, respectively, with high coefficients of determination $R^2$. This is expected for the exact solver, which simply evaluates the formulas of Sec. 2.2. The quadratic polynomial for $T_N\left(I\right)$ makes sense, since the algorithm of Sec. 2.1.2 must necessarily solve for all the nodes on the spacetime mesh to provide the values at the only relevant location $x_{det}$. The transformed spatial meshes all depend on the temporal mesh linearly $\Delta y_n=v_n\Delta t$, implying that the spacetime mesh sizes must depend quadratically on $I$.

The exact solution formulas for simplified NECEEM (Sec. 2.2) provide a standard of comparison for the numerical solution values, since the formulas for $C,L_{eq},T_{eq}$ are exact, while the formulas for $L_{eq},T_{eq}$ can be approximated accurately with some quadrature method (we used the integral function of MATLAB, with a RelTol value of ${10}^{-6}$).

It is important to observe that the ratio of runtimes $r\left(I\right)$ decreases for increasing $I$, and follows a power law with an exponent of $-0.761$ (though the ratio of a linear and a quadratic polynomial should asymptotically follow a hyperbola with an exponent of $-1$). This implies that the gain in runtime using the numerical solver is much larger for lower mesh sizes, the difference ranging between 1-2 orders of magnitude below $I=2000$. For our end goal of an efficient inverse solver, a mesh size of $1000\le I\le 2000$ is sufficient.

Thus we have demonstrated that surprisingly the numerical solver has a significantly lower runtime than even an optimally low runtime of the exact solution formulas. For this to be relevant, the error between the solutions must be reasonable.

To show that the numerical solutions generated by Algorithm 1 are “reasonable”, we demonstrate that the $L^2$-error from the exact values vanishes with an increasing temporal mesh size $I$. Fig. 2 shows that the error decreases for increasing temporal mesh size according to similar power laws for the three concentration components. Taking the weighted average of the three exponents according to the coefficients of determination ($R^2$ values) gives $-0.9798$. Thus we conjecture that all three error components are inversely proportional to the temporal mesh size for large $I$, i.e. precisely of order $O\left(1/I\right)=O\left(\Delta t\right)$ (Sec. 2.1.2).

The convergence of Algorithm 1 for increasing mesh size is self-evident from its design. Clearly, the error in the solution – arising from the derivative approximations, and the multimesh interpolation – vanishes as the spacetime mesh size increases. For an infinitesimally fine mesh, the error is necessarily zero. Thus the above comparison on Fig. 2 – between the exact and simulated solutions in the simplified case – is not intended to be a proof of convergence, but merely a rate analysis of the de facto vanishing of the error with respect to mesh size. It gives an idea of how “reasonable” this convergence is.

Nevertheless, this computational error analysis of the simplified system bears relevance to the original equations. The two differ in the reaction terms: (1) original: $-k_{on}LT+k_{off}C$; (2) simplified: $k_{off}C$, but the solutions differ negligibly. As illustrated on Fig. 3 for the simplified solution, the peaks are separated in a way that the following hold

Consequently, the following approximations hold

So we may conclude heuristically, that the above rate analysis of the error in the simplified case, is closely aligned with the original.

A stable algorithm has been introduced for the efficient generation of accurate solutions to a set of convection–diffusion equations, with constant convection velocities. Stability has been demonstrated algebraically, due to the change of variables to new spatial meshes, for each line of this system of partial differential equations. Both efficiency and accuracy have been demonstrated as well.

The direct solver package corresponding to this work [29] may be extended in two potential directions: the implementation of new initial and boundary conditions, perhaps with a generic classification system; or the derivation of new exact solutions for either other simplifications of the original equations, or perhaps the explicit solution of the complete equations. According to our runtime analysis, however, the computation of closed form solutions is not likely to be more efficient, even in an extremely simple Gaussian case which gives a lower bound on computational complexity, implying the permanent significance of our solver.

The relevance of this algorithm to our further work on an experimental inverse problem has also been highlighted. Indeed, this is the direction our solver package is being extended in, now enabled by the fast computation of parametrized solutions with an error measured from an experimental signal, which is minimized via an optimization algorithm.

The Kinetic Size-Exclusion Chromatography (KSEC) model [3, 2] can be written under certain experimental conditions in a form quite similar to the KCE equations [7, 6]. However, our program package in its current form [29] cannot handle the KSEC equations, mainly due to additional factors in the reaction term and the initial conditions, even in a simplified case [7]. Furthermore, due to their more special form, the KSEC equations can be reduced to an equivalent form, for which a more efficient solver should be developed.

Algorithm 1 can also be adopted to the ACE model [14] and a separation-based approach [22], though it is no doubt applicable to the simulation of many other physical processes.

Both the KCE [15] and KSEC [3] equations describe kinetic separation methods, with concentration components as functions of time and space, where a reversible binding reaction is paralleled by the separation of reactants, convected at constant velocities in a capillary.

The presented discretization method was developed primarily for the KCE equations, which describe the migration of and interaction between three species: the ligand $L$, the target $T$, and the complex $C$. The KCE equations do not incorporate a variety of other “background” species and reactions – such as components of the acid–base equilibration in the background electrolyte – which occur during electrophoresis of these three species.

In practice, the KCE method utilizes experimentally-measured velocities of the three species, and plug lengths corresponding to them. Such parameters are sufficient for the description of the relevant migration–reaction phenomena, only if the electrophoretic experiment is planned and carried out properly.

Most importantly, the concentrations of the analytes in KCE should be well-below the concentrations of the ions in the background electrolyte. In turn, in KSEC experiments, the equilibration between the pores and the free volume should be fast. Beneficially, these conditions are easily satisfied for both KCE and KSEC. Analyte concentrations are typically at least four orders of magnitude below the background electrolyte concentration, while the sub-nanometer size of pores allows fast equilibration between them and the free volume.

This work was supported by the Natural Sciences and Engineering Research Council of Canada (grant: CRDPJ 485321-15). We are grateful to Mirzo Kanoatov for noticing implementation bugs in the direct solver package, and for his suggestions towards improving its user-friendliness. We also appreciate the valuable suggestions made by the referees.