A new mathematical model for monitoring the temporal evolution of the ice crystal size distribution during freezing in pharmaceutical solutions

Abstract

The freezing step plays a key role in the overall economy of the vacuum freeze-drying of pharmaceuticals, since the nucleation and crystal growth kinetics determine the number and size distribution of the crystals formed. In this work, a new mathematical model of the freezing step of a (bio)pharmaceutical solution is developed and validated. Both nucleation and crystal growth kinetics are modeled and included in a one-dimensional population balance (1D-PBM) that describes, given the product temperature measurement, the evolution of the pore size distribution during freezing. The developed model is coupled with the real-time measurements obtained from an infrared video camera. The ending 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, resulted in good agreement with experimental values. The resulting Process Analytical Technology (PAT) has the potential to boost the development and optimization of a freeze-drying cycle and the implementation of a physically grounded Quality-by-Design approach in the manufacturing of pharmaceuticals. A more general mathematical model, including the aforementioned population balance, of a vial filled with a solution of sucrose was also developed and used to further validate the approach.

Introduction

Freeze drying is a low-temperature drying process, particularly appreciated in comparison to other drying processes, as it avoids damage to thermally sensitive Active Pharmaceutical Ingredients (API). In a typical freeze-drying cycle, the liquid formulation is first poured into single-dose containers and frozen, by cooling the shelves upon which the vials are placed. Cooling a sample below its equilibration freezing point induces crystallization of the solvent, typically water, in the liquid formulation. Some of the original solvent remains in a liquid state, bound to the solid structure. Solvent removal is carried out in two phases, namely, primary and secondary drying. In the former, the pressure is reduced below the triple point of water, while the temperature of the shelves is increased to allow the sublimation of the ice crystals formed; in the latter, since the amount of bound water is usually greater than the allowed residual moisture content in the final product, the temperature is further raised to enhance its desorption [1], [2], [3], [4].

In recent years, the number of drugs and pharmaceutical products that require a freeze-drying step in the manufacturing process, for solvent removal and drug stabilization, has increased. Also the requirements in terms of efficiency, quality, and safety have become more stringent, which has motivated efforts to improve the understanding and design of the freezing process [5], [6]. During the subsequent drying, ice crystals will sublime to leave voids within a porous network through which the water vapor will flow on its way from the sublimation interface to the drying chamber. A structure characterized by bigger pores enables a higher sublimation flux [7], due to the lower transport resistance of the dried product, while having an opposite effect on the rate of the secondary drying, because of a lower surface-area-to-volume ratio [8]. As the porous structure directly influences the mass transfer rate, 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 of the product is a critical variable that must be carefully monitored to avoid jeopardizing product quality or impairing the APIs. The formation of small ice crystals, increasing the interfacial area, might also affect the unfolding and aggregation of proteins [9].

The prediction of the crystal size distribution resulting from the freezing of the liquid formulation is of utmost importance for process optimization, not only in the freeze drying of pharmaceutical solutions, and has been deeply investigated in the past years. In the literature, many empirical relationships can be found to correlate the thermal evolution of the product to the average pore size. Bald [10] first proposed to relate mean pore diameter to the cooling rate of the product by

$${\overline{D}}_{p,i}=\alpha {\left(\frac{\mathit{dT}}{\mathit{dt}}\right)}^{-\beta }$$

where α and β are empirical constants that must be fitted to experimental results.

A whole family of empirical models is also available in the literature which relates the average pore diameter to the velocity of the solidification front, P, and the temperature gradient in the solid layer, Q:

$${\overline{D}}_{p,i}={\mathrm{\alpha }}^{\text{'}}{P}^{-{\lambda }_1}{Q}^{-{\lambda }_2}$$

where α’, λ₁ and λ₂ have to be determined experimentally. In fact, the parameter Q is representative of the rate of removal of the latent heat released during crystal formation; the resulting crystal structures are dependent from this. The freezing front velocity instead is representative of the salt rejection rate at the crystals interface which is one of the phenomena governing crystal growth. A non-exhaustive list of applications of Eq. (2) includes freezing of apples [11], starch gels [12], gelatin gels [13], solidification of metal samples at low rates [14] and solidification of alloys at high rate [15]. Eq. (2), with λ₁ = λ₂ = 0.5, has been also used to describe the freezing of aqueous solutions in vials [16], [17], [18], [19]. In the latter case, P and Q are calculated from

$$P\left(t\right)=\frac{d{L}_{\mathit{frozen}}\left(t\right)}{\mathit{dt}}$$

$$Q\left(t\right)=\frac{T_f-T_B\left(t\right)}{L_{\mathit{frozen}}\left(t\right)}$$

The strong dependence of the many parameters on the kind of product, the specific application, and operating conditions that appear in these empirical models limits the predictive reliability and accuracy of the model. Besides, the parameters P and Q provide a macroscopic description of some of the many aspects related to the dynamics of the system. They are consistent with the physics of the system, but they do not provide a characterization of the system which can be given only in terms of thermodynamic state functions.

The state of the art in the modeling of the dependency of the average pore diameter on the thermal evolution of the product is the mechanistic model proposed by Arsiccio et al. [20]. The basic idea underlying their approach is that crystal growth is an exothermic process: the incorporation of new water molecules to the existing crystals, together with the creation of new interfaces, will release energy. The cooling medium and surrounding environment remove some of this energy and, thus, an energy balance in an infinitesimal slice of the material is related to the associated crystal growth dynamics. This approach neglects the nucleation kinetics and, as all the other models previously discussed, can predict the average pore diameter but not the pore size distribution.

In this work, a different approach to modeling the freeze-drying process is proposed and validated. The nucleation kinetics is described by a stochastic model in the form of a chemical Master equation. The parameters of the nucleation kinetics are fit to achieve the best agreement between the simulated and experimental induction times. Given the number and probability distribution function of the ice seed crystals formed, their evolution over time is described by a 1D population balance. The crystal growth rate is modelled as a function of the thermodynamically correct driving force, which is the supersaturation. This model does not require any additional variable to be estimated or measured; whenever a single or multiple measurement of the temperature inside the material is provided [21], the model can be solved and both the whole pore size distribution and the amount of bound water can be determined. An analytical solution of the proposed model is derived using the Method of Characteristics [22], making this approach suitable for on-line monitoring and optimization of the whole freeze-drying process.

The paper is thus organized in five more sections. Section 2 describes the stochastic model and the experimental studies performed to derive the nucleation kinetic parameters, while Section 3 presents the 1D population balance model, the assumptions made, and the derivation of the analytical solution. Chapter 4 presents the experimental setup and the case study discussed. Chapter 5 describes the methodology for the validation of the proposed model, including the mathematical simulation of the whole freezing step of a single dose container (the 2D model used is presented in Appendix A), while Chapter 6 presents the main results of the study. Chapter 7 discusses the main conclusions and future development of the proposed approach.