Next Article in Journal / Special Issue
Dynamic Modeling for the Design and Cyclic Operation of an Atomic Layer Deposition (ALD) Reactor
Previous Article in Journal / Special Issue
Scale-Up Design Analysis and Modelling of Cobalt Oxide Silica Membrane Module for Hydrogen Processing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Modeling of Particulate Processes for the Continuous Manufacture of Solid-Based Pharmaceutical Dosage Forms

Department of Chemical and Biochemical Engineering, Rutgers University, Piscataway, NJ 08854, USA
*
Author to whom correspondence should be addressed.
Submission received: 31 May 2013 / Revised: 11 July 2013 / Accepted: 26 July 2013 / Published: 19 August 2013
(This article belongs to the Special Issue Feature Papers)

Abstract

:
The objective of this work is to present a review of computational tools and models for pharmaceutical processes, specifically those for the continuous manufacture of solid dosage forms. Relevant mathematical methods and simulation techniques are discussed, as is the development of process models for solids-handling unit operations. Continuous processing is of particular interest in the current study because it has the potential to improve the efficiency and robustness of pharmaceutical manufacturing processes.

1. Introduction

In recent years, the pharmaceutical industry has experienced significant changes in the prevailing economic and regulatory environments. Increased global competition, particularly from manufacturers of generic products, has resulted in decreasing competition-free lifespan of products and reduced profit margins as drugs come off-patent [1,2]. Meanwhile regulatory agencies worldwide, including the US Food and Drug Administration (FDA) and the European Medicines Agency (EMA), have begun to adopt the quality by design (QbD) paradigm introduced by the ICH Q8 guidance on pharmaceutical development. QbD requires that companies demonstrate understanding of the way in which variability in raw materials as well as process design and operating conditions affect product quality and use this understanding to implement effective quality control strategies [3,4,5].
Despite the challenges faced by the industry, pharmaceutical manufacturing processes remain relatively inefficient and poorly understood as compared with other those in other chemical process industries [4]. This can be attributed in part to the unique challenges associated with pharmaceutical process design. Throughout the drug development process, material is required for clinical trials. The corresponding timeline for delivery of clinical materials may limit the available resources for process development [6]. In addition, each active pharmaceutical ingredient (API) has its own set of physical and chemical properties which can affect the success of various drug product formulations and manufacturing routes [4]. Perhaps in part as a result of the aforementioned challenges, sequential scale-up of batch processes remains the predominant process development trajectory within the pharmaceutical industry [6]. Unfortunately this does not necessarily result in the most efficient or robust manufacturing processes. Under the current process development paradigm, manufacturing costs consume a large portion of revenue for many pharmaceutical companies, as much as 27% by some estimates [7]. In many cases more is spent on manufacturing than on research and development [6,7]. Insufficient process understanding and lack of robust process development can also result in variability in product quality [8]. In order meet the challenges associated with current economic and regulatory realities, the pharmaceutical industry will need to invest in efficient and reliable manufacturing technologies [8,9]. Continuous processing has a great deal of potential to address issues of cost and robustness in the development of pharmaceutical manufacturing processes. From an economic perspective, continuous processes tend to involve smaller equipment than batch processes. This corresponds to decreased capital investment in equipment and plant space as well as decreased utility requirements. Continuous processes scale readily through increases in operating time, total flow rate or parallelization, reducing the need for scale-up studies throughout the development process. This can reduce time to market, which in turn may increase competition-free lifespan. It can also decrease the amount of potentially expensive API required for process development, as continuous processes can generate a large quantity of data relatively quickly [2,8,10,11]. Finally, continuous processing can mitigate issues of product variability through implementation of on-line process control [8,9,10].
Process systems engineering tools have the potential to play a significant role in the transition from batch to continuous processing in the pharmaceutical industry. Predictive process models can be used as a supplement to experiments throughout process development, enhancing understanding of process variability and contributing to design space exploration [12,13]. Modeling and computational tools such as flowsheet simulations and global sensitivity analysis can also contribute to identification of critical process parameters to support quality risk assessment (QRA) [14,15,16]. Predictive process models can be used to develop and assess control strategies for continuous processes and set the stage for the implementation of model predictive control (MPC) for pharmaceutical processes. These advanced control strategies can help to ensure consistent product quality [10,17,18,19]. Process models can also be used to optimize manufacturing processes and suggest optimal design and operating conditions [20,21].
The objective of this review is to describe the application of process modeling tools to the study of continuous manufacturing processes for pharmaceutical solid dosage forms, specifically tablets. Three common manufacturing routes, direct compaction, wet granulation and dry granulation, are considered. Available unit operation models and relevant computational methods for the study of these processes are discussed. References are provided to specific applications of these tools for modeling and simulation of solids-based pharmaceutical processes. The remainder of this paper is organized as follows. Section 2 provides an overview of continuous tablet manufacturing processes and the unit operations involved. Processing equipment and key design and operating parameters for each unit operation are described. In Section 3 various computational tools for process modeling and simulation are discussed. Section 4 provides an overview of existing process models for the unit operations involved in continuous tablet manufacture. An emphasis is placed on equation-oriented modeling approaches that can be used for process simulation. In Section 5 several techniques for model validation and verification are discussed. Section 6 summarizes the state of process modeling and identifies areas for future work.

2. Continuous Tablet Manufacturing

2.1. Process Overview

Tablets are among the most common oral solid dosage forms for drugs [22]. The exact manufacturing process for solid dosage forms vary from compound to compound, as the properties of the active pharmaceutical ingredient (API) molecule play a significant role in the development of an appropriate formulation. This work will describe a general process for transforming raw materials (API and excipients) into tablets and provide a summary of several manufacturing routes available for tablet production.
Downstream pharmaceutical processes typically begin with the feeding of raw materials to the process. These materials include active pharmaceutical ingredients, excipients such as microcrystalline cellulose or lactose, and lubricants like magnesium stearate. In most cases the concentration of excipient is significantly greater than that of API or lubricant. The material may then be passed through a comill to eliminate any large, soft lumps within the powder. Thereafter the material is mixed (blended) to ensure uniform distribution of active ingredient. Optionally the blend may be granulated via wet or dry granulation. The use of wet granulation necessitates a granule drying step prior to further processing. If a granulation process has been implemented, a milling step is typically required to reduce the granule size to the desired level before tableting. In the absence of a granulation step material may be sent directly to the tablet press after blending. In a continuous tablet press the powder is typically fed via a hopper and a rotary feed frame. The powder blend fills a die and is subsequently compressed to create a tablet. An example flow sheet for a continuous tableting process showing several design alternatives (e.g. wet granulation, dry granulation, direct compaction) is shown in Figure 1. Certain units including feeders, hoppers and tablet presses are common to all manufacturing routes. Mills are present in both wet and dry granulation. Roller compactors are unique to dry granulation while wet granulation equipment such as twin screw extruders and granule dryers are present only in wet granulation based processes. A summary of relevant processing equipment along with design and operating parameters is provided in Section 2.2 and in Table 1.
Figure 1. Three manufacturing routes for the continuous production of pharmaceutical tablets are shown. Common processing steps for all three routes include feeding, blending and tableting. The dry granulation route involves roller compaction followed by milling of the produced ribbons while the wet granulation route involves wet granulation followed by drying and then by granule milling. In the case of direct compaction there is no granulation step and therefore milling is not required prior to tableting.
Figure 1. Three manufacturing routes for the continuous production of pharmaceutical tablets are shown. Common processing steps for all three routes include feeding, blending and tableting. The dry granulation route involves roller compaction followed by milling of the produced ribbons while the wet granulation route involves wet granulation followed by drying and then by granule milling. In the case of direct compaction there is no granulation step and therefore milling is not required prior to tableting.
Processes 01 00067 g001
Table 1. Processing equipment for continuous tablet manufacturing including the adjustable design and operating parameters.
Table 1. Processing equipment for continuous tablet manufacturing including the adjustable design and operating parameters.
Unit nameSymbolDesign parametersOperating parameters
Hopper Processes 01 00067 i001Shape (conical, wedge)
Width
Outlet diameter
Wall angle
Material of construction
Powder flow rate
LIW Feeder Processes 01 00067 i002Tooling (screw, screen)
Hopper size
Operating mode
Screw speed
Flow rate set point
Continuous Mixer Processes 01 00067 i003Vessel length and diameter
Agitator size and configuration
Agitator rpm
Mixer fill level
Twin Screw Extruder Processes 01 00067 i004Number of screws
Screw geometry
Barrel length
Binding solution properties and addition location
Screw speed
Granulation temperature
Liquid to solid ratio
Powder Flow Rate
Binder content
Roller Compactor Processes 01 00067 i005Roll configuration
Roll diameter
Roll surface
Powder feed
Powder feed rate
Roll speed
Compaction pressure
Roll gap
Mill Processes 01 00067 i006Mill type
Mill configuration
Geometry
Screen/selector size
Equipment size
Air nozzle arrangement
Solids feed rate
Rotor speed
Grinding pressure
Tablet Press Processes 01 00067 i007Die and punch size and geometry
Die feeding method
Number of compression stations
Die filling method
Lubrication method
Powder feed rate
Compression force
Tableting speed
The success of a tableting process is assessed based on the quality attributes of the tablets produced. Properties of interest include tablet strength (hardness, friability), tablet API composition and relative standard deviation, tablet weight, weight variability and tablet dissolution.

2.2. Processing Equipment

2.2.1. Hoppers

Hoppers are commonly used in solids processing as a means of holding materials and conveying them gravimetrically. In pharmaceutical tablet production, hoppers are typically found in conjunction with other operations like feeders and tablet presses. In a well-design hopper, powder should flow consistently towards the outlet of the unit at an approximately constant flow rate. This is described as the “mass flow” regime. Hoppers come in a variety of geometries including conical and wedge shaped, both of which are used in pharmaceutical applications. In addition to geometry, hoppers can vary in width, outlet size and wall angle as well as material of construction. The hold-up in the hopper is an important operating parameter that should be monitored to avoid overfilling or, in the case of feed hoppers, running empty while the process is still operational [23,24]. Segregation and poor flow are potential concerns in the operation of industrial hoppers. Flow problems include arching and bridging, which can stop flow from the hopper, and ratholing—the stagnation of material at the hopper walls which effectively reduces capacity. These issues can be addressed via changes in hopper design or modification of the flow path using inserts [23,25,26]. Segregation is a significant concern in pharmaceutical applications, as it can result in variability in composition over time. In a hopper attached to a tablet press, segregation could affect content uniformity. Powder properties play a significant role in segregation behavior. Segregation is more likely to occur in polydisperse materials with wide particle size distributions and a large degree of density variability among the particles [27,28].
Hopper refilling dynamics and their effect on downstream operations are also a concern in continuous processes. The impact of hopper refilling on discharge rate and on the performance of a continuous tableting process has been studied both computationally and experimentally [14,29,30,31]. Flow and discharge behavior of pharmaceutical powders from several different hopper geometries has also been studied experimentally by Faqih et al. [26] who developed a quantitative flow index that was shown to correlate with hopper flow behavior based on data from a gravitational displacement Rheometer.

2.2.2. Loss-in-Weight Feeders

Loss-in-weight feeders, as the name implies, feed powders into a process by mass. Feeding by mass provides the desired accuracy for pharmaceutical processes, where low feed rates may be needed, particularly for APIs and lubricants [30,32]. Loss-in-weight feeders generally consist of a hopper, a load cell that is integrated with a gravimetric controller and a conveying mechanism such as a screw feeder. A secondary conveying mechanism at the base of the hopper may be used to transfer powder to the screws. A discharge screen at the exit of the feeder can be used to break-up any lumps in the powder entering the process. The preferred tooling (screw, discharge screen) for a loss-in-weight is a function of the material properties of the feed powder and the intended flowrate. Engisch and Muzzio [30] have described a characterization method for loss-in-weight feeder equipment that can aid in determining appropriate tooling for a given application. The key operating parameter for loss-in-weight feeders is the screw speed, which dictates the powder flow rate. The hopper fill level can also affect feeder performance. Specifically, feeder accuracy may decrease at low hopper fill levels [33].

2.2.3. Continuous Mixers

The purpose of mixing in downstream pharmaceutical applications is to reduce composition variability in multi component powder blends. In continuous mixing spatial and temporal variations in the composition of the material exiting the mixer are of concern. Therefore a well-designed continuous mixer should be able to produce evenly distributed (not segregated) blends with good control of composition over time.
A summary of available equipment for continuous mixing along with studies used to characterize these mixers is given in Pernenkil and Cooney [34]. Of particular interest in pharmaceutical applications are continuous convective mixers, in which the primary mixing mechanism is convection induced by rotating blades [35]. Key design parameters for these types of mixers include the vessel length and diameter and the agitator size, configuration and geometry [36]. Important operating parameters include mixing angle, agitator rpm and powder flowrate. Powder properties, particularly cohesiveness and flow properties, can also affect mixing performance [35,37,38,39]. As is the case for hoppers, if components in the blend have very different properties (particle size and density) segregation may occur [40]. Several metrics can be used to monitor mixing efficiency. Relative standard deviation (RSD) and variance reduction ratio (VRR) are based on composition variability at the outlet and, in the case of VRR, inlet of the mixer. These metrics are discussed at greater length in Section 4.3, as is the residence time distribution (RTD). These characteristics can be calculated from experimental data such as concentration measured at the outlet of the mixer. Composition can be assessed using near infrared spectroscopy (NIR) [34,35,41,42,43,44]. In addition, flow trajectories within the mixer can be studied using positron emission particle tracking (PEPT) [39].
The effect of operating parameters and powder properties on mixing performance within continuous convective blenders has been studied extensively [35,36,37,38,39,41,45]. Marikh et al. [36,45] have studied the influence of mass flow rate, impeller speed and bulk powder properties on the hold up within a pilot-scale continuous mixer for different agitator types. Based on experimental studies a correlation between agitator speed and the bulk mean residence time was developed which could potentially be used for scale-up. Portillo et al. [39] have studied flow behavior in continuous convective mixers as a function of impeller speed, flow rate and powder cohesiveness. Flow behavior is characterized by particle trajectories, axial dispersion coefficient and residence time. In a separate study Portillo et al. [37] have examined the effect of mixing angle in addition to impeller speed and cohesion on residence time and blend content uniformity. Two mixers that differed according to geometric parameters such as vessel length and diameter and blade size, shape and configuration were compared in this study. Finally Vanarase and Muzzio [38] have studied the impact of impeller rotation rate, flow rate and blade configuration on mixer performance as characterized by RSD and VRR.
Extensive modeling efforts have also been conducted for continuous mixers. These will be discussed in Section 4.3.

2.2.4. Wet Granulation

The purpose of granulation is to create agglomerates, or granules, of powder blends. These granules are of higher bulk density than the bulk powder and thus tend to have better flow properties. The granulation process also contributes to improved control of content uniformity and compactibility [46]. Wet granulation involves combining a dry powder blend with a liquid, known as a binding solution, in the presence of some sort of agitation. When the powder is wetted granule nuclei consisting of binder droplets that have taken up powder begin to form. Thereafter the granule grows and densifies through coalescence and consolidation. Finally breakage may occur as the granules collide at high speed or experience significant shear. The properties of the powder blend, choice of binder, binding solution content relative to the amount of powder, and granulation method all have the potential to affect granule properties [47,48].
Wet granulation can be performed using a variety of techniques and equipment, including high shear mixers, fluid bed granulators and single or twin screw extruders. High shear mixers are primarily used in batch applications although so-called “instant agglomerators”, which are similar to high shear mixers that process very small volumes of material and therefore have an extremely short residence time, can be used for continuous processing [46,49]. Fluid bed granulators, the most common of which are horizontal moving bed granulators, operate by spraying binding solution from above onto a powder bed that is fluidized with air (or an inert gas such as nitrogen) from below. The granules are partially dried by the air flow and granule drying is readily completed within the fluid bed granulators when the addition of binding solution is stopped. Fluid bed granulators tend to operate at relatively high production rates, e.g., 20 kg/h or more [46,50]. Therefore they are not necessarily suited to pharmaceutical processes which may operate at lower throughput, particularly during the formulation development process. Extrusion equipment consists of one (single) or two (twin) screws within a barrel. Often pharmaceutical granulation is carried out in twin-screw extruders with co-rotating screws. Mixing and agglomeration occur as the material is conveyed along the length of the barrel. Binding solution is injected into the process at a certain point or points along the length of the barrel [47,51]. Depending on the degree of densification that occurs within the extruder it may be necessary to mill the granules resulting from an extrusion process down to an appropriate size prior to using them in subsequent processing steps. For all processes except fluid bed granulation, a granule drying step is also required [46,50].
Extrusion processes have been studied most for continuous pharmaceutical processes. This is due in part to the fact that extrusion is readily scaled-up via increased throughput and therefore lends itself to the development environment within the pharmaceutical industry. In these processes, the main design parameters include the length of the barrel, the number of screws (one or two), the geometry of the screws, and the location of a point or points along the barrel where binding solution is added [46,50]. Operating parameters of interest include screw speed, granulation temperature and liquid to solid ratio [51]. Powder properties and binder properties including viscosity and surface tension have the potential to influence the process as well [52]. Performance indicators for wet granulation are related to the granule properties, including the granule density and porosity, mean granule size and granule size distribution and the granule composition with respect to active ingredient [47,52]. In addition, the yield or recovery from a granulation process is important. The loss of material to fines during the granulation process impacts the overall yield and could affect granule composition if the active ingredient forms fines more readily during the granulation process [47,51].
Lee et al. [47] have compared twin screw extrusion (TSE) with high shear mixing (HSM) by measuring the properties of the resulting granules and have found that TSE tends to elicit multimodal granule size distributions with greater porosity. This is consistent with the need for a size reduction step following some extrusion-based granulation processes. Dhenge et al. [52] have studied the effect of binder solution properties including viscosity and surface tension on residence time and torque in a TSE granulation process as well as on granule properties including size distribution and granule strength. Tu et al. [51] have experimentally studied granule properties as a function of operating conditions including screw speed and liquid to solid ratio in two different twin screw extruders. This enabled the authors to develop regime maps for the process, describing the extruder geometry in terms of granulation, extrusion and breakage regimes. Regime maps are useful in the characterization of granulation equipment like twin screw extruders. Several authors have studied wet granulation in the context of other unit operations. Cartwright et al. [53] have addressed the importance of accurate powder feeding for successful granulation performance in a twin screw extruder by comparing two types of loss in weight feeders. The impact of granulation conditions, including throughput, screw speed, screw configuration, angle of kneading elements and barrel temperature, on both granule and tablet properties have been examined by Vercruysse et al. [54]. Finally the continuous drying of wet agglomerates in a fluid-bed dryer has been studied by Palzer [55], who found, among other things, that average residence time in the dryer significantly affects the propensity for undesired secondary agglomeration to occur during the drying process.
Extensive research has also been conducted into the modeling of granulation processes, as will be discussed in Section 4.4.

2.2.5. Roller Compactors

Granulation can also be accomplished via dry granulation. Dry granulation is carried out continuously through roller compaction. Dry granulation differs from wet granulation in that granules are formed only through compression. Thus the powder properties of the raw materials are important in determining whether dry granulation via roller compaction is an acceptable approach. Powders that have low bulk density, small particle size (e.g., micronized materials) or are very cohesive may not perform well in roller compaction [56]. Powders with good flowability and compressibility tend to be better candidates for roller compaction [57,58]. The development of a robust granulation process can be more challenging via dry granulation than wet granulation, but dry granulation offers advantages including the elimination of the granule drying step, shorter processing time and lower capital investment and utility costs. Additionally, roller compaction can be used for moisture-sensitive APIs [50,59,60].
In the roller compaction process, powder is fed into a set of counter rotating rolls, conveyed forward with the motion of the rolls and as it reaches the point where the rolls are closest together undergoes compression and forms a compact ribbon. The ribbon is then conveyed forward and released from the rolls. When the particles are fed to the rolls they are initially considered to be in the slip region of the process, characterized by the particles slipping at the surface of the rolls. Relatively little pressure is exerted on the powder in this region. As the powder is drawn towards the point where the rolls are closest together, the wall velocity of the powder begins to match that of the rolls and the pressure exerted on the powder increases substantially. This is known as the nip region, and it is in this part of the process that powder is compacted. The nip angle, typically denoted as α, defines the angle at which the powder transitions from the slip region to the no-slip (nip) region [61]. Once the roll gap begins to increase again, after the point where the rolls are closest, the ribbon produced in the nip region is said to be in the release region [62]. The ribbons can subsequently be milled to obtain appropriately sized granules [60].
The different equipment available for roller compaction varies mostly in terms of configuration. Rolls can be arranged either horizontally (with ribbons coming out parallel to the floor) vertically (with ribbons coming out perpendicular to the floor) or at some intermediate angle. In addition, it is possible that both rolls will be fixed, resulting in a constant roll gap throughout the process, or that one roll can be adjustable allowing for a change in roll gap during processing. The roll diameter and width vary from one piece of equipment to another. The surface of the rolls can also vary from smooth to knurled (rough/with grooves) to pocketed design. Finally the powder can be fed either gravimetrically or using a screw feeder [50,56]. In addition to the previously discussed powder properties and equipment design parameters, several operating parameters affect the performance of roller compaction processes. These include the feed rate of powder to the rolls, the roll speed and the compaction pressure applied to the powder, which is a function of the roll gap [60,63]. The nip angle, which is related to the length of the compaction zone and therefore to the degree of compression, is an important indicator of the performance of a roller compaction process, as is the peak or maximum pressure [57,62]. The performance of roller compaction processes is assessed according to metrics similar to wet granulation. The propensity for fines can be greater in roller compaction depending on the properties of the powder being formulated [50]. In addition, ribbon density and ribbon density variability can also be used to evaluate the process performance.
Many studies of roller compaction focus on understanding the relationship between design and operating parameters. Bindhumadhavan et al. [62] have studied the influence of roll speed and roll gap on pressure profiles during roller compaction. The influence of rotation angle, feed rate and roll speed on pressure distribution and drive torque applied to the rolls has been experimentally evaluated by Lecompte et al. [58]. In addition to varying the angle at which powders are fed to the rolls Miguélez-Morán et al. [64] have examined the relationship between lubrication and ribbon density, which is an important performance indicator for roller compaction processes. In general lubrication was found to reduce variability in ribbon density. Yu et al. [57] have studied the relationship between powder properties and compaction performance and have shown that powders with better flowability achieve higher peak pressures during roller compaction, resulting in improved compaction behavior.

2.2.6. Milling

Milling is a broad term that describes a wide variety of size reduction techniques. The current work will focus on dry milling techniques that are relevant for downstream pharmaceutical processes. Milling can be used to delump powders, to reduce the particle size of raw materials, particularly APIs, or to reduce the size of granules generated through either wet or dry granulation [65].
Conical screen mills, such as comills, can be used for delumping of powders, or for coarse to fine control of wet or dry granule size. Conical screen mills consist of a cone shaped screen with an impeller inserted into the center. The impeller rotates and material is ground between the impeller and the screen until it is small enough to pass through the holes in the screen and leave the mill. The key design parameters for a conical mill include the screen mesh size, the size of the cone, the impeller shape and the impeller to screen distance, which can be adjusted using spacers. Important operating parameters include powder feed rate and agitator speed [60,66]. The properties of granules entering the process also play a significant role in dictating the degree of size reduction that can be achieved [67]. Granule size can also be reduced using oscillating granulators, which are essentially screens integrated with the roller compaction process. The compact ribbon is forced through a screen using an oscillating rotor, so called because rotor speed and direction can vary in time. Screen size and rotor speed and rotation angle dictate particle size [68].
Comills are often used for granule size reduction but other milling technologies may be used if a smaller final particle size is desired. These include air jet mills and impact mills. For each type of mill a variety of configurations is available. Air jet mills generally consist of a grinding chamber, where particles are broken, and a classification chamber, where particles are separated according to size. Sufficiently small particles are passed to the next unit operation while fines can be collected in a dust filter. For these types of mills the main design parameters of interest include the geometry, number and configuration of the air nozzles. Significant operating parameters include the solid feed rate and the grinding pressure [65]. Impact mills include pin mills and hammer mills. In both types of mills, material is fed through the center of the mill and exits at the outer edge of the chamber. A selector grid can be used to allow only sufficiently small particles to exit the milling chamber. Pin mills grind the product between two disks to achieve size reduction while hammer mills rely on high impact particle wall, particle blade and particle-particle collisions. For pin mills, the mill size is the main design variable. For hammer mills the equipment size, blade configuration and selector grid size are all relevant design variables. Variable operating parameters for pin and hammer mills include solids feed rate and the rotor speed [65,67].
Regardless of the equipment used, the objective of milling is to achieve the desired particle size. The performance of a milling process can thus be assessed in terms of its ability to achieve the desired mean size or size distribution [16]. Specific surface area, which is related to particle size, can also be used to determine milling performance. Yield is also a concern in milling operations, as the production of a large quantity of fines that are recovered in the dust filter could result in significant losses of potentially expensive raw material (e.g., API) [67].
Comills have been studied more extensively for granule size reduction than have impact or air jet mills. Experimental studies have shown that screen size, impeller speed and impeller shape can be varied to affect granule size distribution [66,69]. It has also been demonstrated that decreasing screen size coupled with increased impeller speed results in smaller granule size [70]. Some correlation between granule properties and milling performance has also been demonstrated. Inghelbrecht and Remon [71] have found that low friability corresponded to reduced dust formation during milling. Verheezen et al. [67] have studied impact milling for granule size reduction and found that granule strength has a significant effect on the final particle size achieved. In addition, fines formation was found to be related to the total degree of size reduction.

2.2.7. Tablet Press (with Integrated Hopper and Feed Frame)

The tableting process involves the compaction of powder blends to form a hard compact. This is achieved in a tablet press, which contains several components integrated as a single processing unit. A hopper conveys material into the tablet press. A feed frame is then used to move the powder or granular material into the die, a cavity that defines the tablet size and shape. A punch compresses the material within the die to form a tablet. Cam tracks guide the continuous movement of the dies so that they can be filled, compacted, and discharged [14].
Material properties can significantly affect tableting performance. Wide particle size distributions or variability in density can result in segregation during die filling, causing non-uniform tablet composition. Low bulk density, poor compressibility or flow properties can affect die filling and the pressure profile during compaction, resulting in tablet weight variability and insufficient tablet hardness [72,73]. Experimentally a correlation between the particle size and specific surface area of the material to be compacted and the hardness and dissolution properties of tablets has been demonstrated [74,75]. Compression problems such as capping can occur due to the presence of excessive fines within the particle size distribution [76,77]. Particle or granule moisture content can impact drug product stability, compressibility and tablet physical properties [78]. In addition, the physiochemical properties of the active ingredient could affect the propensity for form change due to temperature increase during the compaction process [79]. Drug and lubricant properties can also affect tableting performance, particularly as drug loading or lubricant content increases [74,80,81,82].
The variable design parameters for a tablet press are related to tooling and method of operation. The tooling includes the die and punch size and geometry, which can be changed according to the desired tablet weight and shape. The number of compression stations also differs between tablet presses. The die can be filled via force feeding or suction filling [83]. Finally lubrication can be internal, mixed with the powder or granule blend, or external, applied directly to the punch and die assembly. Internal lubrication is common in the pharmaceutical industry, though some studies have indicated that external lubrication can mitigate the effect of lubricant on tablet hardness [82]. Operating parameters of interest for a tablet press include the powder feed rate and the compression force applied to the tablets as well as the rate of tablet production. Tableting speed and powder feed rate are related to tablet weight and weight variability while compression force affects tablet properties like hardness and density [22,76]. The performance of tableting processes can be assessed based on tablet active ingredient content, content uniformity, weight variability, and physical properties such as friability, hardness and dissolution performance. Drug content can be measured continuously using near infrared (NIR) spectroscopy [42,78]. Hardness, friability and dissolution must be measured offline, but models can be implemented to predict hardness and dissolution performance based on operating conditions or spectroscopic measurements [73,74,84,85].

3. Computational Tools and Mathematical Modeling Approaches

3.1. DEM Simulation

Particulate systems are of tremendous importance within the pharmaceutical and numerous other industries, yet they remain relatively poorly understood. In order to enhance understanding of macroscopic behavior in solids processes, it is helpful to first understand the particle-particle and particle-environment interactions that give rise to it. Discrete Element and Finite Element Method (DEM/FEM) models are particle level computational tools that can be used to develop this understanding [86]. An extensive review of theoretical developments and applications in discrete particle simulation has been provided by Zhu et al. [87,88]. This review will focus on the application of DEM to pharmaceutically relevant processes.
DEM can be used to understand particle packing, which is important for problems involving powder bed densification and compaction like the manufacture of pharmaceutical tablets. The ability of DEM to accurately model particle packing has been well documented in the literature through studies comparing simulated with experimental behavior [89,90]. Specifically it has been shown that DEM can adequately represent packing density (porosity), coordination number, radial distribution function and the force network within a packed bed [88,91,92]. For instance, Yi et al. [92] have studied systems of multi sized spheres and shown that the porosity and coordination number obtained via DEM agree well with those obtained experimentally.
Particle and particle-fluid flow behavior has also been studied extensively using DEM. Particle flows, and specifically confined flows, are relevant in pharmaceutical operations such as continuous mixing, fluidized bed drying, and powder feeding and conveying. It has been shown that DEM can accurately model confined flows like direct and annular shear, vertical flow and biaxial or triaxial compression [88]. McCarthy et al. [93] have compared solids fraction and particle velocity profiles as well as granular temperature distributions for simulated and experimental studies of a horizontally aligned annular shear cell and shown that as long as the system and particle geometries are well represented the DEM and experimental profiles agree well. In addition, Wu et al. [94] have modeled the flow of powder in confined space in order to study a die filling process and found that appropriately calibrated DEM models can reproduce the powder flow behavior observed experimentally using high speed cameras.
DEM has been used to model several aspects of the tableting process, including die-filling, compression of materials within a die and crushing of particles during compression [95,96]. Mehrotra et al. [72] have used DEM to explore the impact of powder flow properties, specifically cohesiveness, on the die filling and compression process. It was found that as material cohesion increases the time it takes to fill the die also increases, which is consistent with the experimental observation that the tablet weight variability decreases at lower tableting speeds for cohesive powders. Several authors have also coupled DEM with computational fluid dynamics (CFD) to investigate suction filling of dies and to evaluate the effect of air or lack thereof on segregation during die filling [97,98,99,100]. One significant challenge associated with the use of DEM to model compression processes is the issue of deformation at high relative densities. In many computational studies of compression and compaction processes, FEM has been combined with DEM to improve the representation of particle deformation during compaction [76,95,96,101,102,103]. Specifically Frenning [102] has shown that finite and discrete element methods may be combined to simulate the behavior of granules in a densely packed bed and to understand the relationship between individual granules and the behavior of the granule bed under compression. Gethin et al. [101] have demonstrated that a combined DEM/FEM method can be used to determine particle packing in a die as a function of particle shape.
DEM has also been used to study powder mixing processes. Several authors have demonstrated qualitative agreement between DEM simulations and experimentally obtained mixing behavior [104]. For instance, Marigo et al. [105] have obtained qualitative agreement between experimentally obtained axial and radial dispersion trends and those obtained through DEM simulation for a tubular mixer. Remy et al. [106,107] have compared experimental and DEM results for granular flows in a bladed mixer. This work has demonstrated qualitative agreement between the segregation profiles obtained via DEM and experimental studies for polydisperse, cohesionless spheres and has also shown that DEM can reproduce observed surface velocities, granular temperature profiles and mixing kinetics at varying surface and wall roughness. Of particular interest for continuous pharmaceutical manufacturing applications is the use of DEM to study continuous convective mixers. DEM has been used to validate the periodic section approach to modeling continuous convective mixers through demonstrating that the velocity profiles obtained from periodic section simulation are in agreement with those obtained from full blender simulation [108]. DEM simulations have been used to calculate the residence time distribution (RTD) in various blending equipment, including continuous convective mixers [41,109]. Data from DEM has also been used to inform reduced-order mixing models for process simulation [110].
Due to the importance of hoppers in a variety of industries, including agriculture, food processing and pharmaceuticals, numerous detailed studies of hopper flow have been conducted using DEM as well. DEM has been used to study the impact of powder properties, hopper geometry and operating parameters on segregation and flow patterns in hoppers [23,24,25,27,28,29,84,88,95].

3.2. Population Balance Models

The population balance equation can be used to describe the development of a set of properties of interest for a group of entities over time. Population balance equations have been used to model particulate processes including crystallization, mixing, milling, granulation, drying and dissolution, all of which are of great interest in the manufacture of solid dosage forms for pharmaceutical products [66,111,112,113,114,115] The general form of the population balance equation in 2 dimensions is given below.
Processes 01 00067 i008
F(x, z, t) is referred to as the population distribution function, and describes the state of particles with respect to internal and external coordinates in time. x reflects the internal coordinate(s) such as particle size, mass or volume. z reflects the external coordinate(s) such as axial or radial position and t indicates time. The right hand side of the population balance equation (PBE) indicates the rate of formation and depletion of particles, which can occur through a variety of mechanisms including nucleation, aggregation and breakage. Analytical expressions for these mechanisms can be used to express the formation and depletion terms on the right hand side of the population balance equation [116,117].
Population balance models are often discretized with respect to the internal and external coordinates. In this case the differential terms with respect to x and z in Equation (1) are replaced by finite differences. Discretized population balance models can also be parameterized, as show in Equation (2).
Processes 01 00067 i058
y(t,p) are the discretized states and p are the model parameters, which can be estimated based on experimental or simulated data. The parameter estimation can be formulated as an optimization problem with a least squares objective as outlined in Ramachandran et al [118]. Sen et al. [119] have demonstrated the use of a parameterized population balance equation to model a continuous mixer. The parameters were fit by minimizing error between measured and predicted relative standard deviation (RSD) of the product composition and the API mass fraction at the mixer outlet.
High (three or four) dimensional population balance equations, such as those often encountered in modeling wet granulation processes, can be computationally expensive to evaluate [113]. The implementation of hierarchal solution techniques has been discussed extensively in the literature as a means of reducing the computational time associated with solving these equations [120,121,122,123]. In order to implement a hierarchal solution strategy the multidimensional PBE is discretized with respect to both internal and external coordinates. The partial differential equations in the PBE can then be rewritten as ordinary differential equations which can be integrated over time using a first order Euler predictor/corrector method [122].
The aforementioned methods are summarized in Table 2 with respect to their applications, the level of detail with which they model the process and the relative computational cost.
Table 2. Comparison of modeling techniques discussed in Section 3 of the current work.
Table 2. Comparison of modeling techniques discussed in Section 3 of the current work.
MethodDescriptionPharmaceutically relevant applicationsLevel of detailComputational expense
DEMparticle level simulation of powder behaviorpowder flow, powder mixing, and compactionParticle level informationhigh
PBMdescribes the evolution of populations of entities (particles, granules, droplets) over timemixing, crystallization, granulation, millingDescription of population of particlesmoderate to high—depending on problem dimensionality
ROMapproximation of high fidelity models using a variety of estimation and interpolation techniquesvarious unit operations, simulation-based optimizationUnit operation level descriptionlow

3.3. Reduced Order Models

Reduced order models (ROM) are a class of models that represent high fidelity or full scale models in a lower dimensional space. The order reduction can be achieved through multivariate analysis techniques as in principal component analysis (PCA) or proper orthogonal decomposition. [124,125,126,127,128]. Alternatively a computationally expensive model can be replaced by lower dimensional surrogate model obtained through fitting of experimental or simulated data using techniques such as kriging, response surface methodology (RSM), artificial neural networks (ANN) or high dimensional model representation (HDMR) [12,129,130,131,132,133,134]. The motivation for using ROMs is that they are less computationally expensive than the original models and are therefore suitable for process simulation and optimization purposes.

3.3.1. Kriging

Kriging is a black-box interpolating technique that can be used to generate metamodels or response surfaces from input-output data for a process [135,136]. Originally developed for the purpose of predicting mineral deposit distributions, in recent years Kriging has been increasingly used for modeling in a variety of fields [136,137]. Its popularity can be attributed in part to its ability to model complex nonlinear and dynamic processes [138]. In addition, Kriging generates error estimates associated with each predicted point that can be used to assess prediction accuracy and direct future sampling [12,139,140]. This combined with the fact that Kriging models do not rely on a pre-defined closed form allows for the development of accurate model representations from relatively sparse datasets as compared with those obtained through traditional experimental designs [131,139]. The primary disadvantage associated with Kriging is that it does not provide a simple closed-form expression for the relationship between inputs and responses the way that other response surface methods do (see Section 3.3.5) [131].
In Kriging models, the predicted response (xk) associated with a new input point xk is determined as a weighted sum of the known function values f (xi) associated with previously sampled inputs xi as shown in Equation (3). The weight, wi, attributed to each f (xi) decreases with increasing Euclidian distance between the points xk and xi. Thus Kriging is considered an inverse distance weighting method. In general some maximum distance r is defined such that only points within some distance r of xk are considered [131,136,140].
Processes 01 00067 i010
The weights wi in Equation (3) are unknown and their determination is a critical step in the development of the Kriging model. The objective is to select the weights such that the mean square error of prediction is minimized. In practice these weights can be obtained using a fitted variogram model based on N points that are sufficiently close to xk. The weights must sum to unity, a constraint which arises in part from the condition that the Kriging predictor should be unbiased [138,140]. In addition, if input points are closely clustered together they are given lower weight in order to prevent biased estimation [12,131]. For purposes of the subsequent discussion the Euclidian distance between points (h) and the variogram corresponding to a dataset x consisting of NT sample points γ(h) will be define as in Equation (4).
Processes 01 00067 i011
As can be seen in Equation (4) the variogram is calculated for individual pairs of points. This type of variogram is sometimes referred to as semi-variance while the variogram is described as (h), but for the purposes of the subsequent discussion γ(h) will be referred to as the variogram [138]. For a dataset containing NT sample points there are a total of NT(NT − 1)/2 Euclidian distances and corresponding γ values to be determined. The resulting γ vs. h data is usually smoothed and is then fitted to one of five basic models; sphereical, gaussian, exponential, linear or power-law. These are summarized in Jia et al. [131] and Boukouvala and Ierapetritou [139]. The form of the variogram model is chosen such that it minimizes prediction error, though computational efficiency may also be considered in making the selection. If necessary combinations of the various model types may be used to obtain appropriately low error. [131,138] Once an expression for the variogram has been obtained, it can be used to obtain a complementary function known as the covariance function shown in Equation (5).
Processes 01 00067 i012
The term σ2max corresponds to the maximum variance of the variogram function and is also referred to as the sill. The kriging weights wi for a given test point xk can then be obtained from the covariance function by solving the system of Equation (6).
Processes 01 00067 i013
di,j represents the distance between two points xi and xj while di,k represents the distance between a point xi and the test point xk. Likewise Cov(di,j) represents the covariance between data corresponding to input vectors that are distance di,j apart, obtained from Equation (5). λ are the Lagrange multipliers associated with the constraint that the weights must sum to unity. The variance associated with the test point xk is calculated from Equation (7), in which the term Cov(di,k) is the right hand side of Equation (6).
Processes 01 00067 i014
Figure 2. Algorithm for Kriging with the option for updating the sampling space in order to achieve sufficiently low prediction variance.
Figure 2. Algorithm for Kriging with the option for updating the sampling space in order to achieve sufficiently low prediction variance.
Processes 01 00067 g002
Briefly, the Kriging algorithm, depicted in Figure 2 can be described as follows:
  • Select an initial sample set x consisting of NT sample points and evaluate the process or model at these points to obtain the corresponding function evaluations f(x).
  • Using the data obtained in step 1, calculate the Euclidian distances h and the corresponding semi-variances γ(h) using Equation (4) for all NT(NT − 1)/2 sampling pairs.
  • Smooth the γ(h) vs. h data and fit it to an appropriate variogram model according to a least squares error minimization criterion and/or a secondary criterion for computational efficiency.
  • Based on the variogram model, determine the covariance function as in Equation (5).
  • For a test point xk calculate the weights from Equation (6). Calculate the predicted response (xk) from Equation (3) and the associated variance from Equation (7).
  • Optional—If the predicted variance is larger than desired, collect additional sample points in the region of the test point xk and add those to the set NT to develop an updated Kriging model.
Kriging has been used in a variety of applications related to pharmaceutical process modeling. Jia et al. [131] have demonstrated Kriging models to represent the flow variability in a loss-in-weight feeder as a function of flowability and feed rate. Boukouvala et al. [129] have proposed an extension of this model that includes a Kriging-based imputation of missing data. Boukouvala et al. [138] have demonstrated the use dynamic Kriging to study the roll gap and ribbon density in a roller compaction process with respect to variable feed speed, roll speed and hydraulic pressure. Boukouvala et al. [138] have also demonstrated the use of Kriging to map the design space for a continuous convective mixer and a loss-in-weight feeder and Boukouvala and Ierapetritou [139] have described the use of a Kriging-based method for feasibility analysis of a roller compaction process.

3.3.2. Response Surface Methodology (RSM)

Response surface methodology describes a set of computational tools that can be used to establish relationships between a system or process response and multiple input variables. Box and Wilson first proposed RSM in 1951 as a means of optimizing operating conditions for a chemical processing [130]. In recent years RSM has been increasingly used for pharmaceutical applications [12,129,131,134,141,142]. The goal of RSM is to develop a functional representation for the response as a function of input variables. The functional representation, or response surface, is an approximation because the true form of the variable-response relationship is not known explicitly. The response surface is generated through a local sampling and optimization process within a region of interest. This region may be defined as the full operating space of the process or a subset of this space near an optimum or other point of interest [133].
Figure 3. Algorithm for response surface modeling with optimal model development around a region of interest.
Figure 3. Algorithm for response surface modeling with optimal model development around a region of interest.
Processes 01 00067 g003
The main steps in implementing RSM are given below for a system with a single response, y, and m input variables x1, x2,…, xm. The algorithm for developing a response surface model is also depicted in Figure 3.
  • Establishment of an experimental design
    A design, D, consisting of n samples can be generated using a design of experiments (DOE) or other appropriate statistical approach. Processes 01 00067 i015.
    The proper selection of D is critical to ensure that the generated response surface will be an accurate predictor for the response of interest. The number of sampling points should be greater than the number of coefficients to be fitted for the response surface model. For noisy data, the number of sampling points needed may be greater. Further discussion of appropriate designs is provided in the literature [131,132,133].
  • Development of a response surface model in the region of interest
    The initial response surface model is developed around the nominal sampling point. The form of the model is defined by the modeler. Typically second-order polynomial functions as in Equation (8) are selected for the response functions. Justification for the selection of second order polynomials is provided in the literature [132,133]. The sampled data can then be regressed to the specified model using least squares or other appropriate fitting techniques.
    Processes 01 00067 i016
    ŷ is the estimated response. β0, βi, βij and βii are the model coefficients. xi and xj are the input variables.
  • Local model optimization
    Model optimization is performed in order to determine the region where expected process improvement can be maximized. The optimization can be completed using a steepest descent search over the sampling region. In this case the local optimum is found iteratively. An initial model is built based on the first sampling point. The optimum of this model then becomes the nominal point for the next iteration and a new response surface is built and optimized, with the addition of new sample points. As the algorithm converges, the nominal and optimal points become one and the same [131]. Other optimization techniques such as ridge analysis can also be applied [132]. If different process designs are to be considered, binary variables can be introduced to indicate the design configuration. This results in a mixed integer nonlinear program (MINLP) optimization problem formulation [12].
    RSM has been used to model pharmaceutically relevant unit operations, including loss-in-weight feeding [129,131] and granulation [134,142]. It has also been used to explore design space and aid in the identification of critical process parameters for pharmaceutical applications [12,141].

3.3.3. High Dimensional Model Representation (HDMR)

High dimensional model representation (HDMR) describes a set of techniques that can be efficiently used to describe input-response relationships in high dimensional systems. HDMR expresses the system output as a finite hierarchal correlated function expansion of the various inputs. This functional representation accounts for the effect of individual inputs and the effect of interaction between inputs on a particular output [143,144,145].
Processes 01 00067 i017
in Equation (9) f0 is a zero order term which indicates the average value of the output over the domain of x in random sampling HDMR or the value of x at a specified reference point in cut-HDMR. fi (xi) indicates the individual contribution of variable i to the output while fij(xi,xj) indicates the contribution from the interaction of variables i and j.
If Equation (9) is expanded to include all potential parameter interactions the number of component functions to be calculated can become intractable as the number of variables increases. However in most applications it is sufficient to include terms up to second order interactions only. Most higher order terms represent only a small contribution to the overall response [143,145]. Thus Equation (9) can be approximated as a second order expansion.
Processes 01 00067 i018
The form of each of the component functions in Equation (10) can be optimally selected in order to accurately represent the available data. Objective functions and algorithms for selecting optimal component functions are discussed extensively in the literature [143,145,146].
Component functions can be determined differently depending on whether or not the data is randomly sampled. If it is, Random sampling HDMR (RS-HDMR) can be used. The HDMR expansion is determined from the average value of f(x) over the entire domain of the input space. Cut-HDMR can be used when ordered sampling is used for data collection. In this case f(x) is defined relative to a specified reference point x which is used to determine the component functions. Multiple cut-HDMR is an extension of cut-HDMR in which the expansions are determined relative to several different reference points and then combined to represent f(x). Cut-HDMR generally requires fewer samples than RS-HDMR, as the latter relies on the evaluation of high dimensional integrals to obtain the constant terms [12,147]. However the use of analytical basis functions, including orthonormal polynomials or spline functions can reduce the computational expense of fitting RS-HDMR component functions. In terms of basis functions Equation (10) can be written as:
Processes 01 00067 i019
where αr and Processes 01 00067 i059 are coefficients for the basis functions φr(xi) and φpq(xi,xj) and k, l and l' are integers [143]. The use of variance reduction methods has also been shown to reduce the sampling required for accurate Monte-Carlo integration of coefficients for RS-HDMR expansions [148,149]. An overview of the algorithm for developing high dimensional model representations is provided in Figure 4.
Figure 4. Algorithm for high dimensional model representation based on a minimum prediction error criteria.
Figure 4. Algorithm for high dimensional model representation based on a minimum prediction error criteria.
Processes 01 00067 g004
In pharmaceutical applications, Banjeree et al. [150,151] have used HDMR to generate input-output maps for use in feasibility analysis of black-box processes. Boukouvala et al. [12,123] have extended this methodology to determine the design space for a continuous mixing process using cut-HDMR. RS-HDMR can also be used as a tool for variance based global sensitivity and uncertainty analysis. The total and partial variances can be calculated based on the HDMR component functions where D represents the total variance, Di represents the contribution of variable i to the total variance and Dij represents the contribution of the interaction between variables i and j to the total variance.
Processes 01 00067 i020
Based on the variances and partial variances the sensitivities can be determined as:
Processes 01 00067 i021
The use of HDMR for sensitivity and uncertainty analysis has been demonstrated for several applications in environmental and atmospheric chemistry [152,153]. Ziehn and Tomlin have developed a Matlab® based tool for the implementation of HDMR and HDMR-based sensitivity analysis called GUI-HDMR. This tool has been used for sensitivity analysis of air flow models, complex reaction models of pollutants in air and cyclohexane oxidation [154]. An example of HDMR component functions generated using the GUI-HDMR tool developed by Ziehn and Tomlin is given in Figure 5.
Figure 5. Example of fit for of a first order HDMR metamodel for tablet API concentration in a direct compaction process obtained using the GUI-HDMR software of Ziehn and Tomlin [154].
Figure 5. Example of fit for of a first order HDMR metamodel for tablet API concentration in a direct compaction process obtained using the GUI-HDMR software of Ziehn and Tomlin [154].
Processes 01 00067 g005

3.3.4. Artificial Neural Networks

Artificial Neural Networks (ANN) have become quite popular in a variety of fields as a tool for addressing complex science and engineering problems and for developing empirical process models [155]. In particular they are useful in modeling pharmaceutical processes due to their ability to accurately represent nonlinear system behavior [138,156]. Neural networks can be defined in terms of the transfer functions defined by their neurons, the learning rule applied and the connectivity of the system. The neurons are arranged in layers: an input layer which consists of data entering the network, hidden layers containing the neurons that transform the input data, and the output layer which contains the network output corresponding to a specific input or set of inputs. Each hidden layer can contain a single neuron or a group of neurons operating in parallel. The neurons, or processing elements, consist of weights and transfer functions. The weights are coefficients by which the inputs to the neuron are multiplied and the transfer functions are simple linear or nonlinear functions that transform the weighted inputs. The architecture of a network can be described as feedback or feedforward. Feedback networks include connections from output neurons back to input neurons and therefore have some memory of prior states while feedforward networks do not. The architecture of the network can be adjusted through a learning process which involves the use of a training dataset for which the true values corresponding to each set of inputs are known. The learning process is iterative, with weights being modified in order to minimize prediction error, most often via back propagation of the error. Training involves a trade-off between the ability of the network to accurately predict the training set and its ability to generalize to data not included in the training set [138,155,156].
The general procedure for developing an artificial neural network, depicted in Figure 6, involves
  • Obtain a training data set e.g., via DOE;
  • Define the network: the number of hidden layers, the number of neurons to include in each layer and the type of transfer functions to be implemented;
  • Use a training procedure to optimize the weights in such a way that prediction error is minimized. The number of neurons in each layer can also be determined based on the training set, via cross validation;
  • Test the developed network against data that was not contained in the original training set to verify that the network has not been over fitted.
Artificial neural networks (ANNs) have been used to model a variety of pharmaceutically relevant processes. They can be used to recognize patterns for analytical chemistry purposes, to evaluate the influence of molecular structure on material properties, to evaluate pharmacokinetic and pharmacodynamics profiles, [156] and to model pharmaceutical unit operations [16,138].
Figure 6. Algorithm for development of an Artificial Neural Network with the option to add points to the training set or to redefine the network structure in order to improve model predictive ability.
Figure 6. Algorithm for development of an Artificial Neural Network with the option to add points to the training set or to redefine the network structure in order to improve model predictive ability.
Processes 01 00067 g006

3.3.5. Comparison of HDMR, RSM, Kriging and Neural Networks

Previous works have compared Kriging, RSM, HDMR and ANN as black box modeling techniques for pharmaceutical operations [12,131,138]. In these studies it has been found that Kriging tends to outperform other methods with respect to prediction accuracy. For instance, Kriging resulted in a lower prediction error than RSM for predicting flow rate variability in a loss-in-weight feeder [131]. It also provided a more accurate prediction of design space for this unit operation [12]. This may be due in part to the fact that, unlike RSM and HDMR, Kriging provides a measure of prediction uncertainty through the calculation of prediction variance at each point. This uncertainty prediction can be used to suggest regions where additional sampling would be beneficial [139]. In addition, because Kriging does not assume a closed form for the fitted process model it may perform better in modeling complex nonlinear processes [138]. However the selection of a reduced-order modeling methodology is application dependent. In some cases it may be preferable to have a closed-form representation of the process model, in which case RSM would be appropriate. For global sensitivity analysis HDMR provides a convenient framework and thus this modeling methodology might be preferred. Artificial neural networks have the capacity to model pharmaceutical processes accurately, but may require large training sets. In addition, defining the network requires some trial and error which can be quite time consuming [138,155]. A brief comparison of these reduced-order modeling techniques in terms of the number of and type of fitted parameters and basis functions is provided in Table 3.
Table 3. Comparison of reduced order modeling methodologies, Kriging, response surface models, high dimensional model representation and artificial neural networks.
Table 3. Comparison of reduced order modeling methodologies, Kriging, response surface models, high dimensional model representation and artificial neural networks.
MethodFitted parametersNumber of fitted parameters*Common basis functions
Krigingvariogram coefficients, regression coefficients21correlation models: exponential, gaussian, linear, spherical, cubic, spline regression models: polynomial
RSMpolynomial coefficients15Polynomial
HDMRcomponent function coefficients20Analytical basis functions: orthonormal polynomials, spline functions
ANNneuron weights40Transfer functions: linear, threshold, sigmoid
* The number of fitted parameters for second order model with 4 inputs and a single output. For kriging model assume that the correlation model is exponential and the regression model is a second order polynomial. For neural networks look at a feedforward model with 3 layers (1input 1output and 1 hidden) and 4 nodes per layer with sigmoid units.

3.3.6. PCA Based ROM

High fidelity process models like those obtained through CFD, DEM or FEM provide detailed information about distributed parameters like fluid and particle velocities. However these models are computationally expensive to evaluate and thus may not be useful for simulation and optimization purposes. Reduced order modeling based on principal component analysis was introduced by Lang et al. [128] for the co-simulation of CFD models with unit operation models for process equipment. Boukouvala et al. [110] have shown that the same approach can be applied to the reduction of DEM data for use in solids process models.
Principal Component Analysis (PCA) is a statistical tool that can be used to reduce the dimensionality of a dataset through orthogonal transformation. The original dataset, denoted as X, consists of n observations of m variables, which may or may not be correlated. The principal components obtained through PCA are orthogonal to one another and are arranged in order of decreasing variance such that the first component explains the greatest amount of variance in the original dataset and the last component explains the least. The number of components in the model (a) can be selected to achieve the desired percent variance explained and must be less than or equal to the number of variables in the original dataset.
The results of PCA can be expressed in terms of the component scores (T), which are the variables transformed into the latent space, and loadings (P), by which the original variables can be multiplied to obtain the scores. The number of columns in the scores and loadings matrices is given by the number of principal components in the model (a). The original dataset can be reconstructed from the scores and loadings as shown in Equation (14).
X = TP' + ∊
where is a matrix of residuals.
An advantageous feature of PCA from a modeling perspective is that it can tolerate missing elements in the original dataset X reasonably well. Several methods for handling missing data have been proposed in the literature including Expectation-Maximization PCA, the iterative NIPALS algorithm and a nonlinear programming based strategy [157,158,159,160].
The objective of PCA based ROM is to develop a mapping between input variables and the scores or loadings obtained from PCA. The scores and loadings can then be used to approximate the high dimensional data as in Equation (14). In the case of DEM this information can include particle velocities, energies and forces at a number of discrete locations within a processing unit. The basic algorithm for developing a PCA based ROM is outlined below and is also depicted in Figure 7 [110,124].
  • Identification of the input space (U), the state space (X) and the output space (Y)
    The input space consists of operating parameters that can be controlled. The output space contains measured responses at the end or outlet of the process. The state space contains variables monitored within the processing unit at various discrete points. The dimensionality of the state space depends on the spatial discretization of the unit. In the case of a continuous mixing operation the inputs might include blade rpms and configurations as well as fill level and total feed rate while API concentration and relative standard deviation (RSD) at the mixer outlet would make up the output space. The state space could include average particle velocities and energies at discrete positions within the geometry of the mixer.
  • Determination of the domain of the input space and implementation of an experimental design to define the input sampling space. Performing the computer simulations at the defined sampling points.
    The levels of input variables to be investigated can be defined based on the operating regime for the process of interest. A design of experiments (DOE) can be used to sample the input space appropriately, resulting in a total of N distinct sampling points within the input space. Simulating the process at each of the sampling points provides the corresponding state space and output space data.
  • Definition of the discretization of the process geometry in order to extract the state space data
    Boukouvala et al. [110] have indicated that the choice of discretization is critical to successful ROM development. Therefore care should be taken in selecting the mesh for the geometry.
  • Performing PCA on the state space data
    PCA can be performed as discussed above. Note that PCA is conducted separately for each state. Thus if particle velocity data in the axial and radial directions is extracted from a simulation PCA must be conducted on each velocity component separately.
  • Mapping the input space to the output space and to the scores or loadings for the state space
    The functional form of the input-output mapping U → Y and the input-scores or input-loadings mapping U → P is determined by the modeler. Lang et al. [124,128] have described the use of Neural Networks for this mapping, while Boukouvala et al. [110] have described the use of Kriging. Regardless of the type of mapping used, it is important verify the model accuracy e.g., via cross validation.
Figure 7. Algorithm for reduced order modeling based on principal component analysis.
Figure 7. Algorithm for reduced order modeling based on principal component analysis.
Processes 01 00067 g007
Once a ROM has been developed, it can be used to predict output and state space data for new input points that were not part of the original experimental design. The mapping developed is applied to the new point uk to predict a new vector of scores or loadings for the state space tk and a new vector of responses in the output space yk. The full dimensional state space data can then be reconstructed from the loadings using Equation (14).
Boukouvala et al. [110] have demonstrated the use of PCA based ROM to predictively model distributed particle properties including total force, kinetic energy and velocity within a continuous mixer. The developed model was also used to optimize the blending efficiency of the mixer. It would have been impractical to use the original DEM model for optimization purposes due to its high computational cost. Thus the benefit of PCA based ROM is clearly demonstrated by this work.

3.4. Integrated Flowsheet Modeling Tools

Individual unit operation models can be combined into integrated flowsheets using process simulators, which allow for the collection of various types of process models in a single programming environment. This facilitates the development of integrated flowsheet models that can be used for process simulation and optimization as well as design space exploration and sensitivity and uncertainty analysis [14,161,162]. Flowsheet models can also be used to evaluate control strategies for continuous processes [17,19]. The use of process models can reduce the experimental burden associated with process development and help to anticipate and resolve challenges in process scale up [15]. An overview of commercially available for software packages for process simulation is provided in reference [163]. The current work will focus on simulators that can be used to model continuous solids processes, specifically gPROMS™ and ASPEN™.

3.4.1. gPROMS™

gPROMS™ is a modeling platform developed by Process Systems Enterprise, Ltd. It consists of a model development environment (ModelBuilder) that comes with several built-in model libraries. Related PSE applications for solids processing include gSOLIDS and gCRYSTAL, which contain their own model development environment and unit operation models for many solids processes including but not limited to batch, semibatch and continuous crystallization, roller compaction, granulation, milling and fluidized bed drying. Users can also take advantage of custom modeling capabilities. This facilitates the implementation of data-based or hybrid models in conjunction with existing process model library components. gPROMS™ ModelBuilder supports complex models including algebraic, differential, partial differential and integral equations [164,165].
The use of gPROMS™ for the simulation of continuous pharmaceutical processes is documented extensively in the literature. The simulation of continuous tablet manufacturing via wet granulation has been described by Boukouvala et al. [161]. The simulation of manufacturing and control of tableting processes including roller compaction and direct compaction has been documented by multiple authors [14,17,18,19]. A gPROMS™ flowsheet describing a direct compaction process along with corresponding simulation results is shown in Figure 8. gPROMS™ has also been used for combined simulation of upstream and downstream pharmaceutical processes [162].
In addition to process simulation, the development environment has parameter estimation and optimization capabilities. gPROMS™ can be used to perform both steady state and dynamic optimization. It can handle both equality and inequality constraints and can optimize over both continuous and discrete variables. gPROMS™ parameter estimation tool uses a maximum likelihood framework to estimate parameters within process models [166,167]. The use of both of these functionalities has been documented in the literature for industrial processes [168,169]. Both functions have also been used extensively in the simulation and optimization of batch crystallization processes for pharmaceutical applications [170,171,172].
Figure 8. Direct compaction flowsheet simulated in gPROMS™ showing simulated tablet hardness as a function of lubricant concentration.
Figure 8. Direct compaction flowsheet simulated in gPROMS™ showing simulated tablet hardness as a function of lubricant concentration.
Processes 01 00067 g008
While a large number of models are available in the gPROMS™ model libraries, the connection of these models and their integration into processes that will run without error is not always straightforward. Models may take different connection types, and while units akin to stream converters can be used to address this issue it is not always immediately apparent when and why these are required. In addition, some models have a large number of parameters which require specification. The proper values to select may not always be intuitive, though they can usually be determined via heuristics if direct measurement is not possible. Default parameter values are provided by the software, and these can be used if appropriate values cannot be selected any other way. Additional training, provided by Process Systems Enterprise, can supplement the extensive user manual associated with the development environment to familiarize users with the proper implementation of unit operations.

3.4.2. Aspentech’s AspenOne® utilizing AspenPlus® (V7 and V8) and AspenCustomModeler®

AspenTech developed a suite of engineering tools by focusing on chemical analysis in areas such as process control, process engineering, optimization, and supply chain management. Aspen Plus® is the spearhead of AspenTech process modeling products which simulates industrial chemical processes. With the recent acquisition of SolidSim Engineering GmbH, Aspen Plus® Version 8 can model particle behavior with solids handling and separating units [173]. Aspen describes the particulate systems in discrete size classes by a user set number of intervals; each interval varies equally from a lower to upper bound with units selected from angstroms to kilometers. For solids processing, Aspen can handle conventional and non-conventional solid streams alongside mixed liquid/vapor streams with or without particle distribution. These streams can then be manipulated under 19 different types of physical property methods; the most common method as ideal property method including Raoult’s Law and Henry’s Law [174]. The process models in the current Aspen Plus® package consist of convective drying, granulation and agglomeration, crushing and milling, as well as classification and separation units [173,174]. Within varying units, the operating mode and type of unit can be altered; for example, a crystallizer’s operating mode can vary from “Crystallizing”, “Dissolving or melting,” or “Either” while a crusher can vary from cage mill, to cone, to 8 other options; this trend continues for a variety of the available process units [174].
In addition, users can add their own models that utilize unique differential and algebraic equations to simulate behavior of their unit operations. Custom units can be programmed into Aspen Plus® through Fortran, C++, and visual basic; AspenTech also includes a Custom Modeler (ACM) to program user-defined units that can handle any system described by logical, nonlinear, differential and algebraic equations [174,175]. These user defined units can be exported from ACM to the drag and drop interface to interact with the Aspen Plus® simulation [175]. After a user programmed unit becomes exported into Aspen Plus, the unit can be manipulated without programming knowledge; this allows for quicker manipulation and testing of processes. Also, since the unit can be simply drag, dropped, and defined in contrast to reprogramming variables, an additional accessibility encompasses the scientist that’s knowledgeable in chemical processing but not coding. Units exported into Aspen Plus gain the advantage of interacting with the default set of Aspen unit models, the Aspen Plus optimization and analysis utilities, as well as the ease of access that comes with a drag and drop interface.
Within Aspentech software, manipulating the required set of input variables either through the calculation block or equation-oriented modeling becomes available. The user can set up the calculation block to determine certain input parameters using a predefined algorithm and variables within the simulation [176]. The equation oriented modeling can switch inputs and calculated outputs; for example, one can switch from defining the reflux ratio input to the product purity output in a distillation column. AspenTech software can handle both steady state and dynamic process simulation and contains its own optimization software for model and equation oriented simulation based optimization. In order to model time varying conditions, Aspen Plus Dynamics® can add time dependent factors into steady state dynamics developed in Aspen Plus® [175]. Although many Aspen Plus® applications involve continuous processes, the software also can handle batch processes.
One advantage to using AspenTech software is the availability of resources such as guides and forums to aid the user. These can be found both through AspenTech and from independent sources [173,174,176]. The software also comes with an extensive manual along with an integrated guide within the software to ensure proper inputs are given in the varying screens. Along with the included manuals and guidelines, AspenTech software includes an extensive physical property and thermodynamic property library [175].
Use of the AspenTech suite of tools for modeling of continuous, solids-based processes for pharmaceutical applications is not extensively reported in the literature. However pharmaceutically relevant solids-based processes have been successfully modeled using this software. For instance, Aspen Plus® has been used to model batch crystallization processes, which include both liquid and solid phases [177,178]. Given the recent acquisition of SolidSim Engineering by Aspen Technology, Inc. it is possible that publications related to modeling of continuous solids-based pharmaceutical processes in Aspen Plus® will increase in the coming years.
The challenges accompanying model integration in AspenPlus® are comparable to those experienced within the gPROMS™ platform. In addition, while custom built gPROMS™ models can be readily integrated with built-in model library elements (provided inputs and outputs are properly specified) the incorporation of custom models with pre-existing models in AspenPlus® requires a model export process from Aspen Custom Modeler to AspenPlus®. This model export process relies on third party software like Microsoft Visual Studio®. Alternatively custom models can be developed within AspenPlus using Fortran. It would be preferable to be able to integrate existing model library elements with custom models within AspenPlus® itself using the Aspen Custom Modeler language, which is more straightforward than Fortran for unit operation modeling.

4. Unit Operation Models

The focus of this section will be equation-oriented modeling for unit operations that are used in downstream pharmaceutical processes. Many of the operations discussed have also been studied extensively via discrete element or finite element method models. While reference may be made to these studies so that the reader can pursue them further, they will not be considered at length in the current work because such models are not readily used in process simulation and optimization on account of their computational cost.

4.1. Hopper

Several authors have proposed continuum models for the flow of cohesionless granular materials in hoppers. Early efforts focused on empirically describing mass flux out of hoppers for free flowing, frictionless, cohesionless materials under steady flow. For instance Brown [179] proposed that for a cone geometry flow can be described as in 15, where β is the hopper wall half angle (θw)when the hopper is operating in the mass flow regime.
Processes 01 00067 i022
These models tended to assume that behavior near the hopper outlet determines flow rate. Savage [180] and Savage and Sayed [181] introduced an approach based on the laws of motion to predict flow rates for frictional materials in conical and wedge shaped hoppers. For conical hoppers that are sufficiently full mass outflow can be obtained from 16 where φ is the internal friction angle of the granular material.
Processes 01 00067 i023
Equation (16) was found to significantly over estimate flow rates. Numerous extensions to the approach based on laws of motion have been proposed. The introduction of a perturbation to account for the effect of wall friction was found to significantly improve predictive quality [180,181]. Subsequent authors have modified the perturbation method to improve model convergence and predictive ability [181,182,183]. Models of granular flows have continued to increase in complexity to describe additional process physics. Plastic potential models, which consider rigid-plastic or elastic-plastic materials, and the double-shearing model, which can be used to extend models from incompressible to compressible flows, are models that have been used to this end [184]. Recently Weir et al. [24] have proposed a continuum model for variable-density, fully plastic, granular flow. The model was applied to study steady radial flow of a cohesionless granular material in steep-walled wedge and conical hoppers. This model extends previous continuum models by considering compressible flows, allowing for the consideration of pressure effects on density within the hopper. The variation in bulk density as a function of pressure is based on the relationship between pressure and voidage empirically determined by Weir [185].
dP0 = −1.57 × 1011 ρsgdpe−33.1ε
In Equation (17) dP0 is the change in pressure and is the corresponding change in voidage. ρs is the solid density, dp is the particle diameter and g is gravitational acceleration. The flow rate out of the hopper (J), is determined using a set of linear and nonlinear algebraic and partial differential equations given in [24]. The proposed model has shown a qualitative agreement with the observed density variation in hoppers and approximate agreement with experimentally obtained discharge rates for both conical and wedge hoppers.
Sun and Sundaresan [184,186] propose a constitutive model for rate independent granular flows with microstructure evolution, which can be used to determine internal properties of the granular material. This model has been applied to the flow of an incompressible, cohesionless granular material in a conical hopper. The model has been validated against DEM simulations, both with respect to its ability to predict hopper flow rate and its ability to qualitatively describe microstructure within the hopper.
The previously discussed models are particularly relevant in cases where the pressure exerted on material at the hopper outlet by the powder bed above is significant. However in many continuous pharmaceutical processing applications the hoppers used are relatively small and the residence time therein may be quite short. Boukouvala et al. [14] propose a straightforward equation oriented approach that can be applied for such processes.
Processes 01 00067 i024
Processes 01 00067 i057 (t) = H (t)bulk (t)
Equation (18) describes a basic mass balance around the hopper while Equation (19) describes the mass holdup in the hopper, Processes 01 00067 i057, as a function of the height of the powder within the hopper (H), the cross sectional area of the hopper, A, which may vary as a function of the height and the bulk density of the powder. Equation (19) assumes that bulk density is constant throughout the hopper, which is a reasonable assumption for small hoppers with relatively low holdup.
There is also a mean residence time associated with the hopper. This is accounted for by applying a time delay to the propogation of material properties through the hopper as shown in Equation (20).
Processes 01 00067 i025
xi is any relevant property for component i and θrt is the mean residence time for the hopper, which can be determined experimentally.
The above model is appealing due to its relative simplicity. However it should be noted that the model makes several significant assumptions, including that the hopper is operating in the mass flow regime, that the bulk density is constant throughout the hopper and that the experimentally determined mean residence time does not vary significantly over the course of time. These assumptions should be validated experimentally prior to implementing this model.

4.2. Loss-in-Weight Feeders

Feeders are not typically considered in the absence of other processing equipment, so model development for loss-in-weight feeders is not as extensive as for other unit operations. Engisch and Muzzio [30] have described a combined experimental and statistical approach to characterize the performance of loss-in-weight feeders as a function of screw speed, discharge screen size and screw rpms. Analysis of variance (ANOVA) was used to determine the effect of these parameters and their binary interactions on feeder performance. Fast fourier transforms (FFT) were also used to obtain analyze the feed rate data from the equipment. The resulting power spectra could be used to help determine the feeder’s characteristic time.
Boukouvala et al. [12] have demonstrated the use of response surface models, high dimensional model resolution and Kriging to optimize feeder design and operating parameters for a given material. The model is built based on experimental data. Its inputs include screw speed (rpm), screw size, screw configuration (open or closed helix) and powder flow index. The response chosen to indicate feeder performance is relative standard deviation of the outlet flowrate–an indication of variability.
Boukouvala et al. [14] have also described a semi-empirical equation oriented approach to dynamically model a loss-in-weight feeder. Equation parameters that are fitted from experimental data include the process gain Processes 01 00067 i060, a time constant (τ f) and a time delay factor Processes 01 00067 i068. For the purposes of this model, material bulk density and mean particle size is assumed to be constant throughout the feeder.
Processes 01 00067 i026
Processes 01 00067 i062
xin(t) = xout(t)
In Equation (21) ωf is the screw speed (rpms). Equation (22) represents the time delay associated with the feeder, where z is the time delay domain, Processes 01 00067 i063 is the experimentally determined time delay and Processes 01 00067 i064 is the output feed rate. Equation (23) indicates that for invariant properties, like material bulk density, the feeding process has no effect. This should be experimentally verified prior to using this model for a particular process.

4.3. Continuous Mixer

The objective of a powder mixing process is to produce a blend within minimal composition variability. For continuous mixing the blend composition variability is considered with respect to both space and time. Continuous blending processes can be characterized based on several metrics related to product homogeneity including the relative standard deviation (RSD) of the composition at the blender exit, the variance reduction ratio (VRR) and the residence time distribution (RTD) [34]. The relative standard deviation can be calculated based on the exit concentration of a particular component, the API in most pharmaceutical applications, over time.
Processes 01 00067 i027
The variance reduction ratio can be calculated based on the composition variance at the input and output of the blender.
Processes 01 00067 i028
The residence time distribution is a probability distribution that describes the amount of time a solid or fluid element is likely to remain in a particular unit or process. The residence time distribution can be calculated experimentally by injecting pulses of a tracer molecule into the blender and then monitoring the concentration of these particles at the mixer exit over time [41,187]. RTD can be expressed mathematically as an integral in terms of the concentration of tracer particles over time.
Processes 01 00067 i029
It has also been shown that the residence time distribution can be used as a modeling tool for the design of mixing processes. The performance of a mixing process can be viewed as the result of local mixing rate and total mixing time, which is function of residence time. Thus the RTD can be used as a means of determining appropriate throughput such that material remains in the unit long enough to become well mixed [109].
Continuum models for continuous mixing processes have consisted largely of population balance models [14,115,119,123]. A general form of a multidimensional balance that could be used to describe a mixing process is given in Equation (27), where n is the number of components in the mixture, x and y are external (axial and radial) spatial coordinates that describe the position in the blender and r is an internal coordinate that describes the particle size. For pharmaceutical applications the value of n is often three; API, excipient and lubricant.
Processes 01 00067 i030
The number density function F(n,x,y,r,t) indicates the number of particles of type n with size r at position x,y at time t. Thus the particle velocities in the x and y directions are given by Processes 01 00067 i031 and Processes 01 00067 i032. For implementation, the process geometry is usually discretized in space so that the differentials with respect to y and x can be replaced with finite differences. The same can be done for the particle size. It is typically broken up into size categories or bins defined by upper and lower limits on particle size. In order to calculate the right hand size of Equation (27) the material balance for the mixer is needed. This is provided in Equation (28). Equation (29) provides a means of determining composition-dependent properties at the mixer outlet.
Processes 01 00067 i033
Processes 01 00067 i034
Finally Equation (30) describes the calculation of active ingredient concentration in terms of mole fraction at the mixer outlet. From this RSD and VRR can be calculated [14,115].
Processes 01 00067 i035
Predictive models for particle velocities Processes 01 00067 i031 and Processes 01 00067 i032 that appear in the PBM can be developed based on DEM data. Boukouvala et al. [110] have demonstrated the use of a technique called DE-ROM to develop predictive models based on data obtained from DEM and have shown that it can be used to develop a predictive model for particle velocities throughout the mixer as a function of input parameters like blade rpm. The developed reduced order model can then be implemented to provide velocity information for the population balance model.
Periodic section modeling is a method whereby a continuous mixer is modeled as a series of mixer segments in which transverse mixing occurs. Each segment represents a cross sectional slice of the mixer with a limited axial size. Periodic section modeling operates on the principle that continuous mixing can be viewed of as a combination of two processes; powder flow along the length of the blender (axial movement) and transverse powder mixing within the blender. Axial movement is characterized by the axial velocity vx and the dispersion coefficient Ex. Within a periodic section, transverse mixing is understood to occur in a manner similar to batch mixing. The batch-like mixing can be characterized in terms of the variance profile. As variance decreases, the homogeneity of the blend increases. For a periodic section the decay of variance over time has been modeled as an exponential function as shown in Equation (31), where kb is the variance decay rate for the batch mixing process. Processes 01 00067 i065 and Processes 01 00067 i066 describe the composition variance within the periodic section as a function of time, at steady state, and initially (before mixing begins).
Processes 01 00067 i036
Variance can also be expressed as a function of position within the mixer. This empirical relation can be used to determine the mixer length required to ensure homogeneity of the mixture at the outlet.
Processes 01 00067 i037
In Equation (32) E(t, x) is the residence time distribution at position x. Processes 01 00067 i067 is the variance at position x and can be calculated as in Equation (33) where kc is the variance decay rate for a continuous mixing process
Processes 01 00067 i038
This allows continuous mixing processes to be readily compared with batch processes. The relationship between the batch and continuous mixing efficiency can be defined in terms of the variance decay rate as in Equation (34). This equation also allows for estimation of kb based on the values of kc and vx which are readily obtained through DEM simulation.
Processes 01 00067 i039
Gao et al. [188] have demonstrated the application of a periodic section model to mixing case studies with segregating and non-segregating materials. This type of model is shown to accurately characterize the continuous mixing process with respect to variance decay and relative standard deviation of the mixture under steady state conditions. In a subsequent publication the authors have described the use of a periodic section model to optimize a continuous mixing process [20]. In this instance the periodic section model is applied to cases with varying particle properties (diameter, density and cohesiveness) and operating parameters (fill level and blade speed). The results indicate that mixing performance cannot be improved by simply increasing blade speed due to the corresponding increase in axial velocity. The authors suggest strategies which increase impeller speed without significantly increasing axial velocity, such as reducing the blade angle or increasing the weir height at the end of the periodic section. Both of these strategies have been shown to significantly improve the mixing performance. Thus the periodic section modeling method can be used to suggest improved mixer design and operating parameters. [20] An example of a DEM simulation of a periodic section of a continuous convective mixer is shown in Figure 9.
Figure 9. Discrete Element Method (DEM) simulation of particles mixing in a periodic section.
Figure 9. Discrete Element Method (DEM) simulation of particles mixing in a periodic section.
Processes 01 00067 g009
The periodic section modeling approach has several significant advantages. It can be used to draw parallels between batch and continuous mixing. It can also be used to elucidate mechanisms that contribute to improved mixing that might not be apparent otherwise. Finally, the periodic section is significantly smaller than the entire blender. Therefore it is much less computationally expensive to simulate using the discrete element method.
Reduced order modeling approaches including Kriging, HDMR and response surface modeling can also be used to model the entire continuous blending processes, eliminating the need to evaluate the population balance equation. Boukouvala et al. [12,129] have modeled blender performance (RSD) as a function of impeller rpm, powder flow rate and blade configuration using both Kriging and HDMR. These reduced-order models were used to develop a design space for the blending process [12]. One downside of the use of input-output reduced-order models is that they do not provide information about the state of the material within the mixer. Material fluxes, concentration and RSD cannot be determined as a function of position within the mixer using these models, as they typically map input parameters to conditions at the outlet of the mixer and not to the states within the mixer. However, using the DE-ROM approach described above a mapping between input parameters and states within the mixer can be developed [110].

4.4. Wet Granulation

Wet granulation processes are frequently modeled using the population balance equation, the general 2 dimensional form of which is described in Equation (1). Three or four dimensional population balances can be used to model wet granulation processes. In three dimensional models the population distribution is considered with respect to granule solid content (powder or powder blend), liquid content (binding solution), and gas content (porosity). Four dimensional models consider also granule composition with respect to the specific components in the system (e.g., API and excipient) as a fourth dimension. General three and four dimensional population balance equations are shown in Equations (36,37) respectively. s, l and g indicate the solid, liquid and gas volumes in the granule, where the assumption of constant granule density is made. In the 4-D equation s1 and s2 are two different solid phases whose compositions in the granule are defined [46,113,161,189].
Processes 01 00067 i040
Processes 01 00067 i041
The right hand side of the population balance equation for wet granulation typically includes a nucleation kernel, a coalescence or aggregation kernel and may also include a breakage kernel. These describe granule birth (nucleation), granule growth and attrition and are indicated by Rnuc(s,l,g,t), Ragg(s,l,g,t) and Rbreak(s,l,g,t) respectively in Equations (35,36) [120,190].
In practice it is preferable to solve lower dimensional PBMs, as 4-D models can become prohibitively computationally expensive to solve. In order to reduce the dimensionality of a population balance model, one or more of the granule properties can be lumped into the other distributions [113]. For instance, if it is assumed that all granules of the same solid and liquid content also have the same gas content then it is possible to describe the 4-D population balance in 36 as a set of two 3 dimensional equations; the reduced PBE (Equation (37)) and a gas balance equation (Equation (38)). The gas balance equation relates the total gas volume G(s1,s2,l,t) in a particular discretization bin to the volume of gas in each particle g(s1,s2,l,t) and the inflow and outflow of gas for the entire bin [113,191,192].
Processes 01 00067 i042
Processes 01 00067 i043
Barrasso and Ramachandran [113] have investigated a variety of model reduction options and have determined that the best performance is obtained when the gas volume is the lumped parameter, as in Equations (37,38). This reduction gave performance closest to that of the full 4-D population balance model. This is consistent with expectations because gas volume is known to have little effect on aggregation and prior work has indicated that lumped parameters should be chosen such that they have little effect on aggregation rates in order to minimize the error resulting from order reduction [191].
The development of kernels for the nucleation, aggregation and breakage represents a major research area in population balance modeling for wet granulation. Of these the aggregation kernel has been studied most extensively [193]. Kernels can be developed mechanistically or empirically. Often a combination of the two approaches is used and kernels are determined semi-empirically, with parameters fitted based on experimental or simulated data [189,194,195].
Nucleation accounts for the creation of new granule nuclei from primary powder particles. Poon et al. [195] and have described a mechanistic kernel for nucleation that accounts for both wetting kinetics and nucleation thermodynamics. The kernel can be expressed generally as Equation (39) for zero order and Equation (40) for first order dependency on primary particle concentrations. A droplet-controlled nucleation regime is assumed.
Processes 01 00067 i044
Processes 01 00067 i045
N is the number of primary particles in the nucleus and x(s,l,g) is the fraction of particles with properties s, l, g within the nucleus. The nucleation coefficient knuc depends on the binder solution flow rate (Q) and process temperature (T) as well as the ideal gas constant (R). A0 is an adjustable parameter and λ is described as a spreading coefficient, which relates the nucleation coefficient to the particle properties (s, l, g).
Processes 01 00067 i046
N and x(s,l,g) can be determined based on the volume of the nucleus, the effective porosity and the liquid and solid granule volume.
Aggregation or coalescence kernels (Ragg) describe the rate of granule growth, which can occur through consolidation of wetted particles or through layering of fines onto existing granules [48]. The rate of aggregation depends on the probability of wetted particles colliding as well as on the likelihood that they will adhere to one another upon making contact. This in turn is influenced by liquid binder content and granule size as well as the particle velocities [190]. Numerous aggregation or coalescence kernels have been proposed in the literature. While it is beyond the scope of this work to consider them all here, interested readers are referred to Cameron et al. [190] and Liu et al. [193], both of which provide tabular summaries of coalescence kernels discussed in the literature with associated references. The general form of a coalescence kernel contains two terms; a rate constant that relates growth rate to operating conditions and size dependence term that shows the relationship between granule size and growth rate. It is the latter term which dictates the size distribution of the granules and it is thus of great interest to be able to model it accurately [193]. Theoretical models that have been developed to predict the likelihood of aggregation based on powder properties and binder properties often rely on force or energy balances. For instance, several authors [122,193,196] have implemented a coalescence kernel in which granules combine if the kinetic energy of the granules is dissipated through viscous interactions or through plastic deformation of the particles.
Despite the development of mechanistic aggregation kernels, empirical or semi-empirical kernels are often used in simulations due to the difficulty associated with developing mechanistic kernels that are detailed enough to accurately model the process without becoming too computationally expensive to evaluate [193]. For instance, Immanuel et al. [122] describe the theoretical development of a rigorous, mechanistic aggregation kernel, but in order to demonstrate a robust solution strategy for 3D population balances use a semi-empirical kernel.
Many of the aforementioned aggregation kernels apply to single component systems. However in cases where binder is added continuously it may be appropriate to consider the binder as a separate component. Marshall et al. [197] have described a 2 component aggregation kernel for continuous fluid bed granulation where the binder and powder are treated as separate components. Matsoukas et al. [198] have proposed a multi-dimensional aggregation kernel to consider the effect of granule composition on aggregation rate. This is particularly useful for pharmaceutical applications, as multi component powder blends are typically the material being granulated in the pharmaceutical industry. Composition dependence is particularly important in cases where blends are not well mixed prior to granulation or where the materials being blended have disparate properties that affect their propensity for coalescence.
Breakage kernels are relevant for both granulation and milling population balance models and are therefore discussed in Section 4.6.
Granule drying is an inherent part of the wet granulation process. It typically carried out in a fluidized bed, so the models discussed will focus on fluidized bed drying. Mathematical models for granule drying vary in their level of detail from lumped parameter models to more detailed distributed parameter models. Distributed parameter models include continuum models, which use macroscopic laws to describe phenomena like heat and mass transfer, and discrete models like the pore network model which explicitly describe microscopic behavior [199]. Discrete models based on CFD can be used to describe granule fluidization within the dryer. Distributed models based on the population balance equation can be used to represent the distribution of properties like liquid content [112,161,200]. These models can be coupled with experimental information or data from more detailed mechanistic drying models to develop semi-empirical relations for drying rate or granule size as a function of dryer parameters like gas temperature and velocity [112,199]. For instance, Mortier et al. [112] have described a one-dimensional population balance model for granule drying that incorporates information obtained from more detailed mechanistic model. The PBM, shown in Equation (42), contains a single internal coordinate, the wetted radius Rw, and a negative growth term, Gr, associated with the reduction in particle size due to granule drying.
Processes 01 00067 i047
The growth term in Equation (42), can be determined empirically as a function of the drying gas temperature based on a more detailed mechanistic model as described in Mortier, Van Daele et al [201]. A comprehensive review of various mechanistic models for fluidized bed drying is provided in Mortier, De Beer et al [200].

4.5. Roller Compaction

One of the first models to describe behavior of roll compaction processes was the rolling theory of granular solids proposed by Johanson [202]. Bindhumadhava et al. [62] have shown that this theory can accurately predict pressure profiles in the nip region of a roller compaction process, as well as the influence of material properties on nip angle and peak pressure. Johanson’s model makes several simplifying assumptions, including isotropic, cohesive material, no-slip between the powder and the roll surface in the nip region (the material is frictional) and that all material within the nip region is compressed to a ribbon with a thickness equal to the exact gap between the rolls [203]. The required model inputs include effective angle of internal friction and angle of friction, which can be found experimentally, and the relationship between pressure and density for the powder of interest, which can be determined experimentally using a punch-die system similar to that found in a tablet press. The assumptions associated with Johanson’s model are reasonable for gravity fed roller compaction systems with smooth rollers of relatively large diameter (>500 mm) and powders with a high enough friction coefficient that the no-slip condition holds. In such cases this model agrees well with experimental data [63].
Johanson’s model describes the pressure gradients for the slip and nip regions according to Equations (43,44). In these equations, θ denotes the roll angle, S denotes the roll gap, D indicates the roller diameter, μ is the friction coefficient between the powder and the roll, δE is the effective angle of internal friction, the constant K indicates material compressibility and σ is the normal stress applied to the powder. The parameter A can be calculated from Equation (45), where W is the angle of wall friction [61].
Processes 01 00067 i048
Processes 01 00067 i049
Processes 01 00067 i050
The nip angle (α) is determined based on the assumption that the pressure gradients at the boundary between the no-slip and nip regions are equal. Thus setting Equations (43,44) equal to one another it is possible to calculate α [62].
Based on Johanson’s model, the relative density exiting the process can be calculated from Equation (46), where γ0 is the pre consolidation relative density, Rf is the roll force, and W is the roll width. The roll force can be calculated based on the roll diameter and width, the roll gap, the peak pressure, the angle θ and the compressibility K [61].
Processes 01 00067 i051
The relevant equipment parameters in Equations (43–46) can generally be determined from the equipment geometry [63]. In some cases the roll force is determined experimentally as a function of the hydraulic pressure set point for the equipment. The material properties such as the effective angle of internal friction, angle of wall friction, compressibility and preconsolidation relative density must be determined experimentally [61]. The use of instrumented roller compaction equipment to empirically obtain pressure profiles and torque information has been demonstrated by several authors [58,62]. This information can be used to empirically determine friction coefficients and angles as well as roll force. In addition, the obtained pressure profiles can be used to verify the pressures predicted by Johanson’s model.
One of the drawbacks of Johanson’s model is its inability to accurately predict ribbon density for incompressible materials. This is due to the use of the Jenike-Shield yield criterion. An alternative model has been proposed by Marhsall [204] which facilitates modeling of both compressible and incompressible materials by using the Coulomb-Mohr criterion. Although it is more generally applicable, the resulting model is far more complicated than Johanson’s model. Another alternative is the so called “slab method”, which uses a force balance on a slab of material to predict pressure distribution in powder under compression. The yield criterion for this model can also be adjusted to consider incompressible material such as metals [63,203]. Since most pharmaceutically relevant powders are compressible it is not typically necessary to consider alternate yield criterion, but other assumptions within Johanson’s model (e.g., that of smooth, relatively large rolls) are not always accurate.
Finite element models, which are akin to the discrete element method models discussed in Section 3.2, have been shown to describe roller compaction processes with a high degree of accuracy [63,205]. The accuracy of these models can be attributed to the fact that they consider particle-particle and particle-roll interactions at a fundamental level. Information about the friction conditions, contact angles, process geometry, roll surface texture can thus be considered. However the utility of these models is limited by computational expense [203].
For simulation and optimization purposes, Hsu et al. [203] have proposed a dynamic model that is based on Johanson’s theory of rolling granular solids but incorporates the effect of changing roll gap on ribbon density. Johanson’s model assumes a constant roll gap. The ability to consider the process response to dynamic variability in roll gap is useful, as the roll gap typically varies with changes in the feed rate to the roller compactor. The model in Hsu et al. [203] combines Johanson’s model with a material balance to account for the variability in roll gap as a function of powder feed rate. In addition, delay differential equations are incorporated to account for the time dependent response of the feed rate, roll speed, and roll pressure to changes in the set point. This model has been implemented in gPROMS™ by Boukouvala et al. [14] for the simulation of a continuous tablet manufacturing process that incorporates dry granulation.

4.6. Milling

Size reduction by milling is common at both the drug substance and the drug product step of pharmaceutical manufacturing. Milling can be used to reach the desired particle size for an active pharmaceutical ingredient or granule size for a pharmaceutical blend to be used in tableting. Milling processes are often modeled using population balance equations similar to those discussed in Section 3.1 and ‎4.4. The PBEs used to describe a milling process can be discretized with respect to time or spatial coordinates depending on the mill design and residence time [206,207]. Both one-dimensional and multi-dimensional population balance models for milling have been described in the literature [16,206,208]. However for pharmaceutical applications it is preferred to implementa multi-dimensional model as the materials being milled are typically granules, which may consist of multiple solid phases as well as a gas phases. A multi-dimensional population balance facilitates consideration of composition variability with respect to each internal coordinate [207].
A general form of a multi-dimensional population balance that incorporates breakage is given by Equation (47) where s1 and s2 represent two different solid phases (e.g., API and excipient) and g represents a gas phase and t is time. For a wet milling process a liquid phase could also be included.
Processes 01 00067 i052
The breakage term, Rbreak(s1,s2,g,t), in Equation (47) accounts for the likelihood that a particle in a particular size class will break into smaller particles (breakage probability) and the distribution of sizes into which the particle will fragment (breakage distribution) [206,208]. A general form of the breakage term is given in Equation (48), where the integral portion indicates the formation of particles in smaller size classes due to breakage of larger particles into smaller fragments and the second term indicates the loss of particles from a larger size class due to breakage. The breakage kernel, kbreak and the breakage function, b, account for the rate at which particles break and the size distribution into which they fragment [14,66,198,209].
Processes 01 00067 i053
If it is assumed that the breakage rate is constant for a given size class, a simpler linear breakage kernel may be used. For short milling durations this assumption may be valid [16]. However, if this assumption is not made then the associated breakage rate expression should be formulated to incorporate nonlinear effects [206,208]. Breakage kernels can vary in complexity depending on the level of detail accounted for in the model. For multi-component systems a composition specific breakage kernel may be used to account for variable breakage rates and size distribution as a function of granule API or excipient content [14,198]. Semi-empirical breakage functions in which parameters are estimated from experimental data are common in the literature [66,208]. Heuristic approaches are also used in which breakage distributions of a particular form, often corresponding to the limiting case for an assumed fragmentation pattern, are applied [209,210].

4.7. Tablet Press

Modeling of the tablet press unit operation includes modeling of powder flow into the dies and modeling the compaction process for powder blends. The modeling of flow into dies includes consideration of the feed frame, which supplies powder to the dies. Relatively little work has been done in the literature to independently model the feed frame aspect of the tablet press, but Boukouvala et al. [14] have described a response surface model for the residence time in the feed frame as a function of feed frame rotation rate and disk rotation rate.
The compaction of powders in confined geometries has been modeled extensively. Detailed models often rely on discrete or finite element methodology as discussed in Section 3.1 [72,95,96,98,99,100]. Semi-empirical models that describe the compaction behavior of granular materials by relating compaction pressure to relative density are also well established in the literature. A review of models describing compression behavior is provided in Patel et al [75]. Of particular interest are the Heckel [211,212] and Kawakita [213,214] equations, variations of which have been extensively used to model the compaction behavior of powder blends within tablet presses [74].
The Heckel equation (Equation (49)) assumes that the relationship between powder densification and applied pressure is first-order.
Processes 01 00067 i054
In Equation (49) D is the relative density, P is applied pressure and K is a material-specific constant. A is a constant densification term, where Processes 01 00067 i055 relates to the initial die filling process and B corresponds to densification from slippage and rearrangement of particles [75,211,212]. The Kawakita equation (Equation (50)) assumes that the relationship between applied pressure and powder volume is constant because particles are at equilibrium throughout the compression process.
Processes 01 00067 i056
a is the initial porosity of the material and 1/b is the pressure that would result in a compression of the powder by one half of the total possible volume reduction [213,214]. The parameters for the Heckel and Kawakita equations can be found experimentally from force-displacement data [75].
Kuentz and Leuenberger [215] have described a modified Heckel equation that accounts for the effect of relative density on susceptibility to pressure, defined as the change in relative density in response to applied pressure. This model also incorporates a critical density, beyond which pressure susceptibility cannot be defined due to a lack of rigid structure. Singh et al. [216] describe a detailed model that can be used for process modeling and control of a tablet press which is based on the Kawakita equation. This model also incorporates an equation for tablet hardness as a function of compression force, proposed by Kuentz and Luenberger [73].
Modifications to the Kawakita equation to consider binary powders have been proposed to predict compaction behavior for powder blends. This is applicable in pharmaceutical case studies, as the material being compressed in a tablet press typically contains at least two components, an API and excipient. Frenning et al. [217] suggest the use of effective Kawakita parameters for mixtures, which are calculated based on the parameters for the individual components and a mixing rule that assumes the component volumes are additive. Mazel et al. [218] propose a method that does not use effective Kawakita parameters, but applies the Kawakita equation to each component in the blend separately and assumes that volumes are additive. The implementation of these mixture models is particularly useful for blends in which the various components have very different physical properties.

5. Model Validation and Verification

Model verification and validation is a necessary part of the development process. Specifically verification describes the process of determining whether the desired conceptual or mathematical model has been properly implemented while validation describes the process of ensuring that the developed model is sufficiently accurate for its intended purpose [219,220]. Verification and validation can be used to establish the predictive ability of the model and justify the use of the underlying theories and assumptions associated with the model equations [221]. Model validation is also an inherent part of the quality by design paradigm outlined in the ICH Q8 guidance for pharmaceutical development [3]. Model verification and validation should be carried out in parallel with process and model development, so that at each stage the existing model can be assessed and potential gaps in process knowledge and the associated model equations can be identified [222]. Issues with model implementation can also be addressed [219,221,223].
In order to validate a model, it is first necessary to identify its intended purpose, as this will dictate the validation criteria in terms of the model outputs to be considered, the domain over which the model should be valid and the level of prediction accuracy required [221,223]. The approach to model verification also depends on its planned use. Sargent [221] outlines four basic approaches: evaluation by the model development team, evaluation by the end users of the model in collaboration with the model development team, independent verification and validation (IV&V) by a third-party and evaluation using a scoring model. For the types of models discussed in the current work, the first two approaches are most likely to apply. These models are typically used in a process development context, and often the development team creates the model or collaborates closely with the group that creates the model. However in cases where models are included in a regulatory filing the regulatory agency may become involved in assessing the model as an independent third party.
In order to carry out model validation for pharmaceutical unit operations it is necessary to conduct experiments to which model predictions can be compared. These experiments should be carried out in processing equipment comparable to the modeled units. If data-based models have been implemented, the validation data must be distinct from the data used in model development. In some cases it may be appropriate to verify a model against a more detailed simulation such as a discrete or finite element simulation [221]. For instance, in the case of a reduced-order model developed to approximate a detailed simulation using a semi-empirical or black box method such a validation process is appropriate [110].
Model validation can include both qualitative and quantitative assessments. Qualitative evaluations can be done through data visualization and sensitivity analysis. Preliminary evaluation of model predictive ability may be conducted through visualizing the model output in conjunction with experimental data to determine how similar the predicted and observed behaviors are. For instance, parity plots may be used to qualitatively assess the discrepancy between predicted and experimentally obtained model outputs. Alternatively the predicted and observed data may be plotted on the same set of axes for visual comparison [119,221]. Model results may also be qualitatively validated by observing process trends or trajectories and confirming that they correspond with expected process behavior based on prior knowledge of process behavior [113,120]. Sensitivity analysis is the process of attributing variability in model outputs to uncertain model parameters and input variables [224]. Sensitivity metrics, which quantify the sensitivity of a model output to a particular input or collection of inputs, can be used as an indicator of qualitative model agreement with observed process behavior [221,225]. Sensitivity analysis also indicates which model parameters the process outputs are most sensitive to. Thus additional experimental effort can be directed towards fitting these parameters accurately [120,226].
Qualitative model assessment can be helpful in visualizing model performance, but is insufficient to justify the use of a model for quantitative prediction of process outputs. A variety of metrics can be used to quantitatively assess model prediction accuracy [221,227]. In linear regression, the coefficient of determination (R2) can be used as an indicator of model fit. The coefficient of determination describes the strength of the linear relationship between the model inputs and responses, indicating the proportion of variance in the response that can be accounted for by a particular input [228]. In pharmaceutical modeling applications, more frequently used metrics of model performance are based on prediction error, including the sum of square error (SSE), the mean square error of prediction (MSEP) and root-mean square error of prediction (RMSEP) [229,230]. The mean square error of prediction can be defined as the expected squared distance between the predicted and true value of a model output and the root-mean square error of prediction is simply its square root. Unlike correlation-based metrics, such as the coefficient of determination, the mean square error of prediction and root mean square error of prediction are sensitive to the scale of the data [231].
Metrics like the coefficient of determination and the sum of square error and can be calculated based on direct comparison between experimental and simulated data that correspond to the same design and operating configuration. These metrics are relatively easy and computationally inexpensive to calculate and can characterize model accuracy under the specific conditions studied. However they do not provide any information about the error associated with model predictions of points not included in the original dataset. To obtain an estimate of the prediction error (MSEP, RMSEP) alternative validation procedures such as cross-validation and bootstrapping are needed.
Cross-validation, in which a subset of available data is excluded from the model building process and subsequently used to compare with model predictions, can be used to evaluate model performance metrics. In leave-one-out cross validation, a sample of size N is partitioned into a calibration (model-building) set of size N-1 and a validation set of size 1. Cross-validation calculations are then carried out N times with a different sample left out of the calibration set at each iteration. An average error over the N cross validation calculations is then used to indicate model performance. Similar methods wherein a somewhat larger subset (n) of the N samples is used for validation and the calibration set is of size N-n can also be implemented [231]. An advantage of cross-validation is that it can provide a nearly unbiased estimate of prediction error. However, for relatively small sample sizes the variability in the estimates can be high, resulting in unstable error estimation [232,233]. Cross-validation, particularly as it pertains to multivariate regression and latent variable methods, is discussed extensively in the literature [234,235]. An alternative to cross validation is the bootstrap procedure. In this process, a bootstrap size (n) is selected and the input space is randomly sampled n times with replacement. Each set of n samples is referred to as a bootstrap sample. The model evaluation metric of interest (e.g., SSE) is determined for the bootstrap sample. This process is repeated with a large number (B) of bootstrap samples and the resulting metric is averaged over the B bootstraps to provide an estimation of model error. The bootstrapping procedure can provide accurate error estimation for large values of B. In several cases bootstrapping has been shown to provide a more stable estimator of the error than cross validation procedures. This is particularly true for instances with small sample sizes [232,233,236,237].
Quantitative metrics like the mean square error of prediction are straightforward to calculate and useful in determining the validity of model predictions [230]. However the conclusions drawn from these metrics depend on the intended purpose of the model. In early process development a relatively large prediction error may be acceptable since the objective is simply to obtain a general understanding of process behavior. For applications like model predictive control (MPC) large model errors can prove problematic and are generally not acceptable [238,239].
Beyond metrics associated with prediction error, a number of statistical techniques can be used to objectively decide whether a model is sufficiently predictive for a given application. Statistical model validation involves the formulation of a null hypothesis H0 (“the API concentration of a blend predicted from model X is within ±5% of the experimentally determined concentration”) and an alternative that may be the direct opposite of the null hypothesis. Subsequently a statistical metric that will be used to evaluate the hypothesis must be defined. Finally, a rejection criterion for the developed hypothesis should be determined [227,240]. A variety of statistical approaches for model validation are discussed in the literature. The most well-known is classical hypothesis testing in which the t-test, z-test or F-test statistics may be used as a criterion for rejecting or accepting the null hypothesis [227,241]. Bayesian hypothesis testing uses the Bayes factor as a validation metric. The Bayes factor expresses the ratio between the conditional probability of observing the experimentally obtained results given that the null hypothesis is true and the conditional probability of observing the experimental data given that the alternative hypothesis is true [241]. An advantage of using the Bayes factor is that it can be applied to prediction accuracy for a single point or for an entire probability distribution function as described in Ling et al [227]. Statistical tests can also be used to determine confidence intervals for the model predictions relative to experimental observations. These can be used to define the range of model accuracy. In defining confidence intervals care should be taken to set appropriate confidence levels and consider the associated sampling requirements [221].
This section provides only a brief overview of model validation techniques. Model validation is an extensive topic in itself has been the subject of numerous publications. The validation methods discussed in the current work are summarized in Table 4 along with references in which they have been applied to validate various types of models. Any of the aforementioned validation methods can be used alone or in conjunction with other techniques. The choice of model verification and validation approach as well as the validation criteria selected should be consistent with the intended use of the model.
Table 4. Summary of validation methods discussed in the current work, including references in which the proposed validation methodology is implemented.
Table 4. Summary of validation methods discussed in the current work, including references in which the proposed validation methodology is implemented.
Validation MethodMetricsAdvantagesDisadvantagesReferences
Data visualizationQualitativeStraightforward to implement and interpret
-
Relies on visual assessment to determine model quality;
-
Preferable to use in conjunction with quantitative model assessment [242]
[62,113,120,189,192,193]
Sensitivity analysisQualitativeIdentifes important sources of process variability
-
Can require many model runs to calculate some sensitivity indices;
-
Requires knowledge of factors that contribute significantly to sensitivity in practice
[14,225]
Direct comparisonQuantitative (R2, SSE)
-
Straightforward and inexpensive to calculate;
-
Indicates model accuracy for a specific design and operating configuration
Does not provide an estimate of the prediction error [191,242]
Cross-ValidationQuantitative (MSEP, RMSEP)Provides a nearly unbiased estimate of prediction errorCan be an unstable error estimator, particularly for small datasets[42,78,85,110,234,235,243]
BootstrappingQuantitative (MSEP, RMSEP, etc.)
-
Provides a nearly unbiased estimate of prediction error;
-
More stable error estimator than cross-validation, particularly for small sample sizes
Can become computationally expensive as the number of bootstraps increases[244,245]
Hypothesis TestingQualitative result (reject or accept model) based on a Quantitative decision making criteria ( t-test, z-test, or F-test statistic or Bayes factor)
-
Provides objective decision-making criterion;
-
Can be used to provide confidence intervals for model predictions
-
Does not provide an exact indication of model prediction error;
-
Can be difficult to implement relative to other methods
[41,241]

6. Conclusions

Economic and regulatory pressures have led to increased interest in continuous processing for pharmaceutical applications in recent years [4,8,9]. Process systems engineering tools have a significant role to play in the transition from batch to continuous processing. Predictive models are of particular interest, as they can contribute to cost effective and robust process development while enhancing process understanding in a fashion consistent with the ICH Q8 guidance on quality by design [3,10,15]. The use of process systems engineering tools for modeling and simulation of pharmaceutical processes has been increasingly reported in the literature in recent years [13,14,15,95,161,222]. This review has focused specifically on process modeling tools related to the production of tablets via continuous processing. An overview of continuous tableting processes including direct compaction, wet granulation and dry granulation manufacturing routes has been provided, along with a discussion of relevant processing equipment and equation-oriented unit operation models. Computational tools, including those related to detailed process modeling, reduced-order modeling and integrated flowsheet simulation have been discussed in some detail. Finally approaches for model validation have been discussed.
Significant work remains to be done in the area of solids-based process modeling, specifically as it pertains to pharmaceutical manufacturing. Understanding how variability in material properties of APIs and pharmaceutical blends that include these active ingredients affects the quality attributes of drug products remains a challenge in pharmaceutical process development [4,246]. Many existing process models do not sufficiently account for the impact of API physical properties on the performance of unit operations such as powder blending and conveying and the compaction of pharmaceutical blends. DEM simulation has contributed to understanding of the influence of particle-level properties on flow and compaction behavior as described in Section 3.1, but connecting particle level properties in DEM to the bulk powder properties that are typically measured throughout process development remains a challenge. Another area where continued development is needed is in the area of tablet properties prediction. While models exist in the literature to predict tablet properties like hardness and dissolution performance as a function of operating conditions, these are relatively simplistic and may not sufficiently account for the effect of all relevant sources of variability in the tablet manufacturing process [73,161,216,247]. An ongoing challenge in the development of mathematical models for particulate processes is the inherent trade-off between the level of detail included in a model and computational efficiency. Detailed models may be more predictive and provide additional insight into the underlying system behavior, but their computational expense can become prohibitive for process simulation, optimization and MPC applications. Thus the development of reduced order or constitutive models to supplement detailed process models for applications where computational efficiency is important is an area where continued efforts will be needed.

Acknowledgments

The authors would like to acknowledge financial support from Bristol-Myers Squibb as well as from the Engineering Research Center for Structured Organic Particulate Systems at Rutgers University (NSF-0504497, NSF-ECC 0540855).

Conflict of Interest

The authors declare no conflict of interest.

References

  1. Suresh, P.; Basu, P.K. Improving pharmaceutical product development and manufacturing: Impact on cost of drug development and cost of goods sold of pharmaceuticals. J. Pharm. Innov. 2008, 3, 175–187. [Google Scholar] [CrossRef]
  2. Shah, N. Pharmaceutical supply chains: Key issues and strategies for optimisation. Comput. Chem. Eng. 2004, 28, 929–941. [Google Scholar] [CrossRef]
  3. ICH Harmonised Tripartite Guideline: Pharmaceutical Development Q8(r2), Current Step 4 Version. In The International Conference on Harmonisation, Geneva, Switzerland, 2009.
  4. McKenzie, P.; Kiang, S.; Tom, J.; Rubin, A.E.; Futran, M. Can pharmaceutical process development become high tech? AIChE J. 2006, 52, 3990–3994. [Google Scholar] [CrossRef]
  5. Food and Drug Administration, Guidance for Industry, Q8 Pharmaceutical Development; Food and Drug Administration: Silver Spring, MD, USA, 2006.
  6. Reinhardt, U.E. Perspectives on the pharmaceutical industry. Health Aff. (Millwood) 2001, 20, 136–149. [Google Scholar] [CrossRef]
  7. Basu, P.; Joglekar, G.; Rai, S.; Suresh, P.; Vernon, J. Analysis of manufacturing costs in pharmaceutical companies. J. Pharm. Innov. 2008, 3, 30–40. [Google Scholar] [CrossRef]
  8. Plumb, K. Continuous processing in the pharmaceutical industry—Changing the mind set. Chem. Eng. Res. Des. 2005, 83, 730–738. [Google Scholar] [CrossRef]
  9. Buchholz, S. Future manufacturing approaches in the chemical and pharmaceutical industry. Chem. Eng. Process. 2010, 49, 993–995. [Google Scholar] [CrossRef]
  10. Aksu, B.; de Beer, T.; Folestad, S.; Ketolainen, J.; Linden, H.; Lopes, J.A.; de Matas, M.; Oostra, W.; Rantanen, J.; Weimer, M. Strategic funding priorities in the pharmaceutical sciences allied to quality by design (QBD) and process analytical technology (PAT). Eur. J. Pharm. Sci. 2012, 47, 402–405. [Google Scholar] [CrossRef]
  11. Schaber, S.D.; Gerogiorgis, D.I.; Ramachandran, R.; Evans, J.M.B.; Barton, P.I.; Trout, B.L. Economic analysis of integrated continuous and batch pharmaceutical manufacturing: A case study. Ind. Eng. Chem. Res. 2011, 50, 10083–10092. [Google Scholar]
  12. Boukouvala, F.; Muzzio, F.J.; Ierapetritou, M.G. Design space of pharmaceutical processes using data-driven-based methods. J. Pharm. Innov. 2010, 5, 119–137. [Google Scholar] [CrossRef]
  13. Gernaey, K.V.; Gani, R. A model-based systems approach to pharmaceutical product-process design and analysis. Chem. Eng. Sci. 2010, 65, 5757–5769. [Google Scholar] [CrossRef]
  14. Boukouvala, F.; Niotis, V.; Ramachandran, R.; Muzzio, F.J.; Ierapetritou, M.G. An integrated approach for dynamic flowsheet modeling and sensitivity analysis of a continuous tablet manufacturing process. Comput. Chem. Eng. 2012, 42, 30–47. [Google Scholar] [CrossRef]
  15. Gernaey, K.V.; Cervera-Padrell, A.E.; Woodley, J.M. A perspective on PSE in pharmaceutical process development and innovation. Comput. Chem. Eng. 2012, 42, 15–29. [Google Scholar] [CrossRef]
  16. Akkisetty, P.K.; Lee, U.; Reklaitis, G.V.; Venkatasubramanian, V. Population balance model-based hybrid neural network for a pharmaceutical milling process. J. Pharm. Innov. 2010, 5, 161–168. [Google Scholar] [CrossRef]
  17. Ramachandran, R.; Arjunan, J.; Chaudhury, A.; Ierapetritou, M.G. Model-based control-loop performance of a continuous direct compaction process. J. Pharm. Innov. 2011, 6, 249–263. [Google Scholar] [CrossRef]
  18. Singh, R.; Ierapetritou, M.; Ramachandran, R. System-wide hybrid MPC-PID control of a continous pharmaceutical tablet manufacturing process via direct compaction. Eur. J. Pharm. Biopharm. 2013. [Google Scholar] [CrossRef]
  19. Singh, R.; Ierapetritou, M.; Ramachandran, R. An engineering study on the enhanced control and operation of continuous manufacturing of pharmaceutical tablets via roller compaction. Int. J. Pharm. 2012, 438, 307–326. [Google Scholar] [CrossRef]
  20. Gao, Y.J.; Muzzio, F.J.; Ierapetritou, M.G. Optimizing continuous powder mixing processes using periodic section modeling. Chem. Eng. Sci. 2012, 80, 70–80. [Google Scholar] [CrossRef]
  21. Troup, G.M.; Georgakis, C. Process systems engineering tools in the pharmaceutical industry. Comput. Chem. Eng. 2013, 51, 157–171. [Google Scholar] [CrossRef]
  22. Järvinen, M.A.; Paaso, J.; Paavola, M.; Leivisk, K.; Juuti, M.; Muzzio, F.; Järvinen, K. Continuous direct tablet compression: Effects of impeller rotation rate, total feed rate and drug content on the tablet properties and drug release. Drug Dev. Ind. Pharm. 2012. [Google Scholar] [CrossRef]
  23. Anand, A.; Curtis, J.S.; Wassgren, C.R.; Hancock, B.C.; Ketterhagen, W.R. Predicting discharge dynamics from a rectangular hopper using the discrete element method (DEM). Chem. Eng. Sci. 2008, 63, 5821–5830. [Google Scholar] [CrossRef]
  24. Weir, G.J. A mathematical model for dilating, non-cohesive granular flows in steep-walled hoppers. Chem. Eng. Sci. 2004, 59, 149–161. [Google Scholar] [CrossRef]
  25. Gremaud, P.A.; Matthews, J.V.; Schaeffer, D.G. On the computation of steady hopper flows III: Model comparisons. J. Comput. Phys. 2006, 219, 443–454. [Google Scholar] [CrossRef]
  26. Faqih, A.M.N.; Alexander, A.W.; Muzzio, F.J.; Tomassone, M.S. A method for predicting hopper flow characteristics of pharmaceutical powders. Chem. Eng. Sci. 2007, 62, 1536–1542. [Google Scholar] [CrossRef]
  27. Ketterhagen, W.R.; Curtis, J.S.; Wassgren, C.R.; Kong, A.; Narayan, P.J.; Hancock, B.C. Granular segregation in discharging cylindrical hoppers: A discrete element and experimental study. Chem. Eng. Sci. 2007, 62, 6423–6439. [Google Scholar] [CrossRef]
  28. Ketterhagen, W.R.; Hancock, B.C. Optimizing the design of eccentric feed hoppers for tablet presses using DEM. Comput. Chem. Eng. 2010, 34, 1072–1081. [Google Scholar] [CrossRef]
  29. Balevičius, R.; Kačianauskas, R.; Mróz, Z.; Sielamowicz, I. Discrete-particle investigation of friction effect in filling and unsteady/steady discharge in three-dimensional wedge-shaped hopper. Powder Technol. 2008, 187, 159–174. [Google Scholar] [CrossRef]
  30. Engisch, W.E.; Muzzio, F.J. Method for characterization of loss-in-weight feeder equipment. Powder Technol. 2012, 228, 395–403. [Google Scholar] [CrossRef]
  31. Gao, Y.J.; Muzzio, F.; Ierapetritou, M. Characterization of feeder effects on continuous solid mixing using fourier series analysis. AIChE J. 2011, 57, 1144–1153. [Google Scholar] [CrossRef]
  32. Yang, S.; Evans, J.R.G. Metering and dispensing of powder; the quest for new solid freeforming techniques. Powder Technol. 2007, 178, 56–72. [Google Scholar] [CrossRef]
  33. Berthiaux, H.; Marikh, K.; Gatumel, C. Continuous mixing of powder mixtures with pharmaceutical process constraints. Chem. Eng. Process. 2008, 47, 2315–2322. [Google Scholar] [CrossRef]
  34. Pernenkil, L.; Cooney, C.L. A review on the continuous blending of powders. Chem. Eng. Sci. 2006, 61, 720–742. [Google Scholar] [CrossRef]
  35. Portillo, P.M.; Ierapetritou, M.G.; Muzzio, F.J. Characterization of continuous convective powder mixing processes. Powder Technol. 2008, 182, 368–378. [Google Scholar] [CrossRef]
  36. Marikh, K.; Berthiaux, H.; Gatumel, C.; Mizonov, V.; Barantseva, E. Influence of stirrer type on mixture homogeneity in continuous powder mixing: A model case and a pharmaceutical case. Chem. Eng. Res. Des. 2008, 86, 1027–1037. [Google Scholar] [CrossRef]
  37. Portillo, P.M.; Ierapetritou, M.G.; Muzzio, F.J. Effects of rotation rate, mixing angle, and cohesion in two continuous mixers—A statistical approach. Powder Technol. 2009, 194, 217–227. [Google Scholar] [CrossRef]
  38. Vanarase, A.U.; Muzzio, F.J. Effect of operating conditions and design parameters in a continuous powder mixer. Powder Technol. 2011, 208, 26–36. [Google Scholar] [CrossRef]
  39. Portillo, P.M.; Vanarase, A.U.; Ingram, A.; Seville, J.K. Investigation of the effect of impeller rotation rate, powder flowrate, and cohesion on powder flow behavior in a continuous blender using pept. Chem. Eng. Sci. 2010, 65, 5685–5668. [Google Scholar]
  40. Koller, D.M.; Posch, A.; Horl, G.; Voura, C.; Radl, S.; Urbanetz, N.; Fraser, S.D.; Tritthart, W.; Reiter, F.; Schlingmann, M.; et al. Continuous quantitative monitoring of powder mixing dynamics by near-infrared spectroscopy. Powder Technol. 2011, 205, 87–96. [Google Scholar] [CrossRef]
  41. Gao, Y.J.; Vanarase, A.; Muzzio, F.; Ierapetritou, M. Characterizing continuous powder mixing using residence time distribution. Chem. Eng. Sci. 2011, 66, 417–425. [Google Scholar] [CrossRef]
  42. Jarvinen, K.; Hoehe, W.; Jarvinen, M.; Poutiainen, S.; Juuti, M.; Borchert, S. In-line monitoring of the drug content of powder mixtures and tablets by near-infrared spectroscopy during the continuous direct compression tableting process. Eur. J. Pharm. Sci. 2013, 48, 680–688. [Google Scholar] [CrossRef]
  43. Vanarase, A.U.; Alcala, M.; Rozo, J.I.J.; Muzzio, F.J.; Romanach, R.J. Real-time monitoring of drug concentration in a continuous powder mixing process using NIR spectroscopy. Chem. Eng. Sci. 2010, 65, 5728–5733. [Google Scholar] [CrossRef]
  44. Martínez, L.; Peinado, A.; Liesum, L.; Betz, G. Use of near-infrared spectroscopy to quantify drug content on a continuous blending process: Influence of mass flow and rotation speed variations. Eur. J. Pharm. Biopharm. 2013, 84, 606–615. [Google Scholar] [CrossRef]
  45. Marikh, K.; Berthiaux, H.; Mizonov, V.; Barantsev, E. Experimental study of the stirring conditions taking place in a pilot plant continuous mixer of particulate solids. Powder Technol. 2005, 157, 138–143. [Google Scholar] [CrossRef]
  46. Faure, A.; York, P.; Rowe, R.C. Process control and scale-up of pharmaceutical wet granulation processes: A review. Eur. J. Pharm. Biopharm. 2001, 52, 269–277. [Google Scholar] [CrossRef]
  47. Lee, K.T.; Ingram, A.; Rowson, N.A. Comparison of granule properties produced using twin screw extruder and high shear mixer: A step towards understanding the mechanism of twin screw wet granulation. Powder Technol. 2013, 238, 91–98. [Google Scholar] [CrossRef]
  48. Iveson, S.M.; Litster, J.D.; Hapgood, K.; Ennis, B.J. Nucleation, growth and breakage phenomena in agitated wet granulation processes: A review. Powder Technol. 2001, 117, 3–39. [Google Scholar] [CrossRef]
  49. Betz, G.; Junker-Burgin, P.; Leuenberger, H. Batch and continuous processing in the production of pharmaceutical granules. Pharm. Dev. Technol. 2003, 8, 289–297. [Google Scholar] [CrossRef]
  50. Vervaet, C.; Remon, J.P. Continuous granulation in the pharmaceutical industry. Chem. Eng. Sci. 2005, 60, 3949–3957. [Google Scholar] [CrossRef]
  51. Tu, W.D.; Ingram, A.; Seville, J. Regime map development for continuous twin screw granulation. Chem. Eng. Sci. 2013, 87, 315–326. [Google Scholar] [CrossRef]
  52. Dhenge, R.M.; Cartwright, J.J.; Hounslow, M.J.; Salman, A.D. Twin screw wet granulation: Effects of properties of granulation liquid. Powder Technol. 2012, 229, 126–136. [Google Scholar] [CrossRef]
  53. Cartwright, J.J.; Robertson, J.; D’Haene, D.; Burke, M.D.; Hennenkamp, J.R. Twin screw wet granulation: Loss in weight feeding of a poorly flowing active pharmaceutical ingredient. Powder Technol. 2013, 238, 116–121. [Google Scholar] [CrossRef]
  54. Vercruysse, J.; Cordoba Diaz, D.; Peeters, E.; Fonteyne, M.; Delaet, U.; van Assche, I.; de Beer, T.; Remon, J.P.; Vervaet, C. Continuous twin screw granulation: Influence of process variables on granule and tablet quality. Eur. J. Pharm. Biopharm. 2012, 82, 205–211. [Google Scholar] [CrossRef] [Green Version]
  55. Paltzer, S. Drying of wet agglomerates in a continuous fluid bed: Influence of residence time, air temperature and air-flowrate on the drying kinetics and the amount of oversize particles. Chem. Eng. Sci. 2007, 62, 463. [Google Scholar] [CrossRef]
  56. Kleinebudde, P. Roll compaction/dry granulation: Pharmaceutical applications. Eur. J. Pharm. Biopharm. 2004, 58, 317–326. [Google Scholar] [CrossRef]
  57. Yu, S.; Gururajan, B.; Reynolds, G.; Roberts, R.; Adams, M.J.; Wu, C.Y. A comparative study of roll compaction of free-flowing and cohesive pharmaceutical powders. Int. J. Pharm. 2012, 428, 39–47. [Google Scholar] [CrossRef]
  58. Lecompte, T.; Doremus, P.; Thomas, G.; Perier-Camby, L.; Le Thiesse, J.-C.; Masteau, J.-C.; Debove, L. Dry granulation of organic powders—Dependence of pressure 2d-distribution on different process parameters. Chem. Eng. Sci. 2005, 60, 3933–3940. [Google Scholar] [CrossRef]
  59. Boswell, S.; Smith, G. Improving solid dosage forms with dry granulation. Pharm. Technol. Eur. 2011, 23, 31. [Google Scholar]
  60. Sharma, A.; Sharma, A.; Chauhan, C.S. Roller compaction: Imperative process for tablet manufacturing: A review. Int. J. Pharm. Res. Dev. 2012, 4, 40–47. [Google Scholar]
  61. Reynolds, G.; Ingale, R.; Roberts, R.; Kothari, S.; Gururajan, B. Practical application of roller compaction process modeling. Comput. Chem. Eng. 2010, 34, 1049–1057. [Google Scholar] [CrossRef]
  62. Bindhumadhavan, G.; Seville, J.P.K.; Adams, N.; Greenwood, R.W.; Fitzpatrick, S. Roll compaction of a pharmaceutical excipient: Experimental validation of rolling theory for granular solids. Chem. Eng. Sci. 2005, 60, 3891–3897. [Google Scholar] [CrossRef]
  63. Dec, R.T.; Zavaliangos, A.; Cunningham, J.C. Comparison of various modeling methods for analysis of powder compaction in roller press. Powder Technol. 2003, 130, 265–271. [Google Scholar] [CrossRef]
  64. Miguélez-Morán, A.M.; Wu, C.-Y.; Seville, J.P.K. The effect of lubrication on density distributions of roller compacted ribbons. Int. J. Pharm. 2008, 362, 52–59. [Google Scholar] [CrossRef]
  65. Nakach, M.; Authelin, J.R.; Chamayou, A.; Dodds, J. Comparison of various milling technologies for grinding pharmaceutical powders. Int. J. Miner. Process. 2004, 74, S173–S181. [Google Scholar] [CrossRef]
  66. Reynolds, G.K. Modelling of pharmaceutical granule size reduction in a conical screen mill. Chem. Eng. J. 2010, 164, 383–392. [Google Scholar] [CrossRef]
  67. Verheezen, J.J.; van der Voort Maarschalk, K.; Faassen, F.; Vromans, H. Milling of agglomerates in an impact mill. Int. J. Pharm. 2004, 278, 165–172. [Google Scholar] [CrossRef]
  68. Vendola, T.A.; Hancock, B.C. The effect of mill type on two dry-granulated placebo formulations. Pharm. Technol. 2008, 32, 72–86. [Google Scholar]
  69. Samanta, A.K.; Ng, K.Y.; Heng, P.W. Cone milling of compacted flakes: Process parameter selection by adopting the minimal fines approach. Int. J. Pharm. 2012, 422, 17–23. [Google Scholar] [CrossRef]
  70. Motzi, J.J.; Anderson, N.R. The quantitative evaluation of a granulation milling process ii. Effect of output screen size, mill speed and impeller shape. Drug Dev. Ind. Pharm. 1984, 10, 713–728. [Google Scholar] [CrossRef]
  71. Inghelbracht, S.; Remon, J.P. Reducing dust and improving granule and tablet quality in the roller compaction process. Int. J. Pharm. 1998, 171, 195–206. [Google Scholar] [CrossRef]
  72. Mehrotra, A.; Chaudhuri, B.; Faqih, A.; Tomassone, M.S.; Muzzio, F.J. A modeling approach for understanding effects of powder flow properties on tablet weight variability. Powder Technol. 2009, 188, 295–300. [Google Scholar] [CrossRef]
  73. Kuentz, M.; Lunenberger, H. A new model for the hardness of a compacted particle systems, applied to tablets of pharmaceutical polymers. Powder Technol. 2000, 111, 145–153. [Google Scholar] [CrossRef]
  74. Gentis, N.D.; Betz, G. Compressibility of binary powder formulations: Investigation and evaluation with compaction equations. J. Pharm. Sci. 2012, 101, 777–793. [Google Scholar] [CrossRef]
  75. Patel, S.; Kaushal, A.M.; Bansal, A.K. Effect of particle size and compression force on compaction behavior and derived mathematical parameters of compressibility. Pharm. Res. 2007, 24, 111–124. [Google Scholar] [CrossRef]
  76. Wu, C.Y.; Ruddy, O.M.; Bentham, A.C.; Hancock, B.C.; Best, S.M.; Elliott, J.A. Modelling the mechanical behaviour of pharmaceutical powders during compaction. Powder Technol. 2005, 152, 107–117. [Google Scholar] [CrossRef]
  77. Podczeck, F. Methods for the practical determination of the mechanical strength of tablets—From empiricism to science. Int. J. Pharm. 2012, 436, 214–232. [Google Scholar] [CrossRef]
  78. Corredor, C.C.; Bu, D.; Both, D. Comparison of near infrared and microwave resonance sensors for at-line moisture determination in powders and tablets. Anal. Chim. Acta 2011, 696, 84–93. [Google Scholar] [CrossRef]
  79. Zavaliangos, A.; Galen, S.; Cunningham, J.; Winstead, D. Temperature evolution during compaction of pharmaceutical powders. J. Pharm. Sci. 2008, 97, 3291–3304. [Google Scholar] [CrossRef]
  80. Onuki, Y.; Kawai, S.; Arai, H.; Maeda, J.; Takagaki, K.; Takayama, K. Contribution of the physicochemical properties of active pharmaceutical ingredients to tablet properties identified by ensemble artificial neural networks and Kohonen’s self-organizing maps. J. Pharm. Sci. 2012, 101, 2372–2381. [Google Scholar] [CrossRef]
  81. Govedarica, B.; Ilic, I.; Sibanc, R.; Dreu, R.; Srcic, S. The use of single particle mechanical properties for predicting the compressibility of pharmaceutical materials. Powder Technol. 2012, 225, 43–51. [Google Scholar] [CrossRef]
  82. Wang, J.; Wen, H.; Desai, D. Lubrication in tablet formulations. Eur. J. Pharm. Biopharm. 2010, 75, 1–15. [Google Scholar] [CrossRef]
  83. Jackson, S.; Sinka, I.C.; Cocks, A.C. The effect of suction during die fill on a rotary tablet press. Eur. J. Pharm. Biopharm. 2007, 65, 253–256. [Google Scholar] [CrossRef]
  84. Hancock, B.C.; Ketterhagen, W.R. Discrete element method (DEM) simulations of stratified sampling during solid dosage form manufacturing. Int. J. Pharm. 2011, 418, 265–272. [Google Scholar] [CrossRef]
  85. Otsuka, M.; Yamane, I. Prediction of tablet properties based on near infrared spectra of raw mixed powders by chemometrics: Scale-up factor of blending and tableting processes. J. Pharm. Sci. 2009, 98, 4296–4305. [Google Scholar] [CrossRef]
  86. Yu, A.B. Discrete element method—An effective way for particle scale research of particulate matter. Eng. Comput. 2004, 21, 205–214. [Google Scholar] [CrossRef]
  87. Zhu, H.P.; Zhou, Z.Y.; Yang, R.Y.; Yu, A.B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62, 3378–3396. [Google Scholar] [CrossRef]
  88. Zhu, H.P.; Zhou, Z.Y.; Yang, R.Y.; Yu, A.B. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 2008, 63, 5728–5770. [Google Scholar] [CrossRef]
  89. Persson, A.S.; Frenning, G. An experimental evaluation of the accuracy to simulate granule bed compression using the discrete element method. Powder Technol. 2012, 219, 249–256. [Google Scholar] [CrossRef]
  90. Siiria, S.M.; Antikainen, O.; Heinamaki, J.; Yliruusi, J. 3d simulation of internal tablet strength during tableting. AAPS PharmSciTech 2011, 12, 593–603. [Google Scholar] [CrossRef]
  91. Jerier, J.F.; Hathong, B.; Richefeu, V.; Chareyre, B.; Imbault, D.; Donze, F.V.; Doremus, P. Study of cold powder compaction by using the discrete element method. Powder Technol. 2011, 208, 537–541. [Google Scholar] [CrossRef]
  92. Yi, L.Y.; Dong, K.J.; Zou, R.P.; Yu, A.B. Coordination number of the packing of ternary mixtures of spheres: DEM simulations versus measurements. Ind. Eng. Chem. Res. 2011, 50, 8773–8785. [Google Scholar] [CrossRef]
  93. McCarthy, J.J.; Jasti, V.; Marinack, M.; Higgs, C.F. Quantitative validation of the discrete element method using an annular shear cell. Powder Technol. 2010, 203, 70–77. [Google Scholar] [CrossRef]
  94. Wu, C.Y.; Cocks, A.C.F. Numerical and experimental investigations of the flow of powder into a confined space. Mech. Mater. 2006, 38, 304–324. [Google Scholar] [CrossRef]
  95. Ketterhagen, W.R.; Ende, M.T.A.; Hancock, B.C. Process modeling in the pharmaceutical industry using the discrete element method. J. Pharm. Sci. 2009, 98, 442–470. [Google Scholar] [CrossRef]
  96. Wu, C.Y. DEM simulations of die filling during pharmaceutical tabletting. Particuology 2008, 6, 412–418. [Google Scholar] [CrossRef]
  97. Wu, C.-Y.; Guo, Y. Numerical modelling of suction filling using DEM/CFD. Chem. Eng. Sci. 2012, 73, 231–238. [Google Scholar] [CrossRef]
  98. Guo, Y.; Kafui, K.D.; Wu, C.Y.; Thornton, C.; Seville, J.P.K. A coupled DEM/CFD analysis of the effect of air on powder flow during die filling. AIChE J. 2008, 55, 49–62. [Google Scholar]
  99. Guo, Y.; Wu, C.-Y.; Kafui, K.D.; Thornton, C. Numerical analysis of density-induced segregation during die filling. Powder Technol. 2009, 197, 111–119. [Google Scholar]
  100. Guo, Y.; Wu, C.Y.; Kafui, K.D.; Thornton, C. 3D DEM/CFD analysis of size-induced segregation during die filling. Powder Technol. 2011, 206, 177–188. [Google Scholar] [CrossRef]
  101. Gethin, D.T.; Yang, X.S.; Lewis, R.W. A two dimensional combined discrete and finite element scheme for simulating the flow and compaction of systems comprising irregular particles. Comput. Methods Appl. Mech. Eng. 2006, 195, 5552–5565. [Google Scholar] [CrossRef]
  102. Frenning, G. Compression mechanics of granule beds: A combined finite/discrete element study. Chem. Eng. Sci. 2010, 65, 2464–2471. [Google Scholar] [CrossRef]
  103. Nwose, E.N.; Pei, C.L.; Wu, C.Y. Modelling die filling with charged particles using DEM/CFD. Particuology 2012, 10, 229–235. [Google Scholar] [CrossRef]
  104. Bertrand, F.; Leclaire, L.A.; Levecque, G. DEM-based models for the mixing of granular materials. Chem. Eng. Sci. 2005, 60, 2517–2531. [Google Scholar] [CrossRef]
  105. Marigo, M.; Davies, M.; Leadbeater, T.; Cairns, D.L.; Ingram, A.; Stitt, E.H. Application of positron emission particle tracking (PEPT) to validation a discrete element method (DEM) model of granular flow and mixing in the turbula mixer. Int. J. Pharm. 2013, 66, 1811–1824. [Google Scholar]
  106. Remy, B.; Khinast, J.G.; Glasser, B.J. Polydisperse granular flows in a bladed mixer: Experiments and simulations of cohesionless spheres. Chem. Eng. Sci. 2011, 66, 1811–1824. [Google Scholar] [CrossRef]
  107. Remy, B.; Canty, T.M.; Khinast, J.G.; Glasser, B.J. Experiments and simulations of cohesionless particles with varying roughness in a bladed mixer. Chem. Eng. Sci. 2010, 65, 4557–4571. [Google Scholar] [CrossRef]
  108. Sarkar, A.; Wassgren, C.R. Comparison of flow microdynamics for a continuous granular mixer with predictions from periodic slice DEM simulations. Powder Technol. 2012, 221, 325–336. [Google Scholar] [CrossRef]
  109. Gao, Y.J.; Muzzio, F.J.; Ierapetritou, M.G. A review of the residence time distribution (RTD) applications in solid unit operations. Powder Technol. 2012, 228, 416–423. [Google Scholar] [CrossRef]
  110. Boukouvala, F.; Gao, Y.J.; Muzzio, F.; Ierapetritou, M.G. Reduced-order discrete element method modeling. Chem. Eng. Sci. 2013, 95, 12–26. [Google Scholar] [CrossRef]
  111. Griffin, D.W.; Mellichamp, D.A.; Doherty, M.F. Reducing the mean size of api crystals by continuous manufacturing with product classification and recycle. Chem. Eng. Sci. 2010, 65, 5770–5780. [Google Scholar] [CrossRef]
  112. Mortier, S.T.F.C.; Gernaey, K.V.; de Beer, T.; Nopens, I. Development of a population balance model of a pharmaceutical drying process and testing of solution methods. Comput. Chem. Eng. 2013, 50, 39–53. [Google Scholar] [CrossRef]
  113. Barasso, D.; Ramachandran, R. A comparison of model order reduction techniques for a four-dimensional population balance model describing multi-component wet granulation process. Chem. Eng. Sci. 2012, 80, 380–392. [Google Scholar] [CrossRef]
  114. Langham, Z.A.; Booth, J.; Hughes, L.P.; Reynolds, G.K.; Wren, S.A.C. Mechanistic insights into the dissolution of spray-dried amorphous solid dispersions. J. Pharm. Sci. 2012, 101, 2798–2810. [Google Scholar] [CrossRef]
  115. Sen, M.; Ramachandran, R. A multi-dimensional population balance model approach to continuous powder mixing processes. Adv. Powder Technol. 2013, 24, 51–59. [Google Scholar]
  116. Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems Engineering; Academic Press: London, UK, 2000. [Google Scholar]
  117. Majumbder, A.; Kariwala, V.; Ansumali, S.; Rajendran, A. Lattice boltzmann method for population balance equations with simultaneous growth, nucleation, aggregation and breakage. Chem. Eng. Sci. 2012, 65, 4884–4893. [Google Scholar]
  118. Ramachandran, R.; Barton, P.I. Effective parameter estimation within a multi-dimensional population balance model framework. Chem. Eng. Sci. 2010, 65, 4884–4893. [Google Scholar] [CrossRef]
  119. Sen, M.; Singh, R.; Vanarase, A.; John, J.; Ramachandran, R. Multi-dimensional population balance modeling and experimental validation of continuous powder mixing processes. Chem. Eng. Sci. 2012, 80, 349–360. [Google Scholar] [CrossRef]
  120. Ramachandran, R.; Immanuel, C.D.; Stepanek, F.; Litster, J.D.; Doyle, F.J. A mechanistic model for breakage in population balances of granulation: Theoretical kernel development and experimental validation. Chem. Eng. Res. Des. 2009, 87, 598–614. [Google Scholar] [CrossRef]
  121. Immanuel, C.D.; Doyle, F.J. Computationally efficient solution of population balance models incorporating nucleation, growth and coagulation: Application to emulsion polymerization. Chem. Eng. Sci. 2003, 58, 3681–3698. [Google Scholar] [CrossRef]
  122. Immanuel, C.D.; Doyle, F.J. Solution technique for a multi-dimensional population balance model describing granulation processes. Powder Technol. 2005, 156, 213–225. [Google Scholar] [CrossRef]
  123. Boukouvala, F.; Dubey, A.; Vanarase, A.; Ramachandran, R.; Muzzio, F.J.; Ierapetritou, M. Computational approaches for studying the granular dynamics of continuous blending processes, 2—Population balance and data-based methods. Macromol. Mater. Eng. 2012, 297, 9–19. [Google Scholar] [CrossRef]
  124. Lang, Y.D.; Zitney, S.E.; Biegler, L.T. Optimization of IGCC processes with reduced order CFD models. Comput. Chem. Eng. 2011, 35, 1705–1717. [Google Scholar] [CrossRef]
  125. Zitney, S.E. Process/equipment co-simulation for design and analysis of advanced energy systems. Comput. Chem. Eng. 2010, 34, 1532–1542. [Google Scholar] [CrossRef]
  126. Brenner, T.A.; Fontenot, R.L.; Cizmas, P.G.A.; O’Brien, T.J.; Breault, R.W. A reduced-order model for heat transfer in multiphase flow and practical aspects of the proper orthogonal decomposition. Comput. Chem. Eng. 2012, 43, 68–80. [Google Scholar] [CrossRef]
  127. Liberge, E.; Hamdouni, A. Reduced order modelling method via proper orthogonal decomposition (POD) for flow around an oscillating cylinder. J. Fluids Struct. 2010, 26, 292–311. [Google Scholar] [CrossRef]
  128. Lang, Y.D.; Malacina, A.; Biegler, L.T.; Munteanu, S.; Madsen, J.I.; Zitney, S.E. Reduced order model based on principal component analysis for process simulation and optimization. Energy Fuels 2009, 23, 1695–1706. [Google Scholar] [CrossRef]
  129. Boukouvala, F.; Muzzio, F.J.; Ierapetritou, M.G. Predictive modeling of pharmaceutical processes with missing and noisy data. AIChE J. 2010, 56, 2860–2872. [Google Scholar] [CrossRef]
  130. Box, G.E.P. On the experimental attainment of optimum conditions. J. R. Stat. Soc. Ser. B 1951, 13, 1–35. [Google Scholar]
  131. Jia, Z.; Davis, E.; Muzzio, F.J.; Ierapetritou, M.G. Predictive modeling for pharmaceutical processes using kriging and response surface. J. Pharm. Innov. 2009, 4, 174–186. [Google Scholar] [CrossRef]
  132. Khuri, A.I.; Mukhopadhyay, S. Response surface methodology. WIREs Comput. Stat. 2010, 2, 128–149. [Google Scholar] [CrossRef]
  133. Myers, R.H.; Montgomery, D.C. Response Surface Methodology Process and Product Optimization Using Designed Experiments; John Wiley & Sons, Inc.: New York, NY, USA, 2002. [Google Scholar]
  134. Ranjbarian, S.; Farhadi, F. Evaluation of the effects of process parameters on granule mean size in a conical high shear granulator using response surface methodology. Powder Technol. 2013, 237, 186–190. [Google Scholar] [CrossRef]
  135. Metheron, G. Principles of geostatistics. Econ. Geol. 1963, 58, 1246–1266. [Google Scholar] [CrossRef]
  136. Calder, C.A.; Cressie, N. Kriging and Variogram Models. In International Encyclopedia of Human Geography; Kitchin, R., Thrift, N., Eds.; Elsevier Science: Oxford, UK, 2009; Volume 1, pp. 49–55. [Google Scholar]
  137. Krige, D.G. A Statistical Approach to Some Mine Valuations and Allied Problems at the Witwatersrand. Master Thesis, University of Witwatersrand, Johannesburg, North Africa, 1951. [Google Scholar]
  138. Boukouvala, F.; Muzzio, F.J.; Ierapetritou, M.G. Dynamic data-driven modeling of pharmaceutical processes. Ind. Eng. Chem. Res. 2011, 50, 6743–6754. [Google Scholar] [CrossRef]
  139. Boukouvala, F.; Ierapetritou, M.G. Feasibility analysis of black-box processes using adaptive sampling kriging based method. Comput. Chem. Eng. 2012, 36, 358–368. [Google Scholar] [CrossRef]
  140. Kleijnen, J.P.C. Kriging metamodeling in simulation: A review. Eur. J. Oper. Res. 2009, 192, 707–716. [Google Scholar] [CrossRef]
  141. Huang, J.; Kaul, G.; Cai, C.; Chatlapalli, R.; Hernandez-Abad, P.; Gosh, K.; Nagi, A. Quality by design case study: An integrated multivariate approach to drug product and process development. Int. J. Pharm. 2009, 382, 23–32. [Google Scholar] [CrossRef]
  142. Braumann, A.; Kraft, M.; Mort, P.R. Parameter estimation in a multidimensional granulation model. Powder Technol. 2010, 197, 196–210. [Google Scholar] [CrossRef]
  143. Li, G.Y.; Wang, S.W.; Rabitz, H. Practical approaches to construct RS-HDMR component functions. J. Phys. Chem. A 2002, 106, 8721–8733. [Google Scholar] [CrossRef]
  144. Li, G.Y.; Rabitz, H.; Hu, J.S.; Chen, Z.; Ju, Y.G. Regularized random-sampling high dimensional model representation (RS-HDMR). J. Math. Chem. 2008, 43, 1207–1232. [Google Scholar] [CrossRef]
  145. Li, G.; Rosenthal, C.; Rabitz, H. High dimensional model representations. J. Phys. Chem. A 2001, 105, 7765–7777. [Google Scholar] [CrossRef]
  146. Ziehn, T.; Tomlin, A.S. Global sensitivity analysis of a 3d street canyon model—Part I: The development of high dimensional model representations. Atmos. Environ. 2008, 42, 1857–1873. [Google Scholar] [CrossRef]
  147. Li, G.; Wang, S.-W.; Rabitz, H. High dimensional model representations (HDMR): Concepts and applications. Available online: http:// www.ima.umn.edu/talks/workshops/3-15-19.2000/li/hdmr.pdf (accessed on 10 April 2013).
  148. Li, G.; Rabitz, H.; Wang, S.-W.; Georgopoulos, P.G. Correlation method for variance reduction of Monte Carlo integration in RS-HDMR. J. Comput. Chem. 2003, 24, 277–283. [Google Scholar] [CrossRef]
  149. Li, G.; Rabitz, H. Ratio control variate method for efficiently determining high-dimensional model representations. J. Comput. Chem. 2006, 27, 1112–1118. [Google Scholar] [CrossRef]
  150. Banerjee, I.; Pal, S.; Maiti, S. Computationally efficient black-box modeling for feasibility analysis. Comput. Chem. Eng. 2010, 34, 1515–1521. [Google Scholar] [CrossRef]
  151. Banarjee, I.; Ierapetritou, M.G. Model independent parametric decision making. Ann. Oper. Res. 2004, 132, 135. [Google Scholar] [CrossRef]
  152. Li, G.Y.; Wang, S.W.; Rabitz, H.; Wang, S.Y.; Jaffe, P. Global uncertainty assessments by high dimensional model representations (HDMR). Chem. Eng. Sci. 2002, 57, 4445–4460. [Google Scholar] [CrossRef]
  153. Ziehn, T.; Tomlin, A.S. A global sensitivity study of sulfur chemistry in a premixed methane flame model using HDMR. Int. J. Chem. Kinet. 2008, 40, 742–753. [Google Scholar] [CrossRef]
  154. Ziehn, T.; Tomlin, A.S. GUI–HDMR—A software tool for global sensitivity analysis of complex models. Environ. Model. Softw. 2009, 24, 775–785. [Google Scholar] [CrossRef]
  155. Basheer, L.A.; Hajmeer, M. Artificial neural networks: Fundamentals, computing, design, and application. J. Microbiol. Methods 2000, 43, 3–31. [Google Scholar] [CrossRef]
  156. Agatonovic-Kustrin, S.; Beresford, R. Basic concepts of artificial neural network (ANN) modeling and its application in pharmaceutical research. J. Pharm. Biomed. Anal. 2000, 22, 717–727. [Google Scholar] [CrossRef]
  157. De la Fuente, R.L.N.; Garcia-Munoz, S.; Biegler, L.T. An efficient nonlinear programming strategy for PCA models with incomplete data sets. J. Chemom. 2010, 24, 301–311. [Google Scholar]
  158. Walczak, B.; Massart, D.L. Dealing with missing data: Part I. Chemom. Intell. Lab. Syst. 2001, 58, 15–27. [Google Scholar] [CrossRef]
  159. Walczak, B.; Massart, D.L. Dealing with missing data: Part II. Chemom. Intell. Lab. Syst. 2001, 58, 28–42. [Google Scholar]
  160. Noonan, R.; Wold, H. NIPALS path modelling with latent variables. Scand. J. Educ. Res. 1977, 21, 33–61. [Google Scholar] [CrossRef]
  161. Boukouvala, F.; Chaudhury, A.; Sen, M.; Zhou, R.J.; Mioduszewski, L.; Ierapetritou, M.G.; Ramachandran, R. Computer-aided flowsheet simulation of a pharmaceutical tablet manufacturing process incorporating wet granulation. J. Pharm. Innov. 2013, 8, 11–27. [Google Scholar] [CrossRef]
  162. Sen, M.; Chaudhury, A.; Singh, R.; John, J.; Ramachandran, R. Multi-scale flowsheet simulation of an integrated continuous purification-downstream pharmaceutical manufacturing process. Int. J. Pharm. 2013, 445, 29–38. [Google Scholar] [CrossRef]
  163. Papavasileiou, V.; Koulouris, A.; Siletti, C.; Petrides, D. Optimize manufacturing of pharmaceutical products with process simulation and production scheduling tools. Chem. Eng. Res. Des. 2007, 85, 1086–1097. [Google Scholar] [CrossRef]
  164. gPROMS ModelBuilder Guide Release v3.6.; Process Systems Enterprise, Ltd.: London, UK, October; 2012.
  165. Model Developer Guide Release v3.6.; Process Systems Enterprise, Ltd.: London, UK, October; 2012.
  166. Model Validation Guide Release v3.6.; Process Systems Enterprise, Ltd.: London, UK, October; 2012.
  167. gPROMS Advanced User Guide Release 2.3.; Process System Enterprise, Ltd.: London, UK, February; 2004.
  168. Minceva, M.; Rodrigues, A.E. Two-level optimization of an existing SMB for p-xylene separation. Comput. Chem. Eng. 2005, 29, 2215–2228. [Google Scholar] [CrossRef]
  169. Asteasuain, M.; Brandolin, A. Modeling and optimization of a high-pressure ethylene polymerization reactor using gPROMS. Comput. Chem. Eng. 2008, 32, 396–408. [Google Scholar] [CrossRef]
  170. Nowee, S.M.; Abbas, A.; Romagnoli, J.A. Optimization in seeded cooling crystallization: A parameter estimation and dynamic optimization study. Chem. Eng. Process. 2007, 46, 1096–1106. [Google Scholar] [CrossRef]
  171. Nowee, S.M.; Abbas, A.; Romagnoli, J.A. Antisolvent crystallization: Model identification, experimental validation and dynamic simulation. Chem. Eng. Sci. 2008, 63, 5457–5467. [Google Scholar] [CrossRef]
  172. Bermingham, S.K.; Verheijen, P.J.T.; Kramer, H.J.M. Optimal design of solution crystallization processes with rigorous models. Chem. Eng. Res. Des. 2003, 81, 894–903. [Google Scholar]
  173. Beck, R. Aspen Plus v8.0 Solids Modeling: A Brief Introduction; Aspen Technology, Inc.: Burlington, MA, USA, 2012. [Google Scholar]
  174. Levine, J. Jump Start: Solids Process Modeling in Aspen Plus® v8; Aspen Technology, Inc.: Burlington, MA, USA, 2012. [Google Scholar]
  175. Chemical engineering: Software. Available online: https://uwaterloo.ca/chemical-engineering/resources-services/computing-facilities/software#AspenCustomModeler (accessed on 5 January 2013).
  176. Aspen Plus, Getting Started Modeling Processes with Solids; Version v7.1; Aspen Technology, Inc.: Burlington, MA, USA, 2009.
  177. Wei, H.-Y. Computer-aided design and scale-up of crystallization processes: Integrating approaches and case studies. Chem. Eng. Res. Des. 2010, 88, 1377–1380. [Google Scholar] [CrossRef]
  178. Lau, S.-Y.; Gonawan, F.N.; Bhatia, S.; Kamaruddin, A.H.; Uzir, M.H. Conceptual design and simulation of a plant for the production of high purity (S)-ibuprofen acid using innovative enzymatic membrane technology. Chem. Eng. J. 2011, 166, 726–737. [Google Scholar] [CrossRef]
  179. Brown, R.L. Minimum energy theorem for flow of dry granular materials through apertures. Nature 1961, 191, 458–461. [Google Scholar] [CrossRef]
  180. Savage, S.B. The mass flow of granular materials derived from coupled velocity-stress fields. Br. J. Appl. Phys. 1965, 16, 1885–1888. [Google Scholar] [CrossRef]
  181. Savage, S.B.; Sayed, M.S. Gravity flow of coarse cohesionless granular materials in conical hoppers. J. Appl. Math. Phys. 1981, 2, 125–143. [Google Scholar]
  182. Brennen, C.; Pearce, J.C. Granular material flow in two dimensional hoppers. J. Appl. Mech. 1978, 45, 43–50. [Google Scholar] [CrossRef]
  183. Nguyen, T.V.; Brennen, C.; Sabersky, R.H. Gravity flow of granular materials in conical hoppers. J. Appl. Mech. 1979, 46, 529–535. [Google Scholar] [CrossRef]
  184. Sun, J.; Sundaresan, S. A constitutive model with microstructure evolution for flow of rate-independent granular materials. J. Fluid Mech. 2011, 682, 590–616. [Google Scholar] [CrossRef]
  185. Weir, G.J. Sound speed and attenuation in dense, non-cohesive air-granulator systems. Chem. Eng. Sci. 2001, 56, 3699–3717. [Google Scholar] [CrossRef]
  186. Sun, J.; Sundaresan, S. Radial hopper flow prediction using a constitutive model with microstructure evolution. Powder Technol. 2013, 242, 81–85. [Google Scholar] [CrossRef]
  187. Danckwerts, P.V. Continuous flow systems: Distribution of residence times. Chem. Eng. Sci. 1953, 2, 1–13. [Google Scholar] [CrossRef]
  188. Gao, Y.; Ierapetritou, M.; Muzzio, F. Periodic section modeling of convective continuous powder mixing processes. AIChE J. 2012, 58, 69–78. [Google Scholar] [CrossRef]
  189. Bouffard, J.; Bertrand, F.; Chaouki, J. A multiscale model for the simulation of granulation in rotor-based equipment. Chem. Eng. Sci. 2012, 81, 106–117. [Google Scholar] [CrossRef]
  190. Cameron, I.T.; Wang, F.Y.; Immanuel, C.D.; Stepanek, F. Process systems modelling and applications in granulation: A review. Chem. Eng. Sci. 2005, 80, 3723–3750. [Google Scholar]
  191. Hounslow, M.J.; Pearson, J.M.K.; Instone, T. Tracer studies of high-shear granulation: II. Population balance modeling. AIChE J. 2001, 47, 1984–1999. [Google Scholar] [CrossRef]
  192. Biggs, C.A.; Sanders, C.; Scott, A.C.; Willemse, A.W.; Hoffman, A.C.; Instone, T.; Salman, A.D.; Hounslow, M.J. Coupling granule properties and granulation rates in high-shear granulation. Powder Technol. 2003, 130, 162–168. [Google Scholar] [CrossRef]
  193. Liu, L.X.; Litster, J.D. Population balance modelling of granulation with a physically based coalescence kernel. Chem. Eng. Sci. 2002, 57, 2183–2191. [Google Scholar] [CrossRef]
  194. Gannt, J.A.; Cameron, I.T.; Litster, J.D.; Gatzke, E.P. Determination of coalescence kernels for high-shear granulation using DEM simulations. Powder Technol. 2006, 170, 53–63. [Google Scholar] [CrossRef]
  195. Poon, J.M.-H.; Immanuel, C.D.; Doyle, F.J.I.; Litster, J.D. A three-dimensional population balance model of granulation with a mechanistic representation of the nucleation and aggregation phenomena. Chem. Eng. Sci. 2008, 63, 1315–1329. [Google Scholar] [CrossRef]
  196. Liu, L.X.; Litster, J.D.; Iveson, S.M.; Ennis, B.J. Coalescence of deformable granules in wet granulation processes. AIChE J. 2000, 46, 529–539. [Google Scholar] [CrossRef]
  197. Marshall, C.L.J.; Rajniak, P.; Matsoukas, T. Multi-component population balance modeling of granulation with continuous addition of binder. Powder Technol. 2013, 235, 211–220. [Google Scholar] [CrossRef]
  198. Matsoukas, T.; Kim, T.; Lee, K. Bicomponent aggregation with composition-dependent rates and the approach to well-mixed state. Chem. Eng. Sci. 2009, 64, 787–799. [Google Scholar] [CrossRef]
  199. Sahni, E.K.; Chaudhuri, B. Contact drying: A review of experimental and mechanistic modeling approaches. Int. J. Pharm. 2012, 434, 334–348. [Google Scholar] [CrossRef]
  200. Mortier, S.T.F.C.; de Beer, T.; Gernaey, K.V.; Remon, J.P.; Vervaet, C.; Nopens, I. Mechanistic modelling of fluidized bed drying processes of wet porous granules: A review. Eur. J. Pharm. Biopharm. 2011, 79, 205–225. [Google Scholar] [CrossRef]
  201. Mortier, S.T.F.C.; van Daele, T.; Gernaey, K.V.; de Beer, T.; Nopens, I. Reduction of a single granule drying model: An essential step in preparation of a population balance model with a continuous growth term. AIChE J. 2013, 59, 1127–1138. [Google Scholar] [CrossRef]
  202. Johanson, J.R. A rolling theory for granular solids. J. Appl. Mech. 1965, 32, 842–848. [Google Scholar] [CrossRef]
  203. Hsu, S.H.; Reklaitis, G.V.; Venkatasubramanian, V. Modeling and control of roller compaction for pharmaceutical manufacturing. Part I: Process dynamics and control framework. J. Pharm. Innov. 2010, 5, 14–23. [Google Scholar] [CrossRef]
  204. Marshall, E.A. A theory for the compaction of incompressible granular materials by rolling. IMA J. Appl. Math. 1973, 12, 21–36. [Google Scholar] [CrossRef]
  205. Muliadi, A.R.; Litster, J.D.; Wassgren, J.R. Modeling the powder roll compaction process: Comparison of 2-D finite element method and the rolling theory for granular solids (Johanson’s model). Powder Technol. 2012, 221, 90–100. [Google Scholar] [CrossRef]
  206. Bilgili, L.; Capece, M. Quantitative analysis of multi-particle interactions during particle breakage: A discrete non-linear population balance framework. Powder Technol. 2011, 213, 162–173. [Google Scholar] [CrossRef]
  207. Verkoeijen, D.; Pouw, G.A.; Meesters, M.H.; Scarlett, B. Population balances for particulate processes—A volume approach. Chem. Eng. Sci. 2002, 57, 2287–2303. [Google Scholar] [CrossRef]
  208. Bilgili, E.; Scarlett, B. Population balance modeling of non-linear effects in milling processes. Powder Technol. 2005, 153, 59–71. [Google Scholar] [CrossRef]
  209. Tsoy, E.N. A modelling approach for derivation of the breakage functions. Chem. Eng. Sci. 2012, 80, 361–364. [Google Scholar] [CrossRef]
  210. Bilgili, E.; Yepes, J.; Scarlett, B. Formulation of a non-linear framework for population balance modeling of batch grinding: Beyond first-order kinetics. Chem. Eng. Sci. 2006, 61, 33–44. [Google Scholar] [CrossRef]
  211. Heckel, R.W. An analysis of powder compaction phenomena. Trans. Metall. Soc. AIME 1961, 221, 1001–1008. [Google Scholar]
  212. Heckel, R.W. Density pressure relationships in powder compaction. Trans. Metall. Soc. AIME 1961, 221, 671–675. [Google Scholar]
  213. Kawakita, K.; Ludde, K.H. Some considerations on powder compression equations. Powder Technol. 1971, 4, 61–68. [Google Scholar] [CrossRef]
  214. Kawakita, K.; Hattori, I.; Kishigami, M. Characteristic constants in kawakita’s powder compression equation. J. Powder Bulk Solids Technol. 1977, 1, 3–8. [Google Scholar]
  215. Kuentz, M.; Leuenberger, H. Pressure susceptibility of polymer tablets as a critical property: A modified Heckel equation. J. Pharm. Sci. 1999, 88, 174–179. [Google Scholar] [CrossRef]
  216. Singh, R.; Gernaey, K.V.; Gani, R. Icas-pat: A software for design, analysis and validation of PAT systems. Comput. Chem. Eng. 2010, 34, 1108–1136. [Google Scholar] [CrossRef]
  217. Frenning, G.; Nordström, J.; Alderborn, G. Effective Kawakita parameters for binary mixtures. Powder Technol. 2009, 189, 270–275. [Google Scholar] [CrossRef]
  218. Mazel, V.; Busignies, V.; Duca, S.; Leclerc, B.; Tchoreloff, P. Original predictive approach to the compressibility of pharmaceutical powder mixtures based on the Kawakita equation. Int. J. Pharm. 2011, 410, 92–98. [Google Scholar] [CrossRef]
  219. Robinson, S.; Brooks, R.J. Independent verification and validation of an industrial simulation model. Simulation 2010, 86, 405–416. [Google Scholar] [CrossRef]
  220. Davies, P.K. Generalizing Concepts and Methods of Verification, Validation and Accreditation (VV&A) for Military Simulation; RAND: Santa Monica, CA, USA, 1992. [Google Scholar]
  221. Sargent, R.G. Verification and Validation of Simulation Models. In Proceedings of the 2010 Winter Simulation Conference, Baltimore, MD, USA, 5–8 December 2010.
  222. Kremer, D.M.; Hancock, B.C. Process simulation in the pharmaceutical industry: A review of some basic physical models. J. Pharm. Sci. 2006, 95, 517–529. [Google Scholar] [CrossRef]
  223. Balci, O. Golden rules of verification, validation, testing, and certification of modeling and simulation applications. SCS M&S Mag. 2010. No. 4. [Google Scholar]
  224. Saltelli, A.; Chan, K.; Scott, E.M. Sensivivity Analysis; John Wiley & Sons Ltd.: Chichester, UK, 2000. [Google Scholar]
  225. Saltelli, A.; Tarantola, S.; Campolongo, F. Sensitivity analysis as an ingredient of modeling. Stat. Sci. 2000, 15, 377–395. [Google Scholar] [CrossRef]
  226. Sin, G.; Gernaey, K.V.; Eliasson Lantz, A. Good modelling practice (gmop) for pat applications: Propagation of input uncertainty and sensitivity analysis. Biotechnol. Prog. 2009, 25, 1043–1053. [Google Scholar] [CrossRef]
  227. Ling, Y.; Mahadevan, S. Quantitative model validation techniques: New insights. Reliab. Eng. Syst. Saf. 2013, 111, 217–231. [Google Scholar] [CrossRef]
  228. Rodgers, J.L.; Nicewander, W.A. Thirteen ways to look at the correlation coefficient. Am. Stat. 1988, 41, 59–66. [Google Scholar]
  229. Min, F.Y.; Yang, M.; Wang, Z.C. Knowledge-based method for the validation of complex simulation models. Simul. Model. Pract. Theory 2010, 18, 500–515. [Google Scholar] [CrossRef]
  230. Tedeschi, L.O. Assessment of the adequacy of mathematical models. Agric. Syst. 2006, 9, 225–247. [Google Scholar] [CrossRef]
  231. Browne, M.W. Cross-validation methods. J. Math. Psychol. 2000, 44, 108–132. [Google Scholar] [CrossRef]
  232. Efron, B.; Gong, G. A leisurely look at the bootstrap, the jackknife and cross-validation. Am. Stat. 1983, 37, 36–48. [Google Scholar]
  233. Efron, B. Estimating the error rate of a prediction rule: Improvement on cross-validation. J. Am. Stat. Assoc. 1983, 78, 316–331. [Google Scholar] [CrossRef]
  234. Tomba, E.; de Martin, M.; Facco, P.; Robertson, J.; Zomer, S.; Bezzo, F.; Barolo, M. General procedure to aid the development of continuous pharmaceutical processes using multivariate statistical modeling—An industrial case study. Int. J. Pharm. 2013, 444, 25–39. [Google Scholar] [CrossRef]
  235. Hassani, S.; Martens, H.; Qannari, E.M.; Hanafi, M.; Kohler, A. Model validation and error estimation in multi-block partial least squares regression. Chemom. Intell. Lab. Syst. 2012, 117, 42–53. [Google Scholar] [CrossRef]
  236. Steyerberg, E.W.; Harrell, F.E.J.; Borsboom, G.J.J.M.; Eijkemans, M.J.C.R.; Vergouwe, Y.; Habbema, J.D.F. Internal validation of predictive models: Efficiency of some procedures for logistic regression analysis. J. Clin. Epidemiol. 2001, 54, 775–781. [Google Scholar]
  237. Gong, G. Cross-validation, the jackknife, and the bootstrap: Excess error estimation in forward logistic regression. J. Am. Stat. Assoc. 1986, 81, 108–113. [Google Scholar] [CrossRef]
  238. Ji, G.; Zhang, K.; Zhu, Y. A method of MPC model error detection. J. Process Control 2012, 22, 635–642. [Google Scholar] [CrossRef]
  239. Zhu, Y.; Patwardan, R.; Wagner, S.B.; Zhao, J. Toward a low cost and high performance MPC: The role of system identification. Comput. Chem. Eng. 2013, 51, 124–135. [Google Scholar] [CrossRef]
  240. Robinson, A.P.; Froese, R.E. Model validation using equivalence tests. Ecol. Model. 2004, 176, 349–358. [Google Scholar] [CrossRef]
  241. Rebba, R.; Mahadevan, S. Validation of models with multivariate output. Reliab. Eng. Syst. Saf. 2006, 91, 861–871. [Google Scholar] [CrossRef]
  242. Capece, M.; Bilgile, E.; Dave, R. Identification of the breakage rate and distribution parameters in a non-linear population balance model for batch milling. Powder Technol. 2011, 208, 195–204. [Google Scholar] [CrossRef]
  243. Vanarase, A.U.; Jarvinen, M.; Passo, J.; Muzzio, F.J. Development of a methodology to estimate error in the on-line measurements of blend uniformity in a continuous powder mixing process. Powder Technol. 2013, 241, 263–271. [Google Scholar] [CrossRef]
  244. Babamoradi, H.; van den berg, F.; Rinnan, A. Bootstrap based confidence limits in principal component analysis—A case study. Chemom. Intell. Lab. Syst. 2013, 120, 97–105. [Google Scholar] [CrossRef]
  245. Afanador, N.L.; Tran, T.N.; Buydens, L.M.C. Use of the bootstrap and permutation methods for a more robust variable importance in the projection metric for partial least squares regression. Anal. Chim. Acta 2013, 768, 49–56. [Google Scholar] [CrossRef]
  246. Ng, K.M. Design and development of solids processes—A process systems engineering perspective. Powder Technol. 2002, 126, 205–210. [Google Scholar] [CrossRef]
  247. Kimber, J.A.; Kazarian, S.G.; Štěpánek, F. Microstructure-based mathematical modelling and spectroscopic imaging of tablet dissolution. Comput. Chem. Eng. 2011, 35, 1326–1339. [Google Scholar]

Share and Cite

MDPI and ACS Style

Rogers, A.J.; Hashemi, A.; Ierapetritou, M.G. Modeling of Particulate Processes for the Continuous Manufacture of Solid-Based Pharmaceutical Dosage Forms. Processes 2013, 1, 67-127. https://0-doi-org.brum.beds.ac.uk/10.3390/pr1020067

AMA Style

Rogers AJ, Hashemi A, Ierapetritou MG. Modeling of Particulate Processes for the Continuous Manufacture of Solid-Based Pharmaceutical Dosage Forms. Processes. 2013; 1(2):67-127. https://0-doi-org.brum.beds.ac.uk/10.3390/pr1020067

Chicago/Turabian Style

Rogers, Amanda J., Amir Hashemi, and Marianthi G. Ierapetritou. 2013. "Modeling of Particulate Processes for the Continuous Manufacture of Solid-Based Pharmaceutical Dosage Forms" Processes 1, no. 2: 67-127. https://0-doi-org.brum.beds.ac.uk/10.3390/pr1020067

Article Metrics

Back to TopTop