A Novel PBM for Nanomilling of Drugs in a Recirculating Wet Stirred Media Mill: Impacts of Batch Size, Flow Rate, and Back-Mixing

We examined the evolution of fenofibrate (FNB, drug) particle size distribution (PSD) during the production of nanosuspensions via wet stirred media milling (WSMM) with a cell-based population balance model (PBM). Our objective was to elucidate the potential impacts of batch size, suspension volumetric flow rate, and imperfect mixing in a recirculating WSMM. Various specific breakage rate functions were fitted to experimental PSD data at baseline conditions assuming perfect mixing. Then, the best function was used to simulate the PSD evolution at various batch sizes and flow rates to validate the model. A novel function, which is a product of power–law and logistic functions, fitted the evolution the best, signifying the existence of a transition particle size commensurate with a grinding limit. Although larger batches yielded coarser and wider PSDs, the suspensions had identical PSDs when milled for the same effective milling time. The flow rate had an insignificant influence on the PSD. Furthermore, the imperfect mixing in the mill chamber was simulated by considering more than one cell and different back-mixing flow ratios. The effects were weak and restricted to the first few turnovers. These insights contribute to our understanding of recirculating WSMM, providing valuable guidance for process development.


Introduction
The enhancement of bioavailability stands as an important focus area in pharmaceutical research and development, dedicated to improving the therapeutic effectiveness of drug formulations.Nanoparticles have emerged as highly promising vehicles in this pursuit because of their large surface area, enhanced saturation solubility, sustained release, and targeted delivery to enhance bioavailability [1][2][3].Nanoparticles, which are referred to here as crystalline particles with diameters between 50 and 500 nm, have a specific surface area that is about two orders of magnitude higher than that of coarse micro-sized crystals.This greater surface area, in conjunction with the higher overall solute mass transfer coefficient resulting from thinner diffusion layer in fluids and the increased saturation solubility, especially of <100 nm particles, enable faster drug dissolution [4].
Wet stirred media milling (WSMM) is the most widely used process in pharmaceutical engineering for producing stable drug nanosuspensions [5], which provide drug bioavailability at higher drug doses per injection volume [5,6].The wide use of WSMM has originated from its robust, reproducible, scalable, solvent-free, and environmentally friendly operation [2,7,8].A typical pharmaceutical WSMM process runs in recirculation mode [6], wherein a drug suspension is recirculated between a holding tank and a milling chamber multiple times.High-speed rotation of a stirrer (rotor) induces turbulent motion, leading to repeated stressing and fracture of particles captured between colliding media (beads).A challenge in the preparation of drug nanosuspensions lies in ensuring their physical stability, a critical factor for their successful integration into drug delivery systems [2,[9][10][11].While achieving stability in drug nanosuspensions is of paramount importance, this challenge has been addressed through careful and strategic selection, as well as screening of stabilizers and their concentrations [2,[12][13][14][15].
WSMM, nonetheless, comes with inherent drawbacks such as being time-consuming, costly, and energy-intensive [16].Moreover, while recognizing issues like potential product contamination due to bead wear, process-induced solid-state changes, and prolonged milling times required for drug nanosuspension preparation, the pharmaceutical nanotechnology literature has focused on the stabilization of drug nanosuspensions [17][18][19], with relatively scant information on process modeling and the resolution of the processing challenges [20].Given that breakage kinetics govern the cycle time and production rate for achieving a desired fineness, the development and optimization of WSMM as well as resolution of the above-mentioned problems necessitate a profound understanding of breakage kinetics and its controlling process-design parameters.In this regard, the size and material of the beads affect breakage kinetics significantly: there exists an optimal bead size, and ceramic beads enable faster milling than crosslinked polystyrene beads [21,22].Regarding the process parameters, an increase in the stirrer speed enhances the breakage rate, evident from smaller median sizes and 90% passing sizes of particles at any given milling time during WSMM [1,8,23,24].While higher bead loading increases the breakage rate [8,23,25,26], an increase in drug loading reduces it albeit with an improvement in operational efficiency [23,27].In contrast to other process parameters, suspension flow rate and batch size of a recirculating WSMM have been rarely examined in terms of their impact on the breakage kinetics and PSD evolution; hence, investigation of their potential impacts is within the scope of this paper.
Population balance models (PBMs) are often used for simulation, control, and optimization of a wide variety of particulate processes [43,44].In addition to their capability to describe the spatial and/or temporal evolution of the particle size distribution (PSD), they have also been utilized to identify particle breakage mechanisms [45,46].While several PBM studies exist for the simulation of batch WSMM operation [34] and multi-pass continuous WSMM operation [47], Annapragada and Adjei [32] offered the first credible PBM structure for simulating a recirculating WSMM in the literature, wherein dextrose particles were milled to yield an aerosol suspension (not a nanosuspension).Their model assumed perfect mixing in both the milling chamber and the holding tank.Unfortunately, they did not measure or model the evolution of the entire PSD to validate their model.In addition, contrary to the authors' assertions, the model prediction deviated significantly from the limited experimental data.Their simulations of a recirculating WSMM, without any experimental verification, suggested that faster breakage occurred when a higher rotor speed, higher suspension flow rate, and smaller beads were used, whereas slightly finer PSDs were obtained when bead loading was increased in a narrow range (80-85% v/v).Some of these findings are in line with previous WSMM studies (refer to the review [2]); however, the significant positive impact of the suspension flow rate is quite surprising, which has not been verified experimentally.Also, they did not investigate the impact of the batch size.Both aspects will be scrutinized in the current study with experimental verification.
Pharmaceutics 2024, 16, 353 3 of 28 In this study, we aim to elucidate potential impacts of batch size, suspension volumetric flow rate, and imperfect mixing in a recirculating stirred mill.A cell-based PBM, which reduces to the PBM in [32] for one cell (perfect mixing), was developed.First, four different specific breakage rate functions were fitted to experimental PSD data at baseline conditions assuming a well-mixed mill.Then, the PSD evolution at various batch sizes and flow rates was simulated using the best function to compare with respective experimental data and validate the PBM.Furthermore, potential impact of imperfect mixing in the mill chamber was simulated by considering more than one cell and different back-mixing flow ratios of the cell-based PBM.The analysis of the simulation results along with experimental verification will allow us to develop a fundamental understanding of recirculating WSMM operation, providing valuable guidance for process development.

Experimental 2.1. Materials
The drug utilized in this study was fenofibrate (FNB), which was obtained in BP grade from Jai Radhe Sales (Ahmedabad, India).Fenofibrate is considered a poorly water-soluble drug, with an aqueous solubility of 0.8 mg/L at room temperature [24].Hydroxypropyl cellulose (HPC) of L grade, donated by Nisso America Inc. (New York, NY, USA), was used as a non-ionic polymeric stabilizer.An anionic surfactant, sodium dodecyl sulfate (SDS) of ACS grade, was purchased from GFS chemicals (Columbus, OH, USA).Zirmil Y grade yttrium-stabilized zirconia (YSZ) beads from Saint Gobain ZirPro (Mountainside, NJ, USA) with a density of 6000 kg/m 3 and a nominal size of 400 µm were used as grinding media.The actual median size of the beads (405 µm) was measured using a laser diffraction technique in dry dispersion mode, employing a Helos/Rodos particle size analyzer (Sympatec, Pennington, NJ, USA).

Preparation and Characterization Methods
To prepare an FNB pre-suspension, hydroxypropylcellulose-L (HPC-L) and sodium dodecyl sulfate (SDS) were dissolved in 200 mL of deionized water first, followed by dispersion of the FNB powder using a shear mixer (Cat.#14-503, Fisher Scientific, Pittsburgh, PA, USA).The mixer ran at a speed of 300 rpm for 2 h.The formulation composition of the suspension was determined in the view of previous research [35]: 10% FNB, 7.5% HPC-L, and 0.05% SDS.After preparation, the pre-suspension was stored at a temperature of 8 • C overnight.Subsequently, the pre-suspension was then transferred to a Microcer mill (Netzsch Fine Particle Size Technology, LLC, Exton, PA, USA) and ground for a duration of t = 180 min at a stirrer speed of 3000 rpm and bead loading of c = 0.5 (v/v).Bead loading c was determined by taking the ratio of the true volume of the beads to the volume of the milling chamber (V m = 80 mL) on a volumetric basis (v/v).The actual volume occupied by the beads was used to determine the bead loading.To ensure continuous recirculation of the suspension (Figure 1), a peristaltic pump (Cole-Palmer, Master Flex, Vermont Hills, IL, USA) was employed.The pump facilitated the flow of the suspension between the holding tank and the milling chamber at a volumetric flow rate of Q = 126 mL/min for the baseline conditions (Table 1).To retain the beads within the milling chamber, a stainless-steel screen with openings half the size of the nominal bead size (200 µm) was utilized.This screen served as a barrier to prevent the beads from exiting the milling chamber while allowing for the suspension to flow through.To remove heat generated during the milling, a chiller (Model M1-.25A-11HFX, Advantage Engineering, Greenwood, IN, USA) was employed.The chiller effectively cooled the setup, ensuring that the temperature remained under control throughout the milling process.In 4 additional experiments, batch sizes and suspension flow rates were doubled and halved from their respective baseline values to examine their impacts (see Table 1).The volume-based PSD of the drug suspensions at different milling times was determined using laser diffraction with an LS 13-320 Beckman Coulter instrument (Brea, CA, USA).Samples were collected from the outlet of the mill at specific time intervals, denoted as 2 s , where s represents the time interval (s = 0, 1, 2, …, 7), along with additional samples taken at 24, 48, 96, and 180 min.The final sample was obtained from the holding tank.Prior to each measurement, approximately 1.0 mL of the suspension sample was diluted with 5.0 mL of the stabilizer solution using a vortex mixer (Fisher Scientific Digital Vortex Mixer, Model No: 945415, Pittsburgh, PA, USA) at 1500 rpm for one minute.During the measurements, polarized intensity differential scattering (PIDS) was maintained between 40% and 50%, while the obscuration remained below 8%.The PSD was determined by the instrument's software, which utilized the Mie scattering theory, considering refractive indices of 1.55 for FNB and 1.33 for water.The measurements were repeated four times, and the average PSD was used in the model fits.

Preliminaries: Residence Time Distribution
In continuous WSMM, the product PSD is affected by the residence time distribution (RTD) of the particles in the mill depending also on the specific mode of operation: single pass, multi-pass, and recirculation (circuit) [48,49].To measure RTD function experimentally, a tracer pulse is introduced at the mill's inlet point, and the tracer's concentration is then recorded at the outlet of mill in single-pass operation mode [50][51][52].Space time, which equals the average residence time if there are no dead zones and bypass in the milling chamber, can be calculated using τs = Vm(1 − c)/Q.For a single-pass continuous mill, an increase in τs results in a finer PSD; a plug-flow regime yields a narrower and somewhat finer PSD in the mill than the perfect-mixing regime [53,54].It should be noted  The volume-based PSD of the drug suspensions at different milling times was determined using laser diffraction with an LS 13-320 Beckman Coulter instrument (Brea, CA, USA).Samples were collected from the outlet of the mill at specific time intervals, denoted as 2 s , where s represents the time interval (s = 0, 1, 2, . .., 7), along with additional samples taken at 24, 48, 96, and 180 min.The final sample was obtained from the holding tank.Prior to each measurement, approximately 1.0 mL of the suspension sample was diluted with 5.0 mL of the stabilizer solution using a vortex mixer (Fisher Scientific Digital Vortex Mixer, Model No: 945415, Pittsburgh, PA, USA) at 1500 rpm for one minute.During the measurements, polarized intensity differential scattering (PIDS) was maintained between 40% and 50%, while the obscuration remained below 8%.The PSD was determined by the instrument's software, which utilized the Mie scattering theory, considering refractive indices of 1.55 for FNB and 1.33 for water.The measurements were repeated four times, and the average PSD was used in the model fits.

Preliminaries: Residence Time Distribution
In continuous WSMM, the product PSD is affected by the residence time distribution (RTD) of the particles in the mill depending also on the specific mode of operation: single pass, multi-pass, and recirculation (circuit) [48,49].To measure RTD function experimentally, a tracer pulse is introduced at the mill's inlet point, and the tracer's concentration is then recorded at the outlet of mill in single-pass operation mode [50][51][52].Space time, which equals the average residence time if there are no dead zones and bypass in the milling chamber, can be calculated using τ s = V m (1 − c)/Q.For a single-pass continuous mill, an increase in τ s results in a finer PSD; a plug-flow regime yields a narrower and somewhat finer PSD in the mill than the perfect-mixing regime [53,54].It should be noted that the average residence time in the mill for a recirculation operation τ c can be regarded as the effective milling time and is described by τ c = N to τ s , where N to is the number of turnovers defined by N to = Qt/V s .
The RTD of particles yields information about the transport behavior and back-mixing extent, which affect the product PSD [49].Previous RTD studies on single-pass WSMM demonstrated that (i) the RTD function is close to that of perfect mixing especially for short mills with small length-to-diameter L/D ratio [33] and (ii) longer mills exhibit deviations from perfect mixing (imperfect mixing), which can be explained by either the convectiondispersion model with the Peclet number or the cell-based RTD model with (internal) recirculation between the cells [33].The Peclet number for a single-pass WSMM is defined by Pe s = u a L/ D, where u a , L, and D stand for the axial interstitial velocity of the drug suspension, length of the mill chamber, and dispersion coefficient; u a relates to the axial superficial velocity u s through u a = u s /(1 − c), where u s = Q/A c and A c denote the crosssectional area of the mill chamber perpendicular to the suspension flow.Two limiting cases reflect ideal mixing regimes: Pe s = 0 (perfect mixing) and Pe s → ∞ (plug flow, no axial back-mixing).Experimental RTD studies [33,49] suggest 0.6 < Pe s < 5, which signifies that perfect mixing is a decent approximation in a single-pass WSMM, albeit with some deviations for long mills.

Formulation of the PBM
Initially introduced by Whiten [55], the cell-based RTD model has been used as an alternative to the convection-dispersion model [33,47,49,[56][57][58][59].The mill is conceptualized as an assembly of n well-mixed cells with recirculation between them (internal recirculation) superposed on the recirculating bulk flow of the suspension.In particular, the cell-based RTD model explicitly incorporates an axial recirculation rate of the suspension .R between adjacent cells in the mill chamber.A dimensionless axial back-mixing ratio R is defined as R = .R/Q, which modulates the extent of back-mixing or imperfect mixing, along with n.Kwade and Schwedes [49] conducted a significant study affirming the accuracy of the cell-based RTD model in fitting the RTD of a WSMM.In their investigation, the number of cells equaled the number of stirrer discs, and a single parameter-the back-mixing ratio R-was fitted to the measured RTD.Notably, this estimation of R through the cellbased model did not include the evolution of the PSD because milling was not considered.Further advancements in the use of the cell-based PBM approach were presented by Fadhel et al. [57] and Frances [33].These researchers employed a modified cell-based PBM to fit experimental steady-state PSDs in WSMM; however, a recirculating WSMM with transient operation has not been modeled.
We adapted the cell-based PBM developed by [54,55] and tailored it to model a recirculating WSMM.The final cell-based PBM for suspension in the mill chamber and the holding tank can be recorded as where z and T refer, respectively, to the cell index that varies from 1 to n and the holding tank; i and j are the size-class indices.Equations ( 1) and (2) pose a set of N × n ordinary differential equations (ODEs) with the following initial conditions: M i,z (0) = M i,T (0) = M i,ini at time t = 0. Size classes 1 and N contain the coarsest particles and finest particles, respectively.M i , S i , and b ij stand for the mass concentration of (drug) particles in size class i, the specific breakage rate function, and the breakage distribution function, respectively.S i describes the rate at which particles in size class i are broken, whereas b ij represents the mass fraction of particles broken from size class j into size class i.The latter function is related to its cumulative counterpart B ij through b ij = B ij − B i+1j .Both functions depend on particle size x i .The parameter τ cell , denoting residence time within the cell, is defined as τ cell = τ s /n.Similarly, τ T , representing the average residence time within the holding tank, is defined as τ T = V T /Q, where V T is the holding tank volume.The first term on the right-hand-side of Equation ( 1) explains the rate of disappearance of particles in size class i due to breakage.The second term explains the rate of appearance of progeny particles in size class i due to breakage of all large particles in size classes j < i (x j > x i ).The last term in curly brace accounts for particle transport rates between any generic cell z and adjacent cells z − 1 and z + 1 due to recirculating bulk flow quantified by Q and internal recirculation flow quantified by .R. R is the dimensionless axial back-mixing ratio.Equation (1) accounts for finite axial mixing/dispersion in a long mill for a corresponding single-pass operation via R and n, while Equation (2) simulates the dynamics in the holding tank.The axial particle mixing/dispersion arises from the random motion of particles-beads along the longitudinal axis of a mill with an effective length L and/or radial velocity variations.The model allows for the investigation of limiting cases, such as zero axial back-mixing (plug flow) and perfect mixing, by considering scenarios with a large number of cells or a single well-mixed cell, respectively.The model structure is akin to that in ref. [32] when perfect mixing is assumed for short mills (z = n = 1).

Functional Forms of S i and B ij
The specific breakage rate function S i and the cumulative breakage distribution function B ij play an important role in fitting and simulating the breakage of particles in WSMM using the PBM.The following normalized power-law function was used to describe B ij : where the exponent a 1 sets the slope of the B ij function in a log-log plot against the normalized particle size.Such a simple function was successfully used in the PBM of WSMM before [60]; it allows us to minimize the number of parameters to fit.It is well-established in the milling literature that the B ij function is relatively insensitive to the process parameters as compared with the S i function, which is a strong function of the process parameters [60,61].
The following S i functions, referred to as Models A-D, with different mathematical complexity and representation of different breakage kinetics were used: Model A : Model B : Model C : In Equations ( 4)-( 7), A is the specific breakage rate constant of particles with x 0 size, m is the breakage rate exponent, x 0 is the normalizing reference particle size and was set to particle size of class 1 x 1 .In Equations ( 6) and ( 7), s f is defined as the shape factor of the respective logistic function.Note that the power-law model, Model A, is one of the most used models in the milling literature and has been successfully used for a multitude of materials in various types of mills including the WSMM [60,62,63].However, other WSMM investigations [8,[64][65][66] experimentally identified the existence of a grinding limiting size, which is typically reported in terms of the cumulant size of the PSD; below this size, no particle size reduction takes place.However, to the best knowledge of the authors, the notion or phenomenon of a grinding limit or limiting size has not been incorporated into an S i function explicitly before.Bilgili et al. [34] observed a transition of the S i during the WSMM of a magenta pigment; however, such transition was not modeled using any of the above S i functions.Hence, in this study, commensurate with the notion of a grinding limit, we hypothesize that there exists a transition particle size x* below which the S i decreases to 0 sharply (Model B) or gradually yet much more drastically than what a typical power-law model predicts through the logistic functions (Models C and D).Models C and D are identical except for the way the difference between x i and x* was expressed: logarithmic vs. linear.

A Full Back-Calculation Method for the Estimation of PBM Parameters
The parameters of the S i and B ij functions were assumed invariant in different runs (Table 1) because all process parameters, except for the suspension flow rate and the batch size, were kept constant.We hypothesize that while the suspension flow rate and the batch size can affect the PSD evolution, they do not affect the breakage functions unlike stirrer speed, bead loading, etc.Therefore, we first attempted to find the most appropriate S i model among Models A-D by fitting the PBM at the baseline experimental conditions (Run 1), as will be detailed in the next section.Since our Microcer mill is a small lab-scale stirred mill with L/D < 1, its RTD is expected to follow perfect mixing closely [33], and thus we set n = 1.Then, using the calibrated PBM, we predicted the impacts of batch size and volumetric flow rate on the PSD evolution (Runs 2-5) and compared them to the experimental PSD evolution, which allowed us to validate the PBM.Note that we did not examine the impacts of the stirrer speed, bead loading/size, and bead type in this study partly because they were experimentally studied previously [21,23,27] and partly their consideration would entail a more elaborate experimental study and a PBM structure with more than 8 parameters, which is beyond the scope of this paper.
To estimate the breakage parameters of the PBM and discriminate Models A-D, a full back-calculation method was adopted here [67].This method employs a dynamic global optimizer-ODE solver to fit the PBM to the experimental PSD data, and therefore estimate the parameters.Specifically, the GlobalSearch function and fmincon solver within the MATLAB 2022a software [68] were utilized to minimize the following sum-of-squared residuals SSR between the experimentally determined PSDs in Run 1 and the modelpredicted PSDs: where m MOD iz,p and m EXP iz,p denote, respectively, the simulated and the experimental mass fraction density distributions; p is the time index, P is the total number of PSD sampling (time points), and K = 94 is the number of size classes in the laser diffraction measurements.In the fitting, we specifically used z = n = 1 (perfect mixing) and fitted the PSDs at all time points (P = 10), except for 1 min and 2 min due to the onset of initial aggregation, as will be discussed in Section 4. The MATLAB function "GlobalSearch" was used to generate the next set of initial guesses for the next trial point using the scatter search method [69].Throughout this paper, the PSD represented by the mass fraction density distribution in any cell z, i.e., m iz , was calculated from the mass fraction m iz as per Equations ( 9) and (10).Assuming the particle density remains invariant with particle size, the experimental volume-based PSD was taken as the mass-based PSD in Equation (8).Readers are referred to Appendix A for the numerical details of optimization-parameter estimation.

Impacts of the Batch Size, Flow Rate, and Imperfect Mixing and PBM Validation
After fitting the PBM to the PSD evolution of Run 1 (baseline conditions), determining its parameters, and discriminating Models A-D, we used the best S i model in various simulations to examine the impacts of batch size (Runs 2 and 3), volumetric flow rate (Runs 4 and 5), and the degree of imperfect mixing as modulated by R and n (Runs 6-17).The detailed design of the simulation study is outlined in Table 2, providing a structured overview of the parameters and conditions examined in the simulations.As the evolution of PSD in Runs 2-5 was also studied experimentally, these runs allow us to test the predictive capability of the model and validate it.Note that Runs 6-17 have identical process-design parameters as those of Run 1, except that they correspond to different theoretical single-pass RTD in the milling chamber.They enable a theoretical investigation of the impacts of imperfect mixing or deviation from perfect mixing in the milling chamber because n = 1 in Runs 1-5 signifies perfect mixing in a theoretical single pass of a continuous mill.In the absence of recirculation from a holding tank, for a single pass of the suspension through a continuous mill, n = 1 corresponds to perfect mixing.An increase in n or a decrease in R for a single-pass continuous mill represents a decrease in axial back-mixing, which implies an approach to plug flow in the milling chamber [54].For R = 0 (traditional n-equal volume tank in series RTD), even when n was increased from 1 to 5, the product PSD was closer to that for plug flow than that for perfect mixing, and for n = 15 and n = 60, plug-flow behavior was approached and reproduced, respectively [54].All these findings are for a theoretical single-pass continuous mill without recirculation from the holding tank; therefore, the simulations in Table 2 will shed light on the impact of (single-pass) imperfect mixing on the evolution of the PSD in a recirculating WSMM.

Estimating the Breakage Parameters and Discrimination of Models A-D
We fitted the PBM to the PSD evolution of Run 1 (baseline conditions) using Models A-D for S i , estimated the breakage parameters, and discriminated Models A-D based on their SSR.As will be elaborated in the sequel, Model C turned out to be the best S i model with the lowest SSR. Figure 2a,b shows the PBM simulation with Model C (with fitted parameters for N T = 1000) and experimental PSD evolution over the initial 32 min and 64-128 min, respectively.Only selected time points are shown for clarity, e.g., the feed PSD (t = 0 min) is excluded for proper scaling of Figure 2. The PSD evolution starting with the feed PSD is illustrated in Figure S1 (Supplementary Materials) for the sake of completeness.Unless otherwise indicated, the PSDs refer to the suspensions at the mill outlet.

Estimating the Breakage Parameters and Discrimination of Models A-D
We fitted the PBM to the PSD evolution of Run 1 (baseline conditions) using Models A-D for Si, estimated the breakage parameters, and discriminated Models A-D based on their SSR.As will be elaborated in the sequel, Model C turned out to be the best Si model with the lowest SSR. Figure 2a,b shows the PBM simulation with Model C (with fitted parameters for NT = 1000) and experimental PSD evolution over the initial 32 min and 64-128 min, respectively.Only selected time points are shown for clarity, e.g., the feed PSD (t = 0 min) is excluded for proper scaling of Figure 2. The PSD evolution starting with the feed PSD is illustrated in Figure S1 (Supplementary Materials) for the sake of completeness.Unless otherwise indicated, the PSDs refer to the suspensions at the mill outlet.In general, the experimental PSD shifted from right to left, toward finer particles, monotonically as milling was continued (Figure 2).The modal (peak) size also shifted to the left.The feed had the 10%, 50%, and 90% passing sizes of the cumulative PSD of x10 = 5.64 ± 0.04 μm, x50 = 12.3 ± 0.05 μm, and x90 = 22.2 ± 0.06 μm.Upon 128 min of milling, these characteristic particle sizes were drastically reduced to x10 = 0.103 ± 0.001 μm, x50 = 0.149 ± 0.001 μm, and x90 = 0.212 ± 0.000 μm, suggesting a size reduction ratio of ~80.The PBM reasonably predicted these changes, albeit with some notable deviations from the experimental observations, especially at earlier milling times.For example, despite the unimodal PSD of the feed suspension, a bimodal PSD was observed, with a minor mode in the 1-3 μm range, during the first 16 min, which completely disappeared with prolonged milling (after 32 min).Figure S1 illustrates both the unimodal feed PSD and evolving bimodality during the first 8 min.Nevertheless, it is worth noting that the disparity between the simulation and experimental results gradually diminished after the 16-min mark.This is why In general, the experimental PSD shifted from right to left, toward finer particles, monotonically as milling was continued (Figure 2).The modal (peak) size also shifted to the left.The feed had the 10%, 50%, and 90% passing sizes of the cumulative PSD of x 10 = 5.64 ± 0.04 µm, x 50 = 12.3 ± 0.05 µm, and x 90 = 22.2 ± 0.06 µm.Upon 128 min of milling, these characteristic particle sizes were drastically reduced to x 10 = 0.103 ± 0.001 µm, x 50 = 0.149 ± 0.001 µm, and x 90 = 0.212 ± 0.000 µm, suggesting a size reduction ratio of ~80.The PBM reasonably predicted these changes, albeit with some notable deviations from the experimental observations, especially at earlier milling times.For example, despite the unimodal PSD of the feed suspension, a bimodal PSD was observed, with a minor mode in the 1-3 µm range, during the first 16 min, which completely disappeared with prolonged milling (after 32 min).Figure S1 illustrates both the unimodal feed PSD and evolving bimodality during the first 8 min.Nevertheless, it is worth noting that the disparity between the simulation and experimental results gradually diminished after the 16-min mark.This is why 1 min and 2 min PSDs were disregarded from the model fitting.This bimodality could be explained by the presence of a mixture of aggregates and primary drug particles due to the early onset of aggregation.The simulation failed to accurately capture the observed bimodality in the experimental data because the PBM did not account for aggregates that formed during the first 16 min (see Figure 2 and Figure S1).Note that as fine particles were formed, they might have aggregated as polymer-surfactant adsorption is a kinetic process that takes time.This mechanism was supported by the fact that the bimodality disappeared after 32 min as sufficient time was given for adsorption.Moreover, the simulation aligned more closely with the experimental data, indicating an improvement in its fitting capability as the milling progressed.Between 64 min and 128 min (Figure 2b), the PSD became narrower, and the extent of size reduction decreased as the grinding limit was approached.The characteristic sizes at 180 min were x 10 = 0.102 ± 0.008 µm, x 50 = 0.145 ± 0.007 µm, and x 90 = 0.210 ± 0.014 µm, which were almost identical to those at 128 min.This fact further supports our hypothesis that with increasing milling time, the stabilizer adsorption progressed while mitigating aggregation, and an apparent grinding limit was attained.In fact, in our earlier study, 1-week stability testing with the wet-milled FNB suspensions demonstrated that the same HPC-SDS combination provided effective stabilization during the storage [35].Readers are referred to [35] for details of this extensive physical stability study.Overall, Model C fitted the PSD data reasonably well especially for sufficiently long milling times or in the sub-micron region (Figure 2b).The slight deviations in the sub-100 nm region could originate from modeling error, but it may also be due to the limitations of the laser diffraction method for accurately measuring sub-100 nm sizes despite the equipment's use of the PIDS technology.
The breakage parameters, estimated by the GlobalSearch algorithm using Model C, are presented in Table 3.The table includes results for varying numbers of trials, i.e., sets of guessed parameters (N T ), alongside the lower and upper boundaries, as well as the initial guesses of parameters.Interestingly, an increase in N T from 200 up to 1000 did not cause any change in the estimated parameters and SSR.Hence, we conclude that a probable global optimum was reached even at low N T .On the other hand, SSR was reduced, and the parameters were altered when N T was increased for Models A, B, and D (refer to Tables S1-S3 of Supplementary Materials).For all models, a probable globally optimal solution was reached at N T = 1000 (MATLAB's default N T ) as both the SSR and the estimated parameters remained invariant when N T was increased from low values to 1000.Specifically, the SSR and the estimated parameters did not change after 200 (Model C), 400 (Model D), and 800 (Models A and B) trials.Regardless, the parameter estimates were taken from N T = 1000 for all models because a higher N T correlates with a higher probability of having identified a global minimum [70].The outcomes of the parameter estimation for the PBM with Models A-D are presented in Table 4, which allow for direct model discrimination.PSD evolution plots are also presented for Models A, B, and D in Figures S6-S8 of Supplementary Materials for their comparison with Model C (Figure 2).Notably, Model C was the best model with the lowest SSR.Models C and D had similar SSR, which was significantly lower than those of Models A and B. This suggests a similarity in the performance of Models C and D, as can also be seen from a comparison of Figure 2 and Figure S8, making them practically interchangeable.This is not surprising because Models C and D have a similar mathematical structure: a product of power-law and logistic functions.Conversely, Model A, the traditional powerlaw model, had 41% higher SSR than Model C.Although Model B significantly improved upon Model A, the fitted PSD evolution with Model B exhibited a discontinuity in the PSD profile that was absent from the experimental data (Figure S7).Hence, besides their lower SSR, Models C and D were preferred over Model B because of the absence of any discontinuity in their fitted PSD profiles (see Figure 2 and Figure S8).In contrast to most S i parameters of the PBM, the exponent a 1 of B ij function was insensitive to S i model structure: a 1 = 2.33 ± 0.09 (3.9% RSD).In particular, the PBMs with different S i models predict similar breakage distribution/mechanisms of particles.While this is intuitively expected, the optimization-based back calculation is known to be notoriously plagued with interactions between S i and B ij functions in parameter estimation, which can lead to inaccurate estimations when a local optimization scheme is used [59].The small variation of the a 1 , and thus B ij , for different S i models used in the PBM could be attributed to the use of a global optimization scheme, giving further credence to the methods used.Models B-D indicated a transition particle size x* of 195 ± 20 nm (10.2% RSD), which commensurate with the notion of a grinding limit size.This finding is in line with [35], where a rather simple nth-order kinetic model of the WSMM of fenofibrate suspensions under similar milling conditions indicated a grinding limiting size of 142 nm.This interesting phenomenon, first time explored via PBM simulations here, is illustrated in Figure 3. Model A (power-law) exhibited a linear variation with slope m on the log-log plot, whereas the models with x* predicted lower S i and a drastic slope change at x*.The decrease in S i from about 200 nm to 20 nm was remarkable for Models C and D as compared to that for Model A. As compared with the transition for Model B (step function), the transition was smooth for Models C and D, which explains the absence of discontinuity in the PSD evolution (refer to Figure 2 and Figure S8 vs. Figure S7).However, Model C shows a sharper drop in S i with a decrease in particle size below x*, which is more aligned with the existence of a grinding limit.The fitting results overall suggest the following: (i) a PBM without consideration of a transition particle size x* will not be able to describe the breakage phenomenon in nanomilling of drugs accurately and (ii) a novel function, which is a product of power-law and logistic functions, signifying the existence of a transition particle size commensurate with a grinding limit, fitted the evolution the best.
of power-law and logistic functions, signifying the existence of a transition particle size commensurate with a grinding limit, fitted the evolution the best.PBM simulations also suggest that the PSD in the holding tank was coarser than the mill outlet PSD within the first 16 min and these two PSDs became the same thereafter until the end (Figure S9).This finding implies that after a certain milling time, whether the suspension sample for PSD analysis is taken from the tank or the outlet is immaterial.

Impact of Batch Size and Flow Rate
As the PBM with Model C fitted the baseline conditions (Run 1) best, we used this model and its estimated parameters in Table 4 to simulate the impacts of the batch size and the suspension flow rates (Runs 2-5).Figures 4 and 5 present the PBM-predicted and experimental PSD evolution when the batch size was halved (Run 2) and doubled (Run 3) from the baseline value of 236 mL, respectively, keeping the suspension flow rate at 126 mL/min.The predictions of the PBM in Runs 2 and 3 had SSR values of 71.3 and 79.1, respectively, surpassing the SSR of the fitting of Run 1 data, i.e., 38.0.This was also evident from a visual comparison of the deviations of the PBM fitting/prediction and the experimental data in Figures 2, 4 and 5. Albeit being undesirable, this increase in SSR was expected as the PBM parameters were kept the same in the predictions (Runs 2 and 3) as in the fitting (Run 1); they were not re-fitted to experimental data in Runs 2 and 3 wherein the batch size was drastically altered.Despite the elevated SSR, Figures 4 and 5 illustrate the following: (i) an increase in batch size led to notably slower evolution of the PSD, i.e., slower milling of the whole batch and (ii) the PBM predicted the general evolution trend reasonably well, except for the bimodality.The initial bimodality became more evident and sustained up to 64 min when the batch size was increased.As usual, the bimodality disappeared upon prolonged milling in the sub-micron range.
Figures 6 and 7 present the PBM-predicted and experimental PSD evolution when the suspension flow rate was halved (Run 4) and doubled (Run 5) from the baseline value of 126 mL/min, respectively, keeping the batch size at 236 mL.The predictions of the PBM in Runs 4 and 5 had SSR values of 40.9 and 40.8, respectively, which is slightly above the PBM simulations also suggest that the PSD in the holding tank was coarser than the mill outlet PSD within the first 16 min and these two PSDs became the same thereafter until the end (Figure S9).This finding implies that after a certain milling time, whether the suspension sample for PSD analysis is taken from the tank or the outlet is immaterial.

Prediction of the Impacts of Batch Size and Flow Rate Using Model C 4.2.1. Impact of Batch Size and Flow Rate
As the PBM with Model C fitted the baseline conditions (Run 1) best, we used this model and its estimated parameters in Table 4 to simulate the impacts of the batch size and the suspension flow rates (Runs 2-5).Figures 4 and 5 present the PBM-predicted and experimental PSD evolution when the batch size was halved (Run 2) and doubled (Run 3) from the baseline value of 236 mL, respectively, keeping the suspension flow rate at 126 mL/min.The predictions of the PBM in Runs 2 and 3 had SSR values of 71.3 and 79.1, respectively, surpassing the SSR of the fitting of Run 1 data, i.e., 38.0.This was also evident from a visual comparison of the deviations of the PBM fitting/prediction and the experimental data in Figures 2, 4 and 5. Albeit being undesirable, this increase in SSR was expected as the PBM parameters were kept the same in the predictions (Runs 2 and 3) as in the fitting (Run 1); they were not re-fitted to experimental data in Runs 2 and 3 wherein the batch size was drastically altered.Despite the elevated SSR, Figures 4 and 5 illustrate the following: (i) an increase in batch size led to notably slower evolution of the PSD, i.e., slower milling of the whole batch and (ii) the PBM predicted the general evolution trend reasonably well, except for the bimodality.The initial bimodality became more evident and sustained up to 64 min when the batch size was increased.As usual, the bimodality disappeared upon prolonged milling in the sub-micron range.
Figures 6 and 7 present the PBM-predicted and experimental PSD evolution when the suspension flow rate was halved (Run 4) and doubled (Run 5) from the baseline value of 126 mL/min, respectively, keeping the batch size at 236 mL.The predictions of the PBM in Runs 4 and 5 had SSR values of 40.9 and 40.8, respectively, which is slightly above the SSR of the fitting of Run 1 data, i.e., 38.0.The suspension flow rate had an insignificant effect on the PSD evolution, as can be seen from a visual comparison in Figures 2, 6 and 7.
SSR of the fitting of Run 1 data, i.e., 38.0.The suspension flow rate had an insignificant effect on the PSD evolution, as can be seen from a visual comparison in Figures 2, 6 and 7.         To elucidate the impacts of the batch size and the suspension flow rate, the predictions and experimental observations were presented for only 32, 64, and 128 min (Figure 8).In general, one weakness of the current PBM was its inability to explain the bi-modal   To elucidate the impacts of the batch size and the suspension flow rate, the predictions and experimental observations were presented for only 32, 64, and 128 min (Figure 8).In general, one weakness of the current PBM was its inability to explain the bi-modal To elucidate the impacts of the batch size and the suspension flow rate, the predictions and experimental observations were presented for only 32, 64, and 128 min (Figure 8).
In general, one weakness of the current PBM was its inability to explain the bi-modal PSD, especially prevalent during the earlier part of the milling.However, it could capture several salient features of the recirculating WSMM process and the impacts of the process parameters.For any given set of process conditions, as intuitively expected, both the PBM simulations and the experimental data showed a finer PSD for longer milling; hence, the desirable product PSD can be tailored by setting the milling time.The batch size had a dramatic impact on the PSD; in general, a smaller batch was associated with smaller particle sizes and a narrower distribution, indicating a more homogeneous product.Only after a prolonged milling (128 min), the PSD of the larger batch approached those of the baseline and the smaller batches.In contrast, the effect of volumetric flow rate on the PSD was found to be rather small, and mainly notable at 32 min at which bi-modality associated with particle aggregation had a confounding effect.The effect of the volumetric flow rate was not examined experimentally in [32]; however, their simulations indicated a significant impact.In our study, both the simulations and the experiments showed a relatively insignificant impact.The differences may stem from the differences in mill design, materials' formulations, and preparation of micron-sized suspensions vs. nanosuspensions.We will not further speculate why the simulations in [32] showed a notable impact of the volumetric flow rate unlike our simulations; however, it suffices to mention that our experiments support the conclusion based on the simulations.The origin of this insignificant impact of the volumetric flow rate becomes clear when the effective milling time is examined in the next section.
Pharmaceutics 2024, 16,353 15 of 27 PSD, especially prevalent during the earlier part of the milling.However, it could capture several salient features of the recirculating WSMM process and the impacts of the process parameters.For any given set of process conditions, as intuitively expected, both the PBM simulations and the experimental data showed a finer PSD for longer milling; hence, the desirable product PSD can be tailored by setting the milling time.The batch size had a dramatic impact on the PSD; in general, a smaller batch was associated with smaller particle sizes and a narrower distribution, indicating a more homogeneous product.Only after a prolonged milling (128 min), the PSD of the larger batch approached those of the baseline and the smaller batches.In contrast, the effect of volumetric flow rate on the PSD was found to be rather small, and mainly notable at 32 min at which bi-modality associated with particle aggregation had a confounding effect.The effect of the volumetric flow rate was not examined experimentally in [32]; however, their simulations indicated a significant impact.In our study, both the simulations and the experiments showed a relatively insignificant impact.The differences may stem from the differences in mill design, materials' formulations, and preparation of micron-sized suspensions vs. nanosuspensions.We will not further speculate why the simulations in [32] showed a notable impact of the volumetric flow rate unlike our simulations; however, it suffices to mention that our experiments support the conclusion based on the simulations.The origin of this insignificant impact of the volumetric flow rate becomes clear when the effective milling time is examined in the next section.

Milling Time vs. Effective Milling Time
In the context of a recirculating WSMM, the effective milling time is different from the actual milling time; it equals the average residence time in the mill during circuit operation τc and can be expressed by τc = Ntoτs = Vm(1 − c)t/Vs.Hence, an increase in batch size Vs causes a lower τc due to lower number of turnovers, which in turn leads to slower breakage of the whole batch overall.This explains the slower milling operation for the larger batches as illustrated in Figure 8.Note that Nto is directly proportional to the volumetric flow rate, whereas τs is inversely proportional to it.Hence, τc becomes independent of the volumetric flowrate, which explains the insignificant impact of the volumetric flow rate on the PSD evolution (refer to Figure 8).A practical implication of these findings is

Milling Time vs. Effective Milling Time
In the context of a recirculating WSMM, the effective milling time is different from the actual milling time; it equals the average residence time in the mill during circuit operation τ c and can be expressed by τ c = N to τ s = V m (1 − c)t/V s .Hence, an increase in batch size V s causes a lower τ c due to lower number of turnovers, which in turn leads to slower breakage of the whole batch overall.This explains the slower milling operation for the larger batches as illustrated in Figure 8.Note that N to is directly proportional to the volumetric flow rate, whereas τ s is inversely proportional to it.Hence, τ c becomes independent of the volumetric flowrate, which explains the insignificant impact of the volumetric flow rate on the PSD evolution (refer to Figure 8).A practical implication of these findings is that τ c must be kept the same so that the PSD remains invariant when the batch size is changed, or the process is scaled-up.Let us verify this principle via simulations and experiments by comparing the PSD at the same effective milling time τ c .To this end, for the same V m and c, the milling time t was selected in such a way to ensure identical t/V s when the batch size V s was varied.Figure 9 depicts the experimental (a) and simulation (b) results for the PSDs in the smaller batch (118 mL), the baseline (236 mL), and the larger batch (512 mL) conditions at the same effective milling time of 10.8 min, which corresponds to the actual milling times of 32 min, 64 min, and 128 min.At the same effective milling time, the PSDs at various batch sizes were identical, which was confirmed by both the experiments and the PBM simulations.Although the PBM simulations did not perfectly match the experimental PSDs, as discussed before, the PSDs were similar, and the PBM predicted the same outcome for different batch sizes, which gives further credence to the PBM. that τc must be kept the same so that the PSD remains invariant when the batch size is changed, or the process is scaled-up.Let us verify this principle via simulations and experiments by comparing the PSD at the same effective milling time τc.To this end, for the same Vm and c, the milling time t was selected in such a way to ensure identical t/Vs when the batch size Vs was varied.Figure 9 depicts the experimental (a) and simulation (b) results for the PSDs in the smaller batch (118 mL), the baseline (236 mL), and the larger batch (512 mL) conditions at the same effective milling time of 10.8 min, which corresponds to the actual milling times of 32 min, 64 min, and 128 min.At the same effective milling time, the PSDs at various batch sizes were identical, which was confirmed by both the experiments and the PBM simulations.Although the PBM simulations did not perfectly match the experimental PSDs, as discussed before, the PSDs were similar, and the PBM predicted the same outcome for different batch sizes, which gives further credence to the PBM.Figures 10 and 11 present the experimental data and simulation results for the evolution of the 50% passing (median) size x50 as a function of the milling time and the effective milling time, respectively.As expected, x50 decreased drastically during the first 32 min, and the size reduction slowed later on with an approach to the grinding limiting size (Figure 10).The insignificant impact of the suspension flow rate and the faster decrease in x50 for smaller batch were evident.What is most impressive is that when the x50 evolution is expressed in terms of the effective milling time, the curves corresponding to different batch sizes coincided within experimental accuracy (Figure 11).The impact of the suspension flow rate was relatively insignificant regardless of whether x50 evolution was expressed in terms of the milling time or the effective milling time.Overall, we conclude that the PSD equivalency established for 10.8 min effective milling time (Figure 9) is generally applicable during the whole milling, thus further supporting the use of effective milling time in process scale-up and/or batch size changes.Figures 10 and 11 present the experimental data and simulation results for the evolution of the 50% passing (median) size x 50 as a function of the milling time and the effective milling time, respectively.As expected, x 50 decreased drastically during the first 32 min, and the size reduction slowed later on with an approach to the grinding limiting size (Figure 10).The insignificant impact of the suspension flow rate and the faster decrease in x 50 for smaller batch were evident.What is most impressive is that when the x 50 evolution is expressed in terms of the effective milling time, the curves corresponding to different batch sizes coincided within experimental accuracy (Figure 11).The impact of the suspension flow rate was relatively insignificant regardless of whether x 50 evolution was expressed in terms of the milling time or the effective milling time.Overall, we conclude that the PSD equivalency established for 10.8 min effective milling time (Figure 9) is generally applicable during the whole milling, thus further supporting the use of effective milling time in process scale-up and/or batch size changes.

Simulating the Impact of Imperfect Mixing in the Milling Chamber
Up to this point, the PBM assumed perfect mixing in the milling chamber for a theoretical single pass, which is the expected RTD response for a small mill with a low L/D ratio [33].However, for larger mills with high L/D ratio (>1), the RTD behavior for a single-pass continuous operation deviates from perfect mixing [33,49].The cell-based PBM accounts for this deviation via n > 1 and fine-tunes the RTD behavior via the back-mixing ratio R. Interestingly, for a recirculating WSMM, the impact of n and R (different RTD behavior) on the PSD evolution has not been examined in the milling literature.Hence, in this study, we investigated the influence of n at two different R using simulations: R = 0 (Runs 6-9, no internal recirculation between the cells in the milling chamber) and R = 2.52 (Runs 10-13, strong internal recirculation between cells in the milling chamber).
As depicted in Figure 12, in the absence of recirculation between the cells (R = 0), an increase in the number of cells n resulted in a narrower and sharper PSD initially (at t = 1 min, within the first turnover).Here, the prediction by the perfect-mixing model (n = 1) was included as a reference.This finding is similar to what was observed in a single-pass continuous mill without a holding tank [54].

Simulating the Impact of Imperfect Mixing in the Milling Chamber
Up to this point, the PBM assumed perfect mixing in the milling chamber for a theoretical single pass, which is the expected RTD response for a small mill with a low L/D ratio [33].However, for larger mills with high L/D ratio (>1), the RTD behavior for a singlepass continuous operation deviates from perfect mixing [33,49].The cell-based PBM accounts for this deviation via n > 1 and fine-tunes the RTD behavior via the back-mixing ratio R. Interestingly, for a recirculating WSMM, the impact of n and R (different RTD behavior) on the PSD evolution has not been examined in the milling literature.Hence, in this study, we investigated the influence of n at two different R using simulations: R = 0 (Runs 6-9, no internal recirculation between the cells in the milling chamber) and R = 2.52 (Runs 10-13, strong internal recirculation between cells in the milling chamber).
As depicted in Figure 12, in the absence of recirculation between the cells (R = 0), an increase in the number of cells n resulted in a narrower and sharper PSD initially (at t = 1 min, within the first turnover).Here, the prediction by the perfect-mixing model (n = 1) was included as a reference.This finding is similar to what was observed in a single-pass continuous mill without a holding tank [54].As the milling was prolonged, as illustrated in Figure 13, the impact of n disappeared as early as 8 min, and at the end of 32 min of milling, all PSDs converged to the same PSD.The bulk flow emanating from the holding tank dominates the RTD of the circuit rapidly and any imperfect-mixing effect within the milling chamber was nullified upon an increase in the number of turnovers Nto.Although there is no internal recirculation between the cells (R = 0), the recirculating bulk flow in the circuit from/to holding tank dominates the mixing behavior and the RTD of the circuit [49].Hence, any RTD difference for a single-pass operation of the same mill, as modulated by n and R in the milling chamber, has a diminishing effect on the PSD upon an increase in milling time and number of turnovers.One practical implication of this theoretical finding for pharmaceutical engineers is that performing rigorous, elaborate RTD studies for the continuous stirred mills, typically carried out in single-pass mode at the steady-state, may not be warranted if the same mills are eventually used in the recirculation (circuit) mode.One can use the perfect-mixing model (n = 1) for the recirculating WSMM even though the single-pass operation may exhibit significant deviation of the RTD from that of the perfect mixing.Any deviations for such modeling appear to be restricted to early milling times; however, this finding is As the milling was prolonged, as illustrated in Figure 13, the impact of n disappeared as early as 8 min, and at the end of 32 min of milling, all PSDs converged to the same PSD.The bulk flow emanating from the holding tank dominates the RTD of the circuit rapidly and any imperfect-mixing effect within the milling chamber was nullified upon an increase in the number of turnovers N to .Although there is no internal recirculation between the cells (R = 0), the recirculating bulk flow in the circuit from/to holding tank dominates the mixing behavior and the RTD of the circuit [49].Hence, any RTD difference for a single-pass operation of the same mill, as modulated by n and R in the milling chamber, has a diminishing effect on the PSD upon an increase in milling time and number of turnovers.One practical implication of this theoretical finding for pharmaceutical engineers is that performing rigorous, elaborate RTD studies for the continuous stirred mills, typically carried out in single-pass mode at the steady-state, may not be warranted if the same mills are eventually used in the recirculation (circuit) mode.One can use the perfect-mixing model (n = 1) for the recirculating WSMM even though the single-pass operation may exhibit significant deviation of the RTD from that of the perfect mixing.Any deviations for such modeling appear to be restricted to early milling times; however, this finding is expected to depend on other system parameters especially V s and V m , which were 236 mL and 80 mL in these simulations (Runs 6-17).
Pharmaceutics 2024, 16,353 19 of 27 expected to depend on other system parameters especially Vs and Vm, which were 236 mL and 80 mL in these simulations (Runs 6-17).Subsequently, we explored the impact of n when the axial back-mixing ratio R was set to 2.52 (taken from [49] for a single-pass continuous WSMM).Figure 14 illustrates that, akin to the condition of R = 0, a higher n initially led to a narrower and sharper PSD at t = 1 min.A cursory look at Figures 12 and 14 suggests that the impact of n was slightly more discernible for R = 0 than for R = 2.52.Nevertheless, Figure 15 reveals that as milling progressed, the n effect was nullified, and identical PSDs were obtained for different n at 32 min.While the PSDs at 1 min were notably different for R = 0 and R = 2.52, a comparison of Figures 13 and 15 reveals that the PSDs were similar or identical for other milling times.Subsequently, we explored the impact of n when the axial back-mixing ratio R was set to 2.52 (taken from [49] for a single-pass continuous WSMM).Figure 14 illustrates that, akin to the condition of R = 0, a higher n initially led to a narrower and sharper PSD at t = 1 min.A cursory look at Figures 12 and 14 suggests that the impact of n was slightly more discernible for R = 0 than for R = 2.52.Nevertheless, Figure 15 reveals that as milling progressed, the n effect was nullified, and identical PSDs were obtained for different n at 32 min.While the PSDs at 1 min were notably different for R = 0 and R = 2.52, a comparison of Figures 13 and 15 reveals that the PSDs were similar or identical for other milling times.expected to depend on other system parameters especially Vs and Vm, which were 236 mL and 80 mL in these simulations (Runs 6-17).Subsequently, we explored the impact of n when the axial back-mixing ratio R was set to 2.52 (taken from [49] for a single-pass continuous WSMM).Figure 14 illustrates that, akin to the condition of R = 0, a higher n initially led to a narrower and sharper PSD at t = 1 min.A cursory look at Figures 12 and 14 suggests that the impact of n was slightly more discernible for R = 0 than for R = 2.52.Nevertheless, Figure 15 reveals that as milling progressed, the n effect was nullified, and identical PSDs were obtained for different n at 32 min.While the PSDs at 1 min were notably different for R = 0 and R = 2.52, a comparison of Figures 13 and 15 reveals that the PSDs were similar or identical for other milling times.For a fixed n, e.g., n = 3, the effects of the axial back-mixing ratio R were simulated via Runs 7 and 14-17 (refer to Table 2).An increase in R, which indicates a greater extent of axial back-mixing emanating from the higher rate of internal recirculation between the cells, led to coarser and wider PSD at t = 1 min (Figure 16).This effect disappeared, similar to the effect of n, upon further milling and concomitantly an increase in the number of turnovers (Figure 17).The PSDs at 32, 64, and 128 min at various n and R (different RTD response of the mill for a theoretical single-pass operation) were either similar or identical.Again, this finding points out that regardless of the differences of the RTD response of the continuous mill, which is usually obtained via single-pass operation at the steady-state, the recirculating bulk flow of the suspension dominated the circuit RTD, and the impacts of n and R were insignificant after a few turnovers.In particular, the back-mixing in the mill becomes insignificant in a recirculating WSMM operation after a few passes of the suspension through the mill.For a fixed n, e.g., n = 3, the effects of the axial back-mixing ratio R were simulated via Runs 7 and 14-17 (refer to Table 2).An increase in R, which indicates a greater extent of axial back-mixing emanating from the higher rate of internal recirculation between the cells, led to coarser and wider PSD at t = 1 min (Figure 16).This effect disappeared, similar to the effect of n, upon further milling and concomitantly an increase in the number of turnovers (Figure 17).The PSDs at 32, 64, and 128 min at various n and R (different RTD response of the mill for a theoretical single-pass operation) were either similar or identical.Again, this finding points out that regardless of the differences of the RTD response of the continuous mill, which is usually obtained via single-pass operation at the steady-state, the recirculating bulk flow of the suspension dominated the circuit RTD, and the impacts of n and R were insignificant after a few turnovers.In particular, the back-mixing in the mill becomes insignificant in a recirculating WSMM operation after a few passes of the suspension through the mill.
Although we did not simulate n > 8, we expect that the invariant PSDs observed after 32 min will not change even when higher n values are used.As established in ref.
[54], for a single-pass continuous mill, the cell-based PBM with n = 60 produces a PSD that is identical to the PSD obtained for a plug-flow RTD.Our findings from the current PBM simulations imply that the PSDs in the recirculating WSMM (circuit operation) are invariant to n and R (after a few turnovers), and that they correspond to an equivalent single-pass continuous mill with plug-flow behavior approximately.This finding is in line with the sharpening of the circuit RTD upon an increase in the number of turnovers, and the approach of the circuit RTD to that of the plug flow for an equivalent single-pass mill [49].
response of the mill for a theoretical single-pass operation) were either similar or identical.Again, this finding points out that regardless of the differences of the RTD response of the continuous mill, which is usually obtained via single-pass operation at the steady-state, the recirculating bulk flow of the suspension dominated the circuit RTD, and the impacts of n and R were insignificant after a few turnovers.In particular, the back-mixing in the mill becomes insignificant in a recirculating WSMM operation after a few passes of the suspension through the mill.Although we did not simulate n > 8, we expect that the invariant PSDs observed after 32 min will not change even when higher n values are used.As established in ref.
[54], for a single-pass continuous mill, the cell-based PBM with n = 60 produces a PSD that is identical to the PSD obtained for a plug-flow RTD.Our findings from the current PBM simulations imply that the PSDs in the recirculating WSMM (circuit operation) are invariant to n and R (after a few turnovers), and that they correspond to an equivalent single-pass continuous mill with plug-flow behavior approximately.This finding is in line with the sharpening of the circuit RTD upon an increase in the number of turnovers, and the approach of the circuit RTD to that of the plug flow for an equivalent single-pass mill [49].

On the Grinding Limit and the Transition Particle Size
There are two major reasons as to why the specific breakage rate decreases significantly as particles become smaller (refer to Figure 3): (i) capturing efficiency of sub-micron particles between the relatively large beads is very low, and it becomes lower as the particles get smaller, as predicted by the microhydrodynamic model (see, e.g., [35]) and (ii) the number of defects, flaws, and pre-cracks in particles decreases with particle size, which causes a lower breakage probability for the smaller particles [71].When a brittleductile transition occurs, the grinding limit appears because below that size the relaxation

On the Grinding Limit and the Transition Particle Size
There are two major reasons as to why the specific breakage rate decreases significantly as particles become smaller (refer to Figure 3): (i) capturing efficiency of sub-micron particles between the relatively large beads is very low, and it becomes lower as the particles get smaller, as predicted by the microhydrodynamic model (see, e.g., [35]) and (ii) the number of defects, flaws, and pre-cracks in particles decreases with particle size, which causes a lower breakage probability for the smaller particles [71].When a brittleductile transition occurs, the grinding limit appears because below that size the relaxation mechanism changes from brittle fracture to ductile (plastic) flow [64,72,73].Our estimate of the grinding limit, typically reported as a cumulant size, for drugs is in the range of 30-90 nm, based on previous studies [8,74].In the current study, a transition particle size x* of 195 ± 20 nm (10.2% RSD) was estimated.It is critical to note that x* of the PBM is not equal to the grinding limiting size obtained by milling for extremely long milling times (4-16 h) and/or with intensified process conditions [8,64].Also, it is worth noting that Model C predicts a specific breakage rate range of ~10 -9 -7 × 10 -5 min -1 for the size range of 30-90 nm particles (refer to Figure 3).It is speculated that with such extremely low specific breakage rates (close to 0), 30-90 nm FNB particles will not break extensively within the time scale of the milling experiments explored here, implying that the grinding limit of FNB may indeed fall within this range.To put this into perspective, 30-90 nm particles break ~10 10 -10 5 times slower than 1 µm particles (Figure 3).We will not further explore this point as finding the true grinding limit was not the primary objective of this study.

Limitations of the Current PBM
The majority of the PBM's deviation from the experimental PSD data during the early milling times originated from particle aggregation.The PBM did not consider particle aggregation as it would significantly increase the modeling complexity and the number of parameters.We did not consider the volume of the suspension in the tubing between the holding tank and the milling chamber, and the lag time associated with the flow in the tubing.Moreover, the effects of rotor speed, bead loading, and bead size were not varied in this study, which would significantly affect the A parameter.In that case, A should be multiplicatively decomposed into several power law terms for these variables, and their parameters need to be estimated first.These aspects will be investigated in a forthcoming study.

Conclusions and Outlook
We studied the impacts of batch size, suspension flow rate, and imperfect mixing during the production of fenofibrate nanosuspensions in a recirculating WSMM within the context of a cell-based PBM.Four specific breakage rate functions were discriminated by fitting the PBM to the PSD evolution at the baseline process conditions.It is concluded that Model C, which consists of a product of power-law and logistic functions, is the best model.The parameter estimation study revealed the criticality of incorporating a transition particle size commensurate with the notion of a grinding limit.The fitted PBM was then used to predict the salient features of the impacts of the batch size and the suspension flow rate.Both the experiments and the predictions demonstrated that an increase in batch size entails a proportionate increase in milling time to keep the PSD invariant, and this can be carried out by keeping the same residence time in the mill for the circuit operation.This circuit residence time was regarded as the effective milling time that must be kept identical upon batch size increase and in process scale-up.The suspension flow rate was found to have insignificant impact on the PSD.An intriguing finding from the PBM simulations was that no matter what the single-pass RTD of the mill, as modulated by various n and R of the cell-based PBM, the PSD is not affected by them after a few turnovers in the recirculation mode.This finding practically implies that detailed RTD studies may not be warranted for modeling a recirculating WSMM.As aggregation was present in the system, future modeling effort is warranted to include aggregation terms in the PBM to improve the fittingprediction capability.Also, the PBM should be expanded to account for the impacts of all process parameters.Despite its limitations, the presented PBM offered significant insights about recirculating WSMM, also providing valuable guidance for process development and scale-up.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/pharmaceutics16030353/s1,Section S1: All supporting data about parameter estimation-optimization and evolution of the PSD; Table S1: Fitted parameters of Model A, lower and upper boundaries of the parameters, initial guess, and sum-of-squared residuals SSR for different numbers of trials N T ; Table S2: Fitted parameters of Model B, lower and upper boundaries of the parameters, initial guess, and sum-of-squared residuals SSR for different numbers of trials N T ; Table S3: Fitted parameters of Model D, lower and upper boundaries of the parameters, initial guess, sum-of-squared residuals SSR for different numbers of trials N T ; Figure S1: Temporal variation of the mass fraction density distribution during (a) the first 8 min of milling,

Figure 1 .
Figure 1.A schematic of the recirculation mode of wet stirred media milling (WSMM).

Figure 1 .
Figure 1.A schematic of the recirculation mode of wet stirred media milling (WSMM).

Figure 2 .
Figure 2. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the baseline experiment: volumetric flow rate of 126 mL/min and batch volume of 236 mL (Run 1).Simulation: PBM with Model C and its fitted parameters.

Figure 2 .
Figure 2. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the baseline experiment: volumetric flow rate of 126 mL/min and batch volume of 236 mL (Run 1).Simulation: PBM with Model C and its fitted parameters.

Figure 3 .
Figure 3.The variation of the specific breakage rate Si with the particle size xi according to the fitted Models A-D.Si equals 0 for xi ≤ x* per Model B, which is not shown on this log-log plot.

Figure 3 .
Figure 3.The variation of the specific breakage rate S i with the particle size x i according to the fitted Models A-D. S i equals 0 for x i ≤ x* per Model B, which is not shown on this log-log plot.

Figure 4 .
Figure 4. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the smaller batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 118 mL (Run 2).

Figure 5 .
Figure 5. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the larger batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 472 mL (Run 3).

Figure 4 .
Figure 4. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the smaller batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 118 mL (Run 2).

Pharmaceutics 2024, 16 , 353 13 of 27 SSR
of the fitting of Run 1 data, i.e., 38.0.The suspension flow rate had an insignificant effect on the PSD evolution, as can be seen from a visual comparison inFigures 2, 6 and 7.

Figure 4 .
Figure 4. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the smaller batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 118 mL (Run 2).

Figure 5 .
Figure 5. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the larger batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 472 mL (Run 3).

Figure 5 .
Figure 5. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the larger batch size experiment: volumetric flow rate of 126 mL/min and batch volume of 472 mL (Run 3).

Figure 6 .
Figure 6.Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the lower flow rate experiment: volumetric flow rate of 63 mL/min and batch volume of 236 mL (Run 4).

Figure 7 .
Figure 7. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the higher flow rate experiment: volumetric flow rate of 250 mL/min and batch volume of 236 mL (Run 5).

Figure 6 .
Figure 6.Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the lower flow rate experiment: volumetric flow rate of 63 mL/min and batch volume of 236 mL (Run 4).

Pharmaceutics 2024, 16 , 353 14 of 27 Figure 6 .
Figure 6.Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the lower flow rate experiment: volumetric flow rate of 63 mL/min and batch volume of 236 mL (Run 4).

Figure 7 .
Figure 7. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the higher flow rate experiment: volumetric flow rate of 250 mL/min and batch volume of 236 mL (Run 5).

Figure 7 .
Figure 7. Temporal variation of the mass fraction density distribution during (a) the first 32 min of milling and (b) thereafter for the higher flow rate experiment: volumetric flow rate of 250 mL/min and batch volume of 236 mL (Run 5).

Figure 8 .
Figure 8. Experimentally measured and simulated impacts of the batch size and volumetric flow rate on the mass fraction density distribution at 32, 64, and 128 min.

Figure 8 .
Figure 8. Experimentally measured and simulated impacts of the batch size and volumetric flow rate on the mass fraction density distribution at 32, 64, and 128 min.

Figure 9 .
Figure 9. (a) The experimentally measured and (b) simulated mass fraction density distribution for various batch sizes (Vs = 118 mL, Vs = 236 mL, and Vs = 472 mL in Runs 2, 1, and 3, respectively) at the same effective milling time of 10.8 min.

Figure 9 .
Figure 9. (a) The experimentally measured and (b) simulated mass fraction density distribution for various batch sizes (Vs = 118 mL, Vs = 236 mL, and Vs = 472 mL in Runs 2, 1, and 3, respectively) at the same effective milling time of 10.8 min.

Figure 10 .
Figure 10.The measured and simulated evolution of the 50% passing (median) size x50 as a function of milling time: (a) influence of the flow rate and (b) influence of the batch size (Runs 1-5).

Figure 11 .
Figure 11.The measured and simulated evolution of the 50% passing (median) size as a function of effective milling time: (a) influence of the flow rate and (b) influence of the batch size (Runs 1-5).

Figure 10 . 27 Figure 10 .
Figure 10.The measured and simulated evolution of the 50% passing (median) size x 50 as a function of milling time: (a) influence of the flow rate and (b) influence of the batch size (Runs 1-5).

Figure 11 .
Figure 11.The measured and simulated evolution of the 50% passing (median) size as a function of effective milling time: (a) influence of the flow rate and (b) influence of the batch size (Runs 1-5).

Figure 11 .
Figure 11.The measured and simulated evolution of the 50% passing (median) size as a function of effective milling time: (a) influence of the flow rate and (b) influence of the batch size (Runs 1-5).

Figure 12 .
Figure 12.The effects of number of cells n at the back-mixing ratio R of 0 on the mass fraction density distribution at t = 1 min.

Figure 12 .
Figure 12.The effects of number of cells n at the back-mixing ratio R of 0 on the mass fraction density distribution at t = 1 min.

Figure 13 .
Figure 13.The effects of number of cells n at the back-mixing ratio R of 0 on the mass fraction density distribution at t = 8 min, t = 32 min, and t = 128 min.

Figure 14 .
Figure 14.The effects of number of cells n at the back-mixing ratio R of 2.52 on the mass fraction density distribution at t = 1 min.

Figure 13 .
Figure 13.The effects of number of cells n at the back-mixing ratio R of 0 on the mass fraction density distribution at t = 8 min, t = 32 min, and t = 128 min.

Figure 13 .
Figure 13.The effects of number of cells n at the back-mixing ratio R of 0 on the mass fraction density distribution at t = 8 min, t = 32 min, and t = 128 min.

Figure 14 .
Figure 14.The effects of number of cells n at the back-mixing ratio R of 2.52 on the mass fraction density distribution at t = 1 min.

Figure 14 .
Figure 14.The effects of number of cells n at the back-mixing ratio R of 2.52 on the mass fraction density distribution at t = 1 min.

Figure 15 .
Figure 15.The effects of number of cells n at the back-mixing ratio R of 2.52 on the mass fraction density distribution at t = 8 min, t = 32 min, and t = 128 min.

Figure 15 .
Figure 15.The effects of number of cells n at the back-mixing ratio R of 2.52 on the mass fraction density distribution at t = 8 min, t = 32 min, and t = 128 min.

Figure 16 . 27 Figure 17 .
Figure 16.The effects of back-mixing ratio R for three cells on the PSD at t = 1 min.Figure 16.The effects of back-mixing ratio R for three cells on the PSD at t = 1 min.Pharmaceutics 2024, 16, 353 21 of 27

Figure 17 .
Figure 17.The effects of back-mixing ratio R for three cells on the PSD at t = 8 min, t = 32 min, and t = 128 min.

Table 1 .
Conditions of the milling experiments.
Run No.Identifier Batch Size, Vs (mL) Volumetric Flow Rate, Q (mL/min)

Table 1 .
Conditions of the milling experiments.
Identifier Batch Size, V s (mL) Volumetric Flow Rate, Q (mL/min)

Table 2 .
Various simulation parameters used in the cell-based PBM.

Table 3 .
Fitted parameters of Model C, lower and upper boundaries of the parameters, initial guess, and sum-of-squared residuals SSR for different numbers of trials N T .

Table 4 .
The estimated breakage parameters for different specific breakage functions and associated SSR values (N T = 1000).