A Statistical Approach to Estimating Adsorption-Isotherm Parameters in Gradient-Elution Preparative Liquid Chromatography

Without loss of generality, and for the convenience of readers who are interested only in our statistical approach that can be used for other types of real-world problems, we consider the following mass-balance equation of a two-component system for a fixed-bed chromatography column with the Danckwerts boundary condition, which is the most commonly used one for column chromatography (Ruthven (1984); Guiochon and Lin (2003); Horvath (1988); Javeed et al. (2011); Lin et al. (2017); Zhang et al. (2016b)).

where $x$ is distance, $t$ is time, and $i=1,\,2$ refers to the two components. $C$ and $q$ are the concentrations in the mobile and stationary phases, respectively, $u$ is the mobile phase velocity, and $F$ is the stationary-to-mobile phase ratio. $D_{a}$ is the diffusion parameter. $L$ is the length of the chromatographic column, and $T$ is an appropriate time point slightly larger than the dead time of chromatographic time $T_{0}=L/u$. In this paper, we set $T=1.5T_{0}$. In addition, $g(x)$ is the initial condition and $h(t)$ is the boundary condition, which describes the injection profile in the experiment. We adopt a simplified model here for ease of illustration; a more detailed version of gradient-elution liquid chromatography, with further discussion, can be found in Section S1 of the supplementary material.

Throughout this paper, we focus on the case in which the adsorption-energy distribution is bimodal. In this case, the bi-Langmuir adsorption isotherm is usually adopted as follows:

where subscripts I and II refer to two adsorption sites with different levels of adsorption energy.

In this paper, the collection of adsorption-isotherm parameters is denoted by

Now, we consider the measurement-data structure. In most laboratory and industry environments, the total response $R(\bm{\xi},t)$ is observed at the column outlet $x=L$ with

where $C(x,t)$ is the solution to problem (1) with the bi-Langmuir adsorption-isotherm model (2), and $C_{i}(L,t)$ represents the concentration of the $i$-th component at the outlet $x=L$. The parameter-to-measurement map $\mathcal{A}$: $\mathbb{R}^{8}\to L^{2}(\mathcal{T})$ can be expressed as

where the model operator $\mathcal{A}$ is defined through (4). To be more precise, for a given parameter $\bm{\xi}$, a bi-Langmuir adsorption-isotherm model can be constructed according to (2). Then, the concentration in mobile, i.e. $C$, can be obtained by solving PDE (1). Finally, the experimental data can be collected by using the designed sensor with the physical law (4). The aim of this paper is to estimate adsorption-isotherm parameters $\bm{\xi}$ from the time series database $R(\bm{\xi},t)$ and the integrated mathematical model (5) via a statistical approach.

In a liquid-chromatography experiment, a sampler brings the sample mixture into the mobile-phase stream, which carries it into a column, and pumps deliver the desired flow and composition of the mobile phase through the column. The detector located at the end of the column records a signal proportional to the amount of sample component emerging from the column, for time period $\mathcal{T}=[0,T]$. The signal recorded is the observation of interest.

To build up a statistical model, let $\bm{r}=(r(t_{1}),\cdots,r(t_{n}))^{T}$ be the observation points of the experiment, collected at discrete time points $\mathcal{T}_{n}=\{t_{1},\cdots,t_{n}\}\subseteq\mathcal{T}$. The liquid-chromatography data measured at time $t$ can be modeled by

where $\epsilon(t)$ represents the measurement noise. The clean liquid-chromatography data are $R(\bm{\xi},t)=\mathbb{E}[r(t)]$, where $\{C_{i}(L,t;\,\bm{\xi}):~{}i=1,\,2\}$ is the solution of a system of differential equations with static parameters $\bm{\xi}$; $\bm{\xi}=(a_{I,1},a_{II,1},b_{I,1},b_{II,1},a_{I,2},a_{II,2},b_{I,2},b_{II,2})$ represents the parameter of interest.

Let $\bm{R}(\bm{\xi})=(R(\bm{\xi},t_{1}),\cdots,R(\bm{\xi},t_{n}))^{T}$ be the collection of the exact chromatography signals with parameter $\bm{\xi}$. Then, the general framework of the chromatography measurement of $\bm{r}$ is represented by

$\bm{\epsilon}=(\epsilon(t_{1}),\cdots,\epsilon(t_{n}))^{T}$ stands for the measurement noise, and for simplicity we assume the noises $\epsilon(t)$ are uncorrelated between every pair, with zero mean and identical variances $\sigma^{2}_{\epsilon}$.

Throughout this paper, the framework is based on a single observation, but it can be generalized to a case with multiple observations or multiple injection groups. Assuming there is a single observation $\bm{r_{\text{obs}}}=\bm{R}(\bm{\xi}^{*})+\bm{\epsilon}^{*}$ with fixed parameter $\bm{\xi}^{*}$, the aim of the paper is to estimate $\bm{\xi}$ through the posterior distribution as described in the following sections, with all the other parameters apart from $\bm{\Theta}=\{\bm{\xi},\sigma^{2}_{\epsilon}\}$ assumed to be known.

A natural idea is to estimate the adsorption-isotherm parameters with optimization methods. For example, one may think the minimizer of the $L^{2}$ distance between $\bm{r_{\text{obs}}}$ and $\bm{R}(\bm{\xi})$ with respect to $\bm{\xi}$ as a good estimator. However, the non-smooth loss function and the non-unique global minimum prevent us from obtaining valid estimates by optimization.

More specifically, let us consider a four-dimensional parameter $\bm{\xi}$ by setting the parameters related to the second component as $0$. Since all the elements of $\bm{\xi}$ are positive, we can decompose $\bm{\xi}$ into a two-dimensional unbounded parameter $\bm{\nu}=(a_{I,1}+a_{II,1},b_{I,1}+b_{II,1})^{T}$ and a two-dimensional bounded parameter $\bm{\eta}=(a_{I,1}/\nu_{1},b_{I,1}/\nu_{2})^{T}$ without losing any information. Let $\bm{\xi}=g(\bm{\eta},\bm{\nu})$, and

be the loss function. While studying the structure of $\mathcal{L}$, we notice that it is convex with respect to $\bm{\nu}$ for any $\bm{\eta}$, as shown in Fig. 1(a). The topography suggests that an optimal estimator $\bm{\hat{\nu}}(\bm{\eta})=\arg\min_{\bm{\nu}}\mathcal{L}(\bm{\eta},\bm{\nu})$ can be easily obtained via optimization algorithms for each $\bm{\eta}$, as presented in Fig. 3(a)(c). However, if we try to optimize both $\bm{\eta}$ and $\bm{\nu}$ at the same time, the non-smooth loss function and potential multiple solutions, as illustrated in Fig. 1(b), will make the estimation invalid.

Similar loss functions can be observed when dealing with mixture models. Let us consider the simplest case. Assume the signal is the density function of a two-component Gaussian-mixture model with weights $\bm{\eta}$ and mean $\bm{\nu}$. The corresponding loss function is visualized in Fig. 2. Similar to that of the chromatography system, we can also easily calculate the corresponding optimal means for any fixed weights, but multiple solutions to the problem and saddle points make the optimization unreliable. Unfortunately, this problem will get worse as the number of components grows. Therefore, we usually use Markov chain-based Bayesian methods to make inferences about such parameters of interest. The traversability of these chains and the possible existence of ”label switches” also reduce the impact of multiple solutions.

In summary, the adsorption-isotherm parameters are not likely to be directly estimated via gradient-based optimization algorithms, but such optimization methods might lead to a partial solution under certain restricted conditions. Meanwhile, since the mixture models imitate the chromatography system well in terms of the loss function and the former is more computationally efficient, they can be used as toy examples to verify our proposed methods.

To investigate the parameter of interest, we aim to draw samples from the posterior distribution of $\bm{\xi}$ given the entire observation $\bm{r_{\text{obs}}}$. In this section, we adopt a Bayesian framework to sample from the posterior distribution $\pi(\bm{\theta}|\bm{r_{\text{obs}}})$. Throughout this paper, $p(\cdot)$ is a generic symbol for continuous probability densities, $\pi(\cdot)$ denotes the prior densities, $\pi(\cdot|\bm{r_{\text{obs}}})$ denotes the posterior densities, and $q(\cdot|\cdot)$ denotes the proposal densities used in Markov chain Monte Carlo (MCMC) algorithms.

For Bayesian inference, we assume that the measurement errors $\{\epsilon(t),t\in\mathcal{T}\}$ are normally distributed, and impose the following prior distributions on $\bm{\xi}$ and the error term variance $\sigma^{2}_{\epsilon}$:

where $N(\mu,\sigma^{2})$ represents the normal distribution with mean $\mu$ and variacne $\sigma^{2}$, $IG(\alpha,\beta)$ represents the inverse gamma distribution with shape parameter $\alpha$ and scale parameter $\beta$, and $\gamma$ is a tuning parameter. We let $\bm{\psi}=(\alpha,\beta,\gamma)$, which includes all hyperparameters.

Given this prior specification, the posterior distribution of $(\bm{\xi},\sigma^{2}_{\epsilon},\bm{R}(\xi))$ can be written explicitly as

where $\phi(x;\mu,\Sigma)$ stands for the multivariate normal density for the argument $x$ with mean $\mu$ and covariance matrix $\Sigma$, $I_{n}$ is the $n\times n$ identity matrix, and $E=\bm{R}(\bm{\xi})-\bm{r_{\text{obs}}}$.

After identifying the posterior distribution, we can use the Metropolis-Hastings algorithm to sample from the posterior. Because $\bm{\eta}$ and $\bm{\nu}$ are highly correlated in the posterior, an independent random-walk proposal may not converge quickly, or even get stuck somewhere. For this reason, we consider the Metropolis-within-Gibbs sampler instead. The method is outlined in Algorithm 3.1, in which the acceptance probability is calculated as

The hyperparameters $\bm{\psi}=(\alpha,\beta,\gamma)$ together with the variances of proposal distributions are tuned to make the chain as stable as possible and to have an acceptance rate of approximately $40\%$ to $80\%$ for each element of $(\bm{\xi},\sigma^{2}_{\epsilon})$. Then the sampler will generate a recurrent Markov chain whose stationary distribution is the target posterior distribution $\pi(\bm{\xi},\sigma^{2}_{\epsilon}~{}|~{}\bm{r_{\text{obs}}},\bm{\psi})$.

Modeling inverse problems as the posterior distributions will introduce correlation to the parameters of interest. When sampling more than one of the parameters directly from the posterior, the strong correlation might disrupt the motion of the Markov chain. To surmount this problem, we propose a dimension-reduction strategy with two modified MCMC algorithms to stabilize the posterior chains.

We first consider the degenerated case, in which the posterior distribution of $\bm{\xi}$ is almost singular and only supported on a lower dimensional set the size of observed data tends to infinity. Specifically, suppose that the parameter of interest $\bm{\xi}\in\mathbb{R}^{D}$ can be written as $\bm{\xi}=g(\bm{\eta},\bm{\nu})$, where $\bm{\eta}\in\mathbb{R}^{d}$ with $d\ll D$, $\bm{\nu}\in\mathbb{R}^{D-d}$, and $g(\cdot):\mathbb{R}^{D}\to\mathbb{R}^{D}$ is a one-to-one mapping determined with pilot study. Roughly speaking, we first set $d$ as the largest dimension that leads to a stable posterior chain of $\bm{\eta}$. Then the dimensionality is reduced by combining the elements of $\bm{\xi}$ according to the pilot study. Please see Section S4 of the supplementary material for a detailed discussion on how to choose $\bm{\eta},\bm{\nu},g$ and some empirical results.

In the extreme case that $\bm{\nu}$ and $\bm{\eta}$ are perfectly dependent in the posterior as the size of observed data $n$ tends to infinity, the conditional posterior $\pi(\bm{\nu}|\bm{\eta},\bm{r_{\text{obs}}},\bm{\psi})$ is degenerate to a Dirac distribution $\delta_{h(\bm{\eta})}$ for some unknown function $h$, i.e. it is almost sure that $\bm{\nu}=h(\bm{\eta})$. Then, the support of the posterior of $\bm{\xi}$ can be expressed as

and the induced posterior distribution of $(\bm{\eta},\bm{\nu})$ is

In such a scenario, we can find the least squares estimator of $\bm{\nu}$ in terms of $\bm{\eta},\bm{r_{\text{obs}}},\bm{\psi}$ from the following equation:

such that $\bm{\hat{\nu}}(\bm{\eta})$ approximates $h(\bm{\eta})$ as $n\to\infty$. In what follows, we assume that the optimization problem in (9) can be settled by the following gradient descent algorithm (GD) without discussing the optimization accuracy.

To incorporate the optimizing step in the MCMC algorithm, we keep the overall structure of Algorithm 3.1 for sampling $\bm{\eta}$ and $\sigma_{\epsilon}^{2}$. In the extreme case of $\bm{\nu}=h(\bm{\eta})$ for some function $h$, instead of sampling $\bm{\nu}$, we can calculate a $\bm{\hat{\nu}}(\bm{\eta})$ from (9) in each MCMC iteration. Together with the function $g$, $\bm{\hat{\xi}}=g(\bm{\eta},\bm{\hat{\nu}}(\bm{\eta}))$ can be obtained for any specific $\bm{\eta}$. The detailed algorithm is provided in Algorithm 3.3.

Our theoretical justification for this dimension-reduction strategy is as follows. We first introduce a series of technical assumptions, and our theoretical results will depend on different subsets of these assumptions. We define the loss function $L_{R}(\bm{\xi})=\int_{\mathcal{T}}[R(\bm{\xi},t)-R(\bm{\xi}^{*},t)]^{2}dt$ associated with the solution function $R$, for any $\bm{\xi}\in\bm{\Xi}$.

• (A.1)

  $(\bm{\nu},\bm{\eta})$ lies in a compact space $\bm{\Theta}\subseteq\mathbb{R}^{D}$, where $\bm{\eta}\in\bm{\Theta}_{\bm{\eta}}\subseteq\mathbb{R}^{d}$ and $\bm{\nu}\in\bm{\Theta}_{\bm{\nu}}\subseteq\mathbb{R}^{D-d}$. There is a continuous one-to-one mapping such that $\bm{\xi}=g(\bm{\eta},\bm{\nu})$ and $\bm{\Xi}=g(\bm{\Theta})$ is also a compact space for $\bm{\xi}$. The true parameter $\bm{\xi}^{*}=g(\bm{\eta}^{*},\bm{\nu}^{*})$ is an interior point of $\bm{\Xi}$. $\mathcal{T}$ is a compact set, and $\{\epsilon(t):t\in\mathcal{T}\}$ are independent and identically distributed as $N(0,\sigma_{\epsilon}^{2*})$ for some $\sigma_{\epsilon}^{2*}<\infty$.

• (A.2)

  Let $|\mathcal{T}|$ be the Lebesgue measure of $\mathcal{T}$. As $n\to\infty$,

  $\displaystyle\sup_{\bm{\xi}\in\bm{\Xi}}\left|\frac{|\mathcal{T}|}{n}\left\|\bm{R}(\bm{\xi})-\bm{R}(\bm{\xi}^{*})\right\|_{2}^{2}-\int_{\mathcal{T}}[R(\bm{\xi},t)-R(\bm{\xi}^{*},t)]^{2}dt\right|\to 0.$

  (10)

• (A.3)

  There exists a function $h(\bm{\eta})$, and positive constants $a_{1R},c_{1R},\kappa_{1R}$, such that for any $(\bm{\eta},\bm{\nu})\in\bm{\Theta}$,

  $\displaystyle L_{R}(g(\bm{\eta},\bm{\nu}))-L_{R}(g(\bm{\eta},h(\bm{\eta})))\geq\min\left(a_{1R}\|\bm{\nu}-h(\bm{\eta})\|_{2}^{\kappa_{1R}},c_{1R}\right).$

  (11)

• (A.4)

  There exist positive constants $a_{2R},c_{2R},\kappa_{2R}$, such that, for all $\bm{\xi}\in\bm{\Xi}$ that satisfies $\|\bm{\xi}-\bm{\xi}^{*}\|_{2}\leq c_{2R}$, $L_{R}(\bm{\xi})\leq a_{2R}\|\bm{\xi}-\bm{\xi}^{*}\|_{2}^{\kappa_{2R}}$.

• (A.5)

  There exist positive constants $a_{3R},c_{3R},\kappa_{3R}$, such that, for all $\bm{\xi}\in\bm{\Xi}$, $L_{R}(\bm{\xi})\geq\min\left(a_{3R}\|\bm{\xi}-\bm{\xi}^{*}\|_{2}^{\kappa_{3R}},c_{3R}\right)$.

Assumption A.1 assumes a compact parameter space and normal errors. Assumption A.2 assumes the uniform convergence of the empirical $L_{2}$ distance from the function $R(\bm{\xi},t)$ to the true function $R(\bm{\xi}^{*},t)$ over the compact parameter space $\bm{\Xi}$. This assumption usually holds when $\mathcal{T}_{n}=\{t_{1},\ldots,t_{n}\}$ are densely distributed in $\mathcal{T}$ and $R(\bm{\xi},t)$ is continuous in both $\bm{\xi}$ and $t$. Assumption A.3 is the identification condition for $\bm{\nu}$, where the lower bound guarantees that $h(\bm{\eta})$ is the unique minimizer of the loss function $L_{R}(g(\bm{\eta},\bm{\nu}))$ over the $\bm{\nu}$ argument for any given $\bm{\eta}$. This assumption is only used for justifying our later Algorithm 3.3 in the degenerate case of $\bm{\nu}=h(\bm{\eta})$. Assumption A.4 imposes a local-continuity condition on the loss function $L_{R}(\bm{\xi})$ in a small neighborhood of $\bm{\xi}^{*}$. Since we have already assumed that $\bm{\xi}^{*}$ is an interior point of $\bm{\Xi}$ in A.1, the radius $c_{2R}$ can be made small such that the whole ball $\{\bm{\xi}:\|\bm{\xi}-\bm{\xi}^{*}\|_{2}\leq c_{2R}\}$ is included in $\bm{\Xi}$. Assumption A.5 imposes an identification condition on the loss function $L_{R}(\cdot)$ which guarantees that the true parameter $\bm{\xi}^{*}$ is uniquely identified in $\bm{\Xi}$. This is similar to the identification condition for moment estimation (see, for example, ZE.2 in Belloni and Chernozhukov 2009). The power constants $\kappa_{1R}$, $\kappa_{2R}$, and $\kappa_{3R}$ in A.3, A.4, and A.5 are not specified to allow flexibility in the local-continuity property of $L_{R}(\cdot)$. In Section 2.2 of the supplementary material, we provide some partial verification for these assumptions using the model of simulation case 1 in Section 4.1.

Theorem 3.1 provides justification for our Algorithm 3.3, in the sense that, in each MCMC iteration, we can numerically solve for $\bm{\nu}$ in terms of $\bm{\eta}$ from (9). Part (i) shows that, as the sample size of observations $n$ increases, the empirical solution $\hat{\bm{\nu}}(\bm{\eta})$ becomes increasingly close to $h(\bm{\eta})$, which is the unique minimizer of the loss function $L_{R}(\cdot)$. Part (ii) shows that, under the additional continuity condition on $L_{R}(\cdot)$, most of the posterior probability mass of $\bm{\nu}$ is concentrated around $h(\bm{\eta})$. Part (iii) gives the stronger result of posterior consistency of $\bm{\xi}$ to the truth $\bm{\xi}^{*}$ (Chapter 6, Ghosal and van der Vaart 2017) under the additional identification condition A.5 on $\bm{\xi}$, which implies that most of the posterior probability mass will concentrate around the true parameter $\bm{\xi}^{*}$ asymptotically. The detailed proof of Theorem 3.1 can be found in Section S2.1 of the supplementary material. We emphasize that even when the strong assumption of A.3 that assumes $\bm{\nu}=h(\bm{\eta})$ does not hold, Part (iii) of Theorem 3.1 on the posterior consistency for the whole parameter vector $\bm{\xi}$ remains valid as the sample size $n\to\infty$.

After identifying the method of dimensionality reduction and restoration, we propose a new modified Metropolis-embedded gradient-descent-within-Gibbs sampler (MGDG). The detailed algorithm is presented below.

The main strategy of this algorithm is to sample only $\bm{\eta}\in\mathbb{R}^{d}$ and the error term variance $\sigma^{2}_{\epsilon}$ from the posterior distribution, and $\bm{\nu}$ is treated as a function of $\bm{\eta}$ in the posterior. We can restore a posterior sample of $\bm{\xi}=g(\bm{\eta},\bm{\nu})$ using the posterior draws of $\bm{\eta}$ and $\bm{\nu}$. The lower dimensionality of $\bm{\eta}$ helps to improve the mixing of Markov chains and prevent the ”label switching” phenomenon that occurs in mixture models (Jasra, Holmes and Stephens 2005). After reaching the stationary distribution, the samples $\{\hat{\bm{\nu}}^{(k)}\}$ should be distributed around the truth $h(\bm{\nu})$, as shown in Theorem 3.1.

Now we consider the general case in which $\bm{\nu}$ is not a function of $\bm{\eta}$, but may be highly correlated with $\bm{\eta}$ in the posterior distribution $\pi(\bm{\xi},\sigma_{\epsilon}^{2}|r_{\text{obs}},\bm{\psi})$. In such a scenario, we use the Metropolis-adjusted Langevin algorithm (MALA) (Roberts and Rosenthal, 1998; Roberts and Tweedie, 1996) to sample $\bm{\nu}$. We still use the Metropolis-within-Gibbs algorithm to draw $\bm{\eta}$ and $\sigma^{2}_{\epsilon}$, and then use the MALA to draw $\bm{\nu}$ using the gradient information from the conditional posterior of $\bm{\nu}$. We call this algorithm the Metropolis-adjusted Langevin-dynamics-within-Gibbs sampler (MALG).

  Algorithm 4 Metropolis-adjusted Langevin-dynamics-within-Gibbs sampler (MALG)

where the values of the proposal densities, $q_{LD}(\bm{\nu}^{\prime}|\bm{\nu})$ and $q_{LD}(\bm{\nu}|\bm{\nu}^{\prime})$, are calculated from

Because the log-likelihood term $\log\pi(\bm{\nu}|\bm{\eta},\sigma_{\epsilon}^{2},\bm{R}(\bm{\xi}),\bm{r_{\text{obs}}})$ is proportional to $\|\bm{R}(g(\bm{\eta},\bm{\nu}))-\bm{r_{\text{obs}}}\|_{2}^{2}$, the gradient term is the same gradient $\nabla_{\bm{\nu}}\|\bm{R}(g(\bm{\eta},\bm{\nu}))-\bm{r_{\text{obs}}}\|_{2}^{2}$ as used in Algorithm 3.3. After choosing suitable values for the stepsize $\tau$ and the sub-chain length $m$, this algorithm will return MCMC samples of all parameters, including $\bm{\nu}$.

For the first case, we considered a Gaussian-mixture model with two components. Let the parameter of interest be $\bm{\xi}\in\mathbb{R}_{+}^{4}$, with a single observation $\bm{r_{\text{obs}}}=\bm{R}(\bm{\xi}^{*})+\bm{\epsilon}$ with $\bm{\xi}^{*}=(\frac{1}{3},\frac{2}{3},\frac{8}{3},\frac{4}{3})$ and noise variance $\sigma^{2*}_{\epsilon}=0.001$ from the following numerical solver:

To perform the dimension reduction, we selected $\bm{\eta}\in\mathbb{R}^{2}$ and $\bm{\nu}\in\mathbb{R}^{2}$ as follows:

These two parameters can be regarded as weights and means in the Gaussian-mixture model. Because $\bm{\eta}\in[0,1]^{2}$, the MGDG algorithm can be initialized by uniformly sampling $\{\bm{\eta}_{i}\}_{i=1}^{1000}\sim\text{Uniform}(0,1)^{2}$ and setting $\bm{\eta}^{(0)}$ as the one minimizing the $L^{2}$ norm between $\bm{R}(g(\bm{\eta}_{i},\bm{\hat{\nu}}(\bm{\eta}_{i})))$ and $\bm{r_{\text{obs}}}$, i.e.