Development of a Semi-Mechanistic Modeling Framework for Wet Bead Milling of Pharmaceutical Nanosuspensions

This study aimed to develop a practical semi-mechanistic modeling framework to predict particle size evolution during wet bead milling of pharmaceutical nanosuspensions over a wide range of process conditions and milling scales. The model incorporates process parameters, formulation parameters, and equipment-specific parameters such as rotor speed, bead type, bead size, bead loading, active pharmaceutical ingredient (API) mass, temperature, API loading, maximum bead volume, blade diameter, distance between blade and wall, and an efficiency parameter. The characteristic particle size quantiles, i.e., x10, x50, and x90, were transformed to obtain a linear relationship with time, while the general functional form of the apparent breakage rate constant of this relationship was derived based on three models with different complexity levels. Model A, the most complex and general model, was derived directly from microhydrodynamics. Model B is a simpler model based on a power-law function of process parameters. Model C is the simplest model, which is the pre-calibrated version of Model B based on data collected from different mills across scales, formulations, and drug products. Being simple and computationally convenient, Model C is expected to reduce the amount of experimentation needed to develop and optimize the wet bead milling process and streamline scale-up and/or scale-out.


Introduction
Wet bead milling (also known as wet stirred media milling) is a unit operation used in the pharmaceutical industry for the preparation of suspension-based products [1].A recent survey indicates that wet bead milling is the preferred approach for the preparation of ultrafine drug suspensions and nanosuspensions compared to other techniques, such as liquid antisolvent precipitation and high-pressure homogenization [2].This is not unexpected since wet bead milling has several advantages as it is a robust, reproducible, scalable, organic solvent-free, and environmentally friendly process [3,4].It enables the preparation of concentrated stable suspensions of drug particles [3], which can have several applications in drug delivery such as modulating drug dissolution and absorption [5][6][7][8][9] and the design of long acting injectables (LAI) [10][11][12][13].
Developing the fundamental mechanistic understanding of manufacturing processes, underpinned by a science-and risk-based approach, is a key element of the Quality-by-Design (QbD) framework of product development, aligned with the expectations outlined in the International Council on Harmonization (ICH) guidance Q8 (R2), Q9, and Q10 [14].In this work, we have focused our efforts on developing an enhanced level of process understanding of wet bead milling to support the development of the predictive design capability for this unit operation based on process modeling.There is a need to enhance the mechanistic understanding of potential process-related challenges encountered during wet bead milling, guided by the QbD principles of development.These failure modes may include poor milling efficiency due to aggregation and Ostwald ripening [15][16][17], long processing times [3,18,19], bead wear and product contamination [20,21], potential chemical degradation [22], and mechanical or thermal-stress-induced solid-state transformations and generation of higher energy metastable phases during milling, such as amorphization and polymorphic transitions [22,23].
The modeling of the wet bead milling process can significantly help process development and optimization [24,25] and offers multiple benefits; at a minimum, mechanistic or first-principle-based models provide a quantitative, fundamental understanding of the impact of operation-design parameters.Since APIs are typically very expensive, there is a need to predict their behavior in the unit operations prior to committing a significant amount of material in the developmental stage.Evaluating the behavior of a formulation in wet bead milling and establishing a design space typically requires significant investment in time and material.Computational models for milling performance would be beneficial not only because they reduce API quantities needed for process development but also because they would help inform product development teams in advance regarding whether a process is likely to achieve reasonable throughputs at production scale.
Interestingly, statistically based models such as empirical regression fits, response surface methodology (RSM), etc., have been overwhelmingly preferred over the mechanisticphenomenological models according to a review of the wet bead milling modeling in the pharmaceutical nanotechnology literature [25].This is not surprising: (i) Except the SI-SN [26,27] and the MHD models [28,29], the mechanistic-phenomenological models are computationally expensive, requiring specialized software and expertise.(ii) The usage of the SI-SN and the MHD models can be limited if reliable information on model inputs such as average power consumption during milling, the apparent shear viscosity, and the density of the suspension are not available.(iii) Except for PBM, no complex models such as CFD [30,31] or DEM [32][33][34] consider the evolving PSD due to particle breakage.And (iv) PBMs must incorporate the operation-design parameters into the kernels for process predictions, which makes model calibration difficult and requires computationally intensive global optimizers for parameter estimation.On the other hand, statistically based models are easier to develop and use and are more accessible to pharmaceutical engineers and scientists.Characteristic drug particle sizes (x 50 and x 90 ) or specific surface area can be described and/or predicted by RSM and regression analysis as a function of formulation and process parameters.Most of the statistically based empirical studies have correlated the particle sizes with the milling time and speed [35][36][37][38]; only a handful of studies have additionally considered bead size [39,40], bead loading [41], or API loading [42].Other studies have considered only bead loading [43,44] or milling time and bead loading [45,46].Finally, empirical breakage kinetic models such as the first-order kinetics, n th -order kinetics, and warped-time model have been used to describe the timewise evolution of the median drug particle size or specific (external) surface area [35,47,48].An empirical correlation for the breakage rate constant of the first-order kinetic model based on stirrer speed and a scaling factor was developed in [35], but no direct scale-up was demonstrated.Moreover, such empirical correlations are only applicable to the specific process, unless its application is demonstrated for different case studies.The breakage rate constant of the n th -order model was successfully correlated with the MHD parameters to a lab-scale mill; however, explicit correlation of the MHD parameters on the process parameters has not been established, and scale-up has not been considered [48].
The take-away of the above analysis of the statistically based empirical models is that despite the obvious need for pharmaceutical manufacturing engineers to have a good model for prediction and optimization of wet bead milling processes across scale, incorporating the impacts of all process-formulation-design parameters, it is not easy to develop such a model.Besides being practical and computationally convenient, it should also predict impacts of changes in process parameters and scale-up or scale-out across equipment.One such model [47] described the evolution of the specific surface area using a first-order kinetic model with a characteristic time constant, which was correlated with the rotation speed, bead loading, and API loading.The time constant is purely empirical with no physical (microhydrodynamic) basis, and it does not include the effects of bead properties, which have a significant impact on breakage kinetics, cycle time, and energy consumption [20,29].
The aim of this study is to develop a practical semi-mechanistic modeling framework that is computationally efficient, easy to use, and convenient to calculate and that can be used as an engineering tool during process development.The scope of this work included an analysis of six internal GSK drug products and two drug products at NJIT as milled by six mills across six different scales, with an expected outcome that a common model can be used to model every scale and drug product.A challenge to creating a widely encompassing model for wet bead milling is the high number of potential parameters that are important.Developing a practical model with a large number of independent variables, including process parameters, formulation parameters, and equipment-specific parameter was an ambitious undertaking.To this end, in this work, we developed three models with different levels of complexity and mechanistic rigor.First, we linearized the reciprocal of the characteristic particle size quantiles (x 10 , x 50 , x 90 ) to describe the timewise particle size evolution.Then, we used microhydrodynamic theory [20,29] to derive the general functional form of the apparent breakage rate constant, which accounts for the impact of all process parameters and bead properties.The apparent breakage rate constant was then derived as a linear function of the frequency of drug particle compressions between the beads to obtain Model A. As the most general and demanding model, Model A requires one to calculate the granular temperature, which can be estimated by solving the power dissipation equation in the microhydrodynamic model [25].Model B was developed by assuming a general power-law expression for the granular temperature.This simplification has general applicability and great practical utility as we can obtain the exponents of the process parameters in Model B by fitting it to experimental data.Here, we fit Model B to several case studies conducted by GSK and NJIT on various stirred media mills of multiple scales, which resulted in common exponents for each process parameter, resulting in the simplest model, i.e., Model C. Overall, we have formulated a modeling framework encompassing models with varying complexity and practical utility, and we hope that such models can be adopted and/or adapted by engineers for optimization, scale-out, and scale-up of wet bead milling processes.

Materials
The APIs used for investigation at NJIT were fenofibrate (BP grade, purchased from Jai Radhe Sales Ahmedabad, India) and griseofulvin (BP/EP grade micronized, Letco Medical, Decatur, AL, USA).The drug product formulations at GSK are proprietary APIs that will not be listed explicitly but simply referred to as a drug product (DP) number, i.e., DP1, DP2, . . .DP8.Both NJIT studies used Hydroxypropyl cellulose (HPC, L grade, Nisso America Inc., New York, NY, USA) as a nonionic polymeric stabilizer, and sodium dodecyl sulfate (SDS, ACS grade, GFS chemicals, Columbus, OH, USA) as an anionic surfactant where the formulation was 10% drug, 7.5% HPC-L and 0.05% SDS with respect to 200 g DI water [49,50].In addition, 400 µm nominal sized Zirmil Y grade zirconia beads and HCC grade polystyrene beads were purchased from Saint Gobain ZirPro (Mountainside, NJ, USA) and Norstone Inc. (Bridgeport, PA, USA), respectively.Additionally, GSK formulations contained commonly used excipients such as surfactants, stabilizers, etc.

Experimental Setup
Each of the experiments described in this manuscript involved wet bead milling in one or more of six different mills, purchased from Netzsch Fine Particle Size Technology, LLC (Exton, PA, USA).A tabulated list of the relevant mill characteristics and parameter ranges explored is shown in Table 1.The experiments conducted on each mill varied substantially in the process parameter operating space.A total of 64 runs were executed in the mills, respectively, across 8 different drug products.

Particle Size Measurement
For particle size measurements performed at GSK, the PSD of the drug suspension at various milling times was determined using laser diffraction by Malvern Mastersizer 3000 particle size analyzer (software v3.81 or validated equivalent) with Hydro MV Dispersion Unit and Temperature Control Unit.At NJIT, the particle size distribution (PSD) of the drug suspensions at various milling times was determined by LS 13-320 Beckman Coulter instrument (Brea, CA, USA).While the sampling interval varies for each experiment depending on the total milling time at GSK, predefined time intervals (2 s, s = 0, 1, 2, . . .7 min) with the addition of 40 s, 24 min, 48 min, 96 min, 128 min and 180 min were used at NJIT.
For measurements performed in GSK, ~2 mL samples of suspension were taken from the bulk holding tank at predefined time intervals.Two drops from each sample were transferred via an 18G needle (or equivalent) to a microcentrifuge tube containing ~1 mL of water for injection (WFI).The sample with WFI dilution was then mixed to homogeneity.During measurements, obscuration was maintained at 4-6%.Measurements were repeated three times, and the average and standard deviation of these measurements were determined.The Malvern instrument was turned on for no less than 30 min prior to use.
For the particle size measurements performed at NJIT, the samples were taken from the mill outlet at predefined time intervals.The final sample was taken from the holding tank and all samples were measured with laser diffraction.Before each measurement, a suspension sample (~2.0 mL for griseofulvin and ~1.0 mL for fenofibrate) was diluted with 5.0 mL of the respective vehicle using a vortex mixer (Fisher Scientific Digital Vortex Mixer, Model No: 945415, Pittsburgh, PA, USA) at 1500 rpm for a minute.During measurements, polarized intensity differential scattering (PIDS) was maintained between 40% and 50%, while the obscuration was maintained below 8%.PSD was provided by the equipment software, which used the Mie scattering theory.The refractive indices of GF, FNB and water were taken as 1.65, 1.55, and 1.33, respectively [49,50].Measurements were repeated four times, and the average and standard deviation of these measurements were determined.

Young's Modulus of Compacts
A compaction simulator (Styl'One Evolution, Medelpharm, Beynost, France) was used to produce compacts of each API (with no excipients) with a round flat-faced 11.28 mm B punch at a target compression speed.The compression cycle included a precompression phase followed by main compression.The methodology reported by Mazel et al. [51] was used in this work to measure the Young's modulus (YM) of the bulk API materials based on compaction analysis.To obtain the Young's modulus as a function of compact porosity, a range of main compression forces were applied to the powders in the die.For each compaction, ~500 mg of API was used.The tablet weight and thickness values were measured, and the out-of-die porosity of each compact was calculated using ε = mc Vc ρ , where m c and V c are the weight and volume of the out-of-die compact, and ρ is the true density of API powder [51].
For measuring YM, the compact must only be undergoing elastic deformation.This corresponds to the linear part of compression pressure-tablet thickness curve.During elastic deformation, the linear stress-strain relationship can be written as , where σ ax and σ rad are axial and radial stresses equal to the pressure levels (P ax is a mean value of lower and upper punch pressures, and P rad is the radial pressure to the die wall), PR is Poisson's ratio, h is the thickness of the powder bed, and h 0 is the initial thickness.Plotting (σ ax − 2PRσ rad ) vs. h provides a linear relationship with YM as its intercept, and YM can be obtained at different pressure levels that yield different compact porosities [51].
An exponential decay function was fit to YM versus porosity data, and the YM at the zero compact porosity was estimated for the model compounds in this study.

Theoretical
The approach to model the timewise particle size evolution during wet bead milling entails first establishing a data transformation to linearize the characteristic particle size quantiles x 10 , x 50 , and x 90 of the cumulative PSD (Section 3.1) and then developing a model to predict the timewise evolution of linearized particle sizes.
The general model training workflow is presented in Figure 1.
tank and all samples were measured with laser diffraction.Before each measurement, a suspension sample (~2.0 mL for griseofulvin and ~1.0 mL for fenofibrate) was diluted with 5.0 mL of the respective vehicle using a vortex mixer (Fisher Scientific Digital Vortex Mixer, Model No: 945415, Pittsburgh, PA, USA) at 1500 rpm for a minute.During measurements, polarized intensity differential scattering (PIDS) was maintained between 40% and 50%, while the obscuration was maintained below 8%.PSD was provided by the equipment software, which used the Mie scattering theory.The refractive indices of GF, FNB and water were taken as 1.65, 1.55, and 1.33, respectively [49,50].Measurements were repeated four times, and the average and standard deviation of these measurements were determined.

Young's Modulus of Compacts
A compaction simulator (Styl'One Evolution, Medelpharm, Beynost, France) was used to produce compacts of each API (with no excipients) with a round flat-faced 11.28 mm B punch at a target compression speed.The compression cycle included a precompression phase followed by main compression.The methodology reported by Mazel et al. [51] was used in this work to measure the Young's modulus (YM) of the bulk API materials based on compaction analysis.To obtain the Young's modulus as a function of compact porosity, a range of main compression forces were applied to the powders in the die.For each compaction, ~500 mg of API was used.The tablet weight and thickness values were measured, and the out-of-die porosity of each compact was calculated using  = , where mc and Vc are the weight and volume of the out-of-die compact, and ρ is the true density of API powder [51].
For measuring YM, the compact must only be undergoing elastic deformation.This corresponds to the linear part of compression pressure-tablet thickness curve.During elastic deformation, the linear stress-strain relationship can be written as ( − 2 ) =  −  , where σax and σrad are axial and radial stresses equal to the pressure levels (Pax is a mean value of lower and upper punch pressures, and Prad is the radial pressure to the die wall), PR is Poisson's ratio, h is the thickness of the powder bed, and h0 is the initial thickness.Plotting ( − 2 ) vs. h provides a linear relationship with YM as its intercept, and YM can be obtained at different pressure levels that yield different compact porosities [51].
An exponential decay function was fit to YM versus porosity data, and the YM at the zero compact porosity was estimated for the model compounds in this study.

Theoretical
The approach to model the timewise particle size evolution during wet bead milling entails first establishing a data transformation to linearize the characteristic particle size quantiles x10, x50, and x90 of the cumulative PSD (Section 3.1) and then developing a model to predict the timewise evolution of linearized particle sizes.
The general model training workflow is presented in Figure 1.

Particle Size Data Transformation for Linearization
One of the challenges of modeling wet bead milling is the highly nonlinear and potentially sigmoidal nature of the timewise variation in the characteristic particle sizes of x 10 , x 50 , and x 90 .Like most types of models, there is a substantial benefit to simplifying the model by linearizing the input data with respect to time.This section describes the data transformation that we employed for data linearization.While other data transformation strategies may also be useful, we justified the transformation here through two different theoretical approaches in the milling literature and derived them based on breakage kinetics models (Appendix A).
For drug products where the drug crystals tend to exist as individual crystals, we applied model fits to the complete dataset including the initial particle sizes.However, we disregarded a few data points in the neighborhood of time t = 0 [52,53] for the drugs that tend to agglomerate, which can be justified by the following considerations.First, the initial specific breakage rate of the coarse particles and any drug clusters (agglomerates) is so high that it, when compared with that of the fine particles and nanoparticles, may not be adequately captured by simple empirical kinetic models, even if they account for the size dependence of the specific breakage rate.If agglomerates are initially present, they are broken down in the mill during the first few turnovers.Here, the number of turnovers N t is defined by N t = Qt/V tb , where Q is the pumping rate (suspension flow rate) and V tb is the total batch volume.Since crystal agglomeration is subject to significant variation even within the same batch from crystallization, these initial particle size data may be subject to noise (see also the discussion in [35]), making these not meaningful for the prediction of the behavior after completion of the first few turnovers.Thus, it follows that initial breakage of coarse particles and deagglomeration of clusters should be omitted since the main focus of wet bead milling process development is the breakage of particles into sub-micron size (i.e., the behavior after the first few turnovers).Also, from a theoretical standpoint, the transform we used is supported by an analytical solution of a self-similar PBM [54] away from the initial condition.
After discarding the particle size data for the first few turnovers, the second step in the data transformation strategy was to take the inverse of x 10 , x 50 , and x 90 .Our experience suggests that this part of the transformation process mostly linearizes the data, though some curvature with respect to time and sigmoidal behavior is expected if the milling process is pushed to near the apparent grinding limit.When the milling process is not nearing the grinding limit and no sigmoidal behavior is observed, then only the third step of the transformation is needed.In this third step (see Equation ( 1)), a shape factor N j was introduced for each inverse particle size quantile, which minimizes the error when a line is fitted through the data with respect to time.The proposed transform Yt of the j th size quantile, which can be derived by a semi-empirical n th -order rate-based model (see Appendix A), is expressed as in Equation (1): where j = 10, 50 and 90, and k j is the apparent breakage rate constant that varies with the process conditions.Here, B j is a parameter that was obtained by fitting to the dynamic data.From a theoretical standpoint, it is related to the initial feed PSD prior to milling.B j equals either the value of (1/x j (0)) N j of the actual feed PSD at t = 0 or its value for a theoretical feed PSD that is self-similar with the asymptotic self-similar size distribution at longer milling times (see the discussion on the self-similar solution of a PBM for milling in Appendix A).However, B j was ultimately obtained from fitting because it is not expected to satisfy the actual initial condition due to the issues of using the initial time points, as discussed in the previous section.It is assumed that both N j and B j are invariant to the processing conditions and that the kinetic influence of all process variation is captured by k j .Though the model will have higher fidelity if N j and B j are separately fitted to the particle size evolution of each batch/experiment with different processing conditions, such a model will have little predictive capability and will not be useful.Hence, the success of our model largely rests upon how well k j captures the process variations across scales.We note, however, that the fitted values of N j for each particle size quantile can and likely will be different from each other when fitted separately.Equation ( 1) is supported by two theoretical approaches which are self-similar solution of PBM and microhydrodynamic model along with Charles' energy-average particle size relationship (see Appendix A).
Equation (1) and the theoretical approaches suggest x j → 0 and 1/x j → ¥ as t → ¥, which may not pose a serious problem unless the milling is prolonged to approach an apparent grinding limit size x j,inf [55].Thus, if the data for 1/x j display sigmoidal behavior at long milling times, the best approach is to fit a transform that is similar to the "Prout-Thompkins" transform utilized in auto-catalyzed solid state kinetics [56] as shown in Equation (2): Equation ( 2) requires two parameters to describe a sigmoid: the exponent N j describes aspects of the nonlinear behavior with respect to time, while 1/x j,inf describes the value of 1/x j that could be achieved if milling were carried on for an infinite time.Away from the asymptote which emerges after prolonged milling, the expression 1/x j,inf − 1/x j at 1/x j,inf is valid, and Equation (2) reduces to Equation (1) following some algebraic manipulations.In Appendix A, we also present a derivation of Equation ( 2) based solely on breakage kinetics without any transforms.

Model Development for the Apparent Breakage Rate Constant
We now develop models for the functional dependence of the apparent breakage rate constant k j on the operational-design parameters of wet bead milling and bead properties.We describe here three models with varying levels of complexity.We start by describing the development of Model A, which is a mechanistic model derived from microhydrodynamic theory (Section 3.2.1)that requires significant effort to parameterize.Following that, the development of a more computationally efficient semi-mechanistic Model B is described (Section 3.2.2).Model B is a practical, generalizable, semi-mechanistic model that can be used to describe the wet bead milling process development across mills and scales.Finally, a semi-mechanistic model, fitted to the milled particle size data covering the model compounds in this study (Model C), is presented in Section 3.2.3.Model C is a simplified variant of Model B requiring minimum number of milling experiments to parameterize the model (Section 3.2.3).This latter approach (Model C) offers a computationally and experimentally efficient modeling framework compared to fully mechanistic approaches, making it attractive for industrial simulations within timescales of interest.A high-level comparison of the modeling approaches described in this work is presented in Table 2.With these three models, we can effectively describe particle size evolution during wet bead milling for multiple drug products across mill scales.

Model A: A Microhydrodynamics-Based Model
As small molecule organic materials are relatively brittle and easy to break compared to inorganic materials such as ores and minerals, their breakage kinetics are expected to be governed by the stressing frequency rather than the stress intensity, unless the latter is extremely low.Several in-depth studies using the microhydrodynamic (MHD) theory have concluded that the breakage kinetics during wet bead milling are governed by the average frequency of drug particle compression a [29,48].Hence, our starting point in developing this model was to assume that the extent of milling k j t for a given drug nanosuspension formulation processed in a recirculation mill is given by k j t = A * j N t aτ m , where A * j is a constant dimensionless parameter that is presumed to correlate with the brittleness of APIs [57,58], N t is the number of turnovers defined earlier, and τ m is (the single-pass) mean residence time of the suspension in the mill with volume V m containing a true volume fraction c of beads in the milling chamber.Inserting the definitions of N t , τ m = V m (1 − c)/Q into Equation (2) and noting γ API = m API /V tb , we arrived at Equation (3): where γ API was the mass concentration of the drug in g/mL and 1000 was a unit conversion factor since m API was in kg.Equation (3) shows that the apparent breakage rate constant k j equals the pre-factor in front of t; one may consider A * j a as the true breakage rate constant because V m (1 − c)t/V tb appears to be the effective milling time in recirculation milling, also known as mean residence time of the circuit (mill and holding tank).According to the microhydrodynamic theory [28,29,59], the average frequency of drug particle compressions between the beads a is calculated by multiplying the probability of a single drug particle to be caught by the beads p and the average oscillation frequency of a single bead ν, and is a function of bead material density ρ b , Poisson's ratio PR b , Young's modulus YM b , and diameter D b , as well as the granular temperature θ, which is dependent on both operation-design variables and drug suspension properties such as apparent shear viscosity (refer to Appendix B for the principal equations of the microhydrodynamic theory).Inserting a (refer to Equation (A13) in Appendix B) into Equation (3), we obtain Equation ( 4), or Model A, which is the most general microhydrodynamics-based model for the prediction of the particle size evolution during wet bead milling.
Determining the granular temperature θ for different operational and design conditions within the context of the microhydrodynamic theory entails significant effort and is outside the scope of this manuscript.The calculation of θ for a lab-scale mill for various drugs, stirrer speeds, bead loadings, bead types and sizes, etc., is described in earlier MHD studies [16,20,29,48].However, this approach may not be practical as it requires accurate measurements of power P and power density P v = P/V m in mills at various scales, as well as density and apparent shear viscosity of the suspensions.Previous studies showed that the granular temperature θ varies with the rotational speed ω of the rotor [20,29,[48][49][50]60], the volume fraction of the beads in the suspension c [20,29,[48][49][50], the density of the beads ρ b [29,48,49], the diameter of the beads D b [20,50,60], and the apparent shear viscosity of the suspension µ s at the reference temperature [16].The density of the suspensions ρ s is another factor that has an impact on θ (Appendix B).While the comparative impacts of the process variables and bead type/size on θ were established in the aforementioned studies, a mathematical expression to calculate θ directly based on process variables has not yet been formulated.

Model B: A Practical Semi-Mechanistic Model for Wet Bead Milling
In this section, we approach the wet bead milling modeling from a semi-mechanistic standpoint describing the development of semi-mechanistic Model B. In Section 4, we will discuss fits to Model B and specific case studies which were used to support the structure of this model.Influenced by Model A and prior knowledge on parameters with potential impact on the granular temperature as discussed in Section 3.2.1,we formulated a semi-mechanistic model (Model B) with a power-law structure as described in Equation ( 5).
Model B incorporates all the important process parameters governing milling.Similar to Model A, it is a function of the API mass m API and concentration γ API , the mill volume V m , rotational speed ω, the bead loading BL (BL = c/c lim ), the density of the beads ρ b , the Poisson's ratio of the beads PR b , the Young's Modulus of the beads YM b , the diameter of the beads D b , the apparent shear viscosity of the suspension at the reference temperature µ s , and the density of the suspensions ρ s .While the terms with parameter groups (m API , γ API , V m ) and (PR b , YM b ) have been transferred directly from Model A, other parameters were raised to power N. Since these quantities can span orders of magnitude, we normalized each of these terms by a standard reference quantity to ensure that each term remains around an order of magnitude of one.Two additional parameters (Kt and E) were included in Model B. We introduced Kt as a factor that accounts for temperature change from a reference temperature (15 • C) during milling, which is an Arrhenius function with a K2 j parameter to be fitted exp −K2 j .To account for different milling efficiencies across different pieces of milling equipment, we also introduced a mill-scale efficiency factor, E, which is further discussed in Section 4.3.2.E accounts for the energy transfer efficiency of any mill with respect to a reference mill, which captures the impact of mill design differences at different scales; it is constant for a given mill and independent of the API formulation.Model B can be fitted to particle size data from different milling experiments with varying process conditions to obtain the corresponding unknown parameters (i.e., A j , N 1 , N 2 , N 3 , N 4 , N 5 , N 6 , N 7 , K2 j , B j ).The reference values for these parameters were set as 1000 rpm, 0.63 (the packing limit of the beads), 6000 kg/m 3 (the density of the YSZ beads), 0.33 mm (a commonly used bead size in our studies), and 1.13 cP (viscosity of water at 15 • C), respectively.The unit for A j in Model B depends on the exponent of the tip speed N 5 (kg 0.4 m −N5−0.4s − 1.8 ).The connectivity between Equations ( 3) and ( 5) is explained in Appendix C.
To ensure that θ remains invariant upon scale-up from a reference mill, the power density P v must remain invariant according to the microhydrodynamic theory.The asymptotic scaling analysis in Appendix B suggests a scaling of the form θ ∝ P 0.7 v for fully (upper) turbulent flow.In the upper turbulent flow regime, for which Reynolds number Re > 2 × 10 5 , the following scaling is applicable [61]: , where L is the effective length of the mill chamber.This scaling led us to θ ∝ ω 2.1 D 2.1 a /L 0.7 and θ 0.9 ∝ ω 1.9 D 1.9 a /L 0.6 ; therefore, as per Equation ( 5) to Equation ( 6), N 1 = N 5 emerged; thus, we fitted only N 1 -N 4 independently.As will be shown below, fitting data across scales yielded N 1 = N 5 ∼ = 2, supporting the microhydrodynamic result of N 1 = N 5 = 1.9.For the lower turbulent flow regime (3.5 × 10 4 < Re < 2 × 10 5 ), using the same scaling θ ∝ P 0.7 v , we found N 1 = 1.8 and N 5 = 1.6, which are not far from the fitted values.It must be noted that in the lower turbulent flow regime, θ may follow P v 0.7 -P v 0.8 scaling as some minor viscous effects emerge (for laminar flow θ ∝ P v scaling applies).When θ ∝ P 0.8 v was used, N 1 = 2.0 and N 5 = 1.9 were obtained for the lower turbulent regime.The conclusion from this analysis, upon consideration of all complexities of the fluid flow in the mill and the assumptions made in the asymptotic scaling analysis, is that the microhydrodynamic model suggests

Model C: Semi-Mechanistic Milling Model Fitted to Milled Particle Size Data
Equation ( 5) was used along with the Yt transform in Equation ( 2) to fit the reciprocal of x 10 , x 50 , and x 90 for various model compounds in this work and estimate the model parameters, i.e., N 1 -N 4 , A j , B j , N j and x j,inf .The resultant equation is given in Equation ( 6), which we refer to as Model C, where exponents N 1 -N 4 were determined as described in Section 4. While Model C can be applied towards future datasets without re-fitting the exponents, unlike Model B, making it the most computationally and experimentally efficient model, it bears the signatures of the specific mills and formulations used in this study.In Model C/Equation (6), A j , B j , N j , x j,inf , and K2 j are the fitting parameters and vary with the size quantiles, API, and formulation.Parameters were obtained through nonlinear fitting.
If multiple mill scales were used during the Design of Experiments (DoE) and the efficiencies of the mills are unknown, Equation ( 6) can be used with a match function for the mill-scale efficiency factor E term, which can be fitted to the experimental data to estimate efficiencies of different mills.This procedure will be discussed in Section 4.3.The efficiencies for DV50, DV150, DV300, DV2000 and DV4000 mills were estimated as 0.70, 1.00, 0.83, 0.24 and 0.19, respectively.If readers are interested in estimating the efficiencies for mills that are different from the specific Netzsch mills used in this study (see Table 1) but have similar chamber and impeller geometry, they can apply the empirical relationship shown in Equation (7) to estimate E (refer to Section 4.3 for more details).Here, log E in Equation ( 7) is a function of the agitator diameter D a and length L a .
In addition, if milling was not run for too long to approach the apparent grinding limit (i.e., a significant fit for x j,inf cannot be obtained), readers can estimate x j,inf via an empirical relation (Equation ( 8)) with coefficient C j as 0.057, 0.088 and 0.125 for x 10 , x 50 and x 90 , respectively (refer to Section 4.3 for more details).
x j,inf = C j BL −0.25 D 0.12 b YM 0.24 API (8) With the additional flexibility provided by Equations ( 7) and ( 8), Model C (Equation ( 6)) is a convenient model that can be parameterized with as few as one experiment for a specific drug product and formulation.It can then predict the particle size and milling time required to reach a target particle size for any process condition across milling scales (from lab to pilot or commercial scales).

Reversing the Transform
When Equation ( 1) was used as the linearizing transform, reversing the transform was straightforward: the Yt function was raised to the power of (1/N j ) to then calculate 1/x j at different time points.When Equation ( 2) was used as the linearizing transform, Equation ( 9) was used to calculate the reciprocal particle size: Using statistical analysis software with nonlinear solving capabilities like JMP version 17, it is possible to substitute Equation (4), Equation ( 5), or Equation ( 6) for the Yt term in Equation ( 9) and directly solve for the reciprocal particle size.In this study, Model B (Equation ( 5)) was fitted in JMP in this way with reasonable guesses for each parameter as starting points.The reverse-transformed Equation (5), with reference values substituted following Equation (9), would be as shown in Equation (10).
The term Kt j in the above equation was derived as an Arrhenius equation in Section 3.2.2,but in the interest of brevity, this substitution is not shown explicitly.In our studies, we fitted Equation (10) directly to the data from our case studies, solving for the product-specific terms A j , B j , N j , x j,inf and the temperature term K2 j for each size quantile and the nonlinear exponents (N 1 − N 4 ) on ω, BL, D b , and ρ b , which were shown to be consistent in value across all studies and thus product independent, as will be illustrated in Sections 4.1 and 4.2.We note that not all power-law exponents can be estimated in every case study, as the relevant process parameter might not have been investigated in the experimental design, in which case the associated exponents were fixed to those in Equation ( 6) (Model C).A model of this complexity required several case studies to establish validity across numerous products and five different scales as described Sections 4.1 and 4.2.

Results and Discussion
In Section 4.1, an overall summary of the Model B (Equation ( 5)) fits and the impact of process parameters on the particle breakage will be discussed.In support of Section 4.1, Section 4.2 will illustrate the specific examples from case studies involving different drug products.Finally, Section 4.3 will show how Model C was used to derive mill efficiency correlations based on mill properties (Equation ( 7)); grinding limit correlations based on process and material properties (Equation ( 8)); and milling rate correlations based on material properties.

Impact of Process Parameters on Particle Size: Journey from Model B to Model C
The exponents of Model B were found by fitting Equation ( 5) to particle size data from a number of wet bead milling studies, but since not all studies had variations on all w, BL, D b , and ρ b process parameters, Model C was obtained using insights from all case studies collectively.Some of the process parameters, such as rotor speed and bead loading, were investigated in multiple case studies.The exponents obtained for the different particle size quantiles (x 10 , x 50 , and x 90 ) have some variation.In each of the subsequent subsections, we summarize the fitted values of the N 1 , N 2 , N 3 and N 4 exponents.When parameterizing Model C, we only recorded the fitted values of the exponents if all parameter estimates were statistically significant in the model.If the coefficient estimates were more than two standard deviations away from zero, the estimates were identified as statistically significant.

Impact of Rotor Speed
The rotational rotor speed ω (in rpm) is one of the most important and impactful process parameters [32,41] for wet bead milling.It increases the breakage rate, as it provides more frequent and stressful collisions [48].Our microhydrodynamic analysis suggests that ω should have an exponent of 2. We tested this hypothesis in the case studies described in Sections 4.2.1-4.2.3 with Model B (Equation ( 5)).A summary of the findings for all studies is presented in Table 3, which shows N 1 ∼ = 2 is the most frequent exponent, which agrees well with an earlier scale-up study in Ref. [47].The value of 1.39 found for DP1 could be due to a tendency to aggregate at certain mill scales, which takes longer than expected to break down.Bead loading BL has been shown to be the most impactful parameter for particle breakage in the literature [29,47,48].An increase in BL increases the number concentration of the beads and the value of the radial distribution function at contact, decreases the inter-bead distance, and dramatically increases the number of bead-bead collisions, which ultimately increases the frequency of drug particle compressions and the apparent breakage rate constant [20,29].The studies described in Section 4.2.1, Section 4.2.2 (Figure S1), Section 4.2.3 (Figure S2) and Section 4.2.4 are summarized in Table 4 with the bead loading range explored in each study, the fitted exponent for bead loading with Model B (Equation ( 5)), and the relative milling rate due to bead loading (BL N 2 ), which accounts for the differences in the milling rate resulting from the different bead loading levels.Since there were only two levels of bead loading in the DP3 dataset, 85% and 99.8%, and a different exponent was observed between the two of only 1.29, the relative milling rate at 99.8% was calculated as a ratio to the rate at 85% using the following equation: 85% 2.89 99.8% 1.29  85% 1.29 , where 1.29 is the exponent obtained in the DP3 study, and 2.89 was obtained in the DP1 study which used 85% bead loading.Only parameter fits that were significant were included in this analysis.When the relative milling rate due to bead loading is plotted versus bead loading (Figure 2), a nearly perfect cubic relationship is observed over all studies carried out using 56-90% bead loading.The DP1 and DP3 studies at DV300 scale added the complexity that the bead loading was 99.8%, near the maximum packing limit in the mill, where some level of inefficiency was observed.In addition, 90% bead loading and 99.8% bead loading did not lead to markedly different observed milling rates.Observed and predicted milling rates for 56-99.8%bead loading following a more complex relationship (Equation ( 11)).
If the process stays within 56-90% bead loading, then the simple cubic term on bead loading is appropriate (i.e., N2 = 3), but in excess of 90% bead loading, a more complex term instead of BL N2 is needed, as shown in Equation (11).
where "BL" is the bead loading expressed as a percentage.Guner et al. [29] showed that the impact of bead loading becomes more significant when bead loading approaches the packing limit since beads become much closer to each other, as shown by the radial distribution function.However, bead loading beyond 95% was not explored in that study, and the mills used in this study seem to have an efficiency loss when bead loading was higher than 95.0%, as implied by Equation (11).Hence, the right-hand side of Equation (11) was substituted in place of BL 3 in Equation ( 6) for studies where BL exceeded 90%.

Impact of Bead Material
Density is an important property of the bead materials as it directly affects the energy input to the process, and therefore collision frequency and stress [29,48].A detailed microhydrodynamic analysis suggests that because these have higher density, Yttriumstabilized zirconia (YSZ) beads induce more frequent and forceful collisions than crosslinked polystyrene (CPS) beads, which favors drug particle breakage as signified by If the process stays within 56-90% bead loading, then the simple cubic term on bead loading is appropriate (i.e., N 2 = 3), but in excess of 90% bead loading, a more complex term instead of BL N2 is needed, as shown in Equation (11).
where "BL" is the bead loading expressed as a percentage.Guner et al. [29] showed that the impact of bead loading becomes more significant when bead loading approaches the packing limit since beads become much closer to each other, as shown by the radial distribution function.However, bead loading beyond 95% was not explored in that study, and the mills used in this study seem to have an efficiency loss when bead loading was higher than 95.0%, as implied by Equation (11).Hence, the right-hand side of Equation (11) was substituted in place of BL 3 in Equation ( 6) for studies where BL exceeded 90%.

Impact of Bead Material
Density is an important property of the bead materials as it directly affects the energy input to the process, and therefore collision frequency and stress [29,48].A detailed microhydrodynamic analysis suggests that because these have higher density, Yttrium-stabilized zirconia (YSZ) beads induce more frequent and forceful collisions than crosslinked polystyrene (CPS) beads, which favors drug particle breakage as signified by the higher apparent breakage rate constant when using YSZ beads [29,48].Bead density in our model was considered a ratio to the density of the YSZ beads (6000 kg/m 3 ), as it is the most preferred bead material [25].Normalized bead density was raised to the 1.4 power, , as per results from a study completed at NJIT comparing milling rates of crosslinked polystyrene (CPS) beads to YSZ beads [29,48,49].CPS beads in this study were noted to have a density of 1040 kg/m 3 .This study is discussed in more detail in Section 4.2.2.The power of 1.4 was obtained by averaging the N 3 fits to x 50 and x 90 size classes with Model B (Equation ( 5)), as the x 10 fit was insignificant due to noise.The Young's modulus YM b and the Poisson's ratio PR b of the beads are 200 GPa and 0.20 for the YSZ beads [62] and 1.5 GPa and 0.33 for the CPS beads [63].

Impact of Bead Size
Bead size can be an important process parameter as it may provide combined advantages of faster breakage [21,50,64], lower heat generation rate [65][66][67], and lower media wear [21,64].However, its impact on the breakage rate may be more subtle because it has counteracting microhydrodynamic trends where smaller beads capture drug particles more frequently but apply less stress upon contact [20,60].Therefore, the exponent for the bead size term was small, at only 0.3 in our model.The normalized bead size with the reference 0.33 mm bead size, corresponding to a nominal size of 300 µm, affects the apparent breakage rate constant k j through D b 0.33 0.3 .The impact of bead size was rather weak compared to other parameters; so, a carefully designed DoE and data with minimal noise are required to detect this relationship with a significant model fit.Among the studies in this paper, only an NJIT study in Section 4.2.3 and a DP3 study in Section 4.2.4 varied the bead sizes.However, since the DP3 data did not have enough degrees of freedom, and the NJIT study had noise in the x 10 and x 90 data, only the x 50 fit with Model B (Equation ( 5)) provided an accurate prediction for the bead size exponent.Even though a 0.3 fit was obtained by only one quantile particle size of one drug, the fixed 0.3 exponent worked well for both drugs and all quantiles.4.1.5.Impact of Temperature Temperature during milling varies over time and is influenced by the process parameters [65][66][67], which may have an impact on the particle size evolution.To investigate this hypothesis, we proposed an Arrhenius function for the j th size quantile that accounts for the temperature dependency, i.e., Kt j (see Section 3.2.2).K2 j is a fitted Arrhenius parameter that accounts for the temperature effect for the j th quantile.At the reference temperature of 15 • C, Kt j will normalize to a value of one.This prevents uncertainties in the estimate of K2 j from affecting the error estimate of A j .The Arrhenius parameter K2 j is analogous to an activation energy divided by the universal gas constant R and is likely affected by several competing temperature-dependent processes that are all exponential with respect to temperature, including changes in viscosity, Young's modulus, and possibly particle growth mechanisms like Ostwald ripening, where higher temperatures and temperature swings can dissolve the finer particles and recrystallize the material on the larger particles.In general, milling should speed up with a reduction in formulation viscosity that would tend to occur with increasing temperature; however, this is often balanced or overshadowed by the faster aggregation-Ostwald ripening rates at higher temperatures [68], or slower breakage rates associated with changes in Young's modulus [69].Temperature does vary somewhat within each run.For this reason, the time-averaged temperature at the i th timepoint (Tavg i ) was used with Model B and Model C fits (Equation ( 12)).
where t i is time at the i th time point, and T in and T out are the inlet and outlet temperatures to the mill.This equation was used to calculate the time-averaged value up to each point in time for model calculation.In general, GSK sees a trend that milling rates speed up at cooler temperatures [70]; so, negative values of K2 j are estimated during model fitting.It is possible that systems with higher viscosity may exhibit a dominant temperature-dependent viscosity effect, and milling becomes faster at increased temperatures, with resulting fits for K2 j being positive.Future studies examining more formulation effects than those considered in this paper might require breaking Kt j up into separate viscosity, Young's Modulus, and particle growth effects.If the formulation does not change, then Kt j alone appears to be a sufficient descriptor for the temperature effect.We demonstrated successful K2 j estimations for DP1 and DP5, as shown in Section 4.2.1 and Figure S15, respectively.The first version of Model B was developed at GSK in support of Drug Product 1 (DP1).After the development of the final version of Model B (Equation ( 5)), we first fitted it to a dataset coming from a Design of Experiments (DoE) exploring the impacts of tip speed and bead loading, as well as supplemental runs where batch size was changed from 4 kg to 2 kg API mass.This eight-run set of experiments was completed using the DV2000 wet bead mill where the tip speed was varied from 5 to 6 m/s, and bead loading was varied from 80% to 90%.The average temperature among runs also varied from 11.8 to 19.6 • C. Pumping rate Q was set based on V m at one mill volume per minute flowrate.The product team sampled several of the batches with a relatively large sample size, taking 18% of the batch mass in a sample point an hour before completion of milling.For this reason, the time-averaged batch mass (Average_Mass i ) was used with Model B and Model C fits for this drug product (Equation ( 13)), similar to what we did with the time-averaged temperature.

Case Studies in
We have seen in this dataset that an increase in tip speed (or rotation speed for a given mill) and bead loading both led to higher Tavg.This finding is in line with earlier heat transfer-generation studies [65][66][67], which demonstrated that a higher power consumption was observed at the higher rotation speed, and that higher bead loading was the origin of higher heat generation and ensuing temperature rise.
YTZ beads of 300 µm nominal size were used throughout the experiments; thus, exponents of the bead density and bead size terms were not explored in this case study.Since only one mill scale had been used here, the mill-scale efficiency factor (E) was not explored either.The exponents of the parameters that were not explored in this case study were set to those from Model C, which was developed based on a combination of case studies as summarized in Section 4.1.Model B, as shown in Equation (5) with N 3 = 1.4 and N 4 = 0.3 and N 5 = N 1 , was fitted to the reciprocal particle size quantiles to find the ω and BL exponents, N 1 and N 2 .Milling was not carried out for long enough to find an x j,inf fit for the x 90 population; so, this was estimated as 0.219 µm using Equation (8) for the DV2000 dataset.The resulting model parameter fits and model goodness-of-fit (R 2 values) as calculated by JMP version 17 using analytic Gauss Newton nonlinear solver for the prediction of 1/x 10 , 1/x 50 , and 1/x 90 for the DV2000 dataset are shown in Figure 3.In the JMP figures, the dots are for the experimental data, the lines are for the model fit and shades are for the error of the fit.The deviations in the parity plots (Figure 3) appear random, which indicates that the term could capture the impact of the change in API mass well.The fits for the exponent N1 for the prediction of 1/x10, 1/x50, and 1/x90 consistently showed values in the range of 1.6 to 1.9.The bead loading BL was noted to have an exponent N2 in the 2.6-3.2range.The N1 and N2 fits slightly varied among different case studies, as discussed in Section 4.1, since those exponent fits were influenced by measurement and sampling errors, and the fitted numbers in this case study have averages of approximately N1 = 2 and N2 = 3 as in Model C. Most interestingly, the negative K2j values suggest that milling proceeded significantly faster at lower temperatures, and this effect was more pronounced for the coarser sizes (x90 vs. x50), with no temperature impact for x10.This drug product is known to exhibit Ostwald ripening, which proceeds more quickly at warmer temperatures [68].This mechanism may explain why x90 is more strongly impacted by temperature than x10 and x50.The overall goodness-of-fit was 96.9%, 97.9%, and 96.3% for the 1/x10, 1/x50, and 1/x90, respectively.This drug product was scaled up from non-DoE experiments at DV300 scale, to DoE experiments at DV2000, and to multiple runs at DV4000 scale where batch size varied between 10 kg and 30 kg of API, and tip speed varied from 5.5 to 6.5 m/s.Flow rates through the mill were also changed from 2 L/min to 4 L/min at each batch size change.This resulted in 26 runs across three scales.Since the DV300 data included batches with 99.8% bead loading, the more complex version of the bead loading term shown in Equation (11) was used instead of BL 3 .In some of the DV4000 experiments, the suspension was diluted, causing a change in the drug loading.In this case, we used time-averaged drug loading ( , ) in Equation ( 14), where  is the current API loading in the system, similar to what we did with the time-averaged temperature.
, =  , *  +  * ( −  ) Model B (Equation ( 5)) was used to model this set of experiments, with N3 = 1.4,N4 = 0.3 and N5 = N1 and fitted mill efficiencies.The "match" function in JMP was used to assign The deviations in the parity plots (Figure 3) appear random, which indicates that the γ API V m m API term could capture the impact of the change in API mass well.The fits for the exponent N 1 for the prediction of 1/x 10 , 1/x 50 , and 1/x 90 consistently showed values in the range of 1.6 to 1.9.The bead loading BL was noted to have an exponent N 2 in the 2.6-3.2range.The N 1 and N 2 fits slightly varied among different case studies, as discussed in Section 4.1, since those exponent fits were influenced by measurement and sampling errors, and the fitted numbers in this case study have averages of approximately N 1 = 2 and N 2 = 3 as in Model C. Most interestingly, the negative K2 j values suggest that milling proceeded significantly faster at lower temperatures, and this effect was more pronounced for the coarser sizes (x 90 vs. x 50 ), with no temperature impact for x 10 .This drug product is known to exhibit Ostwald ripening, which proceeds more quickly at warmer temperatures [68].This mechanism may explain why x 90 is more strongly impacted by temperature than x 10 and x 50 .The overall goodness-of-fit was 96.9%, 97.9%, and 96.3% for the 1/x 10 , 1/x 50 , and 1/x 90 , respectively.This drug product was scaled up from non-DoE experiments at DV300 scale, to DoE experiments at DV2000, and to multiple runs at DV4000 scale where batch size varied between 10 kg and 30 kg of API, and tip speed varied from 5.5 to 6.5 m/s.Flow rates through the mill were also changed from 2 L/min to 4 L/min at each batch size change.This resulted in 26 runs across three scales.Since the DV300 data included batches with 99.8% bead loading, the more complex version of the bead loading term shown in Equation (11) was used instead of BL 3 .In some of the DV4000 experiments, the suspension was diluted, causing a change in the drug loading.In this case, we used time-averaged drug loading (γ API,avg i ) in Equation ( 14), where γ APIi is the current API loading in the system, similar to what we did with the time-averaged temperature.
Model B (Equation ( 5)) was used to model this set of experiments, with N 3 = 1.4,N 4 = 0.3 and N 5 = N 1 and fitted mill efficiencies.The "match" function in JMP was used to assign parameters to each mill scale, and the E2000 and E4000 terms were estimates of the mill-scale efficiency factor E at the DV2000 and DV4000 mill scales.Here, the DV300 was assigned an efficiency value of 0.83, which is the efficiency estimate of DV300 as described in Section 4.3, and the E2000 and E4000 parameters then describe the relative efficiencies for the DV2000 and DV4000 mL mills as compared to the DV300.The resulting parameter fits for the 1/x 10 , 1/x 50 , 1/x 90 models and parity plots at each scale are shown in Figure 4.
Pharmaceutics 2024, 16, 394 18 of 36 parameters to each mill scale, and the E2000 and E4000 terms were estimates of the millscale efficiency factor E at the DV2000 and DV4000 mill scales.Here, the DV300 was assigned an efficiency value of 0.83, which is the efficiency estimate of DV300 as described in Section 4.3, and the E2000 and E4000 parameters then describe the relative efficiencies for the DV2000 and DV4000 mL mills as compared to the DV300.The resulting parameter fits for the 1/x10, 1/x50, 1/x90 models and parity plots at each scale are shown in Figure 4. We found that the model fit well across every scale considered in the DP1 study with R 2 values ranging from 0.95 to 0.98.Mill-scale efficiency factors did vary by the particle size quantile considered with each mill scale.DV2000 and DV4000 mills were found to be less efficient since their tip diameters are larger than the small-scale mills.When scaling across mills, we also found that the exponent of ω (N1) may need to be fixed at a value of 2, as ω changed with scale, and some covariance can occur between the mill-scale efficiency factor estimate and N1.For this reason, the N1 was fixed at 2 and the mill-scale efficiency factors were re-estimated as shown in Figure S3 of the Supplementary Materials.This way, we obtained a better estimate of the mill-scale efficiency factors so that we could use them in Model C and Equation (11) with higher confidence.

Impact of Bead Loading, Bead Material, Rotor Speed (NJIT Study)
This dataset followed a two-level full-factorial design of experiments (DoE) with stirrer speed (3000, 4000 rpm), bead loading (56%, 79%), and bead material (YSZ, CPS) variations.The formulation of 10% FNB, 7.5% HPC-L, and 0.05% SDS was milled with Microcer 80 mL in all runs.The details of the solid-state characterization of the milled suspension, the particle size evolution in each run, and the stability of the formulation can be found in Refs.[29,48,49].YSZ beads have been the most preferred beads in the wet stirred media milling literature, followed by CPS beads.While YSZ beads provide more stressful and frequent collisions due to their higher density ( = 6000 kg/m 3 vs. We found that the model fit well across every scale considered in the DP1 study with R 2 values ranging from 0.95 to 0.98.Mill-scale efficiency factors did vary by the particle size quantile considered with each mill scale.DV2000 and DV4000 mills were found to be less efficient since their tip diameters are larger than the small-scale mills.When scaling across mills, we also found that the exponent of ω (N 1 ) may need to be fixed at a value of 2, as ω changed with scale, and some covariance can occur between the mill-scale efficiency factor estimate and N 1 .For this reason, the N 1 was fixed at 2 and the mill-scale efficiency factors were re-estimated as shown in Figure S3 of the Supplementary Materials.This way, we obtained a better estimate of the mill-scale efficiency factors so that we could use them in Model C and Equation (11) with higher confidence.

Impact of Bead Loading, Bead Material, Rotor Speed (NJIT Study)
This dataset followed a two-level full-factorial design of experiments (DoE) with stirrer speed (3000, 4000 rpm), bead loading (56%, 79%), and bead material (YSZ, CPS) variations.The formulation of 10% FNB, 7.5% HPC-L, and 0.05% SDS was milled with Microcer 80 mL in all runs.The details of the solid-state characterization of the milled suspension, the particle size evolution in each run, and the stability of the formulation can be found in Refs.[29,48,49].YSZ beads have been the most preferred beads in the wet stirred media milling literature, followed by CPS beads.While YSZ beads provide more stressful and frequent collisions due to their higher density (ρ YSZ = 6000 kg/m 3 vs.ρ CPS = 1040 kg/m 3 ), CPS beads can capture more drug particles in-between as they are softer and form a larger contact circle than YSZ beads (YM b for zirconia = 200 GPa [62] vs. YM b for polystyrene = 1.5 GPa [63]).Even though both bead materials have the potential for providing fast breakage kinetics through different mechanisms, the collision stress and frequency outweigh the contact circle diameter, and zirconia beads provide faster breakage overall [29,48].On the other hand, according to their lower collision stress, CPS beads can provide a gentler particle breakage, which, in turn, may reduce the amorphous transformation rates for sensitive materials [71] and contamination from media wear.In addition, they provided better control over temperature increase during milling [49], making CPS a competitive alternative to YSZ.In order to address the differences in milling rates via different bead materials, since zirconia is the most used material and density is the most important feature of the materials in terms of particle breakage rate, the relative bead density term ρ b /ρ YSZ was used in the model.Model B (Equation ( 5)) with N 4 = 0.3 and N 5 = 2 was fitted to the experimental reciprocal particle size data.
Figure S1 shows the fitted parameters and fitting statistics, where x 50 and x 90 fits were statistically significant, and the exponents of ω and BL (N 1 and N 2 ) for x 50 were found to be 2.13 and 3.37, respectively, which compares well to all the exponents from other studies, which led to the use of 2 and 3 in Model C. The scatter in 1/x 10 resulted in an insignificant model.Therefore, only the 1/x 50 and 1/x 90 fits were kept as the bases for N 1 , N 2 , and N 3 exploration in Figure 5.
Pharmaceutics 2024, 16,394 19 of 36 = 1040 kg/m 3 ), CPS beads can capture more drug particles in-between as they are softer and form a larger contact circle than YSZ beads (YMb for zirconia = 200 GPa [62] vs. YMb for polystyrene = 1.5 GPa [63]).Even though both bead materials have the potential for providing fast breakage kinetics through different mechanisms, the collision stress and frequency outweigh the contact circle diameter, and zirconia beads provide faster breakage overall [29,48].On the other hand, according to their lower collision stress, CPS beads can provide a gentler particle breakage, which, in turn, may reduce the amorphous transformation rates for sensitive materials [71] and contamination from media wear.In addition, they provided better control over temperature increase during milling [49], making CPS a competitive alternative to YSZ.In order to address the differences in milling rates via different bead materials, since zirconia is the most used material and density is the most important feature of the materials in terms of particle breakage rate, the relative bead density term ρb/ρYSZ was used in the model.Model B (Equation ( 5)) with N4 = 0.3 and N5 = 2 was fitted to the experimental reciprocal particle size data.
Figure S1 shows the fitted parameters and fitting statistics, where x50 and x90 fits were statistically significant, and the exponents of ω and BL (N1 and N2) for x50 were found to be 2.13 and 3.37, respectively, which compares well to all the exponents from other studies, which led to the use of 2 and 3 in Model C. The scatter in 1/x10 resulted in an insignificant model.Therefore, only the 1/x50 and 1/x90 fits were kept as the bases for N1, N2, and N3 exploration in Figure 5.For model consistency across different case studies, N1 and N2 were fixed at 2 and 3, and the model was recalibrated for x50 and x90 as shown in Figure 5, yielding N3 of 1.46 and 1.40, respectively.Therefore, an average value of N3 = 1.4 with two significant figures was used in Model C (Equation ( 6)) as the power for the bead density (ρb/ρYSZ) term.The final equation using these fixed parameters for N1, N2, and N3 was fitted to x10 data, and the model became significant despite the noise in the data, as can be seen in the left side of Figure 5.The bead density term ρb/ρYSZ = 1 when YSZ beads are used, and ρb/ρYSZ = 0.17 when CPS For model consistency across different case studies, N 1 and N 2 were fixed at 2 and 3, and the model was recalibrated for x 50 and x 90 as shown in Figure 5, yielding N 3 of 1.46 and 1.40, respectively.Therefore, an average value of N 3 = 1.4 with two significant figures was used in Model C (Equation ( 6)) as the power for the bead density (ρ b /ρ YSZ ) term.The final equation using these fixed parameters for N 1 , N 2 , and N 3 was fitted to x 10 data, and the model became significant despite the noise in the data, as can be seen in the left side of Figure 5.The bead density term ρ b /ρ YSZ = 1 when YSZ beads are used, and ρ b /ρ YSZ = 0.17 when CPS beads are used.Hence, with N 3 = 1.4, the model would predict a slower decrease in particle size during milling for polystyrene beads compared to zirconia beads under otherwise identical conditions, which is in agreement with other findings in the literature [72].This experimental study was designed to have three different rotor speeds ω (3000, 3500, and 4000 rpm), three different bead loadings BL (56%, 68%, and 79%) and two different YSZ bead sizes D b (0.2 and 0.4 mm) with a total of 10 runs.Milling was carried out for 3 h using the same mill (Microcer 80 mL located at NJIT) and the same drug formulation (10% GF, 7.5% HPC-L, and 0.05% SDS).Model B (Equation ( 5)) with N 3 = 1.4 and N 5 = 2 was used to fit the reciprocal particle sizes, allowing for the estimation of N 1 , N 2 and N 4 .The impact of process parameters was assessed by fitting 1/x j data; however, since the scatter in 1/x 10 and 1/x 90 data was usually higher compared to 1/x 50 data, the significance of the model fits was not as strong (see Figure S2 and Figure 6), similar to 1/x 10 in Section 4.2.2.This could be attributed to the low level of aggregation affecting the tails of the PSD at the very high HPC-L loading.This high concentration of HPC-L was selected to increase power consumption P so that accurate P values can be used in the microhydrodynamic model in a previous study [49].The 1/x 10 model fails to have a significant rate term (A), whereas the 1/x 90 model fails to have a significant bead size exponent term (N 4 ), as can be seen in Figure S2 in the Supplementary Materials.Therefore, we focused on the 1/x 50 model for the investigation of w, BL, and D b impacts.The positive exponents indicate that the milling rate increased at higher ω and BL with the use of smaller beads.The exponents of ω and BL were observed to be 1.99 and 2.84, respectively, which were again very close to the exponents obtained in other studies; so, this informed the use of exponents of 2 and 3 in Model C. As can be seen in Figure 6, when the N 1 and N 2 values were set to 2 and 3 for x 50 , N 4 was found to be 0.3, which is the justification of the D bead 0.33 mm 0.3 term in the Model C. The relative impact of process parameters, i.e., BL > ω > D B , on the milling rate agreed with previous studies [20,60].
Although the parameter ranges used in this dataset and those in Section 4.2.2 were the same, the drugs were different; regardless, the impact of stirrer speed ω and bead loading BL on the milling rate was found to be consistent.

Impact of Rotor Speed, Bead Size, and Scale from DV150-DV300 (DP3 Study)
In this dataset, there were six experiments in total, five of which were performed using the DV300, and one was performed using the DV150.While the DV150 experiment was carried out with 85% bead loading, 0.3 mm nominal sized YSZ beads, and 5.5 m/s tip speed, the DV300 experiments were performed with two different bead loadings, which were higher compared to previous case studies (85% and 99.8%), and with two different bead sizes (0.3, 0.65 mm), where all smaller bead studies were loaded at 99.8% and larger beads were loaded at 85%.Since the experimental design was not well suited for capturing the impact of both bead size and bead loading at the same time, we used the previously developed D bead 0.33mm 0.3 term and fitted the bead loading exponent to understand the impact of increasing bead loading in the high bead loading range close to packing.A version of Equation ( 5), where N 1 = 2, N 3 = 1.4,N 4 = 0.3 and N 5 = N 1 , with the predicted x j,inf values shown in Table 5, was fitted to 1/x 10 , 1/x 50 , and 1/x 90 data, and the fit statistics together with the parity plots are shown in Figure 7.The fitted N 2 values were different from previously observed values due to a possible efficiency loss at high bead loading close to the packing limit.The average of the fitted N 2 values was found to be 1.29, which was used in the BL N2 analysis described in Section 4.2.1.
higher ω and BL with the use of smaller beads.The exponents of ω and BL were observed to be 1.99 and 2.84, respectively, which were again very close to the exponents obtained in other studies; so, this informed the use of exponents of 2 and 3 in Model C. As can be seen in Figure 6, when the N1 and N2 values were set to 2 and 3 for x50, N4 was found to be 0.3, which is the justification of the . .term in the Model C. The relative impact of process parameters, i.e., BL > ω > DB, on the milling rate agreed with previous studies [20,60].Although the parameter ranges used in this dataset and those in Section 4.2.2 were the same, the drugs were different; regardless, the impact of stirrer speed ω and bead loading BL on the milling rate was found to be consistent.In this dataset, there were six experiments in total, five of which were performed using the DV300, and one was performed using the DV150.While the DV150 experiment was carried out with 85% bead loading, 0.3 mm nominal sized YSZ beads, and 5.5 m/s tip speed, the DV300 experiments were performed with two different bead loadings, which were higher compared to previous case studies (85% and 99.8%), and with two different bead sizes (0.3, 0.65 mm), where all smaller bead studies were loaded at 99.8% and larger beads were loaded at 85%.Since the experimental design was not well suited for capturing the impact of both bead size and bead loading at the same time, we used the previously developed .
. term and fitted the bead loading exponent to understand the impact of increasing bead loading in the high bead loading range close to packing.A version of Equation ( 5), where N1 = 2, N3 = 1.4,N4 = 0.3 and N5 = N1, with the predicted xj,inf values shown in Table 5, was fitted to 1/x10, 1/x50, and 1/x90 data, and the fit statistics together with the parity plots are shown in Figure 7.The fitted N2 values were different from previously observed values due to a possible efficiency loss at high bead loading close to the packing limit.The average of the fitted N2 values was found to be 1.29, which was used in the BL N2 analysis described in Section 4.2.1.

Benefits of Model C: Underlying Trends in the Factors Governing Milling Other Than Process Parameters
While fitting Model B would result in better fitting statistics since it has more degrees of freedom, a consistent model fit to a group of studies, where different APIs, formulations, and mills were used, would better elucidate whether there is any

Benefits of Model C: Underlying Trends in the Factors Governing Milling Other Than Process Parameters
While fitting Model B would result in better fitting statistics since it has more degrees of freedom, a consistent model fit to a group of studies, where different APIs, formulations, and mills were used, would better elucidate whether there is any underlying behavior governing the impact of drug, suspension, and mill properties, and process conditions on apparent grinding limit, x j,inf , mill-scale efficiency factor, E, and milling rate, A j .In this section, Model C (Equation ( 6)) was fitted to multiple datasets so that the fitted parameters on x j,inf , E, and A j .could be compared across studies.Then, driving forces for the difference in x j,inf , E, and A j fits were explored considering drug, process, and mill properties.

The Apparent Grinding Limit
If the particle size does not approach the apparent grinding limit, x j,inf , within the given milling time, the nonlinear fitting routine will yield a lack of parameter fit for x j,inf and potentially lack of significance in the N j and A j fits, as these parameters are all somewhat covariant.Equation ( 1) can always be used as the linearizing transform in these cases, but it would be convenient to have a model to predict x j,inf in these cases, thus allowing for the estimation of the entire sigmoidal curve.This then engables the use of the transform in Equation (2).To estimate x j,inf , we sought a mathematical relation for the fitted x j,inf values as a function of relevant process parameters and Young's modulus of the APIs (YM API ) and arrived at Equation ( 8), which was previously discussed in Section 3.2.
To derive Equation ( 8), we first needed to obtain x j,inf by fitting Model C (Equation ( 6)) to datasets for four drug products where milling was run for long enough for the particle size to approach the apparent grinding limit.We used data from multiple drug products to derive our empirical x j,inf model to ensure that it can generalize across drug products.Since the particle size approached x j,inf in all runs in the griseofulvin case study from Section 4.2.3,Model C (Equation ( 6)) was fitted to experimental 1/x 50 values for each individual run, yielding fitted x 50,inf values.For griseofulvin, since the x 10 and x 90 measurements were noisy compared to x 50 data, instead of having individual run fits, x 10,inf and x 90,inf were found only by fitting Model C to the complete dataset.The other three drug products where we could obtain a significant x j,inf fit to the data were DP1 in Section 4.2.1, fenofibrate in Section 4.2.2 and DP3 in Section 4.2.4.The datasets are summarized in Table 5 with tip speeds U tip , bead loadings BL, bead sizes D b , YM API , and fitted x j,inf values.Following the Model C fits, we sought a power-law correlation between the fitted x 50,inf values and U tip , BL, D b , and YM API , but the U tip contribution was found to be insignificant.The resulting fitted power-law correlation for x 50,inf is shown in Equation ( 15), with the associated fitting statistics and parity plot in Figure 8.  15)) with associated parity plot and statistics for its parameters.
The x10,inf and x90,inf values were highly correlated with x50,inf.Therefore, the models for x10,inf and x90,inf were built by only fitting Cj, the coefficient of Equation ( 15), and fixing the exponents b1, b2, and b3 to −0.25, 0.12 and 0.24, respectively.C10 was found to be 0.057 with R 2 of 0.85, and C90 was found to be 0.125 with R 2 of 1.0, since it was obtained with only two data points.The resulting equation is captured as Equation ( 8) in Section 3.2.
Using Equation ( 8), we were able to predict the apparent grinding limit (xj,inf) for case studies where milling time was relatively short, i.e., DP3, DP5, DP6 and DP4.When fitting models, we estimated model parameters A, B, N and xinf using the complete dataset instead of obtaining model parameters for each experiment.For the sake of consistency, when estimating xinf with Equation (8), we took the average of BL and Db values if multiple bead loadings and sizes were used in the dataset, in order to obtain a single xinf estimation that is representative of the complete dataset.These grinding limit predictions thus enabled us to use Model C (Equation ( 6)) with the transform in Equation (2) to perform further analysis.
A deeper examination of apparent grinding limits and Equation (8) shows that xj,inf values for GF varied in a relatively narrow range as the process parameters were changed (see Table 5), and the exponents in Equation ( 8) were rather small, signifying no impact of Utip (or w) and a weak impact of BL and Db.This is unsurprising because the grinding limit is mainly determined by material properties, not by the process parameters, provided the milling time is sufficiently long.In fact, the true grinding limit of APIs is independent of process parameters, depending solely on material properties such as YMAPI, hardness H, and fracture toughness Kc, which collectively determine the brittleness of the APIs as scored by the brittleness index BI = H/Kc, and the brittle-ductile transition size [57,58].We speculate that the relatively low exponent of YMAPI could be related to the fact that YMAPI was used here as the sole descriptor of material properties due to a lack of data on other API material properties such as H and Kc.Thus, to refine Equation ( 8), we would need to study a broad range of material properties, which is outside the scope of this paper.

Mill-Scale efficiency Factor (E)
We found that most of the scaling between mills is reasonably well predicted by only considering the batch size and mill volume available for beads.Just these two terms could be used without any further consideration for mill-scale efficiency differences with average errors of approximately 20% (but could be as high as 45%) when predicting the required milling time to reach a specific particle size endpoint.In an effort to reduce this error when scaling between mills, we introduced a mill-scale efficiency factor, E, to score the efficiency of particular mills.The mill-scale efficiency factor is derived from parameter fitting and scores mills relative to a reference mill, here DV150.Thus, the DV150 mill-scale efficiency factor is fixed to 100%, whereas the efficiency factor for other mills is determined  15)) with associated parity plot and statistics for its parameters.The x 10,inf and x 90,inf values were highly correlated with x 50,inf .Therefore, the models for x 10,inf and x 90,inf were built by only fitting C j , the coefficient of Equation ( 15), and fixing the exponents b 1 , b 2 , and b 3 to −0.25, 0.12 and 0.24, respectively.C 10 was found to be 0.057 with R 2 of 0.85, and C 90 was found to be 0.125 with R 2 of 1.0, since it was obtained with only two data points.The resulting equation is captured as Equation ( 8) in Section 3.2.
Using Equation ( 8), we were able to predict the apparent grinding limit (x j,inf ) for case studies where milling time was relatively short, i.e., DP3, DP5, DP6 and DP4.When fitting models, we estimated model parameters A, B, N and x inf using the complete dataset instead of obtaining model parameters for each experiment.For the sake of consistency, when estimating x inf with Equation (8), we took the average of BL and D b values if multiple bead loadings and sizes were used in the dataset, in order to obtain a single x inf estimation that is representative of the complete dataset.These grinding limit predictions thus enabled us to use Model C (Equation ( 6)) with the transform in Equation (2) to perform further analysis.
A deeper examination of apparent grinding limits and Equation (8) shows that x j,inf values for GF varied in a relatively narrow range as the process parameters were changed (see Table 5), and the exponents in Equation ( 8) were rather small, signifying no impact of U tip (or w) and a weak impact of BL and D b .This is unsurprising because the grinding limit is mainly determined by material properties, not by the process parameters, provided the milling time is sufficiently long.In fact, the true grinding limit of APIs is independent of process parameters, depending solely on material properties such as YM API , hardness H, and fracture toughness K c , which collectively determine the brittleness of the APIs as scored by the brittleness index BI = H/K c , and the brittle-ductile transition size [57,58].We speculate that the relatively low exponent of YM API could be related to the fact that YM API was used here as the sole descriptor of material properties due to a lack of data on other API material properties such as H and K c .Thus, to refine Equation ( 8), we would need to study a broad range of material properties, which is outside the scope of this paper.

Mill-Scale efficiency Factor (E)
We found that most of the scaling between mills is reasonably well predicted by only considering the batch size and mill volume available for beads.Just these two terms could be used without any further consideration for mill-scale efficiency differences with average errors of approximately 20% (but could be as high as 45%) when predicting the required milling time to reach a specific particle size endpoint.In an effort to reduce this error when scaling between mills, we introduced a mill-scale efficiency factor, E, to score the efficiency of particular mills.The mill-scale efficiency factor is derived from parameter fitting and scores mills relative to a reference mill, here DV150.Thus, the DV150 mill-scale efficiency factor is fixed to 100%, whereas the efficiency factor for other mills is determined as a parameter fit that minimizes the sum of the squared errors between particle sizes from model and data.We can determine mill efficiency scores for each particle size quantile and drug product.We found that these efficiency scores are reasonably consistent across drug products and particle size quantiles.Mill efficiency scores may also be affected by wear and tear as the components in the mill become worn down over long time periods.We can average mill efficiency for each mill over all drug products and particle size quantile, which yields the overall mill-scale efficiency factor, E, for each mill.We used the E terms that were obtained by averaging the individual efficiencies fitted to each drug product and size class.Readers can follow the same procedure for their own mills, or they may use Equation (7) to have a rough estimate of the efficiency of their mills.The mill-scale efficiency scores for different mills and drug products are shown in Table 6.The parity plots for the efficiency fits with Model C are shown for DP1, DP2, DP3, DP4, DP5 and DP6 in Figure S3, Figure S4, Figure S5, Figure S6, Figure S7 and Figure S8, respectively.We explored how mill properties might affect the differences in the mill efficiencies and developed an empirical correlation (see Equation (7) in Section 3.2) describing mill-scale efficiency factor as a function of tip diameter and the length to diameter ratio of the agitator (blade) length to mill chamber diameter.The model fit and parity plots for Equation (7) are shown in Figure 9.If readers are interested in using the Model C with a mill with different dimensions, Equation (7) (Section 3.2.3)can be used to estimate the mill efficiency.Though the parameters have highly significant estimates in this model, we caution readers that this may be only applicable to Netzsch mills, which have many of the mill dimensions scaled to correlate with the agitator diameter.Since there are many mill dimensions that scale similarly to agitator diameter and length to diameter ratio in our study, it is not unexpected to find those parameters and their interaction to be significant predictors of mill efficiencies that could be explained with a multiple linear regression model such as the following.

Milling Rate (Aj)
In this section, we used Model C, with the mill-scale efficiency factors reported in Table 6, to fit Aj, Bj, Nj, K2j, and xj,inf for each case study.To generalize when to fit the xj,inf term and when to predict it via Equation ( 8), the decision criteria were selected as fit if endpoint x10, x50, x90 < 0.1, 0.15, 0.3 µm, respectively, and predict via Equation ( 8) otherwise.The model fits and parity plots can be seen for FNB, GF, DP1, DP2, DP3, DP4, DP5 and DP6 in Figure S9, Figure S10, Figure S11, Figure S12, Figure S13, Figure S14, Figure S15 and Figure S16, respectively.Table 7 summarizes all the fitted model parameters, where the unit of Aj is . . ., Bj and Nj are dimensionless, xinf,j is in µm, and K2j is in °K.In this section, we used Model C, with the mill-scale efficiency factors reported in Table 6, to fit A j , B j , N j , K2 j , and x j,inf for each case study.To generalize when to fit the x j,inf term and when to predict it via Equation ( 8), the decision criteria were selected as fit if endpoint x 10 , x 50 , x 90 < 0.1, 0.15, 0.3 µm, respectively, and predict via Equation ( 8) otherwise.The model fits and parity plots can be seen for FNB, GF, DP1, DP2, DP3, DP4, DP5 and DP6 in Figure S9, Figure S10, Figure S11, Figure S12, Figure S13, Figure S14, Figure S15 and Figure S16, respectively.Table 7 summarizes all the fitted model parameters, where the unit of A j is kg 0.4 m 2.4 s 1.8 , B j and N j are dimensionless, x inf,j is in µm, and K2 j is in • K. Similar to x j,inf , the parameters A j , B j , and N j can be mathematically correlated with drug material properties.We note that these parameters may also be impacted by formulation properties, but since the formulation was kept constant in each case study, and the measurements for suspension properties such as viscosity and density were not available in most case studies, we did not investigate their relationship with the model parameters in the scope of this paper.For such correlations to be accurate, the suspension must be stable in terms of particle size growth and crystallinity.In addition, having fixed exponents in Model C helps in identifying what factors impact the model parameters other than process parameters in a dataset that consists of multiple drug products.As a first attempt to correlate milling characteristics of drugs with drug properties, in Section 4.3.1,we sought correlations for A j , which had a significant fit in all studies.The correlations of the milling rate term A j with YM API are shown in Figure 10.This correlation for A j was stronger for the coarser quantiles, as shown by the higher R 2 values.Additional predictors such as the viscosity of the suspension and material properties like H and K c could strengthen the correlation.The relationships observed in Figure 10 suggest that wet bead milling processes of materials with a YM API of less than approximately 4 GPa may be impractical for drug product development due to the predicted low milling rates of the x 50 and x 90 .The milling rate A j parameter trends toward zero below a YM API of 3 GPa.This tends to manifest as a bi-modal population where a large size class influencing the x 90 is time consuming to eliminate.Generally, low YM materials will require lower milling temperatures, higher bead loading, and larger tip speeds to achieve desired size specifications within reasonable manufacturing times at large scales.While YM correlates with milling rate, H and K c may improve this correlation and similar correlations can be sought for B and N parameters in the future.Similar to xj,inf, the parameters Aj, Bj, and Nj can be mathematically correlated with drug material properties.We note that these parameters may also be impacted by formulation properties, but since the formulation was kept constant in each case study, and the measurements for suspension properties such as viscosity and density were not available in most case studies, we did not investigate their relationship with the model parameters in the scope of this paper.For such correlations to be accurate, the suspension must be stable in terms of particle size growth and crystallinity.In addition, having fixed exponents in Model C helps in identifying what factors impact the model parameters other than process parameters in a dataset that consists of multiple drug products.As a first attempt to correlate milling characteristics of drugs with drug properties, in Section 4.3.1,we sought correlations for Aj, which had a significant fit in all studies.The correlations of the milling rate term Aj with YMAPI are shown in Figure 10.This correlation for Aj was stronger for the coarser quantiles, as shown by the higher R 2 values.Additional predictors such as the viscosity of the suspension and other material properties like H and Kc could strengthen the correlation.The relationships observed in Figure 10 suggest that wet bead milling processes of materials with a YMAPI of less than approximately 4 GPa may be impractical for drug product development due to the predicted low milling rates of the x50 and x90.The milling rate Aj parameter trends toward zero below a YMAPI of 3 GPa.This tends to manifest as a bi-modal population where a large size class influencing the x90 is time consuming to eliminate.Generally, low YM materials will require lower milling temperatures, higher bead loading, and larger tip speeds to achieve desired size specifications within reasonable manufacturing times at large scales.While YM correlates with milling rate, H and Kc may improve this correlation and similar correlations can be sought for B and N parameters in the future.

Typical Model Use and Validation
Thus far, model development has been described without the mention of the order in which the data were gathered or how the model is typically used during product development.Once Model C was built on the first three drug products (DP1-3), the model started to be used for a priori prediction of scale-up from small scale batches.Typically, Model C uses DV50 and/or DV150 data on small volume batches to estimate the A, B and N terms, thus allowing for the prediction of processing times for large batches produced on the DV300.At times, the model has been used to prospectively predict both results of scale changes, as well as tip speed changes to reduce processing times such that a large batch can be milled in a single shift.The model has been used prospectively to predict scale-up for DP4-6.This provides an ongoing validation that the model can be used to predict milling operations in advance after being trained on very limited data, often only one or two batches at a small scale.
An illustrative example for Drug Product 6 (DP6) is shown in Figure 11.These plots show experimental and predicted particle sizes versus milling time for a run at a DV300 scale.The predictions in Figure 11A are for the model in Figure S16, which was trained by

Typical Model Use and Validation
Thus far, model development has been described without the mention of the order in which the data were gathered or how the model is typically used during product development.Once Model C was built on the first three drug products (DP1-3), the model started to be used for a priori prediction of scale-up from small scale batches.Typically, Model C uses DV50 and/or DV150 data on small volume batches to estimate the A, B and N terms, thus allowing for the prediction of processing times for large batches produced on the DV300.At times, the model has been used to prospectively predict both results of scale changes, as well as tip speed changes to reduce processing times such that a large batch can be milled in a single shift.The model has been used prospectively to predict scale-up for DP4-6.This provides an ongoing validation that the model can be used to predict milling operations in advance after being trained on very limited data, often only one or two batches at a small scale.
An illustrative example for Drug Product 6 (DP6) is shown in Figure 11.These plots show experimental and predicted particle sizes versus milling time for a run at a DV300 scale.The predictions in Figure 11A

Conclusions
In this work, we developed a semi-mechanistic modeling framework with elements from microhydrodynamic theory to predict particle size evolution in pharmaceutical wet bead milling over a wide range of process conditions and scales.Models A, B, and C, with different levels of complexity and practical utility, offer pharmaceutical engineers a wide range of capabilities to simulate, scale-up, and understand their processes.Model A, derived directly from microhydrodynamic theory, is the most general model and entails significant effort to parameterize.The full potential of Model A should be tested in a future study.Being simple and computationally convenient, Model B is expected to reduce the amount of experimentation needed to develop, optimize, scale-up, scale-down, and scaleout the wet bead milling process.We demonstrated the wide applicability of Model B in multiple case studies across drug products and mill scales.Model C, though developed by fitting to data for the specific drug products and mills described here, requires even fewer experiments to parametrize and deploy to future milling batches.We demonstrated the power of such a drug-product-and mill-agnostic model to provide insights into the most important mechanisms that govern milling rates outside of process parameters.
In the future, Model B can be further refined by considering formulations in a wider range of viscosity and flow regimes.The impact of residence time distribution may be incorporated through circuit Pe correlations as a function of tip speed, axial mean velocity

Conclusions
In this work, we developed a semi-mechanistic modeling framework with elements from microhydrodynamic theory to predict particle size evolution in pharmaceutical wet bead milling over a wide range of process conditions and scales.Models A, B, and C, with different levels of complexity and practical utility, offer pharmaceutical engineers a wide range of capabilities to simulate, scale-up, and understand their processes.Model A, derived directly from microhydrodynamic theory, is the most general model and entails significant effort to parameterize.The full potential of Model A should be tested in a future study.Being simple and computationally convenient, Model B is expected to reduce the amount of experimentation needed to develop, optimize, scale-up, scale-down, and scale-out the wet bead milling process.We demonstrated the wide applicability of Model B in multiple case studies across drug products and mill scales.Model C, though developed by fitting to data for the specific drug products and mills described here, requires even fewer experiments to parametrize and deploy to future milling batches.We demonstrated the power of such a drug-product-and mill-agnostic model to provide insights into the most important mechanisms that govern milling rates outside of process parameters.
In the future, Model B can be further refined by considering formulations in a wider range of viscosity and flow regimes.The impact of residence time distribution may be incorporated through circuit Pe correlations as a function of tip speed, axial mean velocity of the suspensions, and number of turnovers.Future research should also consider the dependence of apparent and true grinding limit on the hardness and fracture toughness of APIs in addition to Young's modulus.Overall, this comprehensive theoretical and experimental study has provided a semi-mechanistic modeling framework that can be adopted and adapted in a fit-for-purpose manner by pharmaceutical engineers to simulate, optimize, and scale wet bead milling processes.In Equation (A10), ε m is the nondimensional bead-bead gap thickness at which the lubrication force stops increasing and becomes a constant.Parameters k 1 and k 2 were computed using the following: All parameters and variables in Equation (A6) are either known or experimentally measured except for the granular temperature, which can be easily solved by any nonlinear equation solver.Next, granular temperature can be used to calculate the average frequency of drug particle compressions between the beads a, which is the multiplication of the probability of single drug particle to be caught by the beads p and the average oscillation frequency of a single bead ν.To derive a rough scaling approximation for ω and D a dependence of θ when only ω or D a is varied to affect P v , keeping all process-formulation-design parameters fixed, we neglected ε ht , defined z 2 = θ, and rewrote Equation (A6) as a cubic equation of the form, P v − λ 1 z 2 − (λ 2 +λ 3 )z 3 = 0, where λ 1 , λ 2 , and λ 3 correspond to the constant positive-valued factors obtained from plugging Equation (A8) into Equation (A6).Instead of finding the only real, positive-valued root analytically, for which a numerical solution is much more convenient as mentioned above and as implemented in previous microhydrodynamic studies [29,60], we performed an asymptotic analysis: First, λ 1 << λ 2 + λ 3 was assumed, which corresponds to turbulent flow with negligible viscosity effects (formulations with low viscosity) and yields an approximate solution: θ ∼ = P v 2/3 /(λ 2 + λ 3 ) 2/3 .The other asymptotic solution pertains to viscous formulations for which the viscous effects tend to dominate (laminar flow), i.e., λ 1 >> λ 2 + λ 3 .In this case, we found θ ∼ = P v /λ 1 .The exponent of P v is theoretically bound between 2/3 and 1.Note that these two asymptotic solutions were not intended to replace either a numerical solution or a full analytical solution; rather, in the absence of any other information, it allowed us to estimate the w and D a dependence of θ under the limiting condition of turbulent flow (see Section 3.2).

Appendix C. Derivation of Model B from Model A
In Section 3.2.2,we described the development of a semi-mechanistic model for wet bead milling.The connectivity between Equations (3) and ( 5) is explained here.The particle size term in Equation (A13) has been wrapped in the slope term in Equation ( 4).As the time dependence of x p was already modeled by the transformation Yt in Equations ( 1) and ( 2), for the sake of convenience and preservation of dimensional homogeneity, we set x p = x j (0) = x j,ini .A similar approach was also adopted to assess the impacts of various process parameters and bead properties in earlier microhydrodynamic studies (e.g., [29,48]).Plugging Equation (A13) into Equation (3) and lumping all constants into a new parameter A * * j = 0.02328x j,ini A * j / √ π, we obtained Equation (4) of the main text (Model A).To formulate a microhydrodynamically inspired semi-mechanistic model capable of scale-up, we assumed that θ could be described empirically as a power-law model as follows: where ref stands for reference mill and bead conditions and D a is the agitator's tip (blade) diameter.Here, E is a mill scale factor that accounts for the energy transfer efficiency of any mill with respect to a reference mill, which captures the impact of mill design differences at

Figure 8 .
Figure 8.An empirical model for x 50,inf (Equation (15)) with associated parity plot and statistics for its parameters.

Figure 9 .
Figure 9. Model fit and parity plot for Equation (7): natural logarithm of the mill efficiencies as a function of logarithm of agitator diameter, agitator length to diameter ratio and their interaction.* Indicates parameter significance.

Figure 9 .
Figure 9. Model fit and parity plot for Equation (7): natural logarithm of the mill efficiencies as a function of logarithm of agitator diameter, agitator length to diameter ratio and their interaction.* Indicates parameter significance.4.3.3.Milling Rate (A j )

Figure 10 .
Figure 10.Milling rate A and Young's modulus of the drugs in various case studies for x10 (left), x50 (middle), and x90 (right).

Figure 10 .
Figure 10.Milling rate A and Young's modulus of the drugs in various case studies for x 10 (left), x 50 (middle), and x 90 (right).
are for the model in Figure S16, which was trained by all available experiments.The model was trained on DV50 and DV150 data to predict the DV300 run.The predictions in Figure 11B are for the model trained via only DV50 data.The predictions in Figure 11C are for the model trained via only DV150 data.Finally, the predictions in Figure 11D are for the model trained via DV50 and DV150 data.For these cases, the model prediction fidelity was subsequently verified based on experimental data at the DV300 scale.In all cases, the predictions are very well aligned with the experimental observations, confirming the prediction fidelity of the model.Pharmaceutics 2024, 16, 394 27 of 36 all available experiments.The model was trained on DV50 and DV150 data to predict the DV300 run.The predictions in Figure 11B are for the model trained via only DV50 data.The predictions in Figure 11C are for the model trained via only DV150 data.Finally, the predictions in Figure 11D are for the model trained via DV50 and DV150 data.For these cases, the model prediction fidelity was subsequently verified based on experimental data at the DV300 scale.In all cases, the predictions are very well aligned with the experimental observations, confirming the prediction fidelity of the model.

Figure 11 .
Figure 11.Use of DV50 and/or DV150 Data to prospectively predict DV300 outcome for DP6: (A) model trained by all data, (B) model trained by DV50 data only, (C) model trained by DV150 data only), and (D) model trained by DV50 and DV150 data.

Figure 11 .
Figure 11.Use of DV50 and/or DV150 Data to prospectively predict DV300 outcome for DP6: (A) model trained by all data, (B) model trained by DV50 data only, (C) model trained by DV150 data only), and (D) model trained by DV50 and DV150 data.

Table 1 .
Overview of wet bead mills utilized and associated experiments.

Table 2 .
Comparison of the modeling approaches described in this work.

Micro-Hydrodynamic Model, Model A (Section 3.2.1) Semi-Mechanistic Model, Model B (Section 3.2.2) Pre-Calibrated Semi-Mechanistic Model, Model C (Section 3.2.3)
j , N j , x j,inf , (+more parameters to develop a model for Power)

Table 3 .
Fitted exponents of normalized ω (N 1 ) in various case studies.

Table 4 .
Fitted exponent of BL (N 2 ) averaged across x 10 , x 50 and x 90 quantiles in various case studies.

s) Bead Loading (%) Bead Size (mm) Young's Modulus (GPa) x 10,inf (µm) x 50,inf (µm) x 90,inf (µm)
a Model C was fitted by individual run, b Model C was fitted to the complete dataset at once, and the average of the process parameters was determined.

Table 6 .
Summary of the fitted efficiencies.

Table 7 .
Model C fitted parameters for all drug products.

Table 7 .
Model C fitted parameters for all drug products.

a
The average frequency of drug particle compressions, 1/min API Active pharmaceutical ingredient B jIntercept term of the Yt function for each particle size quantile j (10, 50, 90) Shape factor of the Yt transform for each particle size quantile j (10, 50, 90) ht Power spent on shear of milled suspension of the slurry at the same shear rate but calculated (measured) when no beads were present in the flow, W/m 3 ε m Nondimensional bead-bead gap thickness at which the lubrication force stops increasing and becomes a constant,ε visc Energy dissipation rate due to both the liquid-beads viscous friction and lubrication, W/m 3 ε ε powder compact out-of-die porosity γ Mass concentration, g/mL 2 + (135/64)clog c + 11.26c 1 − 5.1c + 16.57c 2 − 21.77c 3 (A11)