(471e) Mathematical Modeling of the Freezing of Pharmaceutical Solutions

Mathematical Modeling of the Freezing of Pharmaceutical Solutions

A convenient formulation for (bio)pharmaceutical drugs is as a freeze-dried product. The product is formed by lowering the temperature of an aqueous solution to below freezing to generate ice, and then a vacuum is applied to sublimate the ice and remove water in the form of vapor. This process is carried out at low temperatures and is suitable for thermally sensitive drugs, and a freeze-dried product is less costly to store and transport than use of freezing alone. Freeze-dried products can be reconstituted quickly and easily.

A typical freeze-drying cycle consists of three phases, namely freezing, primary drying, and secondary drying. The product is firstly poured into single-dose containers, which are then loaded onto the shelves of the chamber of a freeze drier. A technological fluid flows inside the shelves and removes the heat from the product, thus inducing the crystallization of the (free) water; part of the water remains bound to the product molecules without crystallizing. Then, in the primary and secondary drying phases, respectively, the frozen and bounded water is removed. In the primary drying, the pressure is reduced well below the triple point of water, and heat is provided to sublimate the ice crystals. In secondary drying, the bounded water absorbed into the solid structure is desorbed by increasing the temperature until a desired value of the relative humidity is reached.

The freezing step plays a key role in the overall process, as the nucleation and crystal growth kinetics determine the number and size distribution of crystals formed. During the sublimation stage, the water vapor flows from the frozen product to the drying chamber through the porous structure. A highly porous structure enables a higher sublimation flux, due to the lower transport resistance of the dried product, while having an opposite effect on the rate of the secondary drying, as the surface area to volume ratio is lower. As the porous structure directly influences the drying rates, and product temperature results from a thermal balance between the heat transferred to the product and the heat removed through sublimation, the porous structure also affects the spatial temperature distribution during the drying stages. The temperature is a critical variable that must be carefully monitored to avoid jeopardizing product quality.

A full understanding of the freezing step would enable a real-time prediction of the axial distribution of the pores inside the product, which would advance the monitoring, control, and optimization of the freeze-drying process. The nucleation of ice crystals is intrinsically stochastic, which can induce high variability of the porous microstructure of the material inside the vials. The differences in the nucleation temperatures experienced by the different vials are one of the main sources of in-batch variability of a freeze-drying process.

Empirical approaches have been proposed for estimation of the porous structure from the variables that are directly accessible. A series of empirical formulas have been proposed to link the pore diameters to the velocity of the freezing front and the thermal gradients inside the frozen layer. The empirical dependency of the parameters in these models on the application, product, and process conditions – together with the lack of any physical foundation in this approach – makes the results obtained neither reliable nor accurate. Arsiccio et al. (2018) first proposed a mechanistic model, based on an energy balance applied to a slice of the product, for the prediction of the pore average diameters once the information on the thermal evolution of the process is available. This approach, although much more physically grounded than those previously discussed, neglects the nucleation rate and only focuses on the crystal growth kinetics. Furthermore, the model does not provide a full characterization of the pore distribution.

This presentation develops and validates a new mathematical model of the freezing step of a (bio)pharmaceutical solution. Both nucleation and crystal growth kinetics are modeled and included in a one-dimensional population balance that describes, given a product temperature measurement, the evolution of the pore size distribution. The kinetic parameters are estimated from least-square fitting on the experimentally measured induction times obtained on a set of seven batches of ten vials each. The nucleation rate is modelled using a stochastic model in the form of a Master equation, following the approach of Goh et al. (2010), which assumes that the nucleation of crystals is a nonstationary Poisson process.

The experiments were performed inside laboratory-scale equipment (Telstar, Spain) and, in each batch ten vials, 10R (ISO 8362-1) are filled with 5 ml of a 5% b.w. solution of sucrose. Each test was 2.5 h long and the cooling policy was to simply cool to –50°C as fast as possible. The vials were placed inside the chamber and monitored using the sensor presented by Lietta et al. (2019) based on infrared imaging.

The crystal growth rate was experimentally observed to be linearly dependent on the chemical potential difference between the ice and the liquid solution. The chemical potential of the solid phase was modelled according to Corti et al. (2010), while the parameters for the UNIQUAC model of the sucrose solution were taken from Cattè et al. (1994). Experimental evidence indicated fast nucleation, with the material changing its optical properties (turbidity) within less than one second. To model this behavior, nucleation was implemented as a boundary condition on the population balance. An analytical solution of the model was derived using the method of characteristics, making the model suitable for real-time application.

A more general mathematical model, including the aforementioned population balance, of half a vial filled with a solution of sucrose was developed and used to validate the approach. A modified porous media approach was used to model the heat and mass transfer occurring in the system. The freezing rate (that is, the mass rate of conversion of water into ice), which couples the solution depletion in the mass balance and the heat generation, was calculated from the third moment of the number density function obtained from the population balance. The nonlinear set of differential equations was solved via a finite volume numerical scheme. Both the simulated temperature profile and the forecast pore structure were validated. The spatiotemporal evolution of the temperature of the product in different axial position across the vial wall was extracted from real-time infrared video images, while the experimental pore distribution was obtained from the analysis of Scanning Electron Microscope (SEM) images of two vials after being completely freeze-dried. In all cases, a good agreement was observed between model and experiments. When the equilibrium was reached between the technological fluid at –50°C inside the shelves and the product, the model predicts a 18% b.w. water bounded to the solid structure. Experimental measurements ranging between 16% to 20% are reported in literature. The end time of the primary drying stage, and the maximum temperature inside the material, simulated through a simplified model of the process and the pore distribution forecast, were in good agreement with the experimental values.

The purpose of the developed model was for coupling with the real-time measurements obtained from the infrared video camera. The resulting Process Analytical Technology (PAT) has the potential of accelerating the development and optimization of a freeze-drying cycle and the implementation of a physically grounded Quality-by-Design approach for the manufacturing of pharmaceuticals.

Future developments of this work should validate the performance of this model on different case studies, namely active pharmaceutical ingredients, cooling strategies, and controlled nucleation techniques such as ice fog and vacuum- and pressure-induced nucleation. Also of interest is the development of a deeper understanding and modeling of the dynamics of ice crystal nucleation – including its unstable propagation in a short period of time – and of asymmetries observed in the pore distribution of the drug product.

Keywords: Master equation; mathematical modeling; pore distribution.

References:

Arsiccio, A., Barresi, A.A., Pisano, R. (2017) Prediction of ice crystal size distribution after freezing of pharmaceutical solutions. Crystal Growth & Design 17, 4573-4581.

Catté, M., Dussap, C.G., Achard, C., Gros, J.B. (1994) Excess properties and solid-liquid equilibria for aqueous solution of sugar using a UNIQUAC model. Fluid Phase Equilibria 96, 33-50.

Corti, H.R., Angell, C.A., Auffret, T., Levine, H., Buera, M.P., Reid, D.S., Roos, Y.H., Slade, L. (2010) Empirical and theoretical models of equilibrium and non-equilibrium transition temperatures of supplemented phase diagrams in aqueous systems. Pure and Applied Chemistry 80, 1065-1097.

Goh, L., Chen, K., Bhamidi, V., He, G., Kee, N.S., Kenis, P.J.A., Zukoski, C.F., Braatz, R.D. (2010) A Stochastic model for nucleation kinetics determination in droplet-based microfluidic system. Crystal Growth & Design 10, 2515-2521.

Lietta E., Colucci D., Distefano G., Fissore, D. (2019) On the use of IR thermography for monitoring a vial freeze-drying process. Journal of Pharmaceutical Sciences,108, 391-398.