Adjoint optimization for the general rate model of liquid chromatography

Adjoint optimization for the general rate model of liquid chromatography

Journal Pre-proof

Adjoint optimization for the general rate model of liquid chromatography ¨ Sebastian Osterroth, Peter Menstell, Achim Schwammle, Joachim Ohser, Konrad Steiner PII: DOI: Reference:

S0098-1354(19)30396-5 https://doi.org/10.1016/j.compchemeng.2019.106657 CACE 106657

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

10 April 2019 6 November 2019 22 November 2019

¨ Please cite this article as: Sebastian Osterroth, Peter Menstell, Achim Schwammle, Joachim Ohser, Konrad Steiner, Adjoint optimization for the general rate model of liquid chromatography, Computers and Chemical Engineering (2019), doi: https://doi.org/10.1016/j.compchemeng.2019.106657

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights • Formal derivation of adjoint formulation of general rate model for optimization • Simulation of purification process for parameter identification and optimization • Comprehension of confocal laser scanning microscopy (CLSM) data • Comparison to numerical and experimental data

1

Adjoint optimization for the general rate model of liquid chromatography Sebastian Osterroth1,* , Peter Menstell2 , Achim Schw¨ammle2 , Joachim Ohser3 , and Konrad Steiner1 1

Fraunhofer Institute for Industrial Mathematics ITWM, Kaiserslautern, Germany 2 Life Science Division, Merck KGaA, Darmstadt Germany 3 University of Applied Sciences Darmstadt, Darmstadt, Germany * corresponding author November 23, 2019 Abstract

Liquid chromatography is an important step in the purification of biopharmaceuticals. Here simulation is a helpful tool to assist the development and optimization of chromatographic materials and processes. The most general description of the process is the general rate model. This article describes a general approach to identify the model parameters using the adjoint method. The method minimizes the discrepancy between given (experimental) data and the simulation. As additional information, intraparticle concentration profiles can be visualized by confocal laser scanning microscopy. A method to correlate these data to protein concentrations is derived. The examples show that the approach is successful in identifying the model parameters.

Keywords: Liquid chromatography; Protein purification; Modeling and simulation; Adjoint method; Parameter identification; Confocal laser scanning microscopy

1

Introduction

The use of simulation is more and more a helpful tool to assist the design process or to optimize the experimental workflow in the lab. However, simulation can never fully replace experiments since parameter calibration and validation are needed. For determining the performance and efficiency of chromatographic materials, mostly two modes are used, static (batch) and dynamic (column) experiments. In the static mode, chromatographic media are incubated in a protein solution with defined concentration and incubation time. In the dynamic case, a defined protein solution flows through a column packed with chromatographic media. Increasing complexity and detailing of the experiments make it necessary to adjust and modify the models. However, it is often possible to stick to a simplified or lumped model to reproduce the most important characteristics of the chromatographic process. The most versatile model for chromatography is the so-called general rate model (GRM), see e.g. Guiochon et al. (2006); Michel et al. (2005). Many other models can be derived from its rather general formulation, e.g., the equilibrium dispersive (ED) (e.g. Forss´en et al. (2006)) or the transport dispersive (TD) model (e.g. Xu et al. (2013)). The GRM combines several transport mechanisms in the liquid phase and models the transport in the microporous bead, too. The model can be combined with various isotherms or kinetics modeling the adsorption (and desorption) behavior within the bead. The level of detail of the used model also depends on the available experimental data. Since often no intraparticle data is available, it might not be necessary to describe the behavior within the beads in all details. In its basic form, the model consists of a set of partial differential equations (PDEs). A one-dimensional convection-diffusion/dispersion-reaction (CDR) equation in flow direction is coupled to the dynamics within the beads described by two diffusion equations. Those diffusion equations are formulated in spherical coordinates; however, it is assumed that the transport is independent of the two angles. This results in two one-dimensional equations in the radial direction of the beads. One special form of the general rate model we consider in this article is the case of the static mode (batch) (e.g. Carta

2

et al. (2005); Yang et al. (2006)). Here the beads are incubated in a protein solution without any convective or dispersive transport. This means that the CDR equation reduces to a reaction equation (cf. Lewus et al. (1998)). Frequently the equations are made nondimensional (Bowes and Lenhoff, 2011; DiLeo et al., 2009; Gu et al., 2013; Mehay and Gu, 2014; Michel et al., 2005; Qamar et al., 2016) introducing dimensionless variables. This often simplifies the analysis. Optimization or parameter identification can be performed in different ways. One usually starts with the definition of an objective function which is displaying the difference between the actual configuration and some desired result. This is described by the so-called cost functional. Then one tries to adjust the parameters such that the difference is minimized. Thus, it has to be detected, which parameter has to be increased and which to be decreased. Mathematically spoken, a descent direction has to be found. For the solution of the problem, mainly two different categories of methods can be named: gradient-based or gradient-free algorithms. As a basic result, one can show that the negative gradient of the cost functional always points in the direction of the steepest descent. Therefore, this direction is suitable as descent direction. However, it is sometimes computationally not feasible, very difficult or even impossible to compute the gradient. Then gradient-free methods as genetic algorithms or the Nelder-Mead method can be used (Nocedal and Wright, 2006b). Another way is to approximate the exact gradient. Often finite differences are used to compute the so-called sensitivities (Nocedal and Wright, 2006a). Another possibility is to use automatic (also called algorithmic) differentiation (Nocedal and Wright, 2006a). In this article, we compute the exact gradient by means of the so-called adjoint method (Hinze et al., 2008). In the field of PDEs, the method originates from inverse problems and avoids the computation of an inverse of the constraints of the underlying problem. In this method, one computes the exact gradient of the objective function based on the solution of an additional system of equations (the socalled adjoint equations). Since the underlying constraints of the optimization problem are partial differential equations, there are two ways to derive the adjoint equations: First optimize, then discretize or first discretize, then optimize (Hinze et al., 2008; Tr¨oltzsch, 2010). The first way is called the continuous approach and the second the discrete one. In the field of chromatography, various authors discussed parameter identification problems and used various methods. We first summarize the works using a continuous approach. The adjoint method was applied for the equilibrium model of ideal chromatography in James et al. (1999). The authors derive the adjoint state equation and then use a finite difference discretization and a conjugate gradient algorithm to derive a solution. Recently, the method was used for the equilibrium dispersive and the transport dispersive model (Cheng et al., 2018; Zhang et al., 2017). A larger number of authors worked with the discrete approach. The general rate model is considered in Persson and Nilsson (2001) and Persson et al. (2006) and derivative-free methods are used for parameter identification (Gauß-Newton or Nelder-Mead method). In Forss´en et al. (2006), the ED model is discretized using the Rouchon finite difference scheme. To approximate the gradient, numerical differentiation with complex variables is used. The same model is investigated in Kaczmarski (2007) and the discretization with Rouchon finite differences and the orthogonal collocation method on finite elements are compared. The Levenberg-Marquardt method is used to find the desired isotherm parameters. An approach for the parameter identification in the batch mode based on genetic algorithms is presented in Mesa et al. (2011). A similar derivative-free scheme is applied in Xu et al. (2013). The Martin-Synge algorithm is used for the solution of the underlying equations and a genetic algorithm in combination with the Levenberg-Marquardt method finds the desired parameters. In Hahn et al. (2014), the underlying partial differential equations are first linearized and then discretized by finite elements. Based on this discrete formulation, the optimality conditions for computing the exact (discrete) gradient are derived. Some Newton method is used for the solution. A more involved problem is discussed by the group of von Lieres (Leweke and von Lieres, 2018; Puettmann et al., 2013, 2016; von Lieres and Andersson, 2010). They consider the general rate model in combination with steric mass action (SMA) to describe the adsorption behavior of a multi-species mixture. This results in a system of partial differential algebraic equations (PDAEs). They discretize the system of equations using the finite volume method and tackle the obtained DAE. For the solution, they need the Jacobian which is derived by automatic differentiation. To investigate the influence of the parameters, they use the forward sensitivity approach. For each sensitivity (variation with respect to one parameter) of the state variables, a DAE is derived and solved. Finally, all ingredients are combined to derive the gradient. The optimization procedure is performed using MATLAB. In this article, we derive the adjoint formulation for the general rate model in a form that is valid for the dynamic and the static case and also for the equations in dimensioned and dimensionless form. It is shown that the structure of the underlying equations stays in all cases the same. This allows us to derive also a generic form of the adjoint equations. 3

Besides the optimization of the experimental conditions, modeling of chromatographic separation processes can give further insight into the mechanisms that are taking place during the separation process. The simulation can visualize concentration profiles for all species, both within the column and a single chromatographic bead. Therefore, it allows the comparison to concentration data (from column or batch) and concentration profiles within the bead, as they are visualized for example by confocal laser scanning microscopy (CLSM). The method uses a spatial pinhole in order to block out-of-focus light inducing high spatial resolution as well as sufficiently high image contrast, see e.g. Shinya (2006) for more details. The pixel values of the three-dimensional CLSM images can be interpreted as relative fluorophore density inside the beads, where estimations of densities are significantly influenced by attenuation of the excitation beam as well as of the emitted fluorescence light (Bankston, 2008). Various methods have been developed for the correction of these unwanted attenuation effects (Du and Sun, 2009; Kervrann et al., 2004; Lee and Bajcsy, 2006; Mich´alek et al., 2010; Susanto et al., 2006). A review of the techniques is given in Bradley and Thornley (2006). The article is organized as follows: First, the equations of the general rate model are stated and a generalized form is derived. In the sequel, the first-order optimality conditions are derived based on the formulation of a minimization problem and the usage of an adjoint method. Based on this formulation, a descent algorithm is derived to solve the parameter identification problem in order to obtain the model parameters. Having all mathematical details at hand, the derived method is applied to some numerical test cases and to a set of real-world problems. Therefore, the experimental conditions are stated and the approximation quality of the method is shown. For the use of CLSM data, a correlation of fluorescent concentration to protein concentration is stated. The article closes with some conclusions and an outlook.

2

Theory

Consider a column of height H, which is filled with spherical microporous beads of radius rp and it holds rp  H. The random sphere packing of the column has a porosity ε. The axial coordinate of the column is denoted by z, z ∈ [0, H] and the time variable is t. The radial coordinate of the bead is described by r, r ∈ [0, rp ]. The final time of the chromatography process is denoted by Tend . We assume that the inlet of the column is located at z = H and the outlet at z = 0.

2.1

General rate model

The most detailed model to describe the processes of chromatography is the general rate model (GRM). We stick here to its one-dimensional as for example presented in Michel et al. (2005). An extension accounting also for radial dispersion effects is made in Qamar et al. (2017). Other extensions include the definition of two adsorption zones (Kavoosi et al., 2007) or the influence of pH (Gu´elat et al., 2016). Note that an extension to higher spatial dimensions would also require modeling the flow behavior inside the column (cf. Schnittert et al. (2009)). Consider a multicomponent mixture of Nc different species. The concentration of component i in the column is denoted by ci and depends on the axial coordinate z and time t. In the porous chromatographic material (bead), the concentrations cp,i of mobile (dissolved) and qp,i of stationary (adsorbed) phase are considered. Within the bead, the concentrations are assumed to be independent of the two angles and only depends on the radial coordinate r. Nevertheless, these properties also depend on z and t. Introduce the component vectors T T T c = (c1 , . . . , cNc ) , cp = (cp,1 , . . . , cp,Nc ) , and qp = (qp,1 , . . . , qp,Nc ) . Then the general rate model can be stated as ∂ci ∂ 2 ci 3 1−ε ∂ci + uint − Dax 2 = −kf,i (ci − cp,i (·, r = rp , ·)) , z ∈ (0, H), ∂t ∂z ∂z rp ε   ∂cp,i Dc,i ∂ ∂cp,i − 2 r2 = −ψi (cp , qp ) , r ∈ (0, rp ), ∂t r ∂r ∂r   ∂qp,i Dq,i ∂ ∂qp,i − 2 r2 = ψi (cp , qp ) , r ∈ (0, rp ), ∂t r ∂r ∂r

(1a) z ∈ (0, H),

(1b)

z ∈ (0, H),

(1c)

for t ∈ (0, Tend ). The first equation describes the transport of the concentration ci by convection with velocity uint and axial dispersion with coefficient Dax . The effective mass transfer to the adsorbent bead is described by kf,i . With the second and third equation, the mass transport by diffusion within the bead is described. We consider both, pore diffusion (Dc,i ) and surface diffusion (Dq,i ). An adsorption rate (kinetic) is given by ψi (cp , qp ) and will be specified later on. For more details on the considered phenomena and the derivation 4

of the equations, we refer to Michel et al. (2005). Note that we consider the flow in negative z-direction and therefore uint < 0 (take the negative value of the given data). The boundary conditions are given as ∂ci (H, t) = uint cin,i , ∂z ∂ci Dax (0, t) = 0, ∂z

uint ci (H, t) − Dax

(2a) (2b)

∂cp,i (·, r = rp , t) = kf,i (ci (·, t) − cp,i (·, r = rp , t)) , z ∈ (0, H), ∂r ∂cp,i Dc,i (·, r = 0, t) = 0, z ∈ (0, H), ∂r ∂qp,i Dq,i (·, r = rp , t) = 0, z ∈ (0, H), ∂r ∂qp,i Dq,i (·, r = 0, t) = 0, z ∈ (0, H). ∂r Dc,i

(2c) (2d) (2e) (2f)

The first two conditions are related to the inlet and the outlet of the column and the remaining conditions describe the boundary conditions at the surface and center of the bead. Note that compared to the literature, we have not included the particle porosity εp into the equations. This can be motivated in two different ways: Since the model does not account for a change in the particle porosity, either the concentrations cp and qp can directly be defined to include the particle porosity or the measured particle porosity can be incorporated into the value of the coefficients. The initial conditions depend on the actual mode of operation. In the dynamic mode, the column is fed with an inflow concentration cin,i , and therefore, a zero initial condition is described. In the static mode, the beads are incubated in a given initial concentration cinitial,i (and no boundary condition is prescribed). The beads are always considered as free of protein concentration at the initial time. Hence, the conditions read ( 0, column ci (z, 0) = , z ∈ (0, H), (3a) cinitial,i , batch r ∈ (0, rp ),

cp,i (·, r, 0) = 0,

r ∈ (0, rp ),

qp,i (·, r, 0) = 0,

z ∈ (0, H),

z ∈ (0, H).

(3b) (3c)

The equations are simplified to deal with a finite batch (static mode). In this case we set uint = Dax = 0. The convection-diffusion-reaction equation for the column is then replaced by a reaction equation and the boundary conditions for inlet and outlet vanish since none of them is present. 2.1.1

Adsorption kinetics

Sorption isotherms describe the equilibrium between adsorption and desorption at surfaces. More details can be found in Lienqueo et al. (2012) and Schulte and Epping (2005). In this article, three different isotherms in the multi-component form are considered. Since the equilibrium between adsorption and desorption might not be directly reached, we consider the non-equilibrium form of the isotherms (Ma et al., 1996; Michel et al., 2005). Such a form is often called (adsorption) kinetic (e.g. Golshan-Shirazi and Guiochon (1992); Suen (1996)). Linear kinetics The simplest version is the linear kinetic. In case of an equilibrium, it is also called Henry isotherm. It states that in the equilibrium the adsorbed and the mobile concentration within the bead are proportional to each other. It can be written in non-equilibrium form as ψL,i (cp , qp ) = ψL,i (cp,i , qp,i ) = ka,i cp,i − kd,i qp,i ,

(4)

see e.g. Michel et al. (2005). The parameter ka,i describes the adsorption and the parameter kd,i the desorption of species i. Here the kinetics for species i are independent from all other species j 6= i. Langmuir kinetics The single-component Langmuir isotherm was derived for monolayer adsorption on a uniform surface without any interaction of the adsorption sides. It can be generalized to the multi-component

5

case (e.g. Klatt (1999); Schulte and Epping (2005)). In non-equilibrium form, it reads   Nc X q p,j  ψLM,i (cp , qp ) = ka,i qsat,i 1 − cp,i − kd,i qp,i , q sat,j j=1

(5)

see e.g. Ma et al. (1996). In this model, additionally a concentration qsat,i is introduced which describes a maximum saturation concentration for each species. Modified Langmuir kinetics A modified version of the non-equilibrium form of the Langmuir isotherm which arises from a decreasing filtration efficiency model in depth filtration (Litwiniszyn, 1967), is given as   qp,i ψFLM,i (cp , qp ) = ψFLM,i (cp,i , qp,i ) = ka,i 1 − cp,i − kd,i qp,i . (6) qsat,i This is in principle a simplification of the Langmuir kinetics where it is assumed that the adsorption of one species is not depending on the actual adsorbed concentration of all other species (only j = i in the sum in (5)). For a single species, the formulas are equivalent. Note that each of the kinetics above can be transformed to the isotherm formulation by using Eq. (1c). Assume that the temporal change of qp,i is zero and that no additional transport by diffusion is taking place (Dq,i = 0). Then for example in the linear case (Eq. (4)), the Henry isotherm qp,i = kH,i cp,i with the Henry coefficient k can be obtained. kH,i = ka,i d,i The steric mass action (SMA) (Brooks and Cramer, 1992) is not treated here since the inclusion of the salt concentrations changes the mathematical properties of the system of equations. Instead of partial differential equations (PDE), we have a system of partial differential algebraic equations (PDAE) (von Lieres, 2010). Compared to the framework we are considering subsequently, different solution techniques have to be applied for this kind of system. A closer investigation of this adsorption mechanism will be the topic of a future publication.

2.2

Dimensionless and generalized form of the equations

It is possible to derive a general form of the above equations, which holds for both, the static and the dynamic configuration and which is also valid in dimensionless form. A dimensionless form of the equations is quite frequently used in the literature (cf. Michel et al. (2005)). In Bowes and Lenhoff (2011) for both, batch and column, the dimensionless form is derived in a slightly different framework as we consider it here. However, the overall form of the equations is the same. A similar approach for the column is used in Gu et al. (2013). Generic variables are introduced to show that a general variant of the model can be used for both process modes. An overview of the introduced variables is shown in the first column of Table 1. Such variables have to be replaced by the ones in the following four columns for the specific application. For example, the generic time scale α is replaced by either t or τ , in the first case the equations are considered in its dimensioned form, in the second case in dimensionless form, respectively. Thus, when replacing the generic variables, one of the four columns has to be chosen. In the following, wi denotes either the dimensioned or the dimensionless concentration of species i in the column. The variables wp,i and hp,i describe the dimensioned or dimensionless mobile and stationary concentration of T species i in the beads, respectively. Introduce as before the component vectors w = (w1 , . . . , wNc ) , wp = T T (wp,1 , . . . , wp,Nc ) , and hp = (hp,1 , . . . , hp,Nc ) . The coefficients of the model are denoted by ϑ1 , ϑ2 , and ϑ3,i to ϑ5,i and their specific form is specified in Table 1. Introducing new variables β and λ, β ∈ [β0 , β1 ] for the ¯ for the bead, a new time scale α, and performing those transformations or replacing the column and λ ∈ [0, λ] quantities, the following set of state equations can be identified
∂wi ∂wi ∂ 2 wi 1−ε3 ¯ ·) , + ϑ1 − ϑ2 =− ϑ3,i wi − wp,i (·, λ = λ, 2 ¯ ∂α ∂β ∂β ε λ   ∂wp,i 1 ∂ ∂w p,i − ϑ4,i 2 λ2 = −ψ¯i (wp , hp ) , ∂α λ ∂λ ∂λ   1 ∂ ∂hp,i 2 ∂hp,i − ϑ5,i 2 λ = ψ¯i (wp , hp ) , ∂α λ ∂λ ∂λ

6

β ∈ (β0 , β1 ),

(7a)

¯ λ ∈ (0, λ),

β ∈ (β0 , β1 ),

(7b)

¯ λ ∈ (0, λ),

β ∈ (β0 , β1 ).

(7c)

The corresponding boundary conditions read ∂wi (β1 , α) = ϑ1 win,i , ∂β ∂wi ϑ2 (β0 , α) = 0, ∂β

ϑ1 wi (β1 , α) − ϑ2

(8a) (8b)


∂wp,i ¯ α) = ϑ3,i wi (·, α) − wp,i (·, λ = λ, ¯ α) , β ∈ (β0 , β1 ), (·, λ = λ, ∂λ ∂wp,i (·, λ = 0, α) = 0, β ∈ (β0 , β1 ), ϑ4,i ∂λ ∂hp,i ¯ α) = 0, ϑ5,i (·, λ = λ, β ∈ (β0 , β1 ), ∂λ ∂hp,i ϑ5,i (·, λ = 0, α) = 0, β ∈ (β0 , β1 ). ∂λ

ϑ4,i

(8c) (8d) (8e) (8f)

The initial conditions are given as β ∈ (β0 , β1 ), ¯ λ ∈ (0, λ), β ∈ (β0 , β1 ), ¯ λ ∈ (0, λ), β ∈ (β0 , β1 ).

wi (β, 0) = winitial,i , wp,i (·, ·, 0) = 0,

hp,i (·, ·, 0) = 0,

(9a) (9b) (9c)

The actual value of winitial,i depends on the mode of separation (see Table 1). 2.2.1

Adsorption kinetics

In the general case, the form of the kinetics is preserved; however, the coefficients differ also depending on the static or dynamic case. The analogs to the three parameters ka,i , kd,i , and qsat,i are denoted by ϑ6,i to ϑ8,i . The values of the coefficients are given in Table 2. The adsorption kinetics as shown in Eqs. (4) to (6) read: Linear kinetics ψ¯L,i (wp , hp ) = ψ¯L,i (wp,i , hp,i ) = ϑ6,i wp,i − ϑ7,i hp,i

(10)

Langmuir kinetics 

ψ¯LM,i (wp , hp ) = ϑ6,i ϑ8,i 1 −

Nc X hp,j j=1

ϑ8,j



 wp,i − ϑ7,i hp,i

(11)

Modified Langmuir kinetics   hp,i ψ¯FLM,i (wp , hp ) = ψ¯FLM,i (wp,i , hp,i ) = ϑ6,i 1 − wp,i − ϑ7,i hp,i ϑ8,i

(12)

In the following, we work with this general formulation of the GRM. Model simplifications can then be derived from the subsequent formulations. Such simplifications include lumped rate or ideal models, see e.g. Guiochon et al. (2006); Klatt (1999); Michel et al. (2005). However, they are not further discussed in the following.

2.3

Parameter identification and adjoint problem

In this article, we use the principle first optimize then discretize and derive the optimality conditions on the level of partial differential equations. Therefore, we use an adjoint formulation. This is in contrast to Hahn et al. (2014), where first the partial differential equations of the GRM were discretized and afterward the optimality conditions have been derived based on the discrete equations.

7

batch

batch (dimensionless)

α

t

τ=

β

−

λ

r

wi

ci

c¯i =

wp,i

cp,i

c¯p,i =

hp,i

qp,i

q¯p,i =

win,i

t Tend

column t

column (dimensionless) uint t H

τ=

z

ξ=

z H

r rp

r

ρ=

r rp

ci cinitial,i

ci

c¯i =

cp,i cinitial,i

cp,i

c¯p,i =

cp,i cin,i

qp,i cinitial,i

qp,i

q¯p,i =

qp,i cin,i

−

−

cin,i

1

winitial,i

cinitial,i

1

0

0

ϑ1

0

0

uint

1

ϑ2

0

0

Dax

ϑ3,i

kf,i

kf,i Tend k¯f,i = rp

kf,i

ϑ4,i

Dc,i

¯ c,i = Dc,i Tend D rp2

Dc,i

ϑ5,i

Dq,i

¯ q,i = Dq,i Tend D rp2

Dq,i

αend

Tend

1

Tend

β0

−

−

0

0

β1

−

−

H

1

¯ λ

rp

1

rp

1

− ρ=

1 = Pe



ci cin,i

Huint Dax

−1

kf,i H k¯f,i = uint rp  −1 uint 2 1 = r P ecbead,i Dc,i H p  −1 1 uint 2 = r P eqbead,i Dq,i H p τend =

uint Tend H

Table 1: Variables and coefficients for the generalized form of the GRM

8

ϑ6,i

ϑ7,i

ϑ8,i

kinetics

batch

batch (dimensionless)

column

column (dimensionless)

L

ka,i

k¯a,i = ka,i Tend

ka,i

k¯a,i = ka,i uH int

LM

ka,i

k¯a,i = ka,i Tend cinitial,i

ka,i

Hcin,i k¯a,i = ka,i uint

FLM

ka,i

k¯a,i = ka,i Tend

ka,i

k¯a,i = ka,i uH int

L

kd,i

k¯d,i = kd,i Tend

kd,i

k¯d,i = kd,i uH int

LM

kd,i

k¯d,i = kd,i Tend

kd,i

k¯d,i = kd,i uH int

FLM

kd,i

k¯d,i = kd,i Tend

kd,i

k¯d,i = kd,i uH int

L

−

−

−

−

LM

qsat,i

q¯sat,i =

FLM

qsat,i

q¯sat,i =

qsat,i

cinitial,i qsat,i cinitial,i

qsat,i

q¯sat,i =

qsat,i

q¯sat,i =

qsat,i cin,i qsat,i cin,i

Table 2: Coefficients of the adsorption kinetics for the generalized form of the GRM The goal is to minimize the difference between given (experimental) data and the actual simulation, i.e., the numerical solution of the above equations. This is equivalent to the identification of the correct model parameters. In an abstract way, the problem can be formulated as min J(y, µ)

(13a)

s.t. e(y, µ) = 0,

(13b)

(y, µ) ∈ Wad ⊂ Y × U.

(13c)

(y,µ)

The cost functional J(y, µ) describes some desired state of the process and e(y, µ), e : Y × U → Z, are the constraints. In our case, the constraints for the i-th component are given by the set of equations (7), (8), and T (9) combined with one of the kinetics (10) to (12). The considered state variables are y = (w,   wp , hp ) and N

c the control variables µ are defined as the vector of model parameters µ = (ϑ3,i , . . . , ϑ8,i )i=1 ∈ U ⊂ RNp , where Np is the number of parameters and U the space of control variables. For the moment, we omit the parameters ϑ1 and ϑ2 (corresponding to the convective velocity and the axial dispersion) since they are mostly fixed. The space of state variables Y is composed of a triple for the single variables w, wp , and hp . For the single components it holds wi ∈ Wi , wp,i ∈ Wp,i , and hp,i ∈ Hp,i . Depending on the adsorption kinetics, standard spaces for solutions of parabolic PDEs can be chosen, compare e.g. Behrens et al. (2014); Hinze et al. (2008); Tr¨oltzsch (2010). In general, the spaces are Banach spaces. To derive the adjoint system of equations we use the Lagrange formalism (Hinze et al., 2008; Tr¨oltzsch, 2010). We define the Lagrange function L : Y × Z ∗ × U → R as

L (y, z; µ) = J(y, µ) − hz, e(y, µ)iZ ∗ ,Z ,

(14)

and the variable z is the so-called adjoint variable. It can be seen as a generalization of Lagrange multipliers, however, the multipliers are considered as functions and all differential operators are only formulated in a formal way. More details and a more rigorous derivation can be found in Tr¨oltzsch (2010). Here h·, ·iZ ∗ ,Z denotes the duality pairing of the Banach space Z with its dual space Z ∗ . We assume that both, Lagrange function and the cost functional J can be written as the sum of Nc summands for each species i of the mixture L (y, z; µ) =

Nc X i=1

Li (y, z; µ) =

Nc X i=1

  Ji (y, µ) + hz, e(y, µ)iZ ∗ ,Z , i

(15)

where the i-th component of the duality pairing is specified below. For the moment (also for clarity), we stick

9

to the Lagrangian for the i-th component Li . The i-th summand of the Lagrangian function can be written as Li (w,wp , hp , zi ; µ) =Ji (wi , wp,i , hp,i ; µ)  Z αend Z β1 
1 ∂wi 1−ε3 ∂wi ∂ 2 wi ¯ + + ϑ1 − ϑ2 − ¯ ϑ3,i wi − wp,i (·, λ = λ, ·) zi dβ dα ∂α ∂β ∂β 2 ε λ β0 0   Z αend ∂wi zi2 dα − ϑ1 w i − ϑ2 − ϑ1 win,i ∂β 0 β=β1  Z αend  ∂wi zi3 dα − ϑ2 ∂β 0 β=β0 Z β1 {wi (·, 0) − winitial,i } zi4 dβ − β0

− − − − − − −

Z

αend

0

Z

αend

0

Z

αend

β0 β1

Z

β0 Z β1

β1

3 3 ¯ β0 λ Z αend Z 0

Z

αend

0

Z

αend

0

Z

β1

β0

0

Z

Z

Z

¯ λ

0

β1

β0 β1

Z

β0 β1

Z

β0

3 ¯ λ3

Z

¯ λ

0



∂wp,i 1 ∂ − ϑ4,i 2 ∂α λ ∂λ



λ2

∂wp,i ∂λ



 + ψ¯i (wp , hp ) zi5 λ2 dλ dβ dα


∂wp,i ¯ 2 dβ dα ¯ ϑ4,i − ϑ3,i wi − wp,i (·, λ = λ, ·) zi6 λ ∂λ ¯ λ=λ   3 ∂wp,i ϑ zi7 02 dβ dα 4,i ¯3 ∂λ λ λ=0 3 ¯3 λ



(16)

{wp,i (·, ·, 0)} zi8 λ2 dλ dβ

3 ¯ λ3 3 ¯3 λ 3 ¯ λ3

   ∂hp,i 1 ∂ 2 ∂hp,i ¯ − ϑ5,i 2 λ − ψi (wp , hp ) zi9 λ2 dλ dβ dα ∂α λ ∂λ ∂λ 0   ∂hp,i 10 ¯ 2 ϑ5,i ¯ zi λ dβ dα ∂λ λ=λ   ∂hp,i zi11 02 dβ dα ϑ5,i ∂λ λ=0 Z

¯ λ



Z λ¯ 3 12 2 ¯ 3 0 {hp,i (·, ·, 0)} zi λ dλ dβ, β0 λ
T where zi = zi1 , . . . , zi12 is the i-th vector of adjoint variables. Note that adjoint variables zi2 and zi3 are only defined on the boundaries. This choice is made since it is a priori not clear how the adjoint variable zi1 behaves on the boundary (Tr¨oltzsch, 2010). The same holds for the variables on the boundary of the bead and for the initial conditions. The integral terms result from the duality pairing and represent its i-th component. In the following, it is shown that only three of the adjoint variables are further needed and all others are only restrictions (see Appendix A). In the following, we consider the derivatives with respect to the state variables (wi , wp,i , hp,i ) of species i. Note that when choosing, for example, the Langmuir kinetics (11), the i-th summand also depends on all other hp,j , j = 1, . . . , Nc . Thus we will get additional terms from Lj . To retain consistency, we add the terms in the following. Then the optimality conditions are given as −

β1

∂ Li (w, wp , hp , z; µ)w ˜i = 0, ∂wi Nc X ∂ ∂ Li (w, wp , hp , z; µ)w ˜p,i + Lj (w, wp , hp , z; µ)w ˜p,i = 0, ∂wp,i ∂w p,i j=1

∀w ˜i ∈ Wi ,

(17a)

∀w ˜p,i ∈ Wp,i ,

(17b)

˜ p,i ∈ Hp,i . ∀h

(17c)

j6=i

Nc X ∂ ∂ ˜ p,i + ˜ p,i = 0, Li (wi , wp,i , hp,i , zi ; µ)h Lj (w, wp , hp , z; µ)h ∂hp,i ∂h p,i j=1 j6=i

Here the tilde variables denote the directions for the gradients of Li , thus in total, we consider directional derivatives. More details can be found in Bernauer and Herzog (2011). 10

A detailed evaluation of the optimality conditions and the derivation of the adjoint equations can be found in Appendix A. Renaming the adjoint variable such that zi1 = zi , zi5 = zp,i , and zi9 = gp,i , one can identify the adjoint equations as   3 ∂zi ∂zi ∂ 2 zi 1−ε ¯ ·) = ∂Ji (wi , wp,i , hp,i ; µ), − − ϑ1 − ϑ2 2 + ¯ ϑ3,i zi − zp,i (·, λ = λ, (18a) ∂α ∂β ∂β ε ∂wi λ   ∂zp,i 1 ∂ ∂ ψ¯i 2 ∂zp,i − − ϑ4,i 2 λ + (wp , hp ) (zp,i − gp,i ) ∂α λ ∂λ ∂λ ∂wp,i

−

1 ∂ ∂gp,i − ϑ5,i 2 ∂α λ ∂λ



Nc X ∂ ψ¯j ∂Ji + (wp , hp ) (zp,j − gp,j ) = (wi , wp,i , hp,i ; µ), ∂w ∂w p,i p,i j=1

(18b)

i6=j

λ2

∂gp,i ∂λ



+

∂ ψ¯i (wp , hp ) (zp,i − gp,i ) ∂hp,i

Nc X ∂Ji ∂ ψ¯j + (wp , hp ) (zp,j − gp,j ) = (wi , wp,i , hp,i ; µ), ∂h ∂h p,i p,i j=1

(18c)

i6=j

¯ The boundary conditions can be deduced as where β ∈ (β0 , β1 ) and λ ∈ (0, λ). ∂Ji ∂zi (β0 , α) = (wi , wp,i , hp,i ; µ)|β=β0 , ∂β ∂wi ∂zi (β1 , α) = 0, ϑ2 ∂β

−ϑ1 zi (β0 , α) − ϑ2


∂zp,i ¯ α) = ϑ3,i zi (·, α) − zp,i (·, λ = λ, ¯ α) , (·, λ = λ, ∂λ ∂zp,i ϑ4,i (·, λ = 0, α) = 0, ∂λ ∂gp,i ¯ α) = 0, (·, λ = λ, ϑ5,i ∂λ ∂gp,i ϑ5,i (·, λ = 0, α) = 0, ∂λ ϑ4,i

(19a) (19b) β ∈ (β0 , β1 ),

(19c)

β ∈ (β0 , β1 ),

(19d)

β ∈ (β0 , β1 ),

(19e)

β ∈ (β0 , β1 ),

(19f)

and the terminal conditions as ∂Ji (wi , wp,i , hp,i ; µ)|α=αend , β ∈ (β0 , β1 ), ∂wi ¯ zp,i (·, ·, αend ) = 0, λ ∈ (0, λ), β ∈ (β0 , β1 ), ¯ gp,i (·, ·, αend ) = 0, λ ∈ (0, λ), β ∈ (β0 , β1 ). zi (·, αend ) =

(20a) (20b) (20c)

Note that we have derived an end time problem since a condition at final time αend is prescribed. Additionally, the flow direction is reversed since the sign in front of ϑ1 has changed. This also explains why the inlet and outlet condition types are switched. 2.3.1

Adsorption kinetics derivatives

To complete the adjoint system, the derivatives of the adsorption kinetics with respect to the state variables are missing. Based on the formulation of the kinetics in Eqs. (10) to (12), the derivatives can be derived. They are shown in Table 3. The resulting terms can be plugged into Eqs. (18).

2.4

Design derivatives

Looking back to problem (13), we have not taken the parameter dependence into account until now. The  T Nc parameter vector is defined as µ = (ϑ3,i , ϑ4,i , ϑ5,i , ϑ6,i , ϑ7,i , ϑ8,i )i=1 and thus derivatives with respect to these quantities have to be considered. A detailed evaluation of the optimality conditions and the derivation of the design derivatives can be found in Appendix B. Note that the form of the space of control variables U has to be taken into account. Finally, the derivatives can be identified as 11

ψ¯L,i ∂ ∂wp,i ∂ ∂hp,i ∂ ∂wp,j ∂ ∂hp,j

ϑ6,i −ϑ7,i

ψ¯LM,i   Nc X hp,j  ϑ6,i ϑ8,i 1 − ϑ j=1 8,j −ϑ6,i wp,i − ϑ7,i

− −

−ϑ6,i

ψ¯FLM,i   hp,i ϑ6,i 1 − ϑ8,i −ϑ6,i

wp,i − ϑ7,i ϑ8,i

−

−

ϑ8,i wp,i ϑ8,j

−

Table 3: Derivatives of the isotherms with respect to the state variables

∂Li (wi , wp,i , hp,i , zi ; µ) ∂ϑ3,i !    Z αend Z β1 
3 1−ε ∂Ji ¯ ·) ¯ ·) wi − wp,i (·, λ = λ, − zi − zp,i (·, λ = λ, dβ dα , = ¯ ∂ϑ3,i ε λ β0 0 ! Z αend Z β1 Z λ¯ ∂wp,i ∂zp,i 2 ∂Li ∂Ji 3 (wi , wp,i , hp,i , zi ; µ) = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα , ∂ϑ4,i ∂ϑ4,i 0 β0 λ ! Z λ¯ Z αend Z β1 ∂Li ∂hp,i ∂gp,i 2 ∂Ji 3 (wi , wp,i , hp,i , zi ; µ) = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα , ∂ϑ5,i ∂ϑ5,i 0 β0 λ ! Z λ¯ ¯ Z αend Z β1 ∂ ψi ∂Li ∂Ji 3 2 (wi , wp,i , hp,i , zi ; µ) = − ¯ 3 0 ∂ϑ6,i (zp,i − gp,i )λ dλ dβ dα , ∂ϑ6,i ∂ϑ6,i 0 β0 λ ! Z αend Z β1 Z λ¯ ¯ ∂ ψi ∂Ji 3 ∂Li 2 (wi , wp,i , hp,i , zi ; µ) = − ¯ 3 0 ∂ϑ7,i (zp,i − gp,i )λ dλ dβ dα , ∂ϑ7,i ∂ϑ7,i 0 β0 λ ! Z λ¯ ¯ Z αend Z β1 ∂Li ∂ ψi ∂Ji 3 2 (wi , wp,i , hp,i , zi ; µ) = − ¯ 3 0 ∂ϑ8,i (zp,i − gp,i )λ dλ dβ dα . ∂ϑ8,i ∂ϑ8,i 0 β0 λ

(21a) (21b) (21c) (21d) (21e) (21f)

The partial derivatives of the adsorption kinetics with respect to the parameters are given in Table 4. Note that due to the possible coupling in the adsorption kinetics, one has to consider the complete derivative with respect to ϑ8,i . Thus also ∂ϑ∂8,i Lj has to be taken into account.

2.5

Cost functional

For the following applications, we consider the i-th summand of the cost functional in the form Z Z Z θ1 αend β1 θ2 β1 2 2 Ji (wi , wp,i , hp,i ; µ) = (wi − wref,i ) dβ dα + (wi (β, αend ) − wfinal,i ) dβ 2 0 2 β0 β0 Z θ3 αend 2 + (wi (β0 , α) − wout,i ) dα 2 0 Z Z Z λ¯ θ4 αend β1 3 2 2 + ¯ 3 0 (hp,i + ωwp,i − href,i ) λ dλ dβ dα 2 0 β0 λ 2 Z Z ¯ Z hp,i θ5 αend β1 λ ¯ − I λ2 dλ dβ dα + 2 0 hmax,i β0 0 +

8 X γj j=3

2

ϑ2j,i .

12

(22)

ψ¯L,i

ψ¯LM,i 

∂ ∂ϑ6,i

wp,i

∂ ∂ϑ7,i

−hp,i

∂ ∂ϑ8,i

−

ϑ8,i 1 − 

ψ¯FLM,i

Nc X hp,j j=1

ϑ8,j

−hp,i



 wp,i



hp,i 1− ϑ8,i



wp,i

−hp,i



Nc X  hp,j   wp,i ϑ6,i  1 −  ϑ  j=1 8,j

ϑ6,i

wp,i hp,i ϑ28,i

j6=i

∂ ∂ϑ8,j

−

ϑ6,i

ϑ8,i hp,j wp,i ϑ8,j ϑ8,j

−

Table 4: Derivatives of the isotherms with respect to the parameters The cost functional consists of different parts. The first part measures the difference to data wref,i given in the column (or the batch) over the whole time horizon. The second part rates data wfinal,i given at final time αend . The third term compares to the classical breakthrough curve. Note that for the batch case, this term is coinciding with the first part since the interval (β0 , β1 ) degenerates to a single point. The fourth part accounts for concentration profiles given within the bead. It additionally incorporates a switching factor ω. This factor accounts for the fact that when experimentally measuring the concentration within the bead, it can often not distinguished between concentration which is mobile and already stationary. In most of the cases, the stationary concentration hp,i is dominating wp,i (for example in strong binding conditions) and ω can be set to ω = 0. The fifth term is particularly designed for the comparison to concentration profiles in the beads measured from CLSM images. However, it is difficult to relate the resulting intensity value to an actual concentration value. In Yang et al. (2008, 2006) a method based on equilibrium data is proposed to correlate intensity and adsorbed concentration. We suppose that the maximum of the intensity corresponds to the maximum of the concentration and relate the scaled properties. Therefore, we scale to the maximum of the adsorbed concentration and set hmax,i = maxα∈(0,αend ) maxλ∈(0,λ) ¯ hp,i (·, λ, α). It is compared to a scaled intensity I¯ which is relating the original intensity to its maximum. There are also other scaling variants possible. The last term is a regularization term (Tikhonov regularization) with regularization parameters (weights) γj,i . The derivatives of the cost functional, which are used in (18) to (20) and in (21) can be computed as ∂Ji w ˜i =θ1 ∂wi

Z

αend

0

+ θ3 ∂Ji w ˜p,i =θ4 ∂wp,i ∂Ji ˜ hp,i =θ4 ∂hp,i

Z

Z

β1

(wi − wref,i ) w ˜i dβ dα + θ2

β0

αend

0

αend

0

Z

Z

Z

0

+ θ5

Z

Z

β1

β1

β0

αend

0

∂Ji ˜ ϑj,i =γj ϑj,i ϑ˜j,i , ∂ϑj,i

β1

β0

(wi (β, αend ) − wfinal,i (β)) w ˜i (β, αend ) dβ

(wi (β0 , α) − wout,i (α)) w ˜i (β0 , α) dα

β0

αend

Z

Z

3 ¯ λ3 3 ¯3 λ β1

β0

Z

¯ λ

(hp,i + ωwp,i − href,i ) ω w ˜p,i λ2 dλ dβ dα

0

Z

Z

¯ λ

0 ¯ λ

0

(23a)

˜ p,i λ2 dλ dβ dα (hp,i + ωwp,i − href,i ) h !   ˜ p,i hp,i h h p,i ˜ max,i − I¯ −h λ2 dλ dβ dα hmax,i hmax,i h2max,i

j = 3, . . . , 8.

(23b)

(23c) (23d)

The last term in (23c) is a little problematic since the variation of hmax,i is required. For the moment we neglect ˜ max,i . this and assume the variation of hmax,i to be small. So we neglect the term corresponding to h

13

3

Numerics

Since no analytical solutions for the general rate model are available, the system of equations has to be solved numerically. For this purpose, the equations have to be discretized in both, space and time, to derive a numerical solution.

3.1

Discretization

The system of equations in generalized form (7) can be discretized in different ways. In Hahn et al. (2014) finite elements were used. We use the finite volume method. This is similar to the elaborations in von Lieres and Andersson (2010). However, we use for the discretization of the convective part a flux limiter scheme instead of the WENO scheme. We choose the flux limiter superbee (cf. LeVeque (2002)). The same discretization is also used for the adjoint system, however, note that the convective part has the opposite sign and therefore the discretization of this part is slightly different.

3.2

Time integration

The time integration is split in an explicit and an implicit part and the Euler scheme is used. The convective part is due to its discretization with a flux limiter scheme tackled explicitly. If the kinetic is nonlinear, the nonlinear part is also tackled explicitly (for hp,i ) and all other parts are tackled implicitly. Note that the adjoint problem is an end-value problem. To use the established methods, the transformation α ˜ = αend −α is used which reformulates the system of equations (18) to (20) as an initial value problem. Then the same time discretization scheme as for the state system can be used.

3.3

Optimization

For the solution of the optimization problem (13) a descent method is used. In Zhang et al. (2017) a gradient descent method based on Barzilai-Borwein step size is used. The authors in Hahn et al. (2014) claim to use a Newton-Simpson method. Since we derived the optimality conditions from the Lagrange function, we actually investigated the reduced problem ˆ (24a) min J(µ) µ∈U

ˆad = {µ ∈ U : (y(µ), µ) ∈ Wad } . s.t. µ ∈ U

(24b)

Here Jˆ is called the reduced cost functional and is depending on the control variables µ only. Comparing to problem (13) note that the constraints have vanished, however, they are only moved to the set of admissible ˆad . Details on the theory and the assumptions leading to this result can be found in Hinze control variables U et al. (2008); Tr¨ oltzsch (2010). ˆ we are looking for a direction d, such that J(µ ˆ + sd) ≤ J(µ) ˆ To minimize J, for a suitable step size s > 0. The direction d is called descent direction. Thus, a general descent method can be formulated and the pseudo-code is given in Algorithm 1 (Geiger and Kanzow, 1999a). Algorithm 1 General descent method Input: µ0 , k = 0 while stopping criterion not met do determine descent direction dk ˆ k + sk dk ) ≤ J(µ ˆ k) determine step size sk > 0 such that J(µ k+1 k k k Update µ =µ +s d , k =k+1 end while Looking at the algorithm, there are three more points to specify, the stopping criterion, the computation of the descent direction, and the computation of the step size. The most obvious, but not the best choice for a descent direction is the negative gradient, the so-called steepest descent direction !Np ˆ ∂ J k k k ˆ )=− d = −∇J(µ (µ ) . (25) ∂µi i=1

14

It can be shown that this direction always points in the direction of a minimum (Nocedal and Wright, 2006c). This method basically uses information from the first derivative of the reduced cost functional. For the computation of the gradient of our reduced cost functional, it can be shown that the derivatives of L with respect to the parameters given in Eqs. (21) are equal to the gradient ∇Jˆ (Hinze et al., 2008). If additionally information on the second derivatives is available, Newton methods can be used. In this case, the cost functional is approximated by a second-order Taylor series and the exact minimizer of this approximation can be computed. However, the computation of the Hessian (second derivatives) is often very expensive or not possible. Therefore, quasi-Newton methods, where the Hessian is approximated, are frequently used. In the following, we take a closer look at one of those methods, the so-called BFGS (Broyden-Fletcher–Goldfarb-Shanno) method. For the method and the details of the following see Nocedal and Wright (2006d). Approximate the cost functional in the k-th step with a quadratic function mk as ˆ k ) + ∇J(µ ˆ k )T d + 1 dT B k d, mk (d) = J(µ 2

(26)

where B k is a symmetric positive definite matrix approximating the Hessian. It can be shown that the minimizer dk of this function is given as
−1 ˆ k ). dk = − B k ∇J(µ (27) Then as in Algorithm 1, an update of the iterate with a suitable step size sk can be performed µk+1 = µk + sk dk .

(28)

To compute the approximation of the Hessian, one iteratively updates its form by incorporating actual knowledge
−1 on the objective function and its derivative. To avoid the inversion, one usually defines H k = B k . Then the following update formula can be derived 
T  k 
T 
T H k+1 = I − ρk pk f k H I − ρk f k pk + ρk pk pk , (29) ˆ k+1 ) − ∇J(µ ˆ k ), and ρk = where I is the identity matrix, pk = µk+1 − µk = sk dk , f k = ∇J(µ 0



fk


T

pk

−1

.

The algorithm is mostly started with H = I, i.e., one step in steepest descent direction. For computing the step size sk we use a line search method, namely a modified Armijo rule, where both a decrease and an increase of the initial guess for the step size is taken into account (Geiger and Kanzow, 1999b). As pointed out in Nocedal and Wright (2006d), the use of the Wolfe condition in combination with the BFGSalgorithm is recommended, however, its evaluation requires in every step the computation of the gradient and makes therefore the linear search much more computationally expensive since the adjoint system has to be
T solved. The Armijo rule does not guarantee to satisfy the curvature condition f k pk > 0, and therefore, the approximation of the Hessian is not ensured to be positive definite. If in our computations the curvature condition is not satisfied, we reset H k to H 0 and perform again a steepest descent step. For the optimization algorithm and especially for the step size control it is challenging when the variables are of different magnitude. Such problems are called poorly scaled (Nocedal and Wright, 2006c). If the descent direction value is much larger than the actual value of the variable, a rather small step size has to be chosen. This slows down the algorithm. Therefore, we decompose the variables in its magnitude and a number µ ¯ of magnitude 0, i.e., µi = µ ¯i · 10blog10 (µi )c , (30)

where bxc denotes the floor function yielding the largest integer number less than x. This variant is known as diagonal scaling. For physical properties, this is in principle equivalent to the use of different units. Since the steepest descent algorithm is sensitive to poor scaling, as a kind of fallback strategy a coordinate descent direction is incorporated (Nocedal and Wright, 2006b). Here it is easier to adjust the step size to the current magnitude of the variable. To terminate the algorithm, different stopping criteria are considered. First of all, we investigate the magnitude of the gradient and the relative gradient (comparison to the first iterate). If those values are below a given tolerance, the algorithm stops. Furthermore, we prescribe a maximum number of steps and additionally check if the change in the variables and the objective function in the last step is sufficiently small. Note that also the scaling has to be taken into account here.

15

4

Material and methods

To show the possibilities of the parameter identification method, we compare it to a set of experiments. Since our previous considerations are valid for both, batch and column, experiments for the two configurations are performed.

4.1

Batch experiment

R Monoclonal antibody (mAb05) is used to investigate the protein adsorption of Eshmuno A chromatographic media (bead material) of Merck KGaA (Darmstadt, Germany). The spherical beads consist of a rigid polymeric matrix with a large inner pore system in which protein A ligand is immobilized. The mean diameter of the beads is given as 50 µm. The protein solution was prepared as follows: mAb05 was diluted in phosphate-buffered R saline (PBS) to a final concentration of 1.56 mg mL−1 . In each well of a 96 well MultiScreen HTS DV filter plate, 190 µL protein solution and 10 µL of a 50 % (v/v) resin slurry were dispensed. The media was equilibrated with the corresponding buffer beforehand. The filter plate was incubated for different incubation times on an R orbital shaker at 1100 rpm. Liquid was removed by vacuum filtration using MultiScreen HTS vacuum manifold. Part of the antibody was labeled with a fluorescent dye and the degree of labeling (DOL) was 1:1. This enables the visualization of the fluorescent intensity profile across the radius of a single chromatographic bead using CLSM. The protein concentration was calculated based on Lambert-Beer law using flow-through data with UV measurements at 280 nm. Additionally, the protein adsorption was corrected measuring the absorbance at the specific excitation maximum of the conjugated dye. The experimental results are reported in Rammo et al. (2017). Here, fluorescent intensity profiles across the chromatographic beads were measured for six different time instances and the corresponding protein concentrations in the protein solution were calculated. This allows the application of the methods described above to fit either the concentration data, corresponding to θ1 term in (22), or the intensity data, corresponding to θ5 term in (22), or a combination of both.

4.2

Estimation of the concentration profiles

The concentration profiles of the proteins in the spherical beads were estimated based on CLSM. In the attenuation correction for the data used in this article to model the purification processes, we essentially followed the approach of Stanislawski et al. (2010). The idea of this approach is the description of attenuation by the Lambert-Beer law, where scattering was neglected and the light absorption is assumed to be proportional to the fluorophore density. This leads to the numerical solution of an ordinary first-order differential equation. As an extension of the approach of Stanislawski et al. (2010), the beam divergence was taken into account and the excitation beam was modeled as a Gaussian beam, i.e. as the transverse electromagnetic mode of the solution to the paraxial Helmholtz equation. For the emitted beam, we assumed isotropic fluorescence irradiation with a divergence restricted by the pinhole size. The core of the attenuation correction is a linear filtering of the CLSM image with a beam intensity obtained as the convolution of excitation and emission, which can efficiently be done via frequency space (Roerdink and Bakker, 1993). Finally, we also refer to Strasters et al. (1994) involving the correction of scattering effects and presenting a method for fast iterative image restoration. Since strong attenuation leads to a loss of information on the true fluorophore density, the quality of the restored images depends on the attenuation itself and, as a consequence, the quality decreases with depths, see Fahrbach et al. (2019) for further details. For the data acquisition discussed in this paper, we used the Leica CLSM TCS SP3 with 40×/1.25 NA objective of type HCX PL APO CS (oil, refraction index 1.518), which yields a maximum half angle of 55.4◦ of the light cone. The pinhole size was 159 µm and the wavelength of emission beam was about 780 nm. The CLSM images consisted of about 200 slices, each with 512 × 512 pixels with a size of 0.2 µm (xy-resolution of about 0.37 µm). The pixel size in z-direction (i.e. the slice distance) was approximately the same as the xy-resolution. The concentration profiles of the spherical beads were estimated as the mean radial function of the relative fluorophore density, where the mean was taken over all directions on the upper half unit sphere. Examples for the resulting concentration profiles can be found in Figure 3. Further investigations were carried out using a Leica CLSM TCS SP8 with a 40×/1.10 NA objective of type HC PL APO CS2 (water, refraction index 1.33), where the pixel size was again 0.2 µm.

16

1 initial iteration 1 iteration 2 iteration 3 iteration 4 iteration 5 iteration 6 exact solution

0.7 0.6

normalized concentration (-)

normalized concentration (-)

0.8

0.5 0.4 0.3 0.2

0

1000

2000

3000

1 initial iteration 1 iteration 3 iteration 5 iteration 7 iteration 9 exact solution

0.9

4000

time (s)

0.8 0.7 0.6 0.5 0.4 0.3 0.2

0

1000

2000

time (s)

3000

initial iteration 1 iteration 2 iteration 3 iteration 4 iteration 5 exact solution

0.9

normalized concentration (-)

1 0.9

4000

0.8 0.7 0.6 0.5 0.4 0.3 0.2

0

1000

2000

3000

4000

time (s)

(a) Numerical test using Langmuir ki- (b) Numerical test using linear kinetics (c) Numerical test using FLM kinetics netics

Figure 1: Approximation of the concentration in the batch in different iterations of the descent method

4.3

Column experiment

R A cylindrical column with an inner diameter of 10 mm was packed with Eshmuno CPX cation exchanger resin (Merck KGaA, Darmstadt, Germany, art. no. 1.200 83) to a final bed height of 40 mm. The mean diameter of the resin particles was given as 50 µm, with an ionic capacity of 70 µeq mL−1 (Stone et al., 2019), and the estimated total porosity of the column was 74.2 %, according to iSEC data obtained from Merck. The ¨ column was connected to an AKTA Avant chromatographic workstation (GE Healthcare). Polyclonal human immunoglobulin G (Gammanorm, Octapharma, Austria), diluted with binding buffer to a final concentration of 5 mg mL−1 was loaded on the column with a constant flow velocity of 0.67 mm s−1 (feed flow rate 52.22 mm3 s−1 ), which corresponds to a residence time of 1 min. The protein concentration at the outlet of the column was ¨ monitored using the UV detector of the AKTA system.

5

Results

For illustrating the method, we stick to one component (or rather describe the mixture with one concentration), i.e., Nc = 1. Thus, we neglect in the following the index i. The experimental data is given at discrete time instances only. However, the formulation of the cost functional in equation (22) requires the data over the complete time trajectory. Therefore, the data are interpolated by a piecewise linear function. If not stated otherwise, the weights θj and the regularization parameters γj , j = 1, . . . , 6, are set to zero. The discretization is performed as described in section 3. The radius of the bead is discretized in 16 finite shells (one-dimensional finite volumes in spherical coordinates) and the column height is divided into 10 finite volumes. For the time discretization, a time step of 0.12 s is used. For the column experiment, the time step was slightly smaller (0.1 s). The final time Tend is chosen as 3600 s, for the numerical test of the column it is only 1200 s. However, note that due to the use of non-dimensional equations, the time step and the final time is adjusted according to Table 1.

5.1

Numerical test batch

First, we test the proposed method by means of a numerical test. We start with the static mode. For a given set of parameters µtest , the set of equations (7) to (9) is solved and the data is collected. Then, by setting θ1 = θ4 = 0.5 in Eq. (22) with wref and href resulting from the previous simulation with µ = µtest , we test our method. The initial guess is generated by perturbing the parameter configuration µtest using uniformly distributed random numbers η such that µ0 = µtest · (1 + η). In figure 1 the results for the three different adsorption kinetics are shown. We only present the concentration in the batch and skip the profiles of the concentration within the bead. It can be observed that the method converges for all three adsorption kinetics within a few iterations to the desired solution. The last iteration is in all three cases a coordinate descent step.

17

5.2

Batch experiment

The incubation experiment yields the concentration profiles in the batch at different time instances and fluorescence intensity profiles within the bead from CLSM measurements. Since the front at the outer radius of the bead is not sharp, first the radius of the measured bead is determined. For this purpose, the relative volumetric bead capacity Qrel as defined in Hubbuch et al. (2002) is used. The maximum of Qrel is taken as an indicator for the outer radius of the bead. Afterward, the radial direction is normalized to the interval [0, 1]. For the comparison of different beads and of different time instances, these values are mapped to the reference radius R of the bead as specified by the manufacturer, in the case of Eshmuno A, rp = 25 µm. As a next step, the intensity data is normalized using the maximum intensity value of the data series. For the parameter identification, it is important to balance the two ingredients (θ1 and θ5 ) of the cost functional. Therefore, we present in the following three different combinations of the two values and the corresponding results, variant 1: θ1 = 1 and θ5 = 0, variant 2: θ1 = 0 and θ5 = 1, and variant 3: θ1 = 0.5 and θ5 = 0.5. As adsorption kinetic, the Langmuir type (11) for one species is used. The concentration profile in the batch is shown in Figure 2. Note that in the figure only every 2000-th point is marked to avoid overloading the lines. It can be seen that for variant 1 the experimental and the simulated curve are quite close to each other. For variant 2, the deviation is quite large, however, this is no surprise since the concentration data was not used for parameter identification (θ1 = 0). The third variant provides a reasonable fit to the experimental curve.

1

normalized concentration (-)

1 1

0.8

1

=1, =0, =0.5,

5 5

=0 =1 5

=0.5

experiment

0.6

0.4

0.2

0

0

1000

2000

3000

4000

time (s) Figure 2: Normalized concentration over time In Figure 3, the comparison to the intensity data from CLSM measurements at different time instances is shown. For a short time, all three variants coincide quite well with the experimental data. In the course of time, variant 1 deviates more and more from the experimental data, whereas the other two variants are still close. This can be explained in the same way as before for variant 2. Since variant 1 does not know about these intensity curves (θ5 = 0), there is, from a mathematical point of view, no reason to provide a reasonable fit. At later times, variant 2 provides the best fit and the concentration values of variant 3 are too low. As expected, variant 1 and variant 2 provide a good fit for either the concentration or the intensity data. Variant 3 (as a combination of variant 1 and 2) yields a reasonable fit for the complete set of data. By adjusting the weights θ1 and θ5 , the result can be improved in the desired direction. The used parameter values are shown in Table 5. One can see that the dominating transport mechanism is pore diffusion.

18

=1,

5

=0, 1

0.8

1

=0

=1 5

=0.5,

5

=0.5

experiment

0.6

0.4

0.2

0

0.5

1

1.5

2

bead radius (m)

2.5

=1,

0.8

1

=1 5

=0.5,

5

=0.5

0.6

0.4

0.2

0

0

0.5

1

1.5

2

bead radius (m)

10 -5

(a) 1 minute

2.5

1

1 1

0.8

1

=1,

5

=0,

5

=0.5,

1

=1 =0.5

experiment

0.6

0.4

0.2

0

0.5

1

=1,

1.5

2

bead radius (m)

2.5 10

5

=0.5,

=0 =1 5

=0.5

experiment

0.6

0.4

0.2

0

0

0.5

1

1.5

2

bead radius (m)

10 -5

2.5 10 -5

(c) 5 minutes 1

1 1

0.8

1

=1,

5

=0,

5

=0.5,

=0 =1 5

=0.5

experiment

0.6

0.4

0.2

0

5

=0,

(b) 2 minutes

=0

5

1

0.8

1

normalized concentration/intensity (-)

normalized concentration/intensity (-)

=0

experiment

1

0

5

=0, 1

0

-5

(d) 10 minutes

0.5

1

1.5

bead radius (m)

2

2.5 10

normalized concentration/intensity (-)

0

1 1

normalized concentration/intensity (-)

1 1

normalized concentration/intensity (-)

normalized concentration/intensity (-)

1

1 1

0.8

1

=1,

5

=0,

5

=0.5,

=0 =1 5

=0.5

experiment

0.6

0.4

0.2

0

0

-5

(e) 30 minutes

0.5

1

1.5

2

bead radius (m)

2.5 10 -5

(f) 60 minutes

Figure 3: Normalized concentration within the bead in comparison to normalized intensity values for different variants of the cost functional

batch experiment

column experiment

θ1 = 1, θ5 = 0

θ1 = 0, θ5 = 1

θ1 = 0.5, θ5 = 0.5

θ3 = 1

LM

LM

LM

FLM

2.26 · 10−7

1.90 · 10−7

1.85 · 10−7

2.23 · 10−5

isotherm
kf m s−1
Dc m2 s−1
Dq m2 s−1

1.33 · 10−14

3.61 · 10−14

1.24 · 10−13

9.43 · 10−12

90.12 kg m−3 s−1

1.83 kg m−3 s−1

10 kg m−3 s−1

10 s−1

kd s−1

0.55

0.35

3.15

0.065

11.10

2.43

7.16

660.57

164.83 kg m−3

5.21 kg m−3

3.17 kg m−3

153.65

ka




−3

qsat kg m KH




8.35 · 10−17

2.71 · 10−20

1.65 · 10−20

1.05 · 10−12

Table 5: Parameter values for the comparison to experimental data (compare Figure 2 and Figure 5)

19

5.3

Numerical test column

As for the batch, we first test the column formulation with synthetic data. The FLM kinetics (12) is used. Therefore, we fix the parameters at µtest and then start at an arbitrary initial guess. The approximation of the breakthrough curve in different iterations can be found in Figure 4. It can be seen, that the method approximates the given breakthrough curve very well. For an exact fit, more iterations are necessary. However, the descent is slow since the solution is already near its optimum.

normalized concentration (-)

1

0.8

0.6 initial iteration 1 iteration 2 iteration 3 iteration 4 iteration 5 iteration 6 exact solution

0.4

0.2

0

0

200

400

600

800

1000

1200

time (s) Figure 4: Approximation of the breakthrough curve in different iterations of the descent method

5.4

Column experiment

For the comparison to an experimental breakthrough curve, the weight θ3 = 1 is chosen (all other equal to 0). The FLM kinetics are used. In Figure 5 the comparison of the simulation to a breakthrough curve is shown. It is clearly visible that the simulation can provide a very good fit to the given data. The parameter values are shown in Table 5.

6

Conclusions and outlook

We derived an adjoint formulation for the general rate model of chromatography. With this formulation, it is possible to compute the exact gradient with one additional solution of a set of partial differential equations. This enables optimization and parameter identification via a descent method. The choice of a rather general cost functional enables to prescribe different scenarios for parameter identification. Due to the formulation of the optimality conditions on the PDE level, the formulation is not depending on the actual discretization of the PDE and different linearization methods and discretization schemes can be used for the state and the adjoint system. Comparing to the approximation of the gradient by sensitivities, only one additional PDE solve instead of Np is necessary to compute the gradient. Due to the highly different magnitude of the model parameters, scaling of the variables is an essential tool. From the numerical experiments above, it can be seen that the method works well for both, synthetic and experimental data. It is shown that the method can also be applied to data resulting from CLSM measurements by correlating the maximum intensity to the maximum of concentration in the bead. Although we had developed a model based on static affinity data, we think that the adaption of the model to cation exchange resins under 20

normalized concentration (-)

1 experiment simulation

0.8 0.6 0.4 0.2 0 -0.2

0

50

100

150

200

volume (ml) Figure 5: Comparison of experimental and simulated breakthrough curve, only every 2000-th point is indicated with a marker dynamic conditions is acceptable as long as the pore structure of the resins is comparable and the experimental conditions allow for a strong binding in each case. The fact that the column experiment was described precisely within a few iteration processes supports this assumption. As a next step, the model needs to be extended to steric mass action (SMA). From a mathematical point of view, as already mentioned, this includes the shift from partial differential equations (PDE) to partial differential algebraic equations (PDAE). For future CLSM measurements, a newer version of the microscope is used.

Acknowledgments The research leading to these results was funded by BMBF-grant AMSCHA (no. 05M2016). The authors would like to thank the anonymous reviewers for their comments that helped to improve the manuscript and the mathematical formulation.

21

A

Formal derivation of the adjoint equations

Starting with Eq. (16), use partial integration for all variables α, β, and λ, to bring the derivatives to the adjoint variables. Li (w,wp , hp , zi ; µ) =Ji (wi , wp,i , hp,i ; µ)  Z αend Z β1  ∂z 1 1−ε3 ∂z 1 ∂z 1 1 wi − i − ϑ1 i − ϑ2 i2 + − ϑ z ¯ 3,i i dβ dα ∂α ∂β ∂β ε λ β0 0  β Z αend  Z β1 ∂wi 1 ∂zi1 1 1 1 αend ϑ1 wi zi − ϑ2 wi zi 0 dβ − dα − z + ϑ2 w i ∂β i ∂β β0 0 β0 Z αend Z β1 1−ε3 1 ¯ − − ¯ ϑ3,i wp,i (·, λ = λ, ·)zi dβ dα ε λ β0 0  Z αend  ∂wi zi2 dα − ϑ1 wi − ϑ2 − ϑ1 win,i ∂β 0 β=β1  Z αend  ∂wi 3 z dα − ϑ2 ∂β β=β0 i 0 Z β1 − {wi (·, 0) − winitial,i } zi4 dβ β0

− −

Z

αend

0

Z

Z

β1

0

Z

αend

0

Z

αend

Z

β1

0

− −

αend

0

Z

Z

¯ λ

β1

β1

Z

¯ λ

0 β1

β0 β1

Z

β0

β1

β0

3 ¯ λ3

Z

0

3 ¯ λ3 3 ¯ λ3

Z

¯ λ

0

Z

¯ λ

0

   1 ∂ ∂z 5 ∂z 5 λ2 i λ2 dλ dβ dα wp,i − i − ϑ4,i 2 ∂α λ ∂λ ∂λ ψ¯i (wp , hp )zi5 λ2 dλ dβ dα



α wp,i zi5 0 end λ2

3 ¯3 λ



dλ dβ −

Z

αend

0

Z

β1

3 ¯ λ3

β0



∂wp,i 5 zi −ϑ4,i λ2


¯ wi − wp,i (·, λ = λ, ·)

∂wp,i − ϑ3,i ∂λ   3 ∂wp,i zi7 02 dβ dα ϑ 4,i ¯3 ∂λ λ λ=0 ϑ4,i

∂λ

¯ λ=λ

 λ¯ dβ dα + ϑ4,i wp,i λ ∂λ 0 5 2 ∂zi

¯ 2 dβ dα zi6 λ

{wp,i (·, ·, 0)} zi8 λ2 dλ dβ

0

β0

3 − 3 ¯ β0 λ Z αend Z − Z

0 β1

β0

0

Z

¯ λ

β0

β1

0

−

Z

β0 Z β1

3 − 3 ¯ β0 λ Z αend Z − Z

β1

β0

3 − 3 ¯ β0 λ Z αend Z − −

β1

β0

αend

0

Z

Z

¯ λ

3 ¯ λ3 3 ¯ λ3

Z

¯ λ

hp,i

0

Z

0

¯ λ



∂z 9 1 ∂ − i − ϑ5,i 2 ∂α λ ∂λ

3 ¯ λ3

dλ dβ −

Z

αend

0

Z

λ

β1

β0

 ∂hp,i 10 ¯ 2 ϑ5,i ¯ zi λ dβ dα ∂λ λ=λ   ∂hp,i ϑ5,i zi11 02 dβ dα ∂λ λ=0 

9 2 ∂zi

∂λ



λ2 dλ dβ dα

−ψ¯i (wp , hp )zi5 λ2 dλ dβ dα

α hp,i zi9 0 end λ2

3 ¯3 λ



3 ¯ λ3



∂hp,i 9 −ϑ5,i λ2 zi ∂λ

 λ¯ dβ dα + ϑ5,i hp,i λ ∂λ 0 9 2 ∂zi

{hp,i (·, ·, 0)} zi12 λ2 dλ dβ

Note that the boundary conditions at λ = 0 are kept to obtain a valid boundary condition for the adjoint system there. Additionally, new terms arise from the partial integration above. In the further course of the derivation, the resulting conditions are used to connect the adjoint variables at the boundary with the adjoint 22

variables inside the beads. The terms can be reordered. Collect all terms with the same integration range and group them according to the variables wi , wp,i , and hp,i . Li (w,wp , hp , zi ; µ) =Ji (wi , wp,i , hp,i ; µ)   Z αend Z β1  ∂zi1 3 ∂zi1 ∂zi1 1−ε 1 6 wi − − − ϑ1 − ϑ2 2 + ¯ ϑ3,i z − zi dβ dα ∂α ∂β ∂β ε i λ β0 0   Z αend ∂z 1 − wi (β1 , α) ϑ1 zi1 (β1 , α) + ϑ2 i (β1 , α) + ϑ1 zi2 dα ∂β 0   Z αend ∂z 1 − wi (β0 , α) −ϑ1 zi1 (β0 , α) − ϑ2 i (β0 , α) dα ∂β Z0 αend  ∂wi − (β1 , α) −ϑ2 (zi1 (β1 , α) + zi2 ) dα ∂β Z0 αend  ∂wi − (β0 , α) ϑ2 (zi1 (β0 , α) + zi3 ) dα ∂β 0 Z β1 Z β1   − wi (β, αend ) zi1 (β, αend ) dβ − wi (β, 0) −zi1 (β, 0) + zi4 dβ β0 αend

− − − − − − − − − − −

Z

0

Z

αend

0

Z

αend

αend

0

Z

αend

0

Z

αend

0

Z

αend

0

Z

Z

3 ¯3 λ

β1

β0 β1

Z

β0 Z β1 β0 β1

Z

3 ¯3 λ

3 ¯3 λ

Z

0

Z

wp,i

0

Z

0

¯ λ

β0

¯ λ

0

β1

β0 Z β1

winitial,i zi4 dβ

∂z 5 1 ∂ − i − ϑ4,i 2 ∂α λ ∂λ

  5 2 ∂zi λ λ2 dλ dβ dα ∂λ

ψ¯i (wp , hp )zi5 λ2 dλ dβ dα



ϑ4,i λ



5 2 ∂zi

∂λ

+

ϑ3,i (zi6

−



zi1 )

 dβ dα ϑ4,i λ ∂λ λ=0 5 2 ∂zi

dβ dα ¯ λ=λ

∂wp,i 3  5 7 2 dβ dα λ=0 ¯ 3 ϑ4,i (zi + zi )λ ∂λ λ wp,i zi5 α=α

end

β1

β0

β1

∂wp,i 3  −ϑ4,i (zi5 − zi6 ) λ=λ¯ dβ dα ¯ ∂λ λ

¯ λ

β0

αend



¯ λ

3 wp,i ¯ 3 λ

0

β1

0

Z

Z

3 wp,i ¯ λ

β0 β1

Z

3 3 ¯ λ β0 Z αend Z Z

β1

β0

β1

β0

Z

Z

β0

0

Z

−ϑ1 win,i zi2 dα −

Z

λ2 dλ dβ

 wp,i −zi5 + zi8 α=0 λ2 dλ dβ

3 ¯3 λ 3 ¯ λ3

Z

¯ λ

0

Z

0

¯ λ

   ∂z 9 ∂z 9 1 ∂ λ2 i λ2 dλ dβ dα hp,i − i − ϑ5,i 2 ∂α λ ∂λ ∂λ −ψ¯i (wp , hp )zi5 λ2 dλ dβ dα

 ∂zi9 − ϑ5,i dβ dα ∂λ λ=λ¯ 0 β0   Z αend Z β1 3 ∂z 9 − hp,i ¯ 3 ϑ5,i λ2 i dβ dα ∂λ λ 0 β0 λ=0 Z αend Z β1 ∂hp,i 3  − −ϑ5,i (zi9 − zi10 ) λ=λ¯ dβ dα ¯ ∂λ λ 0 β Z

αend

3 hp,i ¯ λ



0

23

− − −

Z

αend

β1

β0

Z

β1

β0

0

Z

Z

β1

β0

3 ¯ λ3 3 ¯ λ3

Z

¯ λ

0

Z

0

¯ λ

∂hp,i 3  ϑ5,i (zi9 + zi11 )λ2 λ=0 dβ dα 3 ¯ ∂λ λ hp,i zi9 α=α

end

λ2 dλ dβ

 hp,i −zi9 + zi12 α=0 λ2 dλ dβ

Computing the optimality condition as in (17a), one gets

∂ ∂ Li (w, wp , hp , zi ; µ)w ˜i = Ji (wi , wp,i , hp,i ; µ)w ˜i ∂wi ∂wi    Z αend Z β1 ∂zi1 3 ∂zi1 ∂zi1 1−ε 1 6 w ˜i − − − ϑ1 − ϑ2 2 + ¯ ϑ3,i z − zi dβ dα ∂α ∂β ∂β ε i λ β0 0   Z αend ∂z 1 − w ˜i (β1 , α) ϑ1 zi1 (β1 , α) + ϑ2 i (β1 , α) + ϑ1 zi2 dα ∂β 0   Z αend ∂z 1 − w ˜i (β0 , α) −ϑ1 zi1 (β0 , α) − ϑ2 i (β0 , α) dα ∂β Z0 αend  ∂w ˜i − (β1 , α) −ϑ2 (zi1 (β1 , α) + zi2 ) dα ∂β Z0 αend  ∂w ˜i (β0 , α) ϑ2 (zi1 (β0 , α) + zi3 ) dα − ∂β 0 Z β1 Z β1   − w ˜i (β, αend ) zi1 (β, αend ) dβ − w ˜i (β, 0) −zi1 (β, 0) + zi4 dβ. β0

β0

By choosing proper directions of variation w ˜i (see e.g. Tr¨oltzsch (2010, p.86 f.)), the equations (18a), (19a), (19b), and (20a) follow. This includes the conditions zi1 (β1 , α) = −zi2 , zi1 (β0 , α) = −zi3 , and zi1 (β, 0) = zi4 . Analogously, the optimality conditions (17b) and (17c) can be evaluated. In the first case this leads to the ¯ α) = z 6 , z 5 (β, λ = 0, α) = −z 7 , and z 5 (β, λ, 0) = z 8 . equations (18b), (19c), (19d), and (20b) with zi5 (β, λ = λ, i i i i i ¯ α) = z 10 , z 9 (β, λ = 0, α) = −z 11 , and In the second case we get (18c), (19e), (19f), and (20c) with zi9 (β, λ = λ, i i i zi9 (β, λ, 0) = zi12 .

B

Formal derivation of the design derivatives

Based on the Lagrangian in Eq. (16), we can derive the derivative

∂L µi . ∂µi (wi , wp,i , hp,i ; µ)˜

∂Li ˜ ∂ϑj,i (wi , wp,i , hp,i , zi ; µ)ϑj,i

In the case of

the above-defined parameter vector, the derivatives for j = 3, . . . , 8 needs to be Np computed. However, note that µ ∈ U ⊂ R . Since the subset U might not be the complete space the classical concept of setting the gradient to zero is not valid and a variational concept has to be used. In general, we have

24

∂ y , z¯; µ ¯)(µ ∂µ L(¯

−µ ¯) ≥ 0 ∀µ ∈ U , where y¯, z¯, and µ ¯ denote the optimal value. This corresponds to

∂Li (wi , wp,i , hp,i , zi ; µ)(ϑ3,i − ϑ¯3,i ) ∂ϑ3,i    Z αend Z β1 
∂Ji ¯ ·) 3 1 − ε zi − zp,i (·, λ = λ, ¯ ·) wi − wp,i (·, λ = λ, = − dβ dα (ϑ3,i − ϑ¯3,i ), ¯ ∂ϑ3,i ε λ β0 0 ! Z λ¯ Z αend Z β1 ∂Li ∂J 3 ∂w ∂z i p,i p,i 2 ¯ (wi , wp,i , hp,i , zi ; µ)(ϑ4,i − ϑ¯4,i ) = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα (ϑ4,i − ϑ4,i ), ∂ϑ4,i ∂ϑ4,i β0 λ 0 ! Z αend Z β1 Z λ¯ 3 ∂Li ∂Ji ∂hp,i ∂gp,i 2 ¯ ¯ (wi , wp,i , hp,i , zi ; µ)(ϑ5,i − ϑ5,i ) = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα (ϑ5,i − ϑ5,i ), ∂ϑ5,i ∂ϑ5,i β0 λ 0 ! Z λ¯ ¯ Z αend Z β1 3 ∂Li ∂J ∂ ψ i i 2 ¯ (wi , wp,i , hp,i , zi ; µ)(ϑ6,i − ϑ¯6,i ) = − ¯ 3 0 ∂ϑ6,i (zp,i − gp,i )λ dλ dβ dα (ϑ6,i − ϑ6,i ), ∂ϑ6,i ∂ϑ6,i β0 λ 0 ! Z λ¯ ¯ Z αend Z β1 3 ∂ ψi ∂Li ∂Ji 2 ¯ ¯ (wi , wp,i , hp,i , zi ; µ)(ϑ7,i − ϑ7,i ) = − ¯ 3 0 ∂ϑ7,i (zp,i − gp,i )λ dλ dβ dα (ϑ7,i − ϑ7,i ), ∂ϑ7,i ∂ϑ7,i β0 λ 0 ! Z λ¯ ¯ Z αend Z β1 ∂ ψ 3 ∂Li ∂J i i 2 ¯ (wi , wp,i , hp,i , zi ; µ)(ϑ8,i − ϑ¯8,i ) = − ¯ 3 0 ∂ϑ8,i (zp,i − gp,i )λ dλ dβ dα (ϑ8,i − ϑ8,i ). ∂ϑ8,i ∂ϑ8,i β0 λ 0 If U = RNp , then µ − µ ¯ may attain every value µ ˜ ∈ RNp and thus the following equations can be obtained ∂Li (wi , wp,i , hp,i , zi ; µ)ϑ˜3,i ∂ϑ3,i    Z αend Z β1 
3 1−ε ∂Ji ¯ ¯ zi − zp,i (·, λ = λ, ·) dβ dα ϑ˜3,i , = − wi − wp,i (·, λ = λ, ·) ¯ ∂ϑ3,i ε λ 0 β0 ! Z αend Z β1 Z λ¯ ∂Li ∂Ji ∂wp,i ∂zp,i 2 3 ˜ ˜ (wi , wp,i , hp,i , zi ; µ)ϑ4,i = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα ϑ4,i , ∂ϑ4,i ∂ϑ4,i 0 β0 λ ! Z αend Z β1 Z λ¯ ∂Li ∂g ∂J ∂h 3 p,i i p,i 2 ˜ (wi , wp,i , hp,i , zi ; µ)ϑ˜5,i = − ¯ 3 0 ∂λ ∂λ λ dλ dβ dα ϑ5,i , ∂ϑ5,i ∂ϑ5,i 0 β0 λ ! Z αend Z β1 Z λ¯ ¯ ∂Li ∂Ji 3 ∂ ψi 2 ˜ ˜ (wi , wp,i , hp,i , zi ; µ)ϑ6,i = − ¯ 3 0 ∂ϑ6,i (zp,i − gp,i )λ dλ dβ dα ϑ6,i , ∂ϑ6,i ∂ϑ6,i 0 β0 λ ! Z λ¯ ¯ Z αend Z β1 ∂ ψ ∂Li ∂J 3 i i 2 ˜ (wi , wp,i , hp,i , zi ; µ)ϑ˜7,i = − ¯ 3 0 ∂ϑ7,i (zp,i − gp,i )λ dλ dβ dα ϑ7,i , ∂ϑ7,i ∂ϑ7,i 0 β0 λ ! Z αend Z β1 Z λ¯ ¯ ∂Ji ∂ ψi ∂Li 3 2 ˜ ˜ (wi , wp,i , hp,i , zi ; µ)ϑ8,i = − ¯ 3 0 ∂ϑ8,i (zp,i − gp,i )λ dλ dβ dα ϑ8,i . ∂ϑ8,i ∂ϑ8,i 0 β0 λ Since then all variations ϑi,j ∈ R, the tilde variables can be skipped. Considering the reduced cost functional ∂ ˆ L(y, z; µ) , if y = y(µ) and z = z(µ). J(µ) = J(y(µ); µ) introduced in (24), it can be shown that Jˆ0 (µ) = ∂µ

References Bankston, T. E. (2008). Theory and applications of refractive index-based optical microscopy to measure protein mass transfer in spherical adsorbent particles. J. Chromatogr., A 1188:242–254. Behrens, M., Bock, H. G., Engell, S., Khobkhun, P., and Potschka, A. (2014). Real-time PDE constrained optimal control of a periodic multicomponent separation process. In Leugering, G., Benner, P., Engell, S., Griewank, A., Harbrecht, H., Hinze, M., Rannacher, R., and Ulbrich, S., editors, Trends in PDE Constrained Optimization, volume 165, pages 521–537. Birkh¨auser, Basel, Switzerland. Bernauer, M. K. and Herzog, R. (2011). Optimal control of the classical two-phase Stefan problem in level set formulation. SIAM J. Sci. Comput., 33(1):342–363. 25

Bowes, B. D. and Lenhoff, A. M. (2011). Protein adsorption and transport in dextran-modified ion-exchange media. II. Intraparticle uptake and column breakthrough. J. Chromatogr., A, 1218(29):4698–4708. Bradley, R. S. and Thornley, M. T. (2006). A review of attenuation correction techniques for tissue fluorescence. J. R. Soc. Interface, 3:1–13. Brooks, C. A. and Cramer, S. M. (1992). Steric mass-action ion exchange: displacement profiles and induced salt gradients. AIChE J., 38(12):1969–1978. Carta, G., Ubiera, A. R., and Pabst, T. M. (2005). Protein mass transfer kinetics in ion exchange media: Measurements and interpretations. Chem. Eng. Technol.: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology, 28(11):1252–1264. Cheng, X., Lin, G., Zhang, Y., Gong, R., and Gulliksson, M. (2018). A modified coupled complex boundary method for an inverse chromatography problem. J. Inverse Ill-posed Probl., 26(1):33–49. DiLeo, M. V., Kellum, J. A., and Federspiel, W. J. (2009). A simple mathematical model of cytokine capture using a hemoadsorption device. Ann. Biomed. Eng., 37(1):222–229. Du, C.-J. and Sun, D.-W. (2009). Retrospective shading correction of confocal laser scanning microscopy beef images for three-dimensional visualization. Food Bioproc. Techn., 2:167–176. Fahrbach, F., Ohser, J., Haas, P., Menstell, P., Schw¨ammle, A., Osterroth, S., and Dobrovolskij, D. (2019). Attenuation correction for confocal laser scanning microscopy and its application in chromatography. submitted to J. Microsc. Forss´en, P., Arnell, R., and Fornstedt, T. (2006). An improved algorithm for solving inverse problems in liquid chromatography. Comput. Chem. Eng., 30(9):1381–1391. Geiger, C. and Kanzow, C. (1999a). Ein allgemeines Abstiegsverfahren. In Numerische Verfahren zur L¨ osung unrestringierter Optimierungsaufgaben, pages 25–33. Springer, Berlin, Heidelberg. Geiger, C. and Kanzow, C. (1999b). Schrittweitenstrategien. In Numerische Verfahren zur L¨ osung unrestringierter Optimierungsaufgaben, pages 35–44. Springer, Berlin, Heidelberg. Golshan-Shirazi, S. and Guiochon, G. (1992). Comparison of the various kinetic models of non-linear chromatography. J. Chromatogr., A, 603(1-2):1–11. Gu, T., Iyer, G., and Cheng, K.-S. C. (2013). Parameter estimation and rate model simulation of partial breakthrough of bovine serum albumin on a column packed with large Q Sepharose anion-exchange particles. Sep. Purif. Technol., 116:319–326. Gu´elat, B., Khalaf, R., Lattuada, M., Costioli, M., and Morbidelli, M. (2016). Protein adsorption on ion exchange resins and monoclonal antibody charge variant modulation. J. Chromatogr., A, 1447:82–91. Guiochon, G., Felinger, A., and Shirazi, D. G. (2006). Fundamentals of preparative and nonlinear chromatography. Elsevier, Amsterdam, the Netherlands. Hahn, T., Sommer, A., Osberghaus, A., Heuveline, V., and Hubbuch, J. (2014). Adjoint-based estimation and optimization for column liquid chromatography models. Comput. Chem. Eng., 64:41–54. Hinze, M., Pinnau, R., Ulbrich, M., and Ulbrich, S. (2008). Optimization with PDE constraints, volume 23. Springer Science & Business Media, Dordrecht, the Netherlands. Hubbuch, J., Linden, T., Knieps, E., Th¨ommes, J., and Kula, M.-R. (2002). Dynamics of protein uptake within the adsorbent particle during packed bed chromatography. Biotechnol. Bioeng., 80(4):359–368. James, F., Sep´ ulveda, M., Charton, F., Quinones, I., and Guiochon, G. (1999). Determination of binary competitive equilibrium isotherms from the individual chromatographic band profiles. Chem. Eng. Sci., 54(11):1677–1696. Kaczmarski, K. (2007). Estimation of adsorption isotherm parameters with inverse method—possible problems. J. Chromatogr., A, 1176(1-2):57–68.

26

Kavoosi, M., Sanaie, N., Dismer, F., Hubbuch, J., Kilburn, D. G., and Haynes, C. A. (2007). A novel two-zone protein uptake model for affinity chromatography and its application to the description of elution band profiles of proteins fused to a family 9 cellulose binding module affinity tag. J. Chromatogr., A, 1160(1-2):137–149. Kervrann, C., Legland, D., and Pardini, L. (2004). Robust incremental compensation of the light attenuation with depth in 3D fluorescence microscopy. J. Microsc., 214:297–314. Klatt, K.-U. (1999). Modellierung und effektive numerische Simulation von chromatographischen Trennprozessen im SMB-Betrieb. Chem. Ing. Tech., 71(6):555–566. Lee, S. C. and Bajcsy, P. (2006). Intensity correction of fluorescent confocal laser scanning microscope images by mean-weight filtering. J. Microsc., 221(2):122–136. LeVeque, R. J. (2002). Finite volume methods for hyperbolic problems, volume 31. Cambridge University Press, Cambridge. Leweke, S. and von Lieres, E. (2018). Chromatography analysis and design toolkit (CADET). Comput. Chem. Eng., 113:274–294. Lewus, R. K., Altan, F. H., and Carta, G. (1998). Protein adsorption and desorption on gel-filled rigid particles for ion exchange. Ind. Eng. Chem. Res., 37(3):1079–1087. Lienqueo, M. E., Mahn, A., Salgado, J., and Shene, C. (2012). Mathematical modeling of protein chromatograms. Chem. Eng. Technol., 35(1):46–57. Litwiniszyn, J. (1967). On some mathematical models of the suspension flow in porous medium. Chem. Eng. Sci., 22(10):1315–1324. Ma, Z., Whitley, R., and Wang, N.-H. (1996). Pore and surface diffusion in multicomponent adsorption and liquid chromatography systems. AIChE J., 42(5):1244–1262. Mehay, A. and Gu, T. (2014). A general rate model of ion-exchange chromatography for investigating ionexchange behavior and scale-up. J. Microb. Biochem. Technol., 6:216–222. Mesa, M. I., Llanes-Santiago, O., Fern´andez, F. H., Rodriguez, D. C., Neto, A. J. D. S., and Cˆamara, L. D. T. (2011). An approach to parameters estimation of a chromatography model using a clustering genetic algorithm based inverse model. Soft Comput., 15(5):963–973. ˇ Mich´ alek, J., Capek, M., Mao, X. W., and Kub´ınova, L. (2010). Application of morphology filters to compensation of lateral illumination inhomogeneities in confocal microscopy images. Anal. Biomed. Sign. Images, 20:49–54. Michel, M., Epping, A., and Jupke, A. (2005). Modeling and determination of model parameters. In SchmidtTraub, H., Schulte, M., and Seidel-Morgenstern, A., editors, Preparative Chromatography: Of Fine Chemicals and Pharmaceutical Agents, pages 215–312. John Wiley & Sons, Weinheim, Germany. Nocedal, J. and Wright, S. J. (2006a). Calculating derivatives. In Numerical Optimization, pages 193–219. Springer, New York, NY. Nocedal, J. and Wright, S. J. (2006b). Derivative-free optimization. In Numerical Optimization, pages 220–244. Springer, New York, NY. Nocedal, J. and Wright, S. J. (2006c). Fundamentals of unconstrained optimization. In Numerical Optimization, pages 10–29. Springer, New York, NY. Nocedal, J. and Wright, S. J. (2006d). Quasi-Newton methods. In Numerical Optimization, pages 135–163. Springer, New York, NY. Persson, P., Gustavsson, P.-E., Zacchi, G., and Nilsson, B. (2006). Aspects of estimating parameter dependencies in a detailed chromatography model based on frontal experiments. Process Biochem., 41(8):1812–1821. Persson, P. and Nilsson, B. (2001). Parameter estimation of protein chromatographic processes based on breakthrough curves. IFAC Proceedings Volumes, 34(5):113–118.

27

Puettmann, A., Schnittert, S., Leweke, S., and von Lieres, E. (2016). Utilizing algorithmic differentiation to efficiently compute chromatograms and parameter sensitivities. Chem. Eng. Sci., 139:152–162. Puettmann, A., Schnittert, S., Naumann, U., and von Lieres, E. (2013). Fast and accurate parameter sensitivities for the general rate model of column liquid chromatography. Comput. Chem. Eng., 56:46–57. Qamar, S., Akram, N., and Seidel-Morgenstern, A. (2016). Analysis of general rate model of linear chromatography considering finite rates of the adsorption and desorption steps. Chem. Eng. Res. Des., 111:13–23. Qamar, S., Uche, D. U., Khan, F. U., and Seidel-Morgenstern, A. (2017). Analysis of linear two-dimensional general rate model for chromatographic columns of cylindrical geometry. J. Chromatogr., A, 1496:92–104. Rammo, O., Menstell, P., Edelmann, B., Skudas, R., and Schulte, M. (2017). Fundamental toolbox for protein A resin characterization and link to application efficiency. Poster at 253rd American Chemical Society National Meeting & Exposition, April 2-6 2017, San Francisco. Roerdink, J. B. T. M. and Bakker, M. (1993). An FFT-based method for attenuation correction in fluorescence confocal microscopy. J. Microscopy, 169:3–14. Schnittert, S., Winz, R., and von Lieres, E. (2009). Development of a 3D model for packed bed liquid chromatography in micro-columns. In Comput. Modeling and Simul., 2009. EMS’09. Third UKSim European Symposium on, pages 193–197. IEEE. Schulte, M. and Epping, A. (2005). Fundamentals and general terminology. In Schmidt-Traub, H., Schulte, M., and Seidel-Morgenstern, A., editors, Preparative Chromatography: Of Fine Chemicals and Pharmaceutical Agents, pages 9–49. John Wiley & Sons, Weinheim, Germany. Shinya, I. (2006). Foundations of confocal scanned imaging in light microscopy. In Pawley, J., editor, Handbook of Biological Confocal Microscopy (3. Ed.), pages 1–19. Springer Science and Business Media LLC, New York. Stanislawski, B., Schmit, E., and Ohser, J. (2010). Imaging of fluorophores in spherical beads and estimation of radial density distributions. Image Anal. Stereol., 29:181–189. Stone, M. T., Cotoni, K. A., and Stoner, J. L. (2019). Cation exchange frontal chromatography for the removal of monoclonal antibody aggregates. J. Chromatogr., A, 1599:152 – 160. Strasters, K. C., van der Voort, H. T. M., Geusebroek, J. M., and Smeulders, A. W. M. (1994). Fast attenuation correction in fluorescence confocal imaging: a recursive approach. Bioimaging, 2:78–92. Suen, S.-Y. (1996). A comparison of isotherm and kinetic models for binary-solute adsorption to affinity membranes. J. Chem. Technol. Biotechnol.: International Research in Process, Environmental AND Clean Technology, 65(3):249–257. Susanto, A., Herrmann, T., and Hubbuch, J. (2006). Short-cut method for the correction of light attenuation influences in the experimental data obtained from confocal laser scanning microscopy. J. Chromatogr., A 1136:29–38. Tr¨oltzsch, F. (2010). Optimal control of partial differential equations: theory, methods, and applications, volume 112. American Mathematical Soc., Providence, Rhode Island. von Lieres, E. (2010). Chromatography models with Langmuir and steric mass action adsorption isotherms are of differential index one. In AIP Conf. Proc., volume 1281, pages 1004–1007. AIP. von Lieres, E. and Andersson, J. (2010). A fast and accurate solver for the general rate model of column liquid chromatography. Comput. Chem. Eng., 34(8):1180–1191. Xu, J., Zhu, L., Xu, G., Yu, W., and Ray, A. K. (2013). Determination of competitive adsorption isotherm of enantiomers on preparative chromatographic columns using inverse method. J. Chromatogr., A, 1273:49–56. Yang, K., Bai, S., and Sun, Y. (2008). Protein adsorption dynamics in cation-exchange chromatography quantitatively studied by confocal laser scanning microscopy. Chem. Eng. Sci., 63(16):4045–4054. Yang, K., Shi, Q.-H., and Sun, Y. (2006). Modeling and simulation of protein uptake in cation exchanger visualized by confocal laser scanning microscopy. J. Chromatogr., A, 1136(1):19–28. Zhang, Y., Lin, G., Gulliksson, M., Forss´en, P., Fornstedt, T., and Cheng, X. (2017). An adjoint method in inverse problems of chromatography. Inverse Probl. Sci. Eng., 25(8):1112–1137. 28

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

29