Dynamic Flowsheet Model Development and Sensitivity Analysis of a Continuous Pharmaceutical Tablet Manufacturing Process Using the Wet Granulation Route

Abstract

In view of growing interest and investment in continuous manufacturing, the development and utilization of mathematical model(s) of the manufacturing line is of prime importance. These models are essential for understanding the complex interplay between process-wide critical process parameters (CPPs) and critical quality attributes (CQAs) beyond the individual process operations. In this work, a flowsheet model that is an approximate representation of the ConsiGma${}^{\mathrm{TM}}$-25 line for continuous tablet manufacturing, including wet granulation, is developed. The manufacturing line involves various unit operations, i.e., feeders, blenders, a twin-screw wet granulator, a fluidized bed dryer, a mill, and a tablet press. The unit operations are simulated using various modeling approaches such as data-driven models, semi-empirical models, population balance models, and mechanistic models. Intermediate feeders, blenders, and transfer lines between the units are also simulated. The continuous process is simulated using the flowsheet model thus developed and case studies are provided to demonstrate its application for dynamic simulation. Finally, the flowsheet model is used to systematically identify critical process parameters (CPPs) that affect process responses of interest using global sensitivity analysis methods. Liquid feed rate to the granulator, and air temperature and drying time in the dryer are identified as CPPs affecting the tablet properties.

1. Introduction

Flowsheet models are approximate mathematical representations of the manufacturing line. The incentives for flowsheet model development for the pharmaceutical industry have been described in the paper of Escotet-Espinoza et al. [1]. Gernaey et al. [2] also wrote extensively on the value of Process Systems Engineering (PSE) for pharmaceutical process development, in which flowsheet models contribute to combine knowledge and models of different unit operations and different scales for a holistic understanding of the process. The incentives for flowsheet modeling boil down to the in-silico achievement of process design and optimization, control system design and optimization, and an accurate risk assessment tool that could be used for regulatory instances. The first step in attaining these benefits is by the development of a flowsheet model that captures the relevant mechanisms for assessing the desired product properties as a function of process settings and material properties. This foundation built in this work, comprised of several diverse unit operation models that capture critical mechanisms, as a function of the process inputs.

The flowsheet model simulation enables the assessment of unit operation outputs downstream in the process. This has the advantage that, instead of applying the unit operation models separately for model-based research, the flowsheet model allows for targeted optimization of unit operation performance as a part of the entire line. In terms of the Quality-by-Design (QbD) paradigm, a flowsheet model allows for investigating the influence of critical process parameters (CPPs) in one unit to the critical quality attributes (CQAs) of material downstream in the process line. Namely, through flowsheet model development, process phenomena are directly linked to the final product quality downstream. Moreover, analysis of the developed model allows assessment of criticality of the various critical process parameters (CPPs) of the process, and subsequently research effort can be targeted towards those most critical areas.

In this work, composing the integrated system requires extensive synchronization of the unit operation models themselves, as these need to run seamlessly in one simulation, regardless of the different time scales, variable magnitudes, or stiffness of the various models. Population balance models, for example, require specific solution methods, and these need to run at par with other less computationally demanding models. This work therefore captures the research into the simultaneous simulation of these diverse models. In addition, process dynamics are included into the flowsheet model. This allows tracing the material properties throughout the entire line. This feature is included as the foundation of applying the model to assess the propagation of process disturbances, with respect to the product quality, i.e., which products needs to be discarded, or how fast can the process recover and return to a position where product critical quality attributes (CQAs) are within specification limits.

Previous work on the development of flowsheet models are restricted to direct compaction [3,4,5,6] and dry granulation routes [7,8] for continuous solid oral dosage manufacturing. Park et al. [7] created a flowsheet model of continuous dry granulation and applied it for optimization. Boukouvala et al. [9] developed a flowsheet model for the wet granulation route. This model served as a proof-of-concept and with no connection to specific experimental data. Boukouvala and Ierapetritou [10] also investigated a methodology for optimization of computationally expensive flowsheet models. Rogers and Ierapetritou [11] showed a flowsheet modeling case with hybrid models incorporating information from both detailed and reduced-order models.

The work presented in this manuscript includes models systematically developed based on experiments on units in the ConsiGma${}^{\mathrm{TM}}$-25 line for continuous tablet manufacturing using the same formulation and relevant materials across all the units. Specifically, the units involved are feeders, blenders, twin-screw wet granulator (TSWG), fluid bed dryer (FBD), comill, and tablet press. Models for these processes are developed [12,13,14,15,16] and included. Besides these, models for intermediate feeding and blending operations are also included. Transfer lines that lead to material holdup in between the units are added to the flowsheet model as well.

Objectives

The specific objectives of this work are:

• Develop a flowsheet model approximating the ConsiGma${}^{\mathrm{TM}}$-25 wet granulation manufacturing line;
• Demonstrate the use of the flowsheet model for simulating effects of disturbances in the continuous process;
• Identify critical process parameters (CPPs) affecting the properties of intermediate and final product.

Section 2.1 details various models used to build the flowsheet model. To enunciate how the flowsheet model can be used for propagation of information and disturbances across the units, a detailed discussion is provided in Section 3.1 along with supporting case studies in Section 3.2. In addition to building the flowsheet model, a detailed analysis of the developed model is provided. The scenario analysis, explained in Section 2.2 and Section 3.3, serves to ensure that the flowsheet model which is a complex set of equations from various modeling approaches, runs successfully at several values of process variables and the process responses thus obtained are aligned with process knowledge. Following this, critical process parameters (CPPs) that affect product quality are identified through implementation of sensitivity analysis as explained in Section 2.3 and Section 3.4.

2. Materials and Methods

The model formulation consisted of two active pharmaceutical ingredients (API), a lubricant and four excipients. Hereafter, the two two active pharmaceutical ingredients (API)s and the four excipients are referred to as API 1, API 2, and Excipient A, B, C, and D respectively. The formulation was processed using demineralized water as granulation liquid. The formulation used in this work is given in

Table 1. Formulation used for model development.

Component Name	Weight %
API 1	75.58
API 2	8.72
Lubricant	0.58
Excipient A	6.05
Excipient B	1.51
Excipient C	6.05
Excipient D	1.51

.

2.1. Flowsheet Modeling

Since flowsheet models are approximate representations of the integrated manufacturing line, developing individual unit operation models aid in the development of a flowsheet model. In this work, the models developed are based on experiments conducted using the ConsiGma™-25 system (GEA Pharma systems, Collette, Wommelgem, Belgium), which is an oral solid dosage manufacturing line based on continuous wet granulation. In Section 2.1.1, Section 2.1.2, Section 2.1.3, Section 2.1.4, Section 2.1.5, Section 2.1.6, Section 2.1.7 and Section 2.1.8 the individual unit models used in this work are briefly described. Figure 1 pictorially shows transfer of information across the unit operations i.e., feeder, blender, granulator, dryer, mill, and tablet press, in that order. Further, intermediate units are added and relevant information is transferred. These models are implemented in the software gPROMS FormulatedProducts v1.2.1 (PSE, London, UK). It is a platform for flowsheet simulations that uses an equation-oriented approach. An overview on the equation-oriented techniques particular to the gPROMS platform is given in Pantelides et al. [17].

2.1.1. Feeder

Loss-in-weight (LIW) feeders are used to feed the required powder components in the continuous manufacturing line. The feeder used in this work has a hopper and a conveying unit, a refill unit, and a PID controller. The hopper is used as a receptacle for the raw materials whereas the conveying system has a rotating screw that is used to move the material out of the feeder. The refill unit is used to feed material into the hopper when the fractional fill level in the hopper drops below a setpoint value. The PID controller enables the feeder to run in gravimetric mode, i.e., the screw speed in the conveying unit is adjusted to maintain a constant mass flow rate out of the unit. The three models (refill, feeding, and PID controller) work in conjunction to represent the overall feeding operation in the continuous line.

The mass flow rate out of the feeder is simulated using a feed factor model [3]. Feed factor is a time-dependent property ($ff\left(t\right)$), defined as maximum mass of powder fitting in a screw flight and has the unit of mass per screw revolution. It is found to be dependent on the amount of material in the hopper ($w\left(t\right)$) as given in Equation (1). The parameters $f{f}_{max}$, $f{f}_{min}$, and $\beta$ are dependent on the powder bulk density, compressibility, cohesion and permeability. More details on the feed factor model and its dependence on material properties are published in [18]. The mass flow rate of the powder out of the feeder can then be obtained as given in Equation (2) where $\omega \left(t\right)$ is the screw speed that is manipulated by the PID controller.

$$ff\left(t\right)=f{f}_{max}+(f{f}_{min}-f{f}_{max})exp(-\beta w\left(t\right))$$

(1)

$${\dot{M}}_{out}\left(t\right)=ff\left(t\right)\omega \left(t\right)$$

(2)

2.1.2. Blender

Continuous blenders that are used to mix the powder components in the continuous line, dampen the flow rate variations from the feeding units. The build-up of mass in the blender $M\left(t\right)$ was found to be following a first order relationship as given in Equation (3), where $M_{ss}$ is the steady state mass holdup and $\tau$ is the time constant. From the mass holdup and flow rate of the material into the blender (${\dot{M}}_{in}$), flow rate out of the unit (${\dot{M}}_{out}$) can be computed as given in Equation (4). An axial dispersion equation [3] is used to model the mixing calculation in the blending unit as a function of time. The equation as given in (5) is subject to initial and boundary conditions (Equation (6)). The coefficients of the axial dispersion model (${\tau }_{ax}$ and $Pe$) are calculated based on their relationship to a CSTR-in-series model constant i.e., number of tanks $n_t$ as given in Equation (7). Experimental data were used to develop regression models that predict the model constants $\tau$, $n_t$, and $M_{ss}$ as a function of flow rate and blade speed. More details on the blender model are available in [18].

$$\tau \frac{dM\left(t\right)}{dt}+M\left(t\right)=M_{SS}$$

(3)

$$\frac{dM\left(t\right)}{dt}={\dot{M}}_{in}-{\dot{M}}_{out}$$

(4)

$${\tau }_{ax}\frac{d{C}_{out}^i}{dt}=\frac{1}{Pe}\frac{{\partial }^2{C}_{out}^i}{\partial {\xi }^2}-\frac{\partial C_{out}^i}{\partial \xi }$$

(5)

$$\begin{array}{l}\hfill I.C:C_{out}^i=0,t=0\\ \hfill B.C:C_{out}^i=C_{in}^i,\xi =0\\ \hfill \frac{d{C}_{out}^i}{dt}=0,\xi =1\end{array}$$

(6)

$$Pe=n_t+{(8n_t+n_t^2)}^{1/2}$$

(7)

2.1.3. Twin-Screw Wet Granulator

The modeling of the change in particle size distribution (PSD) of the material through granulation is executed with the compartmental population balance model of Van Hauwermeiren et al. [13]. The twin-screw wet granulator (TSWG) model mathematically links the aggregation and breakage behavior in the granulator barrel to the granulator process settings of mass flow rate, screw speed, and liquid flow rate. It moreover distinguishes two compartments in the barrel: the wetting zone (i.e., the zone where the liquid is added to the dry powder blend) where only aggregation occurs, followed by the kneading zone (i.e., kneading elements are present in the screws) with different aggregation behavior complimented with breakage. Each compartment is thus modeled by its own population balance model (PBM).

The PBM equation is given in Equation (8). The change in number of particles n of a certain size x over time t is thereby described based on aggregation kernel $\beta (t,x,\epsilon )$, breakage selection function $S\left(\epsilon \right)$ and breakage fragment distribution $b(x,\epsilon )$. The aggregation kernel $\beta (t,x,\epsilon )$ can be modeled as the product of aggregation efficiency ${\beta }_0$ with collision frequency $\beta (x,\epsilon )$, as the relation is in this case independent of time t.

The formula of the collision frequency in the first PBM, describing the wetting zone, is given in Equation (9). It comprises a two-dimensional stepping function and a product kernel in order to reach bimodal granule (PSDs) starting from a monomodal powder (PSD). Kernel parameters ${\beta }_0$, $R_1$, $R_2$, $\mathrm{top}{}_1$, $\mathrm{top}{}_2$, ${\delta }_1$, and ${\delta }_2$ are needed to achieve this mathematical connection [13].

$$\begin{array}{ll}\hfill \frac{\delta n(t,x)}{\delta t}=& \frac{1}{2}{\int }_0^x\beta (t,x-\epsilon ,\epsilon )n(t,x-\epsilon )n(t,\epsilon )d\epsilon \hfill \\ & -n(t,x){\int }_0^{\infty }\beta (t,x,\epsilon )n(t,\epsilon )d\epsilon \hfill \\ & +{\int }_x^{\infty }b(t,x,\epsilon )S(t,\epsilon )n(t,\epsilon )d\epsilon -S(t,x)n(t,x)\hfill \end{array}$$

(8)

$$\begin{array}{ll}\hfill \beta (x,\epsilon )=& \left(\frac{{top}_1}{2}\left(1+tanh\left(\frac{R_1^3-{\left(x^2+{\epsilon }^2\right)}^{1/2}}{{\delta }_1}\right)\right)\right.\hfill \\ & -\frac{{top}_1-{top}_2}{2}\left(1+tanh\left(\frac{R_2^3-{\left(x^2+{\epsilon }^2\right)}^{1/2}}{{\delta }_2}\right)\right))\hfill \\ & ·\left(x^{1/3}·{\epsilon }^{1/3}\right)\hfill \end{array}$$

(9)

The breakage in the kneading zone is modeled by a linear breakage selection function:

$$S\left(\epsilon \right)=S_0{\epsilon }^{1/3},$$

(10)

and a breakage fragment distribution $b(x,\epsilon )$ describing a combination of erosion and uniform breakage:

$$b(x,\epsilon )=f_{prim}\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{{\left(x^{1/3}-\mu \right)}^2}{2{\sigma }^2}}\frac{\epsilon }{{\mu }^3}\frac{1}{3x^{\frac{2}{3}}}+\left(1-f_{prim}\right)\frac{2}{\epsilon },$$

(11)

with $\epsilon$ the volume of the breaking particle, x the volume of the fragment, $S_0$ the breakage rate constant, $\sigma$ and $\mu$ respectively the standard deviation and mean of the Gaussian distribution representing the size distribution of the small eroded particles, and $f_{prim}$ the volume fraction of erosion in the overall breakage (as opposed to a fraction $(1-f_{prim})$ of uniform breakage).

Aggregation in this zone could be described by a sum kernel (Equation (12)).

$$\beta (x,\epsilon )=x+\epsilon$$

(12)

Overall, kernel parameters $R_2$, ${\beta }_0$, and $to{p}_1$ in the wetting zone and ${\beta }_0$, $S_0$, $f_{prim}$, $\mu$, and $\sigma$ in the kneading zone are linearly related to the process setting values in the units given in

Table 2. Process setting ranges of the validated twin-screw wet granulator (TSWG) model [13].

Process Setting	Lower Bound	Upper Bound
Mass flow rate (kg/h)	10	20
Screw speed (RPM)	500	900
Liquid/solid-ratio (kg/kg )	0.08	0.18

.

2.1.4. Dryer

The fluidized bed dryer (fluid bed dryer (FBD)) model consists of prediction of granule batch drying kinetics based on single granule drying kinetics for one dryer cell [15]. The single granule drying kinetics are governed by Stefan diffusion of water vapor through the granule pores, from the source of evaporation to the edge of the granule. The mass transfer rates are corrected with the equilibrium moisture content $X_e$. When the moisture content of the granule is larger than its pore fraction, the remaining liquid is modeled as a layer of water of uniform thickness around the granule, evaporating according to a droplet:

$${\dot{m}}_v=h_D({\rho }_{v,s}-{\rho }_{v,\infty })A_d$$

(13)

with mass transfer rate ${\dot{m}}_v$ ($\mathrm{k}\mathrm{g} {\mathrm{s}}^{-1}$), mass transfer coefficient $h_D$ ($\mathrm{m} {\mathrm{s}}^{-1}$), partial vapor density over the droplet surface ${\rho }_{v,s}$ ($\mathrm{k}\mathrm{g} {\mathrm{m}}^{-3}$), partial vapor density in the ambient air ${\rho }_{v,\infty }$ ($\mathrm{k}\mathrm{g} {\mathrm{m}}^{-3}$), and droplet surface area $A_d$ (${\mathrm{m}}^2$). The energy balance paired with this drying behavior is described by:

$$h_{fg}{\dot{m}}_v+c_{p,w}{m}_d\frac{d{T}_d}{dt}=h(T_g-T_d)4\pi R_d^2$$

(14)

with specific heat of evaporation $h_{fg}$ ($\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}$), specific heat capacity of the liquid $c_{p,w}$ ($\mathrm{J}\mathrm{k}{\mathrm{g}}^{-1} {\mathrm{K}}^{-1}$), droplet mass $m_d$ ($\mathrm{k}\mathrm{g}$), uniform droplet temperature $T_d$ ($\mathrm{K}$), heat transfer coefficient h ($\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$), drying gas temperature $T_g$ ($\mathrm{K}$), and droplet radius $R_d$ ($\mathrm{m}$). After this layer of water is depleted the wet granule enters the subsequent drying phase. Herein the moisture is conceptualized as a sphere with radius $R_i$ ($\mathrm{m}$), also referred to as the wet core, filling up the pore volume of the granule with radius $R_p$ ($\mathrm{m}$). The mass transfer rate ${\dot{m}}_v$ in this stage is given by [19]:

$${\dot{m}}_v=-\frac{8\pi {ϵ}^{\beta }{D}_{v,cr}{M}_w{p}_g}{\Re (T_{cr,s}+T_{wc,s})}\frac{R_p{R}_i}{R_p-R_i}ln[\frac{p_g-p_{v,i}}{p_g-(\frac{\Re }{4\pi M_w{h}_D{R}_p^2}{\dot{m}}_v+\frac{p_{v,\infty }}{T_g})T_{p,s}}]$$

(15)

with $ϵ$ the granule porosity (-), $\beta$ an empirical coefficient, $D_{v,cr}$ the vapor diffusion coefficient (${\mathrm{m}}^2{\mathrm{s}}^{-1}$), $M_w$ the liquid molecular weight ($\mathrm{k}\mathrm{g}{\mathrm{mol}}^{-1}$), $p_g$ the pressure of the drying air ($\mathrm{Pa}$), ℜ the ideal gas constant ($\mathrm{J}{\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$), $T_{cr,s}$ and $T_{wc,s}$ respectively the temperature of solids at the granule surface and at the gas-liquid interface ($\mathrm{K}$), $p_{v,\infty }$ and $p_{v,i}$ respectively the partial vapor pressure in the drying air ($\mathrm{Pa}$) and at the gas-liquid interface, and $T_{p,s}$ the temperature of the particle solids ($\mathrm{K}$). Equation (14) is assumed to apply for the energy balance of the granule during this drying phase. The physical properties of the solids were assumed to be at environment conditions of 25 ${}^{\circ }\mathrm{C}$ and atmospheric air pressure in the model, whereas liquid properties were modeled as those of pure water. A value of 35% was assumed for the granule porosity. Finally, the mass transfer rates ${\dot{m}}_v$ are corrected with the effect of $X_e$:

$${\dot{m}}_{v,res}=\frac{X-X_e}{X}{\dot{m}}_v$$

(16)

Overall this means that the course of moisture content X of a granule amounts to:

$${\dot{X}}_{SPDM}=\frac{{\dot{m}}_{v,res}}{m_p}$$

(17)

with ${\dot{X}}_{SPDM}$ the change in moisture content X over time and $m_p$ the total mass of the granule (particle). This is solved from time $t=0$ until the FBD cell drying time $t=t_{dry}$.

Connecting the single granule drying kinetics to those of the batch, along with the continuous filling of the batch, is done according to the following simplified approach. Drying curves are calculated for several size fractions of the granules, in which the arithmetic mean size is representative in the $X_{SPDM}$. Thus the average moisture content ${\overline{X}}_f$ per size fraction f then equals the average moisture content of different drying curves, starting a certain time $\tau$ later over the cell filling time $t_{fill}$:

$${\overline{X}}_f=\frac{{\sum }_{\tau =0}^{t_{fill}}{X}_{SPDM}(t-\tau )}{n_{\tau }},$$

(18)

with $n_{\tau }$ the amount times $\tau$ that the drying curve was shifted over the cell filling time interval $t_{fill}$.

A constant ideal fluidization behavior, constant relative air humidity in the dryer cell and atmospheric air properties in the drying chamber (with exception of the drying agent temperature) were assumed in the batch approach.

The dryer model discussion so far dealt with the drying behavior of the material, other material properties are directly governed by the dynamic output of the TSWG model. For each dryer cell, these are mass-averaged over the drying cell filling period. This is illustrated in Equation (19) for the concentration of active pharmaceutical ingredients (API) in the dryer cell $C_{API,dryer}$, and is calculated the same way for the material true density, PSD and the mass of the material in the dryer cell. They are weighted by the mass flow rate at which they are flowing from the TSWG at time t ($MF{R}_{TSWG}\left(t\right)$).

$$C_{API,dryer}=\frac{{\int }_t^{t+t_{fill}}{C}_{API,TSWG}\left(u\right)MF{R}_{TSWG}\left(u\right)du}{{\int }_t^{t+t_{fill}}MF{R}_{TSWG}\left(u\right)du}$$

(19)

After drying, the breakage of the material through pneumatic transport through a tube to the evaluation module is calculated by a PBM. As aggregation is assumed not to take place based on the experimental work of De Leersnyder et al. [14], only the last two terms in the right hand side of Equation (8) need to be used. The same breakage kernel as in the kneading zone of the TSWG was found to apply, i.e., breakage rate $S\left(\epsilon \right)$ from Equation (10) and breakage fragment distribution $b(x,\epsilon )$ from Equation (11). Parameter $S_0$ from from Equation (10) is linearly related to the remaining moisture content after drying.

Finally, the six cells are simulated in parallel with the filling–drying–emptying sequences according to the operation in the actual process. A cell is idle until the predecessor cell has been filled, at which point the current cell enters the filling stage. Hereafter, the remainder of the drying time is completed in the drying stage, where the mass flow rate at time t ($MF{R}_{TSWG}\left(t\right)$) in that cell is zero. In the final emptying stage, the mass transfer rates in Equations (13) and (15) are set to zero, and the change in PSD is calculated by the PBM model based on the residual moisture content of the material. These material property values are thus those perceived at the inlet of the evaluation module.

2.1.5. Mill

In the wet granulation continuous manufacturing route, comilling is used to break the granulated product through collisions from a rotating impeller and walls. Granules that are broken to the required size exit the comill through a screen. In this work, the comill model published in Metta et al. [16] is used. Briefly, the mill model is a hybrid model that includes a PBM approach and a partial least squares (PLS) approach. Trajectories of change in mass of particles of various sizes over time is predicted through the PBM as shown in Equation (20) where $M(w,t)$ represents the mass of particles of volume w at time t, $R_{form}$ and $R_{dep}$ represent the rates of formation and depletion of particles respectively. ${\dot{M}}_{in}$ and ${\dot{M}}_{out}$ are the mass flow rates of particles entering and exiting the mill respectively. The rate of depletion $R_{dep}$ is defined in the model (Equation (21)) using a breakage kernel $K\left(w\right)$, which represents the probability that a particle of volume w undergoes breakage. A classification kernel as given in Equation (22) is used in this work, where $v_{imp}$ is the impeller speed, $v_{imp,min}$ is the minimum impeller speed and the parameter $\beta$ is calibrated using data from experiments. The rate of formation $R_{form}$ as shown in Equation (23) uses the breakage kernel and a breakage distribution function. The breakage distribution function $b(w,u)$ represents the distribution of daughter particles formed when a particle of volume w undergoes breakage. The Hill–Ng distribution function given in Equation (24) is used in this work, where the parameters p, q are estimated using experimental data. The mass flow rate out of the mill ${\dot{M}}_{out}(w,t)$ as given in Equation (25) is modeled using the feed particle size distribution ($d_{in}$) and a parameter $\Delta =d_{screen}\delta$, where $d_{screen}$ is the screen size and $\delta$ is referred to as critical screen size ratio. A linear model is used to define the function $f_d$. The parameter $\delta$ is formulated as given in Equation (26) which represents the phenomenon of reduced apparent screen size available for a particle to exit the mill as impeller speed increases.

$$\frac{dM(w,t)}{dt}=R_{form}(w,t)-R_{dep}(w,t)+{\dot{M}}_{in}(w,t)-{\dot{M}}_{out}(w,t)$$

(20)

$$R_{dep}(w,t)=K\left(w\right)M(w,t)$$

(21)

$$K\left(w\right)=\left\{\begin{array}{cc}\beta {\left(\frac{v_{imp}}{v_{imp,min}}\right)}^2\left(\frac{w}{w_{ref}}\right)\hfill & \mathrm{if}w\ge w\left(24\right)\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(22)

$$R_{form}(w,t)={\int }_w^{\infty }K\left(u\right)M(u,t)b(u,w)du$$

(23)

$$b(w,u)=\frac{p{\frac{u}{w}}^{q-1}{(1-\frac{u}{w})}^{r-1}}{wB(q,r)}$$

(24)

$${\dot{M}}_{out}(w,t)=(R_{form}(w,t)-R_{dep}(w,t)+\gamma d_{in}(w,t))(1-f_d)$$

(25)

$$\delta =ϵ{\left(\frac{v_{imp,min}}{v_{imp}}\right)}^{\alpha }$$

(26)

Impeller speed showed little effect on the milled product when the comill feed is obtained from the fluid bed dryer because of breakage that occurs during transport to and from the fluid bed dryer in the horizontal ConsiGma^TM-25 configuration. Hence, $v_{imp}$ is considered equal to $v_{imp,min}$ in the breakage kernel given in Equation (22).

The PBM is thus used to predict milled granule size distribution. The PLS model is an empirical modeling approach used to predict the milled product bulk density and tapped density, using the granule size distribution and moisture content as inputs. To use the mill model in the flowsheet model, batches of material are added to the milling unit each time the FBD completes a drying cycle and initiates an emptying cycle. When a batch of material is added to the existing material in the mill, breakage occurs. Properties of the material exiting the mill are obtained through mass averaging over the milling period $t_{mill}$. This is illustrated in Equation (27) for the bulk density of milled granules and is calculated the same way for the tapped density, span, and true density of the material exiting the mill. The instantaneous bulk density ${\rho }_{bulk,PLS}$ is obtained from the PLS model and weighted by the total mass flow rate of material exiting the mill, ${\dot{M}}_{out,total}$, at time t.

$${\rho }_{bulk,milled}=\frac{{\int }_t^{t+t_{mill}}{\rho }_{bulk,PLS}\left(x\right){\dot{M}}_{out,total}\left(x\right)dx}{{\int }_t^{t+t_{mill}}{\dot{M}}_{out,total}\left(x\right)dx}$$

(27)

2.1.6. Tablet Press

The tablet press model consists of four submodels. Firstly, a residence time distribution (RTD) model is used to describe the propagation of material properties through the powder dosing valve and the tablet press feed frame into the tablet die. The other models work with the material properties modeled as present in the die. The weight model relates material densities to the mean weight of the tablets semi-mechanistically (Equations (28) and (29)). The tablet potency model is a first-principles model (Equation (31)), and the tablet mean hardness model harbors literature empirical correlations related to tablet hardness and tensile strength (Equations (32) to (37)).

Propagation of the material properties through the feed frame is modeled according to a series of a continuously stirred tank reactor (CSTR) and a plug flow RTD model, with respective delay times $t_{cstr}$ and $t_d$ related to the feed frame turret speed. The solution of the feed frame model is performed in the same way as the transfer line models explained in Section 2.1.8. The tablet mean weight model uses the tooling dimensions (cup volume $V_{cup}$, die surface $A_{die}$, and punch cup depth $D_{cup}$), the material densities at the die (true density ${\rho }_{true}$, bulk ${\rho }_{bulk}$ and tapped ${\rho }_{tapped}$ density), a fill density factor $p_{{\rho }_{fill}}$, and the fill depth $d_{fill}$ to calculate the mean tablet weight $\overline{m}$, according to the following equations.

$${\rho }_{fill}={\rho }_{bulk}+p_{{\rho }_{fill}}\left({\rho }_{tapped}-{\rho }_{bulk}\right)$$

(28)

$$\overline{m}=(V_{cup}+A_{die}{d}_{fill}){\rho }_{fill}$$

(29)

The volume of the solids in the die $V_{solid}$ is then:

$$V_{solid}=\overline{m}/{\rho }_{true}$$

(30)

The potency P simply follows from the API concentration $C_{API}$ and $\overline{m}$:

$$P=\overline{m}{C}_{API}$$

(31)

The hardness model in the end uses $V_{solid}$ from the tablet mean weight model, the main compression height process setting $MCH$ and the tablet dimensions (width W, thickness T, and upper punch penetration $upp$) to calculate the tablet hardness through estimation of the tablet tensile strength. W and T are calculated according to:

$$W=MCH-upp$$

(32)

and

$$T=\left(MCH-upp\right)+2D_{cup},$$

(33)

resulting in a tablet volume $V_{tablet}$ and relative density ${\rho }_{rel}$ through:

$$V_{tablet}=\left(W{A}_{die}\right)+2V_{cup}$$

(34)

and

$${\rho }_{rel}=\frac{V_{solid}}{V_{tablet}}.$$

(35)

In order to relate the experimentally obtained hardness N values to the material properties for the biconvex tablets, a tensile strength ${\sigma }_T$ normalization needs to be applied w.r.t. the tablet dimensions [20].

$${\sigma }_T=\frac{10N}{\pi {\mathrm{D}}^2}{\left(2.84\frac{\mathrm{t}}{\mathrm{D}}-0.126\frac{\mathrm{t}}{\mathrm{W}}+3.15\frac{\mathrm{W}}{\mathrm{D}}+0.01\right)}^{-1}$$

(36)

The tensile strength is linked to the critical density for tableting of the material ${\rho }_c$, and maximum achievable tensile strength ${\sigma }_{max}$ [21]:

$$\sigma ={\sigma }_{max}\left[{\rho }_{\mathrm{c}}-{\rho }_{rel}-ln\left(\frac{1-{\rho }_{rel}}{1-{\rho }_{\mathrm{c}}}\right)\right].$$

(37)

The constants of ${\sigma }_{max}$ and ${\rho }_c$ are obtained through empirical linear correlations with the relative size span $R_{span}$ of the material and the moisture content X.

Finally, two modes of using the tablet press model are defined. In the process setting mode, tablet press process setting values (fill depth and MCH) are supplied to the model along with the material properties. Tablet properties such as the mean weight and hardness then result from the combination of the user defined process setting values with the incoming properties. These might however not match and result in unrealistic tablet properties as well as infeasible calculations, such as in the case where Equation (38) is not satisfied. In this mode, this is amended by the inclusion of a boolean variable that is only True when the input variables are such that condition Equation (38) is fulfilled.

$${\rho }_c<{\rho }_{rel}<1$$

(38)

Another mean weight control (MWC) mode was created to allow for the tablet press model to calculate its process setting values in order to achieve a certain mean weight and hardness. These two variables thus serve as an input to the model in the flowsheet model, and the fill depth and MCH are calculated. This allows avoidance of the situation where user-specified tablet press process settings are not compatible with the simulated material properties governed by the process setting values upstream, possibly leading to a long simulation where no tablet weights or hardnesses could be calculated.

2.1.7. Overview

Once the individual unit operation models are developed and parameterized, the flowsheet model is built via connecting the inlet of a unit to the outlet of the preceding unit. Required information (material properties, operating conditions, etc.) are transferred from one unit to the succeeding unit as shown in Figure 1. The interaction between various modeling approaches in the individual unit operations is pictorially shown in Figure 2 to illustrate how empirical, semi-empirical, statistical, and mechanistic models interact to yield a flowsheet model. The system is solved using a backward differentiation solver with variable time steps, which is one of the two built-in solvers in gPROMS.

2.1.8. Intermediate Feeders, Blender, and Transfer Lines

The ConsiGma${}^{\mathrm{TM}}$-25 line includes an intermediate feeder for feeding powder blend from blender to the granulator (hereafter referred to as ’powder blend feeder’), and an intermediate feeder for feeding granulated material from the mill (hereafter referred to as ’granule feeder’). Thus, to achieve a more accurate representation of the ConsiGma${}^{\mathrm{TM}}$-25 line, transfer lines and the intermediate feeder and blender units are also included in the flowsheet model. A key difference in the implementation of the feeder model to the intermediate units is that the powder blend feeder and granule feeder do not have refilling units assigned to them, as these units receive material from their preceding units (powder blend feeder receives material from the blender and granule feeder receives material from the mill).

In addition, a blending unit is also included to add lubricant to the granulated material from the mill. This intermediate blender will hereafter be referred to as ’granule blender’. The granule blender is modeled based on the axial dispersion equation as explained in Section 2.1.2. The blender is assumed to be filled with granulated material and the axial dispersion equation is used to transfer information regarding new or changing material properties from the upstream mill and the lubricant feeder. The true density of lubricated material from the granule blender is taken as weighted average of bulk densities of the milled granule product and lubricant.

In the dynamic flowsheet simulation, an accurate representation of the material properties at all times and at every location in the process would not be complete without accounting for material hold-up between unit operations. Hereto, transfer line models have been implemented to propagate variable values in between unit operation models. These models delay the propagated values according to a plug flow regime, assuming no back-mixing or axial dispersion is taking place during the transfer of materials between unit operations.

The plug flow propagation of these materials is often modeled applying a convolution of the inlet concentration profile at the modeled system, with the residence time distribution function of the material, as for instance described in [22]. This convolution requires information on the inlet concentration values over a range of time, which is not accessible for calculation in gPROMS. Therefore, this plug flow behavior is emulated using an axial dispersion model. A new simulation domain z is created to represent the normalized length of the transfer line under consideration, therefore this domain [0, L] is always equal to [0, 1]. Over this domain, the change in input signal S over time t at $z=0$ is propagated over the domain as in Equation (39) where the plug flow time delay is represented as ${\tau }_{delay}$. The input and output of the transfer lines are hence given by $S(t,0)$ and $S(t,1)$ respectively. The smaller $\partial z$ is chosen, the smoother and more accurate the signal will propagate through the domain z, yet more computational burden is involved with this choice of more grid points. A value of $1/1000$ for $\partial z$ has been found to give a good balance between smoothness and computational burden. Finally, it has been found that normalization of the delayed signal S drastically improves CPU time in the gPROMS solvers.

$${\tau }_{delay}\frac{dS(t,z)}{dt}=\frac{\partial S(t,z)}{\partial z}$$

(39)

The above unit operation models are connected to develop the flowsheet model representing the ConsiGma${}^{\mathrm{TM}}$-25 line as shown in Figure 3. The flowsheet includes intermediate units and transfer lines as well. When the flowsheet model is simulated, it is important to have a clear understanding of the initial states of the model as it impacts the dynamic model state. For example, feeder hoppers could have various amounts of material at the start of the simulation. Similarly, the mill could start empty or have a certain amount of mass held up at the start of the simulation.

Table 3. State of unit models used in flowsheet model development.

Unit	Initial State	Dynamic State
Feeders	Full	Refills when empty
Blender	Empty	Reaches steady state
Granulator	Empty	Output obtained when input is in studied range
		Plug flow delay added to instantaneous response
Dryer	Empty	Releases batches of material at the end of drying time
Mill	Empty	Releases material in semi-continuous mode
Tablet Press	Empty	Output obtained when input is in studied range
		Delay from RTD model in feed frame
Powder blend feeder	Full with powder blend	Continuous feed from blender
Granule feeder	Full with granules	Refill from mill
Granule blender	Always full	Delay from axial dispersion

gives an overview of initial and dynamic states of the various units in the developed flowsheet model. The full flowsheet model thus developed can be used to simulate the continuous process as described in Section 3.1. Two case studies are provided in Section 3.2 where the full flowsheet model is used to understand the effect of step changes in process settings. The flowsheet model can also be used for advanced process analyses as described in the next sections.

2.2. Scenario Analysis

It is important to ensure that the flowsheet model developed successfully runs simulations at various process conditions. This affirms that the models are successfully integrated. In addition, it is important to verify that the process responses from these simulations are aligned with process knowledge. Scenario analysis provides a structured framework to achieve this, as several simulations at various combinations of process settings can be run in parallel and the resulting process responses can be analyzed. In this work, only process responses at the end of simulation, i.e., steady state tablet properties, are analyzed.

However, this exercise is computationally demanding as one simulation takes several hours to run. In this work, a more pragmatic approach is applied and the flowsheet model as given in Figure 3 is adapted in order to implement scenario analysis. Specifically, the intermediate units (powder blend feeder, granule feeder, granule blender) and the transfer lines are not considered. Since only the end of simulation responses are studied, the flowsheet model simplification is valid as the intermediate units and transfer lines do not affect the steady state process responses. The modified model lowers the computational expense and allows implementation of study required to affirm that the process responses from integration of unit operation models are meaningful. The full flowsheet model takes approximately 4 h to run a six-cell drying cycle, whereas the modified model only takes approximately 7 min for the tablet properties to reach steady state. Hence, the adapted model as shown in Figure 4 is used for the implementation of scenario analysis. In addition, the factors used for the study are also chosen judiciously in order to keep the total number of flowsheet model evaluations low. For example, factors such as blender blade speed, mill impeller speed are not considered. The blender blade speed influences the mass holdup in the blender, as explained in Section 2.1.2, but not the steady state flow rate, as this depends on the incoming feed flowrate. The mill impeller speed does not have an effect on granules obtained from a FBD [16]. Overall, the factors considered, their corresponding lower and upper bounds, and number of levels for each factor is listed in

Table 4. Factors and responses for scenario analysis and sensitivity analysis.

Unit	Factor	Bounds	Number of Levels	Response
Blender	flowrate setpoint, kg/h	[10, 20]	3	Mean Residence time
				Number of tanks
Granulator	LS ratio, kg/kg	[0.08, 0.18]	4	PSD: d10, d50, d90
	Screw speed, rpm	[500, 900]	4	Moisture content
Dryer	Air temperature, deg C	[40, 60]	4	PSD: d10, d50, d90
	Drying time, s	[200, 1080]	4	Moisture content
Mill				PSD: d10, d50, d90
				Span
				Bulk density
				Tapped density
Tablet Press				Tablet hardness
				Tablet potency
				MCH
				Fill depth

. In addition, the process responses that are recorded for each simulation run are also listed in Table 4. Three levels for flow rate setpoint, and four levels for LS ratio, granulator screw speed, dryer air temperature, and drying time, are chosen. In Section 3.3, an analysis of process responses from the simulations run is discussed in detail.

While the adapted flowsheet model is useful in implementing further analyses, it is also important to understand its limitations. The simplified flowsheet does not capture the effects of the intermediate units. For example, effects of refilling and propagation of disturbances from the intermediate units are ignored. In addition, a successful scenario analysis on the simplified model does not capture scenarios where the full flowsheet model fails due to the intermediate units. The simplified flowsheet does not support study of dynamic behavior of the line. Analyses such as dynamic sensitivity analysis or identification of dynamic feasible region cannot be accomplished.

2.3. Sensitivity Analysis

Sensitivity analysis is one of the key process systems engineering tools that can be used for quality risk assessment through identification of critical process parameters (CPPs). Sensitivity analysis is the investigation of how variability in the model inputs contributes to variations in model outputs [23]. It is an effective tool to rank and prioritize the process variables based on the effects they have on the output variables of interest. In the context of operation of a continuous manufacturing line, it helps identify the source of issues in meeting product quality or production demands. In the context of process model development, sensitivity analysis has been extensively used to identify the parameters that affect model outputs, thus helping focus experimental and model calibration efforts [23,24]. This helps researchers identify the areas where further model development needs to be focused on. Specifically for a flowsheet model where there is high number of input factors, sensitivity analysis can be used to reduce the number of input factors that need to be studied further. This helps in simplifying a high dimensional problem by filtering out the variables that have negligible effects on the outputs of interest. With this simplification, other tools can be applied for identification of design space of the process and its optimization [3].

Sensitivity analysis can be categorized into local and global methods. Local methods study the effect of input variables around a nominal point (or base case), whereas global methods study the effects over an entire input space. In this work, we focus on global sensitivity analysis as this is more relevant for pharmaceutical processes. For pharmaceutical processes, the input factors may include operating variables such as blender impeller speed, granulator screw speed, etc. The output variables of interest may include product properties such as tablet hardness, granule mean particle sizes, or process variables such as mill mass hold up, total flow rate, etc. There are various global sensitivity analysis methods available. The choice of method usually depends on the computational cost of evaluating the models, sampling budget available, and the detail of sensitivity information desired. In this work, Elementary effects method and Variance based sensitivity analysis methods are used, details of which are described in the next section. These methods are chosen as they are available in gPROMS FormulatedProducts with parallel computing capability.

2.3.1. Morris Method

Morris method also referred to as Elementary effects method is categorized under screening methods for sensitivity analysis. Screening methods are the most effective way to identify the most influential factors with relatively fewer samples [23]. Morris method is based on OAT (one-at-a-time) design where each of the input factors is varied and effects on the model outputs are studied.

For a model with k number of inputs, at a selected base point ($x_1$, $x_2$, …, $x_k$), the elementary effect $E{E}_i$ of the ith factor is given by Equation (40) where ${\Delta }_i$ is the step change in the ith input factor, Y represents the model output and 0≤${\Delta }_i$≤1. In order to represent the sensitivity information accurately, the sample points must be spread in the input space.

$$E{E}_i=\frac{Y(x_1,x_2,\dots ,x_i+{\Delta }_i,\dots ,x_k)-Y(x_1,x_2,\dots ,x_i,\dots ,x_k)}{{\Delta }_i}$$

(40)

Based on the calculation of $E{E}_i$, the sensitivity metrics as given in Equations (41)–(43) can be calculated, where r is the number of trajectories or radial base points for sampling. ${\mu }_i$ represents the average $E{E}_i$, ${\sigma }_i^2$ represents the variance and reflects non linearity or interactions in the ith input. ${\mu }_i\ast$ represents the average elementary effect using absolute $E{E}_i$ to ensure the negative and positive effects do not cancel each other. It is suggested to look at all three metrics together to understand sensitivity information. Total sampling cost for this method is $r(k+1)$ where r can be less than 20 [25]. Hence, it is especially useful for models with a large number of input factors (factor of ten) or when the model is computationally expensive. In this work, the value of r is chosen as 20.

$${\mu }_i=\frac{1}{r}\sum _{j=1}^{r}E{E}_i^j$$

(41)

$${\sigma }_i^2=\frac{1}{r-1}\sum _{j=1}^r{(E{E}_i^j-{\mu }_i)}^2$$

(42)

$${\mu }_i\ast =\frac{1}{r}\sum _{j=1}^r\mid E{E}_i^j\mid$$

(43)

After the metrics are obtained, input factors with large ${\mu }_i$ and/or ${\mu }_i\ast$, ${\sigma }_i^2$ are considered to be significant. Practically, if the metric of an input factor is less than 10% of the largest value of this metric, the input factor is considered insignificant. While the method can be used to rank the factors, it does not quantify how much an input factor is more important than the other factors.

2.3.2. Variance Based Method

In this category of methods, the variance of the output is decomposed into several components including the individual inputs and the interactions between the inputs [26]. For an independent set of input factors, the variance $V\left(y\right)$ is expressed as given in Equation (44), where $V_i$ is the variance term solely due to the input factor $x_i$, and $V_{i,j}$ is the variance term due to the interaction between the input factors $x_i$ and $x_j$. Based on this variance decomposition, sensitivity measures can be defined as given in Equations (45)–(47).

$$V\left(y\right)=\sum _{i=1}^k{V}_i+\sum _{1\le i<j\le k}{V}_{i,j}+\dots +V_{i,j,\dots ,k}$$

(44)

$$S_i=\frac{V_i}{V\left(y\right)}$$

(45)

$$S_{ij}=\frac{V_{i,j}}{V\left(y\right)}$$

(46)

$$S_{Ti}=\frac{V_i+{\displaystyle \sum _{j\ne i}}{V}_{i,j}+V_{1,2,..k}}{V\left(y\right)}=1-\frac{V_{\sim i}}{V\left(y\right)}$$

(47)

For the input $x_i$, $S_i$ represents the ‘first-order sensitivity index’ whereas $S_{ij}$ represents the ‘second-order sensitivity index’ which is the interaction effect of $x_i$ and $x_j$ on the process output. The metric $S_{Ti}$ indicates the ‘total sensitivity index’, which accounts for the main effects as well as all the higher order interaction effects.

Specifically, this method uses Monte–Carlo techniques to compute the sensitivity indices as given in Equations (48) and (49) where $E(.)$ is the expected value and $X_{\sim i}$ represents all possible combinations of input factors with ith input factor $X_i$ fixed. For this method, total number of samples required is $N(k+2)$ where N is recommended to be at least 500 [23]. In this work, the value of N is chosen as 500.

$$S_i=\frac{V_{X_i}\left(E_{X_{\sim i}}(Y\mid X_{\sim i})\right)}{V\left(Y\right)}$$

(48)

$$S_{Ti}=\frac{V_i+{\displaystyle \sum _{j\ne i}}{V}_{i,j}+V_{1,2,..k}}{V\left(y\right)}=1-\frac{V_{X_{\sim i}}(E_{X_i}(Y\mid X_{\sim i})}{V\left(y\right)}=\frac{E_{X_{\sim i}}(V_{X_i}(Y\mid X_{\sim i})}{V\left(y\right)}$$

(49)

For the variance based method, higher values of the metrics $S_i$ and $S_{Ti}$ indicate larger influence of the input factor. Also, $S_i$ is always lower than or the same value as $S_{Ti}$. Hence, the difference between these metrics reflects interaction effects of the input factor with other factors. The adapted flowsheet model as explained in Section 2.2 is used for implementing the sensitivity analysis methods as well. In this work, details of factors and responses listed in Table 4 are also applicable for executing sensitivity analysis. In the next sections, results from simulating the flowsheet model, case studies to demonstrate dynamic simulation capabilities of the model and further analysis are presented and discussed.

3. Results and Discussion

3.1. Simulation Results

The flowsheet model as described in Section 2.1 is simulated using process settings as given in

Table 5. Values of process variables used for simulating the flowsheet model.

Unit	Process Variable	Units	Value
Feeders	API 1 flow rate	kg/h	11.337
	API 2 flow rate	kg/h	1.308
	Lubricant flow rate	kg/h	0.087
	Excipient A flow rate	kg/h	0.907
	Excipient B flow rate	kg/h	0.227
	Excipient C flow rate	kg/h	0.907
	Excipient D flow rate	kg/h	0.227
	Powder blend feeder flow rate	kg/h	14.913
	Granule feeder flow rate	kg/h	14.913
Blenders	Bladespeed	rpm	250
Granulator	Liquid-solid ratio	kg/kg	0.12
	Screw speed	rpm	500
Dryer	Air flow	m${}^3$/h	360
	Air temperature	deg C	40
	Drying time	s	450
	Filling time	s	180
Tablet Press	Turret speed	rpm	29.8
	Mean weight	g	0.43

. The flow rate setpoints for the feeders are based on the formulation given in Table 1 and a total flow rate of 15 kg/h. The operating variables for other units were set within the ranges that were used to develop the individual unit operation models. The flowsheet model is simulated for 1500 s to complete a drying cycle using six dryer cells. Feeder levels in the seven component feeders (two API, one lubricant, and four excipient feeders) decrease until refill occurs at a fractional fill level of 0.1. The fill level in the powder blend feeder reduces until the blender starts feeding powder blend to it. The fill level in the granule feeder reduces until the mill starts feeding granules in a semi-continuous mode. Fill levels of all the component feeders and flow rate are shown in Figure 5a,b, respectively. Figure 5c shows fill levels and flow rates of the two intermediate feeders, i.e., powder blend and granule feeder.

Powder blend flows continuously to the granulator. Wet granules start filling the first cell of the dryer and drying begins. After a filling time of 180 s, the second dryer cell starts filling. Dried granules in the first cell are emptied to the mill after a total drying time of 450 s. Thus, the cycle of filling, drying and emptying continues for the duration of simulation. Batches of dried granules fed to the mill are broken and leave the mill. The left axis in Figure 6a shows flow rate of the blend from the granulator. The right axis in Figure 6a shows mass of the batches of dried granules fed to the mill and corresponding change in holdup in the mill due to granules entering and leaving the mill. Figure 6b shows the evolution of moisture content of wet granules from the granulator and dry granules from the dryer as simulation progresses. As the simulation of the drying behavior requires a moisture content value at the beginning of its simulation, the results for the first dried batch are not representative because the wet granule moisture content at time $t=0$ s equaled zero. This is the moisture content output of the FBD-model at time $t=450$ s, which just indicates that the first dryer cell has emptied. This should be improved towards the future so that the FBD-model has a fully dynamic response towards its input (see Section 4).

As milled granules exit the mill, properties of the batch of milled granules are mass averaged and this information is propagated to the subsequent units. Figure 7a shows profiles of bulk density, tapped density, and true density of milled granules as the simulation progresses. Profiles of LOD, span, and API concentration of the milled product are shown in Figure 7b. It can be observed that the true density, LOD, and API concentration of milled product shows a step change around 600 s. This is because the intermediate powder blend feeder initially contains powder blend with a true density of 1291 kg/m${}^3$, at the start of the simulation. This is eventually replaced when the powder blend of true density of 1344 kg/m${}^3$, coming from the unit operations upstream, reaches that intermediate feeder in the simulation. Similarly, powder blend containing no API is replaced by a powder blend containing API 1 and API 2. This is reflected in the API concentration profile as shown in Figure 7b. In addition, LOD of the first batch of dried granules is lower (as shown in Figure 6b), which reflects in the milled granule LOD profile as well.

Milled granules fed to the granule feeder eventually replace granular material in it. Milled granules exiting the granule feeder are mixed in the granule blender with lubricant. Granules thus lubricated are sent to the feed frame, which is modeled as a PFR and CSTR in series. Propagation of properties (bulk density, tapped density, true density) in the granule feeder is shown in Figure 8a. The density profiles in this figure shows replacement of granules in the granule feeder with granules from the mill. Similarly, Figure 8b shows LOD, span, and API concentration profiles that simulate replacement of material existing in the granule feeder (2% LOD, span of 2, and 0.7 fractional API concentration) with milled granules from the upstream unit. Figure 9a shows density profiles of granules entering and leaving the feed frame. Similarly, Figure 9b shows profiles of LOD, span, and API concentration of granules entering and leaving the feed frame. The profiles changes seen in these figures is self explanatory based on the milled granule profiles (Figure 7a,b) and granule feeder profiles (Figure 8a,b).

The propagation of bulk density, tapped density, true density, API concentration, LOD, and span affect the profiles of tablet properties, namely, tablet hardness, weight, and potency. Figure 10a shows dynamic evolution of tablet properties as powder from component feeders replace material existing in the intermediate units (powder blend feeder and granule feeder). The tablet press hardness model was developed for material with bulk density greater than 300 kg/m${}^3$. Hence, an initial tablet hardness of zero is shown in the hardness profile. In addition, since the tablet press is used in a mean weight control mode, main compression height and fill depth are adjusted as shown in Figure 10b in order to make tablets with a weight of 0.43 g.

3.2. Case Study

To clearly demonstrate the use of the flowsheet model for dynamic simulation purposes, two case studies are presented in this section. In both case studies, the full flowsheet model developed as explained in Section 2.1 is used. The initial simulation process variable settings are the same as explained in Section 3.1, Table 5. Later, in case study 1, a step change from 1.3 kg/h to 4 kg/h is made to the API 2 feeder flow rate setpoint at 200 s. In case study 2, a step change from 40 ${}^{\circ }\mathrm{C}$ to 50 ${}^{\circ }\mathrm{C}$ is made to the dryer air temperature at 455 s. A comparison between simulation results with and without the step changes implemented, and a demonstration of propagation of effects of the changes is discussed in the following sections.

3.2.1. Case Study 1: Step Change in Feeder Flow Rate

A step change in API 2 feeder flow rate setpoint from an initial value of 1.3 kg/h to 4 kg/h is expected to lead to a change in API concentration (API 1 and API 2) in the powder blend, the granules, and subsequently potency of the tablets. Figure 11 shows the step change made in setpoint at 200 s leads to a change in API 2 feeder flow rate.

A change in fractional API concentration at the blender outlet from 0.85 to 0.87 is seen in Figure 12a as a result of the step change. Figure 12a also shows a change in API concentration of the powder blend leaving the powder blend feeder as initial material in the feeder (with no API) is eventually replaced by powder blend with API concentration of 0.87. A comparison with profiles from the simulation explained in Section 3.1 is also shown in this. Since the change is made at 200 s, the first cell in the dryer is filled (dryer filling time = 180 s) with granules from powder blend already present in the powder blend feeder. Hence, profiles of API concentration from the outlet of the granulator and dryer from the simulation explained in Section 3.1 and this case study are the same until the first dryer cell is emptied. This is shown in Figure 12b. Similarly, Figure 12c shows the propagation of change in API concentration at the outlet of respectively the mill and the granule feeder. As a result of this, a change in API concentration profiles at the outlet of granule blender and the feed frame are shown in Figure 12d. Finally, due to the step change in the amount of API 2 in the feed components an eventual deviation in the potency of tablets from 0.364 g to 0.374 g is shown in Figure 13. Thus, the case study demonstrates the use of the flowsheet model developed to track the effects of disturbances in upstream units on the final product quality.

3.2.2. Case Study 2: Step Change in Dryer Air Temperature

In case study 2, a step change from 40 ${}^{\circ }\mathrm{C}$ to 50 ${}^{\circ }\mathrm{C}$ is made to the dryer air temperature. The step change is made at 455 s as shown in Figure 14 (right axis). At 455 s, filling, drying, and emptying of the first cell are completed. In the second dryer cell, filling is completed and drying is in progress. In the third dryer cell, filling, and drying are in progress. Profiles of dried granule moisture content from the simulation explained in Section 3.1, along with this case study, are also plotted in Figure 14.

It can be observed that the dried granule moisture content from the first dryer cell is same in both cases as the step change occurs after emptying the first cell. In the second cell, dried granule moisture content is lower in the case where the step change is imposed. This is because granules in this cell are exposed to a higher air temperature for a duration of 175 s (180 × 1 + 450 − 455), which leads to a lower granule moisture content. In the third cell, dried granule moisture content is further lowered as drying at 50 ${}^{\circ }\mathrm{C}$ occurs for a longer duration of 355 s (180 × 2 + 450 − 455). In the fourth cell, the moisture content is further lowered to a steady state value as all the granules are dried at 50 ${}^{\circ }\mathrm{C}$.

Figure 15 shows propagation of LOD for both cases in the feed frame, which is the result of profile changes in the dryer and subsequent mill, feeder and blender units. It can be seen that the profiles follow the same path until about 750 s, after which a lower steady state value is reached for the case where step change to a higher air temperature occurs. The effect of difference in LOD profiles is also reflected in the tablet hardness as shown in Figure 16. Thus, implementation of the flowsheet model allows analysis of effects of such dynamic changes made to an upstream process variable on the final product quality of interest.

3.3. Scenario Analysis Results

The adapted flowsheet model as shown in Figure 4 was used to implement a scenario analysis as explained in Section 2.2. Scenario analysis entails running the flowsheet model at various process setting values. It serves as a useful step before running a more computationally demanding sensitivity analysis. The flowsheet model is developed and errors are typically debugged at fixed values of process settings. However, this masks errors that may occur at other process setting values and does not provide confidence that the flowsheet model can run seamlessly. Debugging model errors at this stage before running more expensive analyses such as variance based sensitivity analysis serves as an effective modeling practice. Total flow rate, LS ratio, granulator screw speed, dryer air temperature, and drying time as tabulated in Table 4 are the input factors considered for scenario analysis. A total of 768 simulations from three levels of total flow rate, four levels each for LS ratio, granulator screw speed, dryer air temperature, and drying time ($3\times 4\times 4\times 4\times 4$ = 768) were successfully run. Process responses from blender, granulator, mill, dryer, and tablet press models as given Table 4 were recorded at the end of each simulation. For brevity, only few process responses are discussed in this section.

Specifically, wet granule d50, dry granule d50, dry granule LOD, and tablet hardness are plotted and discussed. Since, this is a multivariate analysis (five variables) plotting and analyzing responses from simultaneous change in all the variables is not possible. Hence, a matrix of plots as shown in Figure 17 are used. The matrix consists of 10 plots, each of which shows process response plots from varying two distinct factors with the three other factors fixed at baseline values. Figure 17a can be used to visualize and understand the effect of the five variables on wet granule d50. It can be observed that wet granule d50 increases with LS ratio and screw speed, and decreases with flow rate. This is in accordance with the experimental data used for granulator model development [13]. Figure 17b,c can be used to understand effects on dry granule size and moisture content respectively. Figure 17c shows that moisture content decreases with drying time which is an expected phenomenon. Figure 17b shows that dry granule size increases with LS ratio, screw speed, drying time, and air temperature, and decreases with flow rate. This is in accordance with experiments used for dryer model development [15]. Drying is expected to increase granule strength and lower breakage rate, which leads to a larger size. Hence, increase in drying time and air temperature increases granule size. The effect of flow rate, LS ratio, and screw speed on the dry granule size is due to the propagation of effects of these variables on the size of the wet granule feed to the dryer. Similarly, the effect of changes in the variables on tablet hardness is plotted in Figure 17d. It can be observed that low drying time leads to tablets with very low hardness. In other words, granules with high moisture content cannot be used to make tablets. This is in accordance with process knowledge as gained from experiments. It is also worth noting that development of a flowsheet model has enabled study of the effect of change in a process variable in an upstream unit on the final product quality.

3.4. Sensitivity Analysis Results

Sensitivity analysis was conducted using the adapted flowsheet model as explained in Section 2.2 and using five input factors, flow rate, LS ratio, granulator screw speed, dryer air temperature, and drying time. The effect of these input factors on 20 model responses are studied. The list of factors and responses, ranges for the factors are tabulated in Table 4. The input factors are considered to vary uniformly within these ranges. Morris analysis is implemented using 120 samples (= 20 × (5 + 1)). The results of the analysis are shown in Appendix A.1 that lists the metrics $\mu$, $\mu \ast$, and $\sigma$. From the metrics, we observe that LS ratio influences wet granule moisture content, which is expected. Wet granule size is influenced by all three granulator process variables. This finding conforms with experiments that show effect of these variables on granulation rate [12]. For the dryer model outputs, flow rate and granulator screw speed do not show influence on dry granule moisture content and all five input factors show influence on dry granule size and dry granule properties (bulk density, tapped density, and angle of repose). For the mill model outputs, LS ratio shows the most influence on milled granule size, bulk, and tapped density. This is due to the effect of LS ratio on the size distribution of feed to the mill. The tablet press variables, main compression height and fill depth, are shown to be most influenced by air temperature, drying time and LS ratio. This is due to the effect of all these factors on the granule moisture content as granules with LOD higher than 3% cannot be used to make tablets. Thus, these factors also show an effect on tablet hardness.

From Morris analysis, dryer air temperature, drying time, and LS ratio are identified as the significant factors that influence tablet properties. Variance based analysis as explained in Section 2.3.2 can be applied to this subset of factors to obtain a quantitative understanding of their influence. However, for this work, variance-based sensitivity analysis is performed using all five factors and compared to results obtained from Morris method. The analysis was implemented using 3500 (= 500 × (5 + 2)) samples. First order sensitivity indices ($S_i$), as well as total sensitivity indices ($S_{Ti}$), were computed to quantify effects of the five input factors on 20 output responses. For brevity, the indices are tabulated in Appendix A.2. Here, only the total sensitivity indices, $S_{Ti}$ are pictorially represented in Figure 18 as the first order effects $S_i$ are close to $S_{Ti}$ for responses from granulator, dryer, mill, and blender units. The responses from tablet press show interactions for the factors air temperature, drying time, and LS ratio. Generally speaking, findings from variance based sensitivity analysis agree with the findings from Morris analysis. The variance based method identified LS ratio as a significant factor for granulator. All five factors were identified as significant for dryer. LS ratio showed the greatest influence for the mill. All of these findings align with conclusions obtained from implementing Morris analysis. For the tablet press, air temperature, drying time, and LS ratio were identified as significant factors that influence tablet hardness, main compression height, and fill depth. In addition, all five factors were identified as significant for tablet potency. The influence of these factors on tablet potency was not observed from Morris analysis. This is because, in Morris analysis, the metrics are not dimensionless. For example, the metrics have a unit of g for tablet potency. Hence, any factor that shows an effect less than 0.01 g was identified as insignificant for Morris analysis. However, in variance based analysis, the metrics are dimensionless and the effect of various factors on the responses are normalized. A potency difference in the order of 1 × 10${}^{-6}$ g is also accurately identified in the sensitivity indices. Overall, both Morris and Variance-based methods serve in identifying critical factors. While, the Morris method requires fewer samples and allows ranking the factors by the order of influence, it does not provide quantitative information on relative effects of the factors. On the other hand, variance based analysis requires much higher number of samples but can provide detailed and quantitative information on the relative effects of the factors.

4. Conclusions and Future Direction

In this work, a flowsheet model that approximates the ConsiGma${}^{\mathrm{TM}}$-25 line for continuous tablet manufacturing through the wet granulation route is developed. The flowsheet model is based on models that are developed from experimental runs on units included in the continuous line using the same formulation and materials. For a complete virtual representation of the continuous line, models for intermediate units (powder blend feeder, granule feeder, and granule blender), as well as transfer lines, are also included in the flowsheet model. The developed model successfully demonstrates its ability to simulate the effect of changes in the process variables through case studies where step changes in API flow rate and dryer air temperature are implemented, and their effect on final tablet properties is understood. The robustness of the developed model is established by systematically running the flowsheet model at several combinations of process settings and analyzing the corresponding process responses. The model is also used to identify CPPs that affect intermediate and final product critical quality attributes (CQAs).

Throughout this article, several applications and capabilities of the developed flowsheet model have already been alluded to in the discussion. However, it is also worth noting some gaps in the developed model, which helps throw light on areas where future research efforts can be focused. For instance, the developed flowsheet model is computationally expensive. A simulation of about 1600 s takes approximately 4 h. While a simplified model is adapted in this work to implement steady state sensitivity analysis, it is not feasible to run dynamic sensitivity analysis using this model. Other areas of improvement include further development of the unit operation models. The blender model currently used predicts only 10%, 50%, and 90% percentile diameters for the powder blend. However, a much higher resolution PSD is required in the TSWG model. In addition, the TSWG model predicts only steady state PSD output based on process setting values. For the dryer model used in this work, drying behavior is based on input material properties at the start of the drying cycle of the cell. This could be further improved by incorporating dynamic modeling of drying behavior. For the mill model, the model parameters used currently are not a function of drying time in the fluid bed dryer. In addition, some of the submodels for unit operations use empirical relationships which are formulation specific. Another valuable verification of the model would be to check the mass balance over the entire system. The approximation of the axial dispersion models modeling the material flow propagation could thus be achieved. As research efforts continue on improving the unit operation models, the flowsheet model in its current state has already shown to be a useful tool for enhancing process understanding and enabling better decision making.

Overall, the developed flowsheet model is a prerequisite for identification of design space and optimization of the continuous line. Future research efforts should be focused on reducing computational expense of the model, as well as improving the capability of the unit models to capture dynamics and their applicability for other formulations suited for continuous solid oral dosage manufacturing.