{'fraction_non_alphanumeric': 0.04504300588501584, 'fraction_numerical': 0.020122227252150293, 'mean_word_length': 4.394505494505495, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model-based process design of a ternary protein separation using multi-step gradient ion-exchange SMB chromatography

Model-based process design of ion-exchange simulated moving bed (IEX-SMB) chromatography for center-cut separation of proteins is studied. Use of nonlinear binding models that describe more accurate adsorption behaviours of macro-molecules could make it impossible to utilize triangle theory to obtain operating parameters. Moreover, triangle theory provides no rules to design salt profiles in IEX-SMB. In the modelling study here, proteins (i.e., ribonuclease, cytochrome and lysozyme) on the chromatographic columns packed with strong cation-exchanger SP Sepharose FF is used as an example system. The general rate model with steric mass-action kinetics was used; two closed-loop IEX-SMB network schemes were investigated (i.e., cascade and eight-zone schemes). Performance of the IEX-SMB schemes was examined with respect to multi-objective indicators (i.e., purity and yield) and productivity, and compared to a single column batch system with the same amount of resin utilized. A multi-objective sampling algorithm, Markov Chain Monte Carlo (MCMC), was used to generate samples for constructing the Pareto optimal fronts. MCMC serves on the sampling purpose, which is interested in sampling the Pareto optimal points as well as those near Pareto optimal. Pareto fronts of the three schemes provide the full information of trade-off between the conflicting indicators of purity and yield. The results indicate the cascade IEX-SMB scheme and the integrated eight-zone IEX-SMB scheme have the similar performance that both outperforms the single column batch system.


Introduction
Protein purification and separation are major concerns in the downstream processing of pharmaceutical industries (Carta & Jungbauer, 2010;Scopes, 2013). They require a series of processes aimed at isolating single or multiple types of proteins from complex intermediates of process productions. Choosing a purification method, depending on the purpose for which the protein is needed, is a critical issue. Chromatography is a prevailing separation technology. Simulated moving bed (smb) chromatography (Broughton & Gerhold, 1961) is an excellent alternative to single column batch chromatography, because of its continuous counter-current operation, its potential to enhance productivities and to reduce solvent consumption (Seidel-Morgenstern et al., 2008;Rajendran et al., 2009;Rodrigues et al., 2015). According to the position of columns relative to the ports (i.e., feed, raffinate, desorbent, extract), the process is divided into four zones (i.e., I, II, III, IV), each has a specific functionality in the separations (see Fig. 1).
However, the prominent features of smb processes are based on the prerequisite that optimal operating conditions can be determined, which can be very challenging in practice. Determination of operating conditions is predominantly based on classical and extended triangle theories (Storti et al., 1993), which are originally derived from the ideal chromatographic model with the linear isotherm. Powerful criteria to design processes with nonlinear Langmuir pattern (e.g., Langmuir, modified Langmuir and bi-Langmuir isotherms) have been developed in the past years (Charton & Nicoud, 1995;Mazzotti et al., 1997;Antos & Seidel-Morgenstern, 2001). Versatile results on triangle theory for the design of smb processes have been published (Nowak et al., 2012;Kim et al., 2016;Lim, 2004;Kazi et al., 2012;Bentley & Kawajiri, 2013;Bentley et al., 2014;Sreedhar et al., 2014;Toumi et al., 2007;Silva et al., 2015;Kiwala et al., 2016).
Substances separated by smb processes have recently been evolved from monosaccharides to proteins and other macro-molecules. The adsorption behaviors of macro-molecules, as observed both in experiments (Clark et al., 2007) and molecular dynamic simulations (Dismer & Hubbuch, 2010;Liang et al., 2012;Lang et al., 2015), are much more complex than that described by the linear isotherm, the Langmuir pattern or the steric mass-action (sma) model (Brooks & Cramer, 1992). Macro-molecules can undergo conformation changes in the mobile phase, orientation changes on the functional surfaces of the stationary phase and aggregation in the whole process. Clark et al. (2007) used an acoustically actuated resonant membrane sensor to monitor the changes of surface energy and found that the changes persist long after mass loading of the protein has reached the steady state. Therefore, more sophisticated adsorption models, such as the multi-state sma (Diedrich et al., 2017) based on the spreading model (Ghosh et al., 2013(Ghosh et al., , 2014, have been proposed to describe dynamic protein-ligand interactions. Although these models describe more accurately the adsorption behaviors by taking orientation changes into consideration, the strong nonlinearity of the binding process causes tremendous difficulties in deriving analytical formulas to calculate flowrates, as required in the triangle theory. Ion-exchange (iex) mechanism is commonly used to separate charged components in chromatography. In iex chromatography (both single column and smb processes), an auxiliary component acting as a modifier is used to change electrostatic interaction forces between functional groups and macro-molecules; sodium chloride is most frequently used. The higher the salt concentration is, the lower the electrostatic force is, such that the binding affinity is decreased and bound components are successively eluted (Lang et al., 2015). In the single column iex processes, this principle can be implemented in an isocratic mode or in a linear gradient mode. The linear gradient mode has frequently been applied to achieve better separation capability (Osberghaus et al., 2012b). While in the iex-smb processes, the isocratic mode (i.e., the salt strength is identical in all zones) is originally applied. The gradient mode can be applied, by using a salt strength in the desorbent stream that is stronger than that in the feed (cf. Fig. 1), thus, high salt concentration in zones I and II (high salt region), low salt concentration in zones III and IV (low salt region). As shown in some studies, applying the gradient mode can further improve the separation performance (Antos & Seidel-Morgenstern, 2001;Houwing et al., 2002Houwing et al., , 2003Li et al., 2007). The salt profiles of the gradient mode are referred to as two-step salt gradient in this work; an extension to multi-step gradients for ternary separations is straightforward. However, the construction of the multi-step salt gradient imposes difficulty on the design of iex-smb processes.
This difficulty can be detoured by using open-loop multicolumn chromatography (e.g., jo and mcsgp), where it is possible to apply linear salt gradients. The open-loop feature and possibility of linear gradient implementation provides great application potential for multicolumn continuous chromatography in downstream processing (Faria & Rodrigues, 2015), as dealing with the model complexity induced by linear gradients is much easier than that by multi-step ones; moreover, the open-loop feature renders more degrees of freedom in the process designs.
Conventional four-zone smb processes are tailored for binary separations. Ternary separations are currently requested in many applications; the selection of network configurations is application-specific. Ternary separations, also referred to as center-cut separations, can be performed by the jo process which operates in semi-continuous mode exploiting batch chromatography followed by the smb operation without feed delivery (Borges da Silva & Rodrigues, 2008;Graça et al., 2015). They can also be achieved by cascading two four-zone smb units (Wooley et al., 1998;Nicolaos et al., 2001;Nowak et al., 2012). Cascade schemes require designing and operating two four-zone smb processes, which can be laborious and costly. Hence, integrated schemes have been developed using a single multi-port valve and fewer pumps than cascade schemes (Seidel-Morgenstern et al., 2008;da Silva & Seidel-Morgenstern, 2016). Moreover, process design is simplified due to simultaneous switching of all columns. As the complexity of network configurations, in addition to the nonlinearity of the binding models and multi-step gradient, a comprehensive column model, fast numerical solver and efficient optimizer are required for the process design.
In this study, we shall present a model-based process design of iex-smb units for separating a protein mixture of ribonuclease, cytochrome and lysozyme, using cation-exchange columns packed with SP Sepharose FF beads. Two network configurations for the ternary separation will be used, that is, a cascade scheme and an eight-zone scheme. For comparison, a conventional single column batch system, with the same column geometry and resin utilization, will be studied. In this study, smb processes are modeled by weakly coupling individual models of the involved columns (He et al., 2018).
The operating conditions of the network configurations will be optimized by a stochastic algorithm, Markov Chain Monte Carlo (mcmc), with respect to conflicting objectives of purity and yield. Pareto fronts will be computed for illustrating the best compromises between the two conflicting performance indicators. Unlike multi-objective optimization algorithms (e.g., non-dominated sorted genetic algorithm, strength Pareto evolutionary algorithms) that try to eliminate all the non-dominated points during optimization, mcmc serves on the sampling purpose, which is interested in sampling the Pareto optimal points as well as those near Pareto optimal. For sampling purpose, mcmc not only accepts proposals with better objective values, but also accepts moves heading to non-dominated points with certain probability. Therefore, mcmc could render greater information for process design, such as, uncertainties of parameters that is related to the robustness of process design can be examined. As shown in the Pareto fronts, the performance indicators of the single column are dominated by that of both the cascade and the eight-zone schemes. The performances of the two iex-smb schemes are quite similar, though the cascade scheme has a slight advantage over the eight-zone scheme.

Theory
This section introduces the transport model, the binding model, the load-wash-elution mode of single column chromatography, the smb network connectivity, the performance indicators and the optimization algorithm.

Transport model
The transport of proteins in the column is described by means of the general rate model (grm), which accounts for various levels of mass transfer resistance (Guiochon et al., 2006).
grm considers convection and axial dispersion in the bulk liquid, as well as film mass transfer and pore diffusion in the porous beads: where c j i , c j p,i and q j i denote the interstitial, stagnant and stationary phase concentrations of component i ∈ {1, . . . , M } in column j ∈ {1, . . . , N }. Particularly in the iex chromatography, when the salt component is introduced, i ∈ {0, 1, . . . , M } and i = 0 denotes the salt ions, whose concentration will be required by the steric mass-action (sma) model (cf. Eq. (4)). Furthermore, z ∈ [0, L] denotes the axial position where L is the column length, r ∈ [0, r p ] is the radial position, r p is the particle radius, t time, ε c column porosity, ε p particle porosity, u j int interstitial velocity, D j ax is the axial dispersion coefficient, D j p,i the effective pore diffusion coefficient, and k j f,i the film mass transfer coefficient. At the column inlet and outlet, Danckwerts boundary conditions (Barber et al., 1998) are applied: where c j in,i is the inlet concentration of component i in column j; the calculation of inlet concentration in smb chromatography is deferred to Eq. (7). The boundary conditions at the particle surface and center are given by:

Binding model
The sma model has widely been used for predicting and describing nonlinear ion-exchange adsorption of proteins. The binding of components i 1 is expressed as , the adsorptive fluxes of components i 1 are influenced by the salt concentrations, q j 0 , c j p,0 (i = 0): where k a,i and k d,i denote adsorption and desorption coefficients, and ν i is the characteristic charge of the adsorbing molecules. The concentrationq j 0 of available counter ions has to be algebraically determined from an electro-neutrality condition, as Eq. (4) is only valid for components i 1, Here Λ is the ionic capacity, and σ i the shielding factor. Then, the salt concentration in the stationary phase is given by:

Single column batch system
There are various possibilities to operate a single column in batch manner. The specific one described here is used for comparison with the smb processes; we refer you the protocol to Osberghaus et al. (2012b). Apart from regeneration steps, an operating cycle of this single column process consist of three phases: load, wash and elution, see Fig. 2. In phase I, the column is equilibrated with the running buffer, then the protein solution is injected to the column for time t load . This is followed by a wash step (II) for time t wash . Afterwards, different gradients including a gradient of zero length, i.e., isocratic elution, are set to elute the bound protein for time t elute . The elution phase (III) can be implemented by multi-step solvent gradients. As depicted in Fig. 2 as an example, there are two parts (i.e., IIIa, IIIb) with linear gradients but different slopes (m 1 , m 2 ), and eventually an isocratic part (IIIc). Multi-step gradients are used, because in center-cut separations peak overlaps at the front and at the back of the target peak need to be minimized. The regeneration of the column (IV) completes one operating cycle of the single column batch system.

SMB network connectivity
In smb chromatography two adjacent columns, j and j+1, are connected via a node j, see Fig. 3. A circular smb loop is closed when two column indices point to the same physical column (i.e., by identifying column j = N +1 with column j = 1), see Fig. 1. In this work, node j is located at the downstream side of column j and the upstream side of column j + 1. At each node only one or none of the feed (F), desorbent (D), raffinate (R), or extract (E) streams exists at a time. The inlet concentration of component i ∈ {0, 1, . . . , M } in column j +1 is calculated from mass balance of the node j: Figure 3: Schematic of network connection between two adjacent columns.
where c j out,i = c j i (t, z = L) denotes the outlet concentration of component i in column j, Q j = ε c u j int πd 2 c /4 the volumetric flowrates and d c the column diameter. Meanwhile, δ j determines the current role of node j (i.e., F, D, R, E or none): where c F in and c D in are the component concentrations at feed and desorbent ports, and Q F , Q D , Q R , Q E the volumetric flowrates at the feed, desorbent, raffinate, and extract ports, respectively. δ j = 0 indicates nodes that are currently not connected to a port (i.e., in the interior of a zone); this occurs when more columns than zones are present, such as eight columns in a four-zone scheme. As there might be multiple feed, desorbent, raffinate, and extract ports in smb processes (e.g., the cascade scheme and the integrated eight-zone scheme), an straightforward extension of the indexing scheme is required. Column shifting is implemented by periodically permuting δ each switching time t s .

Performance indicators
Purity, yield and productivity are commonly used for evaluating the performance of chromatographic processes. All indicators can be defined within one collection time t c (Bochenek et al., 2013). In the case of smb systems, t c is the switching time t s , while in single column systems it is the length of pooling time interval t p .
In this study, performance indicators are all defined in terms of component, i ∈ {1, . . . , M }, withdrawn at a smb node, j ∈ {E, R, . . . }, or collected at the outlet of single column systems, j ∈ {B}, within the pooling interval. Note that we do not calculate the performance indicators of the salt component, i = 0.

The definitions of performance indicators are all based on concentration integrals
In smb systems, τ is the starting time of one switching interval and all performance indicators are calculated when the system is upon cyclic steady state (css). In single column systems, τ is the starting point of the pooling time interval. The purity, Pu j i , is the concentration integral of component i relative to the integral sum of all components, Eq. (10). The yield, Y j i , is the ratio of the withdrawn mass and the feed mass, Eq. (11).
In smb systems, when the solution stream is continuously fed via the feed node and the extracts and raffinates are continuously withdrawn, t load = t s = t c and Q F differs from Q j . While in single column systems, the flowrates are the same at the inlet node and outlet node, Q j = Q F ; and t load = t p = t c . The productivity, Pr j i , is the withdrawn mass of the component i per collection time relative to the total volume of the utilized packed bed in all columns, In single column systems, N = 1. According to the definitions of Pr j i and Y j i , both indicators increase with increasing amounts of the product collected, i.e., A j out,i . However, Pr j i can also be improved by reducing the collection time t c .

Optimization
Generally, operating conditions of both smb and single column processes are systematically optimized by numerical algorithms such that the above performance indicators are all maximized. In studies of center-cut separations, the middle component i is of interest. To specifically account for a vector of the conflicting performance indicators of Pu j i , Y j i (referred to as PI hereafter) in this study, multi-objective optimization is applied.
A set of objectives can be combined into a single objective by optimizing a weighted sum of the objectives (Marler & Arora, 2010), or keeping just one of the objectives and with the rest of the objectives constrained (the ε-constraint method (Mavrotas, 2009)). The latter method is used in this work, that is, maximizing the yield of the target component of the processes with the purity constrained to be larger than a threshold ε j i : where θ is a vector of optimized parameters with boundary limitations of [θ min , θ max ]. The inequality c j i (θ) is lumped into the objective function using penalty terms in this study, such that it can be solved as a series of unconstrained minimization problems with increasing penalty factors, σ k , In Eq. (14), the penalty function is chosen as g(θ) = j∈{R,E} min{0, c j (θ)} 2 .
Different types of methods can be applied to solve the minimization of H(θ; σ k ), such as deterministic methods and heuristic methods. Additionally, stochastic methods can be chosen. In order to use a stochastic method, the minimization of H(θ; σ k ) is further formulated to maximum likelihood estimation: where the likelihood is defined as an exponential function of H(θ; σ k ):

Numerical solution
All numerical simulations were computed on an Intel(R) Xeon(R) system with 16 CPU cores (64 threads) running at 2.10 GHz. The mathematical models described above for each column of smb processes are weakly coupled together and then iteratively solved. The open-source code has been published, https://github.com/modsim/ CADET-SMB. cadet-smb repeatedly invokes cadet kernel to solve each individual column model, with default parameter settings. cadet is also an open-source software, https://github.com/modsim/CADET. The axial column dimension is discretized into N z = 40 cells, while the radical bead is discretized into N r = 10 cells. The resulting system of ordinary differential equations is solved using an absolute tolerance of 10 −10 , relative tolerance of 10 −6 , an initial step size of 10 −14 and a maximal step size of 5 · 10 6 .
A stochastic multi-objective sampling algorithm, mcmc, is applied in this study to optimize the operating conditions. A software computing the Metropolis-Hastings algorithm with delayed rejection, adjusted Metropolis and Gibbs sampling has been published as open-source software, https://github.com/modsim/ CADET-MCMC. For sampling purpose, mcmc not only accepts proposals with better objective values, but also accepts moves heading to non-dominated points with certain probability. Samples are collected until the Geweke convergence criterion (Geweke et al., 1991) is smaller than 10 −4 or the desired sample size is reached. Non-dominated stable sort method of Pareto fronts were is applied to generate the frontiers (Duh et al., 2012).

Case
A protein mixture of ribonuclease, cytochrome and lysozyme is of interest in this study, i ∈ {RNase, cyt, lyz}. Sodium chloride, i = 0, acts as the modifier of binding affinities of the proteins. It is a prototype example used in the academic field for modeling purposes (Osberghaus et al., 2012a,b). The center component, cyt, is targeted and the other two components are regarded as impurities. We have modeled separation of these proteins on chromatographic columns packed with strong cation-exchanger SP Sepharose FF (see the column geometry in Tab. 1). Simplified models (e.g., tdm and edm) that assume instance equilibrium of the binding process can not be applied here, because mass transfer dynamics of macro-molecules in ca. 100 µm beads can be rate limiting (Lodi et al., 2017). Hence, the grm model is used to describe the mass transfer in the porous beads. In addition, the sma model inherently considers the impact of salt on binding affinity, the hindrance effect of macro-molecules. The mathematical modeling of the prototype example has been presented in Montesinos et al. (2005); Püttmann et al. (2013); Osberghaus et al. (2012a,b); the transport and binding parameters from experiments are shown in Tab. 1. The transport parameters can also be calculated from correlation formulas. D ax = 5.75 · 10 −8 m 2 s −1 , D p = 6.07 · 10 −11 m 2 s −1 and k f = 6.90 · 10 −6 m s −1 are taken for all proteins because they have almost the same molecular weight. Relatively large values are used for the rate parameters of the binding model, k a,i and k d,i , as the adsorption-desorption kinetics is typically much faster than kinetics of mass transport in pores. The geometry, transport and binding parameters are directly used in this study, while the operating conditions will be optimized.
In pharmaceutical manufacturing typically at least 99% purity of cyt is needed. Therefore, taking expensive computation cost in iex-smb simulations to generate Pareto optimal points with lower purities is not required here. In this study, the ε j cyt in Eq. (13) is set to 99%; thus, only the Pareto optimal fragment with cyt purity higher than 99% is concerned in iex-smb. This implementation reduces the heavy computational burden of multi-objective sampling, as much fewer samples are requested to construct the frontiers. Pareto optimal fragment can be further extended by repeatedly solving Eq. (15) with varying values of ε j cyt . This can be one of the advantages of the ε-constraint method. In single column systems, the whole Pareto optimal front is concerned as the computational cost is not expensive.

Single column
In the single column system, the equilibration, load-wash-elution and regeneration phases are performed. The interstitial velocity is 5.75 · 10 −4 m s −1 , resulting in a retention time for a non-retained component of 24.35 s. A sample with 1 mol m −3 concentration of RNase, cyt and lyz is injected in the load phase. In both the load and wash steps, a salt concentration of 50 mol m −3 is applied. The operating time intervals for the load and wash phases are 10 s and 40 s. The shapes of the multi-step elution phase, that can be characterized by the bilinear gradients (i.e., m 1 and m 2 ), the operating time intervals, and the initial salt concentration (i.e., c init,0 ), are optimized. Fig. 2 shows the elapsed times, from which the operating time intervals are calculated (i.e., ∆t 1 = t 1 elute − 50 and ∆t 2 = t 2 elute − t 1 elute ). Thus, the optimized operating parameters are θ = {∆t 1 , ∆t 2 , m 1 , m 2 , c init,0 }.

SMB
Depending on initial states, smb processes often undergoes a ramp-up phase and eventually enter into a css. In this study, the system is assumed to be upon css when a difference between two iterations (k, k+1) falls below a predefined tolerance error, e t , see Eq. (17).
In this study, the columns of iex-smb are initially empty except for the concentration of bound salt ions q j 0 that is set to the ionic capacity Λ, in order to satisfy the electro-neutrality condition. Additionally, the initial salt concentrations in column j, c j 0 (t = 0, z), c j p,0 (t = 0, z), q j 0 (t = 0, z), are set equal to the salt concentrations in the respective upstream inlet nodes. Taking the four-zone scheme as an example (cf. Fig. 1), the salt concentrations of the mobile and stationary phases (c j 0 and c j p,0 ) in zone III and IV are set to the value of c F 0 ; while c j 0 and c j p,0 in zone I and II are set to the value of c D 0 : In contrast to the single column system, the chosen initial state of smb processes is irrelevant to the performance indicators calculated upon css.
A cascade of two four-zone units (named as U 1 and U 2 in the following) and an integrated eight-zone scheme are applied. The two sub-units of the cascade schemes can be connected via either the raffinate or the extract ports; while in the eight-zone schemes, either the raffinate or the extract streams of the first sub-unit can be internally bypassed. In the cascade scheme, the flow between the raffinate stream of U 1 and feed stream of U 2 is defined as a bypass stream; while in the eight-zone scheme, the connection between the raffinate-I of the first sub-unit and the feed-II of the second sub-unit is defined as a bypass stream. This selection depends on the equilibrium constants of the model proteins; as k eq,i = [7.70, 1.59, 35.5] · 10 −3 , i ∈ {RNase, cyt, lyz}, thus δ (|k eq, RNase − k eq, cyt |) < δ (|k eq, lyz − k eq, cyt |) indicates we should choose the configurations (as shown in Fig. 4) that firstly separate the strongest adsorbed component, lyz, at the extract port of the first sub-unit in both schemes, and the two sub-units are connected via the raffinate of the first sub-unit. Thereafter, the binary separation of the target component, cyt, and the weakest adsorbed component, RNase, are achieved in the second sub-units. In both schemes, three single columns in each zone and 24 columns at all are designed.
In each sub-unit of the cascade scheme, the optimized operating parameters are the switching time (same for both units as they are synchronously switched), recycle flowrate, the inlet and outlet flowrates, and the salt concentrations at the inlet ports, While in the eight-zone scheme, the switch time, the recycle flowrate, the inlet and outlet flowrates, and salt concentrations at inlet ports are optimized, that is, The flowrates Q R and Q R2 are calculated from flowrate balances.

Single column
The gradient shapes of the single column process are first optimized in order to have a reference for the iex-smb processes. Search domain of the parameters, θ, is listed in Tab. 2. The intervals are based on literature values with additional safety margins. The maximal sampling size here is 1.2 · 10 4 and the burn-in length is set to 50% of the samples. Fig. 5 shows the Pareto optimal front between purity and yield. At a rather low purity requirement of 85%, a yield of 0.9 can be achieved. However, at purity of 98%, the yield drops dramatically to 0.1. The Pareto front provides full trade-off information of the single column process; the corresponding operating conditions, that render the Pareto optima on demand in application, can be chosen on purpose.
Three characteristic points (i.e., a, b, c) on the Pareto front are exemplified and compared. In Fig. 6, the corresponding chromatograms and multi-step salt gradients are shown; the gray areas indicate the pooling time intervals, t p , of the target component. The calculated performance indicators of the three characteristic points are as follows: PI B cyt = Pu B cyt , Y B cyt is [90.0%, 0.85], and the productivity is Pr B cyt = 4.94 · 10 −3 mol m −3 s −1 of the point a (see Fig. 6a). While they are [95.59%, 0.71] and 7.98 · 10 −5 mol m −3 s −1 of the point b (see Fig. 6b), [97.02%, 0.49] and 1.36 · 10 −4 mol m −3 s −1 of the point c (see Fig. 6c). The respective operating parameters are listed in Tab. 3.
The pooling time length t p = 5808 s in Fig. 6b is much longer than that in Fig. 6a (i.e., 112 s) and Fig. 6c (i.e., 2352 s). In this study, pooling time intervals were calculated such that the concentrations of other two components (i.e., RNase and lyz) are lower than a predefined threshold, µ = 7.5 · 10 −5 . Based on the same  2.81 · 10 3 3.12 · 10 3 6.26 · 10 2 s ∆t2 3.53 · 10 3 1.56 · 10 3 3.68 · 10 3 s m1 1.28 · 10 −3 4.29 · 10 −3 6.47 · 10 −3 m2 1.11 1.32 · 10 −2 1.62 · 10 −2 cinit,0 77.2 71.7 51.0 mol m −3 operating parameters, if we set µ to a larger value, higher yield and lower purity will be obtained. Therefore, the selection is also a trade-off, which is illustrated in a Pareto front in Fig. 7, where µ was decreased from 10 −4 to 2.5 · 10 −5 with equivalent gap of 2.5 · 10 −5 using operating parameters listed in the b column (i.e., point b) of Tab. 3. Shorter pooling time interval causes the cyt peak to be more concentrated. This can result in high yield but larger overlap areas (thus low purity), see Fig. 6a. When pooling time interval is larger, the cyt peak is more spread and it renders relatively lower yield but higher purity, see Fig. 6b. As observed in Fig. 6, the peak of cyt always overlaps with the other two peaks, base line separation is not achieved. In Fig. 6b and 6c, even when lyz is hardly washed out, it begins to be eluted early at small salt concentrations. Additionally, as lyz is almost not eluted out in Fig. 6b and 6c, a strip step of lyz is needed before the regeneration step, in order to make the column reusable. These overlaps may indicate that cyt can not be collected with 100% purity with the single column batch system.

Cascade scheme
Numerical optimization of the iex-smb schemes require suitable initial values. To this end, empirical rules have been developed with big manual efforts to achieve the center-cut separation first. With the operating conditions listed in the empirical column of Tab. 4, lyz is collected at the R port of U 1 with PI R lyz U1 = [98.23%, 0.99] (see Fig. 8). In U 2 , cyt spreads towards the R port, resulting low yield at the E port (i.e., Y E cyt U2 = 0.33) and low purity at the R port (i.e., Pu R RNase U2 = 58.39%) (the chromatogram is not shown). Based on the initial values with safety margins, search domain as shown in Tab. 5 is used for numerical optimization.
smb processes often undergo a long ramp-up phase first and eventually enter into a css. Elaborately   Zone I II III IV choosing the initial state of columns can also accelerate the convergence (Bentley et al., 2014). In this case, a total of ca. k = 108 switches was required for each iex-smb simulation to fall below the tolerance error of e t = 10 −5 , see Eq. (17). Thus, gaining points for generating Pareto fronts in smb chromatographic processes is computational expensive. Li et al. (2014) have proposed computationally cheap surrogate models for efficient optimization of smb chromatography. In this study, only the operating parameters of U 2 are optimized; the operating conditions of U 1 are taken directly from the empirical design. Moreover, only the Pareto optimal fragment with purity higher than 99% is concerned, such that fewer points are requested. The maximal sampling size of mcmc is 300, and the burn-in length is 50. As seen from the Pareto fronts in Fig. 9, the cascade scheme outperforms the single column system in both purity and yield. Three characteristic points (i.e., a, b, c), ranging yields from high to low, on the Pareto front are compared. The corresponding operating conditions are listed in Tab. 4; and the resulting chromatograms are shown in Fig. 10. Combined axial concentration profiles in the columns are displayed at multiples of the switching time. The performance vector PI E cyt U2 = Pu E cyt , Y E cyt U2 of the point a is [99.27%, 0.96] (cf. Fig. 10a), [99.81%, 0.95] of the point b (cf. Fig. 10b) and [100%, 0.39] of the point c (cf. Fig. 10c). After numerical optimization, the performance indicators tremendously increase from that of the empirical design, [97.92%, 0.33]. The main difference between the chromatograms of a and b is that less RNase keeps retained to the port E in b, resulting in higher purity at the port E. The productivities of cyt, Pr The salt profiles, c 0 (t, z), along the two iex-smb units at time upon css, t = τ , are depicted in Fig. 11. In U 1 , the salt concentrations at inlet ports, c D 0 , c F 0 U1 , used for constructing the two solvent gradients are [290,420] Fig. 11b).
The salt gap of solvent between zones I,II and zones III,IV of U 1 facilitates to firstly separate the most strongly adsorbed lyz from the ternary mixture; RNase and cyt are still not rather separated in U 2 . Because of the salt profile in zones III and IV (i.e., ca. 320 mol m −3 ), RNase and cyt have approximately the same electrostatic interaction forces with the stationary phase, causing rather similar desorption rates off the stationary phase. In order to separate RNase and cyt in U 2 , the salt concentrations at the inlet ports of U 2 must be reduced, leading one component (i.e., cyt) to be more strongly bounded. This can be achieved by diluting the bypass stream with pure buffer: In the numerical optimization, Q F U2 and Q R U1 are optimized, are subsequently changed. As shown in Fig. 11b, the salt concentration in U 2 is decreased to ca. 200 mol m −3 in zones I,II and 160 mol m −3 in zones III,IV. With these salt concentrations, RNase and cyt can be separated in U 2 . With the dilution, the protein concentrations are decreased as a side effect. Therefore, the peaks in U 2 are not as concentrated as observed in U 1 ; they distribute over several zones, which makes the process design of the U 2 challenging.

Eight-zone scheme
As in the design of cascade scheme, the initial values were obtained from empirical designs with manual efforts. With the operating conditions listed in the empirical column of Tab. 6, cyt spreads from zone VI to zone VII, RNase from zone V to VI, such that Pu R2 RNase = 73.23% and Pu E2 cyt = 87.50% (the chromatogram is not shown). The search domain for the numerical optimization is based on the above initial values with safety margins, see Tab. 7.
In this case, a total of ca. k = 148 switches was required for each eight-zone numerical simulation to fall below the tolerance error of e t = 10 −5 . The convergence time from an initial state to the css, thus the total computational time in the eight-zone scheme, is longer than that in the cascade scheme. The maximal sampling size of mcmc is 300, and the burn-in length is 50. The Pareto optimal front with purity larger than 99% is illustrated in Fig. 12; the Pareto fronts of single column and cascade systems are superimposed for comparison. The blue shadow areas show the constraint region.
By using the multi-objective sampling algorithm, mcmc, not only the Pareto front solutions but also the points near Pareto front can be studied. In this study, three points (i.e., a, b, c; the point c is on the Pareto front), ranging purities from low to high, were chosen to illustrate the impact of the composition of R1 stream on the performance of the second sub-unit. The corresponding operating conditions are listed in Tab. 6; and the chromatograms upon css are shown in Fig. 13. Combined axial concentration profiles in the columns are displayed at multiples of the switching time. In Fig. 13a Fig. 13c. As seen from Fig. 13, the higher concentration of the cyt composition is in the R1 stream, the higher of the yield and productivity indicators are at the E2 port. To be specific, Y E2 cyt = 0.95 and Pr E2 cyt = 3.10 mol m −3 s −1 of the point a; Y E2 cyt = 0.90 and Pr E2 cyt = 2.93 mol m −3 s −1 of the point b, and Y E2 cyt = 0.97 and Pr E2 cyt = 3.25 mol m −3 s −1 of the point c. Although in the center-cut separation, the other two components (i.e., RNase and lyz) can be viewed as impurities, they can be withdrawn with high purity and yield.  (Nicolaos et al., 2001).
The salt profile, c 0 (t, z), along the eight-zone iex-smb unit upon css, t = τ , is shown in Fig. 14 443,272,240,217] mol m −3 of the point c. Separating lyz at the port E1 is facilitated with the salt gap between zones I,II and zones III,IV; RNase and cyt are both weekly bounded in zones III,IV and to be separated in the second sub-unit. Meanwhile, the bypass stream was diluted as described by Eq. 20. As seen in Fig. 14 cyt is strongly bounded in the second sub-unit. As seen from Fig. 14, the deviation of the salt profile along the columns of iex-smb is small, but the performance indicators significantly varied. This indicates process design of iex-smb processes might be sensitive to operating conditions. Eight-zone schemes have fewer degrees of freedom than cascade schemes. Hence, Pu E2 cyt decreases from 99.81% of the cascade scheme to 97.18%, based on the same yield, 0.95. Moreover, in the cascade scheme, almost pure cyt can be collected, which can not be achieved with the eight-zone scheme. However, the productivity of the eight-zone scheme can be higher than that of the cascade scheme. The productivity of the cascade scheme, 3.13 · 10 −4 mol m −3 s −1 (PI E2 cyt = [99.27%, 0.96]), is slightly lower than that of the eight-zone scheme, 3.25 · 10 −4 mol m −3 s −1 (PI E2 cyt = [99.27%, 0.97]). This suggests that the stationary phase is more efficiently utilized with integrated smb processes.
As no linear gradients can be applied in the iex-smb schemes, the eight-zone and cascade schemes have fewer degrees of freedom than multi-column continuous chromatography (e.g., jo and mcsgp). Thus, we could conjecture that with proper process design, the jo and mcsgp would render more efficient and accurate elution of the model proteins. However, this contribution is focused on the closed-loop iex-smb systems and the challenge of designing multi-step salt profiles by using modeling and high performance computing.
Three columns in each zone were chosen and used in both the cascade and eight-zone schemes. In the empirical designs, neither one column nor two columns in each zone could lead to good separation performance. Having more columns in each zone and correspondingly shorter switching time can approximate the true moving bed chromatography. However, introducing more columns into iex-smb system is computationally challenging and practically unrealistic. Hence, three columns in each zone are a good compromise.

Sensitivity of operating conditions
As stated in the previous subsection, the css of both the cascade and the eight-zone schemes are sensitive to the operating conditions. By utilizing mcmc sampling (cf. He & Zhao (2019)), uncertainties of operating conditions are characterized by parameter distributions. Posterior predictive checks were carried out to test the robustness of the two investigated processes. Random parameter combinations from respective parameter distributions are used as arguments of forward simulations and to generate chromatograms. The possible deviations of the investigated processes are illustrated with curve clusters, as shown in Fig. 15.
As seen in Fig. 15  are quite significant; the same applies to the zones II and III of the cascade scheme. However, in the outlet streams of the eight-zone configuration (i.e., E1, R2 and E2) and the cascade configuration (i.e., E and R), the deviations are tiny, which means perturbations in the operating conditions characterized by the parameter distributions do not essentially affect the withdrawn streams and consequentially the performance indicators. Fig. 15 also indicates the eight-zone scheme is more robust than the cascade scheme. Therefore, the studied smb configurations are robust provided the processes are performed within the perturbation tolerance.

Conclusions
We have presented model-based process designs of multi-step salt gradient single column batch chromatography and iex-smb chromatographic processes for separating ribonuclease, cytochrome and lysozyme on strong cation-exchanger SP Sepharose FF. The general rate model with steric mass-action binding model was used in the column modeling. Two network configurations of closed-loop iex-smb have been studied, i.e., one cascade scheme of two four-zone iex-smb sub-units and one integrated eight-zone scheme; the single column batch system has been studied for comparison. The multi-objective sampling algorithm, Markov Chain Monte Carlo (mcmc), has been used to generate samples for approximating the Pareto optimal fronts. This study shows that it is possible to design closed-loop iex-smb units for collecting the target, cytochrome, with high purity and yield.
In the numerical optimizations of iex-smb, initial values are required. Empirical designs have been developed in this study to obtain the initial values, and the search domains of numerical optimization are based on the initial values with safety margins. mcmc serves on the sampling purpose, which is interested in sampling Pareto optimal points as well as those near Pareto optimal. With ε-constraint method of the multi-objective optimization, Pareto optimal fragments can be concerned to reduce expensive computational cost in large-scale iex-smb processes.
The Pareto fronts show the information of the trade-off relationship between the conflicting indicators (i.e., purity and yield). The cascade and eight-zone schemes have approximately the same performance, and both outperform the single column system. In addition, the cascade scheme can produce almost pure target component with high yield, but lower robustness; the stationary phase is utilized more efficiently in the integrated eight-zone scheme in terms of productivity. However, single column system does not undergo a time-consuming ramp-up phase before entering into a cyclic steady state, like it is in smb systems. Starting a single column batch system prior to complex multi-column and smb systems is recommended. Dilution steps were applied in the bypass streams of the iex-smb schemes to reduce the salt concentrations, leading cytochrome to be more strongly bounded in the second sub-unit. Without the dilution steps, the salt concentrations are too high such that ribonuclease and cytochrome are both weekly bounded to the stationary phase of the second sub-unit. Though it is beneficial to implement dilution, it consumes more running buffer on one hand, and the cytochrome product needs to be more concentrated in the subsequent downstream processes on the other hand.
There are many options for modeling mass transfer in a column, e.g., the equilibrium-dispersive model, the transport-dispersive model and the general rate model that can be combined with a multitude of binding models. The software package used in this study is flexible in choosing between various modeling options for the columns and hold-up volumes. smb processes do not necessarily need equivalent column numbers in each zone. Different column numbers in each zone can be applied, leading to mixed-integer programming. The computational burden of smb processes in numerical optimizations can be reduced by optimizing initial states of the involved columns, or by using computationally cheap surrogate models. Multi-dimensional Pareto optimal fronts can also be further investigated, taking buffer consumption, productivity or processing time into account.