Predicting the Temperature Evolution during Nanomilling of Drug Suspensions via a Semi-Theoretical Lumped-Parameter Model

Although temperature can significantly affect the stability and degradation of drug nanosuspensions, temperature evolution during the production of drug nanoparticles via wet stirred media milling, also known as nanomilling, has not been studied extensively. This study aims to establish both descriptive and predictive capabilities of a semi-theoretical lumped parameter model (LPM) for temperature evolution. In the experiments, the mill was operated at various stirrer speeds, bead loadings, and bead sizes, while the temperature evolution at the mill outlet was recorded. The LPM was formulated and fitted to the experimental temperature profiles in the training runs, and its parameters, i.e., the apparent heat generation rate Qgen and the apparent overall heat transfer coefficient times surface area UA, were estimated. For the test runs, these parameters were predicted as a function of the process parameters via a power law (PL) model and machine learning (ML) model. The LPM augmented with the PL and ML models was used to predict the temperature evolution in the test runs. The LPM predictions were also compared with those of an enthalpy balance model (EBM) developed recently. The LPM had a fitting capability with a root-mean-squared error (RMSE) lower than 0.9 °C, and a prediction capability, when augmented with the PL and ML models, with an RMSE lower than 4.1 and 2.1 °C, respectively. Overall, the LPM augmented with the PL model had both good descriptive and predictive capability, whereas the one with the ML model had a comparable predictive capability. Despite being simple, with two parameters and obviating the need for sophisticated numerical techniques for its solution, the semi-theoretical LPM generally predicts the temperature evolution similarly or slightly better than the EBM. Hence, this study has provided a validated, simple model for pharmaceutical engineers to simulate the temperature evolution during the nanomilling process, which will help to set proper process controls for thermally labile drugs.


Introduction
Low water solubility is a characteristic of most of the newly discovered active pharmaceutical ingredients (APIs) [1][2][3]. The APIs with low solubility have been regarded as highly risky development candidates [4]; their development into a marketed product has been more challenging, regardless of how desired their pharmacological properties are [5]. Nevertheless, either due to high-throughput screening technologies [6,7] or the need for high lipophilicity and molecular weight for the treatment of bio-complex diseases [8], the low solubility issue remains a challenge [9]. Therefore, there is increasing motivation to study solubility and dissolution enhancement methods by pharmaceutical companies and scientists [10]. Some of the popular methods entail formulation of amorphous solid for the 27 process runs enabled us to train the power law (PL) and the machine learning (ML) model. By using the trained models, we predicted Q gen and UA of the five test runs. Minitab was used for the PL predictions, whereas Google Colab was used for the ML predictions. By inserting the predicted Q gen and UA values in the LPM, the temperature profiles of the test runs were simulated. The advantages and disadvantages of the LPM and the EBM as well as the limitations of the LPM were discussed. Not only will this study reveal the fitting capability of the LPM as compared with the EBM, but it will also enable us to assess their comparative predictive capabilities and usefulness for process development and understanding. Overall, this study provides pharmaceutical engineers with a validated, simple model (LPM), which simulates the temperature evolution during the production of drug nanosuspensions and predicts the impact of process parameters, thereby eventually helping engineers to control and optimize the process.

Materials
A BCS Class II API, fenofibrate (FNB) was used as a model poorly water-soluble drug, which was purchased from Jai Radhe Sales (BP grade, Ahmedabad, India). FNB is a lipophilic compound with a molecular weight of 360 g/mol and an aqueous solubility of 0.8 mg/L at room temperature [55]. Hydroxypropyl cellulose (HPC) was used for stabilizing the drug suspension as a non-ionic polymer, which was generously donated by Nisso America Inc. (L grade, New York, NY, USA). Its L grade has a molecular weight of 140,000 g/mol [56] and is readily water-soluble [57]. HPC is known to adsorb onto FNB particles, thereby imparting steric stability [58]. Moreover, an anionic surfactant, sodium dodecyl sulfate (SDS), was used for wettability enhancement, and it was purchased from GFS chemicals (ACS grade, Columbus, OH, USA). SDS has a molecular weight of 288 g/mol and a critical micelle concentration of 8 mM [55].

Wet Stirred Media Milling
The suspension formulation was the same for all runs where the w/v% of the ingredients was 10% FNB, 8% HPC-L, and 0.05% SDS with respect to 200 mL of deionized water. Based on our prior studies, this formulation is known to be physically stable during milling and storage [30,50]. A pre-suspension was prepared by adding the powders to deionized water gradually under constant mixing with a shear stirrer (Cat#. 14-503, Fisher Scientific, Pittsburgh, PA, USA) operating at 300 rpm for 2 h. The theoretical batch size was fixed for all processing runs: 236 g. Pre-suspensions were stored at 8 • C overnight prior to milling to let it settle and get rid of the foaming that occurred during shear mixing. On the milling days, suspensions were stirred on a magnetic stirrer until they equilibrated close to room temperature.
A Microcer wet-stirred media mill (Netzsch Fine Particle Size Technology, LLC, Exton, PA, USA) was used for milling the pre-suspensions. It has an 80 mL chamber volume lined with zirconia, a zirconia shaft, and stainless-steel screens whose openings are of half the size of the yttrium stabilized zirconia beads (YSZ, Saint Gobain ZirPro, Malvern, PA, USA). A Cole-Palmer peristaltic pump (Master Flex, Vernon Hills, IL, USA) recirculated the suspension between the holding tank and mill chamber at a 126 mL/min flow rate, which was kept the same for all processing runs. The milling conditions are shown in Table 1, where 27 experiments were used as training runs and 5 additional runs were used for testing the model prediction capability. Stirrer speed, bead loading, and bead size were varied at 3 levels for the training runs. The low-high values of the stirrer speed and the bead loading were selected based on our prior wet media milling studies using FNB [36,59], the limitations of our milling equipment, and the objective of preparing drug nanosuspensions with a median particle size below 0.5 µm within 60 min. A stirrer speed lower than 2000 rpm and/or a bead loading lower than 0.4 could result in extremely low breakage rates, and in turn coarser particles with a median size greater than 0.5 µm. Hence, Pharmaceutics 2022, 14, 2840 4 of 19 they were excluded. The design limit of the equipment (4200 rpm) dictates 4000 rpm as the high value with a safety margin. The bead loading above 0.60 could cause significant pressure build-up and mill shut-down, and it is also constrained by the maximum packing limit of the beads (~0.63 for the randomly packed monodispersed spherical beads). The bead size range of 200-800 µm covers the range of the most widely used bead sizes in WSMM (refer to [60,61] and the references cited therein). Whereas the use of beads smaller than 200 µm can be advantageous for the production of sub100 nm particles [60], bead sizes of 50 and 100 µm in size have not been widely used in pharmaceutical manufacturing because of practical clogging issues [62] as well as dust-handling issues. Despite the use of a chiller with an initial temperature of 6.1 • C, the temperature rise during the process was inevitable due to heat generation as the drug suspension with the beads was stirred. Milling was started when the chiller temperature reached 6.1 • C, and the mill outlet temperature was equal to or below 18 • C. Even though keeping the initial temperatures for both the chiller and the mill outlet the same would be a better approach, only the initial chiller temperature could be kept the same; the initial temperature at the mill outlet varied in a narrow range (13-18 • C) because of variation in the ambient temperature, the pre-suspension temperature, and operator practice. During the experiments, the mill outlet temperature was recorded every minute (Runs 1-17) or every 30 s (Runs [18][19][20][21][22][23][24][25][26][27]. The effective milling time was the same for all runs (60 min), whereas the operating time was variable (60-380 min) because of intermitting milling to prevent temperature exceeding the gelation temperature of the polymer (45 • C) [57]. In an intermittent milling cycle, the mill was shut down while cooling continued, and whenever the mill outlet temperature reached 18 • C, the milling continued. Note that we considered only the first milling cycle in the simulations. The average power consumption P was calculated by dividing the total energy consumption read in the mill panel by the effective milling time.

Formulation of the Lumped-Parameter Model (LPM)
During the milling of drug suspensions, heat is generated because of the conversion of mechanical energy input by the stirrer of the mill. The heat generated is removed by a coolant passing through the jacket of the milling chamber. Ignoring the enthalpic effects associated with the suspension recirculation between the holding tank and the milling chamber, we can come up with a simple, low-fidelity model that retains the essential elements of the heat generation-transfer. The difference between the heat generation rate and heat removal rate will cause the internal energy build-up in the mill as milling continues and temperature in the mill rises. Applying the lumped capacity method [63], within the context of a transient enthalpy balance on the mill chamber, we derived the following semi-theoretical lumped-parameter model: where t is milling time, m is the mass in the mill chamber, C p is the specific heat capacity, T is the temperature at the mill outlet, Q gen is the apparent heat generation rate during milling, UA is the apparent overall heat transfer coefficient times surface area, and T ch is the chiller temperature. Strictly speaking, Equation (1) represents a transient enthalpy balance for a perfectly mixed batch process. The perfect mixing implies that the mill outlet temperature is equal to the temperature of the suspension in the mill chamber. The wellmixedness in the milling chamber has been established as a valid approximation to the residence time distribution in small mills (small length-to-diameter ratio) [64]. Hence, for a recirculation mill operating with a fixed batch size and recirculation rate, Q gen and UA may only represent the heat generation rate and overall heat transfer coefficient times surface area, in some approximate, apparent, and statistical manner because they are obtained by fitting to experimental data directly. Although UA can be estimated based on heat transfer correlations for the internal and external convective heat transfer coefficients [33], such correlations are approximate, and none exists for the specific stirrer-mill chamber geometry. We also assumed time-invariant, constant Q gen , C p , UA, and T ch (6.1 • C). Upon separating the variables in Equation (1), integrating both sides, and imposing the initial condition, i.e., t = 0, T = T 0 , the following equation for the time-wise evolution of mill outlet temperature (shortly temperature hereafter) was obtained: Here, m and C p were determined considering the materials in the mill chamber: the beads (zirconia, C p = 0.46 J/g • C [65]), the suspension (10% FNB with respect to water, C p = 3.93 J/g • C), and the stirrer element (zirconia, C p = 0.46 J/g • C). Whereas the stirrer element mass was constant, the bead and suspension mass varied when bead loading was changed in various runs (refer to Table 1). The C p was calculated as the weighted average of the C p of individual materials and the mC p was found to be 465. 6 (2) to the experimental T vs. t data in SigmaPlot 12, Q gen and UA were estimated. Then, these parameters were mathematically expressed as a function of the process parameters for the 27 training runs and predicted as a function of the process parameters for the 5 test runs using a power law (PL) model and a machine learning (ML) model. Minitab was used for the PL predictions, whereas Google Colab was used for the ML predictions. Among the several applied machine learning approaches using Google Colab, as shown in Table S1 of Supplementary Materials, k-nearest neighborhood (KNN) [66] with k = 5 was selected because of its low mean squared error (MSE) and mean absolute error (MAE) compared to other methods for the test runs. Therefore, ML refers to KNN (k = 5) for the rest of this study.

Properties of the Milled Suspensions and Particles
As the scope of this study is the simulation of temperature rise during WSMM via a lumped-parameter model, readers are referred to previous investigations for full characterization of particle sizes, viscosity, crystallinity, and morphology of the particles after milling [33,34]. Here, it suffices to summarize the key findings. All runs yielded nanoparticles upon 60 min milling, where the median particle sizes varied between 149-400 nm (refer to Table S2 of Supplementary Materials). The HPC-SDS combination successfully stabilized the drug nanoparticles by mitigating their aggregation during milling and storage. The nanoparticles were visible in SEM, confirming the laser diffraction results. XRD results of the nanoparticles showed the characteristic peaks of as-received FNB, indicating the crystal structure of the FNB was largely preserved during the milling.

Fitted LPM Parameters and the Origin of Temperature Rise during the Milling
The data on the timewise evolution of the mill outlet temperature was fitted by the LPM, as represented by Equation (2) for each training run (Runs 1-27). The fitted parameters are presented in Table 2 along with the root-mean-squared error (RMSE). The RMSE values ranged between 0.15-0.90 • C. Such low RMSE values suggest that the LPM has excellent fitting or descriptive capability of the temperature profiles despite having only two parameters. Figure 1 demonstrates that the apparent heat generation rate Q gen is linearly and strongly correlated with the average mechanical power consumption P (R 2 = 0.97). The value of the constant slope of the linear correlation in Figure 1 indicates that about 64% of the power consumption (rate of shaft work) dissipates as heat. This is not surprising at all: only a small fraction of the mechanical energy spent on mixing the suspension-bead mixture is used to deform the particles [31]. Most is converted into heat through dissipative processes such as viscous losses, inelastic bead-bead and bead-wall collisions, etc. [67]. Some of the shaft work is also spent on generating new particle surfaces (surface energy), sound, and the elastic parts of bead-bead and bead-wall collisions [33].   Before we delve into the experimental temperature profiles and their fitting by the LPM, let us quickly assess how the temperature rise Trise in the mill at 6 min was affected by the apparent heat generation rate Qgen parameter of the LPM. Based on Equation (1), we expect that Qgen is the driving force for the temperature rise, which was illustrated in Figure 1. Correlation between apparent heat generation rate Q gen and power consumption P. Before we delve into the experimental temperature profiles and their fitting by the LPM, let us quickly assess how the temperature rise T rise in the mill at 6 min was affected by the apparent heat generation rate Q gen parameter of the LPM. Based on Equation (1), we expect that Q gen is the driving force for the temperature rise, which was illustrated in Figure 2. Overall, the temperature rise was more pronounced for higher Q gen or higher P, in view of Figure 1. The Gompertz growth function in Equation (3) fitted the temperature rise well (R 2 = 0.95). As Q gen approaches zero, it predicts a negligibly small temperature rise (~0.7 • C).
Pharmaceutics 2022, 14, x FOR PEER REVIEW 8 of 21 Figure 2. Overall, the temperature rise was more pronounced for higher Qgen or higher P, in view of Figure 1. The Gompertz growth function in Equation (3) fitted the temperature rise well (R 2 = 0.95). As Qgen approaches zero, it predicts a negligibly small temperature rise (~0.7 °C). It is worth mentioning the caveat that the LPM is too simplistic with a multitude of assumptions; therefore, the Qgen values do not reflect the actual heat generation rate; by and large, Qgen is a fitting parameter affected by the accuracy of the experimental measurements and the assumptions made in the model development. On the other hand, Fig-Figure 2. Temperature rise at 6 min as a function of the apparent heat generation rate (Runs 1-25).
It is worth mentioning the caveat that the LPM is too simplistic with a multitude of assumptions; therefore, the Q gen values do not reflect the actual heat generation rate; by and large, Q gen is a fitting parameter affected by the accuracy of the experimental measurements and the assumptions made in the model development. On the other hand, Figures 1 and 2 and the correlations therein strongly associate Q gen with the underlying physics of the conversion of shaft work (power consumption) into heat and ensuing temperature rise. The upshot of these findings is that the LPM differs from a purely empirical model. The latter would fit temperature evolution as a function of time with parameters that have no connection to the physics of the heat generation-transfer phenomena.  (Table 2), the fitted profiles visually corroborate that the LPM has excellent fitting capability despite its simplicity. A cursory look at the experimental temperature profiles suggests that the mill outlet temperature rose during the milling due to the conversion of the shaft work into heat and ensuing heat generation. Although the temperature rise was monotonic, the temperature attained a steady-state value for 2000 and 3000 rpm runs. The heat generation rate was so high at 4000 rpm that the mill was shut down earlier than 60 min and many intermittent milling cycles were conducted. In all profiles, the slope of the temperature profile decreased during the milling. Guner et al. [34] attributed the decreasing rate of temperature rise to a decrease in the instantaneous power consumption during the milling, which originates from the reduction of viscosity at the higher temperatures and particle size reduction during the milling [68,69]. Figures 3-5 also imply that a higher stirrer speed led to higher heat generation rate; stirrer speed is the dominant process parameter, whereas the impact of the bead size is the weakest.    To quantify the impact of the process parameters, motivated by our earlier work [33], a power-law (PL) model and machine learning (ML) model were trained using the directly fitted values of Qgen and UA in Runs 1-27, which were conducted at various stirrer speeds ω, bead loadings c, and the bead sizes Db. For the PL model training, Minitab was used, and Equations (4) and (5) were obtained via fitting. To quantify the impact of the process parameters, motivated by our earlier work [33], a power-law (PL) model and machine learning (ML) model were trained using the directly fitted values of Q gen and UA in Runs 1-27, which were conducted at various stirrer speeds ω, bead loadings c, and the bead sizes D b . For the PL model training, Minitab was used, and Equations (4) and (5)  Equations (4) and (5) signify through their exponents that the stirrer speed and the bead loading had the most significant impact on Q gen and UA, whereas the bead size impact was much weaker. The relative impact of these parameters can be rank-ordered as follows: stirrer speed > bead loading >> bead size. It is well-known that an increase in the stirrer speed and bead loading increases the power consumption in a wet stirred media mill [31,70], which in turn leads to higher Q gen . Similarly, it is well-known that an increase in stirrer speed and bead loading also leads to an increase in the internal convective heat transfer coefficient in the mill chamber, which could lead to a higher UA [71,72].

LPM-Fitted Temperature Profiles and LPM-PL/LPM-ML Predictions in the Training Runs
By augmenting the PL model, Equations (4) and (5), and the ML model (KNN) with the LPM, we predicted the temperature profiles and compared them with the experimental profiles as well as the profiles generated by direct fitting with the LPM alone (see . The associated mean squared error (MSE) and mean absolute error (MAE) are reported in Table S1 of Supplementary Materials. The predictions by the LPM-PL and the LPM-ML deviated from the experimental data more than the direct fits by the LPM alone, which is intuitively expected. These deviations can be reduced by increasing the number of training runs, i.e., additional experiments at other stirred speeds and bead loadings. Both the LPM-PL and the LPM-ML predicted the profiles very well in many runs (e.g., Runs 2, 12, and 21). The LPM-PL predictions were generally closer to the experimental profiles than the LPM-ML predictions, with a few notable exceptions for each stirrer speed. The maximum deviation between the LPM-PL prediction and the experimental data ranged from 1.5 to 6.5 • C, with a mean and standard deviation of 3.1 ± 1.3 • C. Overall, the LPM successfully captured some salient qualitative patterns of the temperature profiles: (i) the monotonic temperature increase with a decreasing rate, (ii) attainment or approach to a steady-state temperature, (iii) the drastic decrease of time to reach 45 • C at 4000 rpm, and (iv) the relative impact of the process parameters.

Comparative Analysis of LPM and EBM Fits and Their Predictions for the Test Runs
We reserved Runs 28-32 data for testing the simple LPM in comparison to the more elaborate EBM. We first fitted the LPM and the EBM to the experimental temperature profiles directly, estimated Q gen and UA (see the RMSE in Table 3), and illustrated the fitted profiles in Figure 6. The EBM data were retrieved from Guner et al. [33] for comparison. Figure 6 and Table 3 data suggest that the LPM fitted the experimental temperature profiles slightly better than the EBM. The average ± standard deviation of the RMSEs are 0.50 ± 0.12 • C for the LPM and 0.96 ± 0.43 • C for the EBM. Moreover, LPM has a lower deviation from the experimental data and is more consistent, according to the lower standard deviation. Table 3. Parameters of the LPM estimated by direct fitting as well as predicted using the PL and ML models coupled to the LPM along with the associated statistics.

Direct Fitting PL Prediction ML Prediction
Q gen (J/min)    The LPM-PL and LPM-ML model predictions are compared with the EBM-PL and the EBM-ML predictions, which were retrieved from Guner et al. [33], as seen in Figures 7 and 8. Whereas all LPM predictions approached an asymptote as the experimental data did, the EBM predictions tended to have a maximum temperature (Figure 7). Moreover, there are better and worse prediction examples when the LPM and EBM predictions are compared. The average ± standard deviation of the RMSEs are 1.57 ± 0.99 • C for all LPM predictions and 1.66 ± 1.00 • C for all EBM predictions. Overall, the LPM is better for predictions but based on the average RMSEs, the difference is not as drastic as it is when the RMSEs associated with direct fittings are compared. The LPM-PL and LPM-ML model predictions are compared with the EBM-PL and the EBM-ML predictions, which were retrieved from Guner et al. [33], as seen in Figures  7 and 8. Whereas all LPM predictions approached an asymptote as the experimental data did, the EBM predictions tended to have a maximum temperature (Figure 7). Moreover, there are better and worse prediction examples when the LPM and EBM predictions are compared. The average ± standard deviation of the RMSEs are 1.57 ± 0.99 °C for all LPM predictions and 1.66 ± 1.00 °C for all EBM predictions. Overall, the LPM is better for predictions but based on the average RMSEs, the difference is not as drastic as it is when the RMSEs associated with direct fittings are compared.    [33].
The average ± standard deviation of the RMSEs are 1.79 ± 1.35 °C and 1.36 ± 0.54 °C, respectively, for the LPM-PL and LPM-ML predictions. The LPM-PL exhibited better or similar predictive capability for most of the test runs, as can be seen from Figure 8, which represents the comparative RMSEs visually. However, the LPM-PL prediction for Run 32 was notably bad. A similarly bad prediction was made by the EBM-PL for the same run. Interestingly, despite its larger deviations for the training runs, when augmented with either the LPM or the EBM, the prediction with the ML model had much lower RMSE than that with the PL model for Run 32. The process conditions of this run were purposefully chosen to be outside the domain of the training runs to test the models under extreme cases. They are not likely used in a lab or industrial setting due to unrealistically low bead loading and very small bead sizes. In general, small beads, e.g., 100 µm, are more difficult to handle compared with 300-600 µm beads during the operations. Nonetheless, these findings suggest that either the EBM or the LPM should be used with caution outside the domain of the training runs.

The LPM and the EBM Comparison and the Limitations of the LPM
The full EBM developed by Guner et al. [33] is a comprehensive model of the heat generation-transfer in WSMM that considers recirculation and batch size as well as the cooling rate provided by the jackets of the milling chamber and the holding tank. It entails using the values of power consumption, physico-chemical properties of the fluid and the beads, and the dimensions of the mill setup, which were needed for the calculation of the overall heat transfer coefficient U and heat transfer surface area of the mill chamber Am (refer to [33] for details of the EBM). The EBM is so versatile that it can be used to investigate the impacts of different coolant types and flow rates as well as the material of construction of the mill on the cooling rate and temperature evolution. Despite all these capabilities and higher fidelity to the actual milling process, its use entails more time, effort, and accurate numerical methods for the simulations/parameter estimation. Moreover, one must obtain appropriate data and correlations for the physical-thermal-heat transfer properties. The EBM consists of five ordinary differential equations (ODEs) with the fraction ξ of the power consumption P that is converted into heat being the sole parameter of the EBM. The estimation of ξ entails using a sophisticated optimizer coupled to the ODE The average ± standard deviation of the RMSEs are 1.79 ± 1.35 • C and 1.36 ± 0.54 • C, respectively, for the LPM-PL and LPM-ML predictions. The LPM-PL exhibited better or similar predictive capability for most of the test runs, as can be seen from Figure 8, which represents the comparative RMSEs visually. However, the LPM-PL prediction for Run 32 was notably bad. A similarly bad prediction was made by the EBM-PL for the same run. Interestingly, despite its larger deviations for the training runs, when augmented with either the LPM or the EBM, the prediction with the ML model had much lower RMSE than that with the PL model for Run 32. The process conditions of this run were purposefully chosen to be outside the domain of the training runs to test the models under extreme cases. They are not likely used in a lab or industrial setting due to unrealistically low bead loading and very small bead sizes. In general, small beads, e.g., 100 µm, are more difficult to handle compared with 300-600 µm beads during the operations. Nonetheless, these findings suggest that either the EBM or the LPM should be used with caution outside the domain of the training runs.

The LPM and the EBM Comparison and the Limitations of the LPM
The full EBM developed by Guner et al. [33] is a comprehensive model of the heat generation-transfer in WSMM that considers recirculation and batch size as well as the cooling rate provided by the jackets of the milling chamber and the holding tank. It entails using the values of power consumption, physico-chemical properties of the fluid and the beads, and the dimensions of the mill setup, which were needed for the calculation of the overall heat transfer coefficient U and heat transfer surface area of the mill chamber A m (refer to [33] for details of the EBM). The EBM is so versatile that it can be used to investigate the impacts of different coolant types and flow rates as well as the material of construction of the mill on the cooling rate and temperature evolution. Despite all these capabilities and higher fidelity to the actual milling process, its use entails more time, effort, and accurate numerical methods for the simulations/parameter estimation. Moreover, one must obtain appropriate data and correlations for the physical-thermal-heat transfer properties. The EBM consists of five ordinary differential equations (ODEs) with the fraction ξ of the power consumption P that is converted into heat being the sole parameter of the EBM. The estimation of ξ entails using a sophisticated optimizer coupled to the ODE solver. Hence, for facile modeling of the temperature profiles as well as effective control and optimization of the WSMM process, development of simpler, low-fidelity models such as the LPM may be warranted.
In contrast to the EBM, the LPM is a simple semi-theoretical model, with two adjustable parameters, which was developed based on various assumptions mentioned in Section 2.2.2. Hence, it has a closed-form analytical solution and can be easily used without much effort and time, while obviating the need for sophisticated ODE solvers/optimizers. The main limitation of the LPM is that it is only valid for a given mill set-up with constant batch size and recirculation rate. This aspect may not be critical because, to the best knowledge of the authors, these parameters are usually kept invariant whereas the impacts of the stirrer speed, the bead loading, the bead size, and milling time have been investigated in most WSMM process studies (see e.g., the reviews [15,16]). Another limitation of the LPM is that although it can predict the temperature of the suspension in the milling chamber, it does not provide any information about the temperatures of the suspension in the holding tank, temperatures of the beads, the stirrer of the mill chamber, and the stirrer of the holding tank. Contrarily, the EBM predicted a lower temperature of the suspension in the holding tank than that in the milling chamber and established a thermal equilibrium among the suspension, the stirrer element, and the beads, owing to fast convective heat transfer and relatively small size of the beads and the mill stirrer [33].
Despite its simplicity, when augmented with the PL and the ML models, the LPM predictions of the temperature profiles in the test runs were better than or similar to the EBM predictions (refer to Figures 7 and 8). Moreover, the LPM predicts a monotone increasing temperature profile with a steady-state temperature approach for sufficiently long milling, which is in excellent qualitative agreement with the profiles in Figures 3-7. In fact, Equation (2) clearly shows that in the limit t→∞, the temperature reaches a steady-state temperature of T ch + Q gen /UA. Although the EBM has higher fidelity to the real WSMM operation, its predictions could be worse, and some predictions exhibited a maximum and ensuing drop in temperature instead of a monotonic approach to a steady-state in the temperature profiles (Figure 7). In general, the temperature drop from the maximum is within a couple of degrees Celsius, and this error was acceptable. Note that even the EBM has its own assumptions and sources of modeling errors; refer to Guner et al. [33] for a detailed discussion of the modeling errors. For example, in the UA m calculations, the mixture correlations for the physical properties and the internal/external convective heat transfer coefficients have some errors. The LPM is free of that source of modeling error because UA was used as a fitting parameter; therefore, having two fitting parameters, as opposed to the one fitting parameter of the EBM, enhanced LPM's fitting and prediction capability.
The main drawback of the LPM is that it does not consider the enthalpic effects associated with the recirculation of the drug suspension and the thermal inertia effects associated with the batch volume of the suspension in the holding tank. Precisely because of this simplification, the LPM is referred to as a semi-theoretical model here, which considers some physical aspects of the process in some time-average, approximate, and statistical sense. Its Q gen and UA are fitting parameters that are not equal to the true heat generation rate and the product of the overall heat transfer coefficient and the heat transfer surface area. In fact, the actual heat generation rate and even the overall heat transfer coefficient vary with time and temperature [33,34]. However, as established in this study, Q gen is strongly and positively correlated with the power consumption and is the driver for temperature rise, whereas the UA correlation with the process parameters revealed a similar qualitative dependence of the convective heat transfer coefficients on the process parameters in their correlations (refer to such correlations in [33]).
The LPM's treatment of the recirculation operation as an equivalent batch operation seems unreasonable for the modeling of temperature profiles on purely theoretical grounds. Being aware of this limitation, we accept the resulting modeling error for the sake of simplicity as LPM is intended to be a low-fidelity model for industrial use. It is also worth mentioning that a similar approach has already been adopted for the modeling of the evolution of the median particle size and even the whole particle size distribution during the recirculation operation of the WSMM by multiple research groups (e.g., [29,73,74]).

A holistic Perspective on the Impact of Process Parameters
This study has focused on the impact of process parameters on heat generation and temperature evolution during the WSMM and a comparative analysis of the newlydeveloped LPM and the recently developed EBM. In a previous study, Guner et al. [34] investigated the impact of the process parameters on the particle sizes of the milled suspensions (refer to Table S2). An increase in the stirrer speed and bead loading led to smaller d 50 and d 90 , whereas an increase in the bead size led to bigger d 50 and d 90 . As the increase in the stirrer speed and bead loading also leads to higher heat generation rate and temperature rise, existence of an optimal set of process conditions is anticipated. Clearly, the milling conditions that are conducive to higher milling efficiency also cause higher heat generation and temperature rise. This is not surprising as both of them are largely determined by the mechanical power consumption during milling. Guner et al. [34] considered a desirable product specification of d 10 < 150 nm, d 50 < 200 nm, and d 90 < 250 nm and max. temperature below 37 • C. Their holistic consideration of a newly-defined thermal desirability score, power (energy) consumption, and total cycle time suggested that Run 13 or 14 had the optimal set of milling conditions: 3000 rpm with 50% loading of 200 or 400 µm beads. Obviously, in general, the optimal conditions depend on the specific pharmaceutical application of the drug nanosuspension with desired particle size specifications as well as the sensitivity of the drug to temperature and stressing in terms of physical stability and chemical degradation.

Conclusions and Future Outlook
This study has developed and implemented a semi-theoretical, lumped-parameter model (LPM) as an alternative to the recently developed enthalpy balance model (EBM) for simulating and predicting the temperature evolution during the nanomilling of drug suspensions. It has two fitting parameters which make it more flexible compared to the enthalpy balance model; therefore, it provides better fitting results. Although the apparent heat generation rate and the apparent heat transfer coefficient multiplied by the heat transfer surface area are not equal to the actual values, they could be predicted by process conditions via the machine learning (ML) and the power-law (PL) approaches. The fittings and predictions of the LPM were found to be slightly better than those of the EBM. Overall, our experimental and modeling results suggest that the LPM has excellent descriptive/fitting capability in addition to its reasonably good predictive capability. Coupled with its simplicity, which obviates the need for using a sophisticated coupled optimizer-ODE solver, it could be selected for facile modeling of the WSMM process. Hence, we provide the pharmaceutical engineering literature with two different models, the LPM and the EBM, which can be used on a fit-for-purpose basis. If a quick analysis and modeling of the temperature profiles are needed without significant effort/time, then the LPM is advantageous and should be used. The LPM can also be easily used for process control. However, if the aim is to intensify and optimize the process that entails a detailed and deep understanding of the impact of the recirculation rate, batch size, cooling type/capacity, and material of construction, then the EBM must be used. In that case, power consumption and the physico-chemical and thermal properties of the suspensions and the beads, as well as reliable correlations for the convective heat transfer coefficient, must be obtained from the literature and/or determined experimentally. Future work will involve applying the LPM and EBM to pharmaceutical wet media milling processes wherein (i) polymeric beads as opposed to ceramic beads are used, (ii) a large batch size (e.g., 7 L) of drug suspensions is to be milled, and (iii) pilot and large-scale commercial wet media mills are used.