Next Article in Journal
Assessing the Effects of VEGF Releasing Microspheres on the Angiogenic and Foreign Body Response to a 3D Printed Silicone-Based Macroencapsulation Device
Next Article in Special Issue
Spin Freezing and Its Impact on Pore Size, Tortuosity and Solid State
Previous Article in Journal
Improving Membrane Activity and Cargo Delivery Efficacy of a Cell-Penetrating Peptide by Loading with Carboranes
Previous Article in Special Issue
Investigating Alternative Container Formats for Lyophilization of Biological Materials Using Diphtheria Antitoxin Monoclonal Antibody as a Model Molecule
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development and Application of a Mechanistic Cooling and Freezing Model of the Spin Freezing Step within the Framework of Continuous Freeze-Drying

1
Laboratory of Pharmaceutical Process Analytical Technology, Department of Pharmaceutical Analysis, Faculty of Pharmaceutical Sciences, Ghent University, 9000 Ghent, Belgium
2
Escuela Profesional de Química, Facultad de Ciencias, Universidad Nacional de Ingeniería, Puerta 5—Av. Tupac Amaru N° 210 Rimac, Lima 15333, Peru
3
Pharmaceutical Engineering Research Unit, Department of Pharmaceutical Analysis, Faculty of Pharmaceutical Sciences, Ghent University, 9000 Ghent, Belgium
4
Laboratory of Pharmaceutical Technology, Department of Pharmaceutics, Faculty of Pharmaceutical Sciences, Ghent University, 9000 Ghent, Belgium
5
Coriolis Pharma, Fraunhoferstraße 18 b, 82152 Martinsried, Germany
6
RheaVita, Frieda Saeysstraat 1, 9052 Zwijnaarde, Belgium
*
Authors to whom correspondence should be addressed.
Pharmaceutics 2021, 13(12), 2076; https://doi.org/10.3390/pharmaceutics13122076
Submission received: 27 October 2021 / Revised: 22 November 2021 / Accepted: 26 November 2021 / Published: 3 December 2021
(This article belongs to the Special Issue New Trends in Freeze-Drying of Pharmaceutical Products)

Abstract

:
During the spin freezing step of a recently developed continuous spin freeze-drying technology, glass vials are rapidly spun along their longitudinal axis. The aqueous drug formulation subsequently spreads over the inner vial wall, while a cold gas flow is used for cooling and freezing the product. In this work, a mechanistic model was developed describing the energy transfer during each phase of spin freezing in order to predict the vial and product temperature change over time. The uncertainty in the model input parameters was included via uncertainty analysis, while global sensitivity analysis was used to assign the uncertainty in the model output to the different sources of uncertainty in the model input. The model was verified, and the prediction interval corresponded to the vial temperature profiles obtained from experimental data, within the limits of the uncertainty interval. The uncertainty in the model prediction was mainly explained (>96% of uncertainty) by the uncertainty in the heat transfer coefficient, the gas temperature measurement, and the equilibrium temperature. The developed model was also applied in order to set and control a desired vial temperature profile during spin freezing. Applying this model in-line to a continuous freeze-drying process may alleviate some of the disadvantages related to batch freeze-drying, where control over the freezing step is generally poor.

Graphical Abstract

1. Introduction

Freeze-drying or lyophilization is a low-temperature drying method under vacuum conditions used for aqueous drug solutions with poor stability [1]. Freeze-drying is divided in three sequential steps: freezing, primary drying, and secondary drying. Freezing is the initial step during which the aqueous solution is solidified, allowing sublimation of the ice and desorption of unfrozen water in the subsequent drying steps. First, the solution is cooled down until ice nucleation occurs, i.e., the moment when the first ice crystal is formed, generally at a temperature several degrees below the equilibrium freezing temperature. The exothermic nature of ice crystallization results in an instantaneous increase in product temperature to the equilibrium freezing temperature (see below, Figure 1). Further cooling of the product leads to further exothermic ice crystal growth, characterized by a slowly decreasing overall product temperature as a result of a growing ice layer. The product temperature decreases sharply once the aqueous solution is completely frozen (see below, Figure 1). The freezing step is considered complete when the product temperature reaches a sufficiently low value that allows for further product processing with a minimal risk of exceeding the eutectic temperature ( T e ) or the glass transition temperature of the maximum freeze-concentrated solution ( T g ).
The ice nucleation temperature and the kinetics of ice crystal growth determine the physical state and morphology of the frozen solution and, consequently, some of the final properties (e.g., the pore size distribution) of the freeze-dried product [2]. Ice nucleation is a stochastic phenomenon that is related to process and product properties such as container surface roughness and the amount of suspended particles present in the solution [3,4]. Due to the stochastic nature of ice nucleation, freezing in a conventional batch freeze-drying process presents significant disadvantages, as it causes uncontrolled product quality variations from vial to vial within the same batch and between batches. Therefore, much effort has been made to develop methods for controlled nucleation in batch freeze-drying, with the aim of achieving a more similar degree of supercooling for all vials [2,5,6,7,8,9,10,11,12,13]. The degree of supercooling is defined as the difference between the equilibrium freezing temperature and the temperature at which the first ice nucleus is formed (i.e., primary nucleation) [2,13]. Secondary nucleation follows primary nucleation, and secondary ice nuclei are formed until the equilibrium temperature is reached [4]. A higher degree of supercooling results in a larger amount of heat freed by exothermic crystallization upon ice nucleation. Hence, a larger fraction of the water is instantaneously frozen, resulting in the formation of more ice nuclei during secondary nucleation compared to ice nucleation at a higher temperature (i.e., a lower degree of supercooling). It is expected that the amount of ice nuclei determine the amount of resulting ice crystals after complete solidification, as during solidification the existing ice nuclei grow, but no new ice nuclei are formed [14]. In this way, a lower degree of supercooling generally leads to a smaller amount of larger ice crystals, resulting in larger pores in the dried product layer during freeze-drying [13]. Additionally, the temperature-controlled shelves most often cool down at a set constant rate until a final temperature during freezing. In this way, more time will have passed before nucleation in the case of a high degree of supercooling compared to a low degree of supercooling. Hence, in the case of high supercooling the shelves will be colder at the moment of nucleation compared to low supercooling. The product equilibrium temperature will remain unchanged (e.g., at 0 °C in the case of pure water) in both cases, resulting in a larger temperature difference between shelf and product in the case of higher supercooling. This increased temperature gradient is associated with faster heat removal during ice-crystal growth, yielding ice crystals with a different morphology (i.e., needle-like ice crystals for fast heat removal and dendritic crystals for slow heat removal) [14]. This difference in ice-crystal morphology may lead to relevant differences in the dried product characteristics (e.g., dried product resistance to water vapour). In essence, for batch freeze-drying, the uncontrolled freezing and the stochastic nature of ice nucleation inherently causes vial-to-vial variability in product characteristics within a batch and between batches [13].
Taking the above into account, it is evident that a good understanding and control of the freezing step and its different phases (i.e., liquid cooling, nucleation, ice crystallization, and subsequent solid cooling) during freeze-drying is of importance. In practice, model-based approaches are often utilized to gain process knowledge and to develop control mechanisms [3,15,16,17,18,19]. Several approaches have been proposed to model the freezing step in batch freeze-drying, such as the one developed by Nakagawa et al. [3,15]. They proposed a two-dimensional axisymmetric cooling and freezing model based on the conductive heat equation (i.e., Fourier’s law of heat conduction). The modelling was divided into two phases: the cooling step of the aqueous solution before ice nucleation, followed by the freezing step including the nucleation event and ice-crystal growth.
The cooling step was simulated using the conductive heat equation expressed in Equation (1) [3,15]. This equation describes the conductive heat flux in relation to the physical properties of the product (i.e., the mass density, the specific heat capacity, and the thermal conductivity), as well as the temperature distribution of the product under study.
ρ c p ( δ T δ t ) = ( k T )
Here, ρ is the solution mass density (kg/m3); c p describes the specific heat capacity of the solution (J/(kg K)); k is the thermal conductivity (W/(m K)); and T is the temperature distribution of the system (K).
Afterwards, the freezing step starts at the moment ice nucleation occurs until ice-crystal growth is completed. The second modelling phase is based on the same base equation but includes two source terms corresponding to the total latent heat released due to ice nucleation, Q n , and the total latent heat released due to ice crystallisation, Q c , as given in Equation (2) [3,15].
ρ c p ( δ T δ t ) = ( k T ) + Q n + Q c
The positive source term Q n was estimated as follows:
Q n = Δ H f k i ( T e q T s )
where Δ H f represents the heat released due to ice crystallisation (J/kg); k i describes the nucleation rate (kg/(m3 s K)); T e q is the equilibrium freezing temperature (K); and T s is the temperature in the supercooled liquid (K).
The nucleation rate k i (kg/(m 3 s K)) was calculated using the freezing front velocity ν (m/s), the thickness of the undercooled zone s (m), and the homogeneous undercooled temperature T (K). [3]:
k i = ρ ν s ( T e q T )
The positive-source term Q c corresponding to the latent heat of crystallization was estimated using the following equation:
Q c = Δ H f δ ( ρ χ i c e ) δ t
where χ i c e represents the ice fraction, which is a suspension of ice in the liquid water phase [3].
The mean ice crystal size L * (m) was estimated using the following equation:
L * = a R 0.5 G 0.5
with a as an empirical constant based on experimental data (ms/K), R as the freezing front rate (m/s), and G as the temperature gradient in the frozen zone (K/m).
Based on the above equations, freezing curves for aqueous solutions were calculated and compared to experimental temperature data, and a good agreement was found. The authors determined the ice-crystal growth rates and the temperature gradients in the aqueous solution, to estimate the ice crystal mean sizes and the resulting water vapour mass transfer permeability. Their results showed that water vapour mass transfer resistance decreases as the nucleation temperature increases and the cooling rate decreases. This was in accordance with the ice crystal size estimation obtained from image analysis. The author concluded that the freezing conditions have a direct impact on the permeability of the dried layer during the sublimation step [3]. This model described the data well and provided important insights into the freezing process. However, in order to apply their model to next-generation technologies such as spin freezing [1], several additions need to be made.
In this work, a mechanistic cooling and freezing model using elements of the model by Nakagawa et al. and applied to the spin freezing step of a continuous freeze-drying technology was developed. During continuous freeze-drying, a constant inflow of vials filled with an aqueous formulation are rapidly rotated along their longitudinal axis, while a flow of a cold, inert, and sterile gas is used for the cooling and freezing of the product. These spin-frozen vials are further processed in the consecutive continuous drying steps, making use of radiative heat [1,20]. This continuous freeze-drying method has several advantages compared to batch freeze-drying. Importantly, similar processing conditions are created for all vials, and the process can be monitored using in-line PAT tools and controlled at a single vial level [21]. The freezing step in particular can be controlled well, since the flow rate and temperature of the cooling gas can be set in order to control the rate of heat transfer during the entire freezing process, as will be shown in this work. The goal of this work was the development and application of a mechanistic model to predict and control the vial temperature during the spin freezing step using thermal imaging. The developed model was evaluated thoroughly by global sensitivity and uncertainty analyses, as well as experimentally.

2. Materials and Methods

2.1. Spin Freezing

All spin freezing experiments were conducted in a single-vial continuous freeze-dryer (RheaVita, Zwijnaarde, Belgium) (Figure 2). Ten mL type I glass vials (Schott, Müllheim, Germany) were filled with 3.0 mL of demineralized water. The vial was placed in a holder supported with grippers inside the stainless-steel processing chamber. The vial was rotated around its longitudinal axis at approximately 2900 rpm. Simultaneously, dehumidified air, cooled using a heat exchanger, was jetted onto the vial. The heat exchanger consisted of polyurethane tubing with an outer diameter of 8 mm and an inner diameter of 5 mm. Three m of this tubing was inserted into a Dewar container that was filled with liquid nitrogen. The cooled air passed through a stainless steel gas diffuser and travelled 1.5 cm before hitting the rotating vial. The cooling gas flow rate was measured and controlled using a Bronkhorst® F-203AV digital mass flow controller (Flowcor, Belgium). The temperature of the cooling gas was measured using a type K thermocouple positioned between the rotating vial and the outlet of the gas diffuser, at a distance of 1 mm from the rotating vial. The temperature of the vial was continuously followed up using a FLIR A655sc IR camera (Thermal Focus, Ravels, Belgium) as described in previous work [21]. The IR camera was installed outside the chamber at a distance of 20 cm from the vial, in front of an IR-transparent germanium window (Figure 3). The transmission of the germanium window was 0.86, and the emissivity of the measured borosilicate vial glass was 0.93. The vial temperature, the gas temperature, and the gas flow-rate data were simultaneously logged every 0.50 s and collected using a custom-made Labview® virtual instrument. The recording of data started as soon as the gas flow rate stabilized, usually within 4 s of starting the gas flow.
It should be noted that when the gas temperature changes, the thermocouple does not instantly report the correct new value, as the thermocouple has to equilibrate with the new gas temperature. This is of particular importance in this work, as the gas temperature has a substantial impact on the cooling process and changes rapidly at the start of the process. Indeed, at the start of spin freezing, the cooling gas is still relatively warm due to heating from the cooling system (e.g., the tubing the gas flows through, which is initially at room temperature). However, the system rapidly cools along with the gas temperature, necessitating a measurement lag correction. A lag error T e , l a g (K) was defined and calculated by the following equation [22]:
T e , l a g = N t d T g a s d t
with d T g a s d t as the rate of change of the gas temperature (K/s) and N t as the time constant (s).
The time constant N t is defined as the ratio of the thermal sensor heat capacity to the thermal conductance of the thermal sensor material. Majdak et al. developed an empirical least-squares regression model that allows calculation of N t based on the type of thermocouple used, which was applied in this research to the used thermocouple type and dimensions. The following equation was used [23]:
N t = 1 N t , a + N t , b w
Here, N t , a (1/s) and N t , b (m 1 2 s 1 2 ) are regression model parameters that depend on the type of thermocouple used, and w is the air flow velocity (m/s).
Using T e , l a g , T g a s was corrected at every timepoint by subtracting the lag error.

2.2. Mechanistic Cooling and Freezing Model

The mechanistic cooling and freezing model developed during this study predicts the temperature profile of the outer-vial wall for every discrete time step of the spin freezing process. From the outer-vial wall temperature, the product temperature and temperature at the inner vial wall can be calculated as described below. In our developed model, the spin freezing step was divided into four phases: liquid cooling, ice nucleation, crystal growth, and solid cooling (Figure 1). Liquid cooling refers to the phase where the vial contents are present in the liquid state and cooling of this system (i.e., the vial including its contents) takes place. The liquid cooling stage ends at the moment ice nucleation occurs and the first ice crystals are formed. The subsequent crystal growth phase refers to the transition from the remaining liquid to the solid state. Finally, solid cooling refers to further cooling of the completely solidified product, down to the final freezing temperature where the model simulation ends.
In the developed spin freezing model, the rate of heat transfer Q ˙ (W) is used to describe the amount of energy per time unit that is removed from the system. Q ˙ is calculated at every time point of the spin freezing process using Newton’s law of cooling [24]:
Q ˙ = π h Ø v i a l h v i a l ( T v , o T g a s )
where h is the heat transfer coefficient (W/(m 2 K)); Ø v i a l is the diameter of the vial (m); h v i a l is the height of the vial (m); T v , o is the temperature at the outside of the vial wall (K); and T g a s (K) is the temperature of the cooling gas.
The heat transfer coefficient h depends on several factors, including the rotation speed of the vial, the distance between the rotating vial and the gas nozzle, the vial type and dimensions, the longitudinal position of the vial, the properties of the freezing chamber, the gas diffuser properties, the cooling gas properties, the heat-exchanger properties, environmental factors (e.g., ambient temperature around the freezing chamber), and the gas flow rate [25]. However, most of these parameters remain constant during spin freezing. Only the gas flow rate is changed during spin freezing, as this parameter is used to control h, and therefore Q ˙ . Hence, a linear regression model (i.e., with the gas flow rate as the independent variable and the heat transfer coefficient as the dependent variable) that is applicable for the certain set of spin freezing conditions mentioned above was created. In this way, h can be calculated from the gas flow rate using the regression model, at any point in time.
In order to create the linear regression model, eight calibration experiments were performed at constant flow rates (i.e., 20, 27, 34, 41, 48, 55, 62, and 69 L/min), whereafter the heat transfer coefficient was calculated at every time point for every experiment by rearranging Equation (9) as follows:
h = Q ¯ π Ø v i a l h v i a l ( T v , o T g a s )
Here, the vial dimensions Ø v i a l and h v i a l were known constants, and T v , o and T g a s were measured at every time point. The mean rate of heat transfer Q ¯ was calculated using the length of the crystal growth phase t c r y s t according to the following equation:
Q ¯ = Q c r y s t t c r y s t
Here, Q c r y s t is the heat released due to the crystallization of the vial contents (J). The slight decrease in temperature of the solid ice layer during this phase was considered to be negligable with regards to energy transfer, and it was assumed that all heat released due to crystallization is removed at a rate Q ¯ . When the crystal growth is complete, the vial temperature profile displays a marked drop, since the exothermic process of ice crystallization stops. This allows determination of t c r y s t using vial temperature profile data.
Q c r y s t was calculated as follows:
Q c r y s t = m w a t e r Δ H f
where m w a t e r is the mass of water contained in the rotating vial (kg), and Δ H f is the latent heat released during ice crystallization (J/kg).
The linear regression model was then constructed using the least-squares method, where an average h was calculated per calibration experiment (i.e., one value for h per gas flow rate) from the heat transfer coefficients calculated at every timepoint (Equation (10)):
h = h a V ˙ + h b
Here, V ˙ is the volumetric gas flow rate (m 3 /s), h a is the regression coefficient (J/(m 5 K)), and h b (W/(m 2 K)) is the regression error term. Using this equation, h, and subsequently Q ˙ , were calculated for every timepoint of the process. The root mean square error of cross-validation (RMSECV) was calculated using the leave-one-out method [26].
As mentioned previously, the developed model predicts T v , o during every phase of spin freezing for every discrete time step. Firstly, during the liquid cooling phase, T v , o decreases every time step with a certain amount, according to the following equation:
T v , o ( i ) = T v , o ( i 1 ) Q ˙ ( i 1 ) C t o t , w d t
Here, T v , o ( i ) refers to the temperature at the outer wall of the vial at the current time step; T v , o ( i 1 ) refers to this temperature at the previous time step; and d t refers to the length of the time step (s). In this way, T v , o decreases every time step with a value of Q ˙ ( i 1 ) C t o t , w d t , where Q ˙ ( i 1 ) refers to the rate of heat transfer of the previous step, and C t o t , w is the total heat capacity of the system before crystal growth (J/K). C t o t , w depends on the specific heat capacity of the used vial glass c p v i a l (J/(kg K)), the weight of the vial m v i a l (kg), the specific heat capacity of water c p w a t e r (J/(kg K)), and the mass of water contained in the vial m w a t e r (kg):
C t o t , w = c p v i a l m v i a l + c p w a t e r m w a t e r
The calculation of the decreasing product temperature continues until the product temperature reaches the nucleation temperature T n u c . As mentioned above, ice nucleation is a stochastic phenomenon, which is difficult to accurately predict, even with controlled nucleation methods. In this model, nucleation is not predicted and, consequently, was an input for the model. However, it is possible to detect nucleation during spin freezing using in-line thermal imaging, as it is associated with a sharp increase in the vial temperature. At high rates of heat transfer (i.e., at high gas flow rates), it is possible that the developed temperature gradient within the liquid layer becomes sufficiently large so that not all contents of the vial are below T e q at the time of nucleation. Therefore, it is assumed that nucleation then only occurs in the zone of the product where the temperature is below T e q . This zone is defined by a volume V n u c l , calculated as follows:
V n u c l = π h ( r v , i 2 r s c 2 )
where r v , i is the inner vial radius (m), and r s c is the radius of ice nucleation (m) outside of which ice nucleation takes place, which is calculated as follows:
r s c = r v , i e 2 π k i c e h v i a l ( T e q T v , i ( i ) ) Q ˙ ( i )
where k i c e is the thermal conductivity of ice at 0 ° C (W/(m K)); T e q (K) is the equilibrium freezing temperature; and T v , i ( i ) refers to the temperature at the inner vial wall (K). Here, T v , i ( i ) must be calculated from T v , o ( i ) , by using Fourier’s law of heat conduction adapted to a hollow cylinder [27]:
T v , i = T v , o + Q ˙ ln ( r v , o r v , i ) 2 π k g l a s s h v i a l
with r v , o as the outer vial radius (m), and k g l a s s as the thermal conductivity of the vial glass (W/(m K)).
It was assumed that nucleation occurs homogeneously in the zone of the product between r s c and r v , i [2,28]. An ice nucleation fraction ( χ n u c ) was calculated at the starting point of crystallization:
χ n u c = C t o t , w ( T e q T n u c ) Δ H f V n u c l ρ w a t e r
where ρ w a t e r (kg/m 3 ) is the water density. From here on, the crystal growth phase of the model starts, where a certain amount of ice crystallizes every discrete time step. Contrary to shelf freezing, during spin freezing, energy is removed from the curved outer surface of the vial, and not from the planar bottom surface. This results in a freezing front that travels centripetally from the edge of the inner glass vial wall towards the centre of the rotating vial (Figure 4). This is in contrast to a circular planar freezing front travelling from the bottom of the vial towards the top, as seen in shelf freezing. The ice crystallizing results in an increase in thickness of the ice layer, as illustrated in Figure 4. Outside of r s c , ice grows where ice nuclei are already present, which has to be taken into account:
m i c e ( i ) = m i c e ( i 1 ) + Q ˙ ( i 1 ) ( 1 + χ n u c ) Δ H f d t
within r s c , no ice was present right after nucleation, and the total mass of ice as the freezing front propagates through this layer was calculated as follows:
m i c e ( i ) = m i c e ( i 1 ) + Q ˙ ( i 1 ) Δ H f d t
In both Equations (20) and (21), m i c e ( i ) refers to the mass of ice present in the current time step, and m i c e ( i 1 ) refers to the mass of ice in the previous time step. Just like in the liquid cooling phase, Q ˙ ( i 1 ) refers to the rate of heat transfer in the previous time step. During crystallization, the temperature of the liquid product is assumed to be the equilibrium temperature.
However, to calculate T v , o , the temperature gradient across the growing ice layer and the glass wall of the vial must be taken into account. The temperature gradient across the ice layer can again be calculated using Fourier’s law of heat conduction (Equation (18)), with some minor adjustments:
Q ˙ = 2 π k i c e h v i a l T e q T v , i l n ( r v , i r v , i t h i c e )
with k i c e as the thermal conductivity of ice (W/(m K)). In order to calculate T v , i , the ice layer thickness t h i c e is required, which was calculated according to the following equations:
V c r y s t = m i c e ρ i c e
t h i c e = r v , i ( π r v , i 2 h v i a l ) V c r y s t π h v i a l
with V c r y s t as the volume of the crystallized ice (m 3 ) and ρ i c e as the density of ice (kg/m 3 ).
Equation (22) can then be rewritten as:
T v , i = T e q Q ˙ ln ( r v , i r v , i t h i c e ) 2 π k i c e h v i a l
T v , o should then be calculated from T v , i as before using Equation (18). The ice layer thickness increases until all liquid water has been converted to ice. At this point, the solid cooling phase starts, which is the final phase of the model. This phase is analogous to the liquid cooling phase, but the specific heat capacity of ice c p i c e (J/(kg K)) instead of water should be taken into account when calculating the total heat capacity:
C t o t , i = c p v i a l m v i a l + c p i c e m i c e
with C t o t , i as the total heat capacity of the frozen system (J/K). In this way, Equation (14) becomes:
T v , o ( i ) = T v , o ( i 1 ) Q ( i 1 ) C t o t , i
The spin freezing process ends when the final vial temperature (i.e., the temperature where the model simulation ends) is reached during the solid cooling phase. In this work, the final temperature was set at −50 °C.

2.3. Uncertainty Analysis and Global Sensitivity Analysis

An uncertainty analysis (UA) and a global sensitivity analysis (GSA) were conducted in order to characterize the developed mechanistic model. The UA allows quantification of the impact of the uncertainty of the model input parameters on the uncertainty in the model outcome, while the GSA allows to apportion the uncertainty in the model output to different sources of uncertainty in the model parameters [29].
A total of seven parameters were considered to be uncertain and therefore were included in the UA and GSA. These parameters, their uncertainty level, and the reason for their inclusion are listed in Table 1. The uncertainty in the heat transfer coefficient h was specified as the root-mean-square error (RMSE) of the linear regression model. The uncertainty of the vial properties was their permissible variation, as obtained from the vial manufacturer. The error on the mass of the water was calculated from the density of water and the accuracy of the pipetted volume, which was obtained from the pipette manufacturer. The uncertainty in the gas flow rate was based on the accuracy of the mass flow controller measurement. The equilibrium temperature was measured using an infrared camera, meaning its uncertainty value was equal to the accuracy of the used infrared camera. Finally, the uncertainty in the gas temperature originates from the accuracy of the used thermocouple.
A sampling-based approach was used for the UA on the freezing model [17,30]. For this method, the model was run for different combinations of input process variables, based on their uncertainty range (i.e., Monte Carlo simulations). The input matrix containing 10,000 samples (i.e., 10,000 different combinations of input parameters) was constructed using the Sobol sampling technique [17,30]. For all of the 10,000 samples, the mechanistic model was applied to predict a vial temperature profile, resulting in 10,000 predicted temperature profiles. In this way, 10,000 temperature values were obtained at every time point (i.e., every 0.5 s) for calculation of the uncertainty limits. At every time point, these 10,000 temperature values were ordered from low to high. The lower temperature limit of the prediction interval was defined as the temperature value corresponding to sample number 250 (i.e., 2.5%), while the upper limit was defined as the value corresponding to sample number 9750 (i.e., 97.5%). This approach was used for every time point of the process, resulting in a lower and an upper uncertainty limit, together forming a 95% prediction interval. Experimental temperature data were subsequently compared to the uncertainty limits.
For the GSA, a variance-based sensitivity analysis was conducted as it involves no prior assumptions about model output and allows quantification of overall interaction effects between factors [16,19]. A total order effect S T i was calculated from the decomposition of the total output variance into the contribution of the input factors [16,19]. S T i represents the total effect of factor i, which is the sum of the first-order effects, higher-order effects, and all interactions with other factors. The design proposed by Saltelli et al. was applied to calculate S T i in this study [29]. The mathematical details regarding the computation of S T i are described by Mortier et al. [16]. The GSA was conducted at different timings, each matching with a specific phase of the spin freezing process (i.e., the liquid cooling phase, the ice-crystal growth phase, and the solid cooling phase), in order to evaluate the difference in impact of the uncertain model input parameters at each process phase. The GSA was also performed for multiple gas flow rates: at a low gas flow rate (20 L/min), an intermediary gas flow rate (50 L/min), and a high gas flow rate (80 L/min).

2.4. Experimental Model Verification of Constant Gas Flow Rate and Imposed Cooling Profile Experiments

The model was verified using the data of 10 spin freezing experiments where the following constant flow rates were used: 20, 26, 32, 38, 44, 50, 56, 62, 68, and 80 L/min. An uncertainty analysis was performed as described above at every flow rate using data (e.g., gas and vial temperature) from the respective experiment.
Additionally, in order to demonstrate its applicability, the model was used to impose a targetted vial temperature profile. This profile was set to have a liquid cooling phase with a cooling rate of 20 ° C/min, a crystal growth phase with a duration of 150 s, and a solid cooling phase with a cooling rate of 20 ° C/min. The required rate of heat transfer Q ˙ to achieve this during liquid and solid cooling was calculated from the cooling rate:
Q ˙ = C r C t o t
where C r is the cooling rate (K/s), and C t o t is C t o t , w for the liquid cooling phase and C t o t , i for the solid cooling phase. The required rate of heat transfer Q ˙ during the crystallization phase was calculated using Equations (11) and (12). From Q ˙ , h was calculated for every timepoint using Equation (9). Finally, the required cooling gas flow rate was calculated from Equation (13) at every time point, which was then an input for the mass flow controller. The same imposed cooling profile was applied to perform 20 freezing cycles in order to assess the robustness of this model application. An uncertainty analysis was performed on the model vial temperature predictions for one representative replicate of the 20 imposed cooling temperature profile runs (i.e., with a cooling rate of 20 ° C/min during the liquid cooling phase, a crystal growth phase with a duration of 150 s, and a solid cooling phase with a cooling rate of 20 ° C/min) and compared to the recorded vial temperature data.
All calculations were performed using the numeric computing platform MATLAB version R2018b (The MathWorks, Natick, MA, USA).

3. Results and Discussion

3.1. Experimentally Obtained Vial Temperature Profiles and Subsequent Uncertainty Analysis

A linear regression model for determining the linear relationship between the gas flow rate and the heat transfer coefficient was constructed as described above, and regression coefficients were determined (Figure 5). h a was calculated to be 71.11 E3 J/(m 5 K), and h b was 32.05 W/(m 2 K). The RMSE was 4.3058 W/(m 2 K) and the RMSECV 5.6293 W/(m 2 K).
It must be noted that the product temperature was not directly predicted by the model developed in this work. The first reason is that model verification is significantly easier when considering the outer vial temperature, since this temperature can be directly verified using infrared camera measurements. Secondly, the product temperature is spatially different across the product layer, depending on the different phases of spin freezing. For example, during the crystal growth phase, the local temperature in the liquid layer will be the equilibrium temperature (i.e., 0 ° C for liquid water). At the same time, the growing layer of solid ice is already cooling down below the equilibrium temperature as there is no longer any exothermic crystallization. Additionally, there is a temperature gradient across the product layer due to the cooling flux during the process. In essence, the product temperature cannot be defined or predicted as a single value, which is why the outer vial temperature was calculated in this model. It should be noted that the outer vial temperature is sufficient to obtain good process control, as it is representative for the average product temperature, and the product temperature can be calculated precisely at any location in the product layer using Equations (18) and (22).
The temperatures of the vial wall measured using the infrared camera and the calculated model prediction, including the model prediction uncertainty interval, are shown for six representative spin freezing runs in Figure 6. During the liquid cooling phase, T v , o decreased until nucleation occurs. It can be seen that the cooling rate during this first phase depends on the gas flow rate (Table 2). As the gas flow rate V ˙ increases, h also increases (Equation (13)). This results in a larger rate of heat transfer Q ˙ (Equation (9)) and subsequently, a higher cooling rate (Equation (28)).
Nucleation occurred at a measured (i.e., by the infrared camera) outer vial wall temperature between −15 °C and −9 °C. By taking the temperature gradient across the glass vial wall into account, this results in a product temperature at the inner vial wall at the point of nucleation between −2 °C and 0 °C (Table 2). The range of nucleation temperatures measured at the outer vial wall was larger than the actual range of nucleation temperatures. The reason for this is that the lower nucleation temperatures were generally the result of experiments at high gas flow rates. Here, the rate of heat transfer across the glass vial wall was high, resulting in a large temperature gradient (Equation (18)). This seemingly results in very low measured nucleation temperatures, while the temperature at the inside of the vial wall at the point of nucleation was still high and similar for all performed experiments.
It should be noted that nucleation is a stochastic event and was not controlled in this study. Interestingly, although uncontrolled, nucleation took place at relatively high temperatures (i.e., ≥−2 ° C for all studied experimental conditions) compared to what is commonly seen in the literature, and the range wherein nucleation occurs appeared small (i.e., 2 ° C) [4,18,31]. Taking into account that the rotating vials were not closed with stoppers, it is possible that a certain amount of “ice fog” was created, which caused early nucleation of the product. Controlled nucleation via ice fog is a known technique whereby cold nitrogen gas is introduced into the freezing environment, creating a fine mist of ice crystals from the moisture in the freezing chamber. These crystals then come into contact with the solution, inducing nucleation if the product temperature is sufficiently low [4,5,32]. In the case of spin freezing, since the gas used to cool the vials is very cold (i.e., −40 to −90 °C), it is possible that a mist of ice crystals is generated as mixing occurs with air containing moisture present in the freezing chamber, which then results in rapid nucleation when the product temperature drops below the equilibrium freezing temperature. Additionally, the mechanical agitation due to the spinning of the vial may be a reason for early nucleation. Indeed, Konstantinidis et al. described mechanical induction of nucleation, suggesting that vibrational disturbances may play a role in this case [33]. Low values and small ranges of supercooling as those seen here may be beneficial with regards to ice crystal size (i.e., large crystals are formed) and vial-to-vial variability, respectively [4]. However, this apparent controlled nucleation phenomenon during spin freezing and its effect on subsequent drying steps must be investigated in more detail. To this end, nucleation behaviour and potential controlled nucleation approaches during spin freezing are the topic of a future investigation.
Once nucleation had taken place, the product temperature increased until T e q for all studied conditions (Figure 6). The duration of the subsequent crystal growth phase depends on the rate of heat transfer. The experiments at higher gas flow rates had shorter crystal growth phases, as expected (Table 2). Again, as the gas flow rate V ˙ increases, Q ˙ also increases (Equation (9)), which subsequently results in faster crystallization of ice (Equations (20) and (21)). During the crystal growth phase, T v , o decreases slightly. As the crystal growth phase progresses, the ice layer forming on the inner vial wall grows. This results in an increase in the temperature gradient between the liquid product (i.e., at equilibrium temperature) and the inner vial wall (Equation (22)). Consequently, this explains the decreasing temperature over time at the outer vial wall (Equation (18)).
After completion of the crystal growth phase, indicated by a sharp temperature drop (i.e., measured as well as predicted by the model), the vial wall cooled down until the final temperature (i.e., −50 ° C), at which point the process was considered complete. It should be noted that, in practice, there is a lower limit of this final temperature, determined by T g a s . In turn, T g a s depends on the characteristics of the heat exchanger (e.g., cooling medium temperature) but also on the gas flow rate. As the gas flow rate decreases, T g a s increases, since the slower-moving cold gas has more time to heat up when travelling from the heat exchanger to the gas nozzle (i.e., a distance of approximately 50 cm) and finally the vial. As the cooling process is driven by the difference between T v , o and T g a s , the vial wall can never cool below T g a s , as Q ˙ reaches a value of zero when T v , o equals T g a s . Similarly to the liquid cooling phase, the solid-phase cooling rate depends on the gas flow rate.
As seen in Figure 6, the model generally described the data well, as the experimental data obtained through infrared measurements consistently lay within the calculated prediction interval limits of the model predictions, for all tested gas flow rates.

3.2. Imposed Cooling Profile

The uncertainty analysis of one representative-imposed cooling profile can be found in Figure 7. Similar to the experiments with constant gas flow rates, the experimental data (i.e., obtained using the infrared camera) for the imposed cooling experiments also fell within the uncertainty limits of the model prediction. The measured cooling rates during the liquid cooling and solid cooling phases were 22.19 ± 2.14 ° C/min and 23.05 ± 2.64 ° C/min, respectively (mean ± SD). The measured crystal growth duration was 151 ± 7.8 s (mean ± SD). Both cooling rates were slightly higher than the set value. As seen further (Section 3.3), the uncertainty in the heat transfer coefficient was very influential with regards to the overall model error. Hence, it is likely that the true values for the regression model parameters were higher than those calculated during the calibration. As a result, the calculated value of the gas flow rate (Equation (13)) was larger than the real required value necessary to obtain a cooling rate of 20 ° C/min. Using a more elaborate calibration procedure (e.g., performing more calibration experiments resulting in more calibration points), it may be possible to reduce this systematic error. Although it does appear that the model slightly overestimated T v , o at higher gas flow rates during the liquid cooling phase, the data were still within the model prediction intervals. It should be noted that alternative control mechanisms such as closed control loops (e.g., proportional integral derivative control) can be used to control the cooling phases with significantly less uncertainty than the open-loop-model-based approach presented here. Using these control loops, measured temperature data would be fed into the control loop, which then outputs the required process settings (e.g., the gas temperature and the gas flow rate).
Contrary to batch freeze-drying, in spin freeze-drying it is possible to tightly control the rate of heat transfer during the freezing phase using the model developed in this work. Indeed, during batch freeze-drying, the rate of heat transfer during freezing is governed by the temperature difference between the vials containing the product and the temperature-controlled shelves. The temperature-controlled shelves are slow to react to changes in the set point temperature, while a change in the gas flow rate during spin freezing is near instantaneous. This allows for rapid changes in the rate of heat transfer at any point (e.g., immediately after nucleation) during the spin freezing process. As mentioned previously, since nucleation is stochastic, vials nucleate at different shelf temperatures during batch freeze-drying, resulting in different rates of heat transfer between vials during ice crystallization. In spin freezing, it is possible to compensate for the stochastic behaviour of nucleation by changing the rate of heat transfer after nucleation to a desired value, which allows for the manipulation of the freezing rate to obtain desired product characteristics such as dried layer resistance [3]. Spin freezing also has a broad range of possible freezing rates, limited only by the temperature and the flow rate of the used cooling gas. On the contrary, batch freezing is limited by the cooling speed of the temperature-controlled shelves (i.e., usually no faster than 1 ° C/min). Many products that are sensitive to the freezing rate with regards to product integrity may benefit from this increased freezing rate range and level of control. In a recent work regarding the distribution of COVID-19 vaccines, the low freezing rate associated with batch freeze-drying was mentioned as a limitation that potentially causes phase separation, impairing product quality [34]. There have been similar findings in studies of other vaccines but also in studies of monoclonal antibodies and other protein-based products [35,36,37,38]. Additionally, the freezing rate is considered very important when freezing cell-based products such as microorganisms or mammalian-cell-based therapeutics, which are classes of products that are expected to become increasingly important in medicine [39,40,41,42,43].

3.3. Global Sensitivity Analysis

In order to explain the results of the uncertainty analysis discussed above, a global sensitivity analysis was performed. In this way, the contribution of the uncertainty of the model input parameters to the overall model output uncertainty can be determined. As seen in Figure 8, the uncertainty in h was consistently the most influential on the model output uncertainty during the liquid cooling phase. After h, T g a s follows as the second most important source of uncertainty of the model during the liquid cooling phase. h and T g a s were responsible for 98.3%, 98.4%, and 98.0% of the prediction uncertainty for gas flow rates of 20, 50, and 80 L/min, respectively. The four remaining parameters (i.e., Ø v i a l , m v i a l , m w a t e r , and V ˙ ) were considered as having a negligible impact on the uncertainty during the liquid cooling phase. In essence, the error on these four parameters was too small to have any noticeable impact on overall model uncertainty.
During the crystal growth phase, the uncertainty in the measurement of the equilibrium temperature had the most important impact for all gas flow rates. The equilibrium temperature is the starting point for T v , o during the crystal growth phase, resulting in a large impact on the uncertainty during this phase. However, T g a s and h still have a relevant impact, especially at higher gas flow rates. Equation (25) shows how the thickness of the ice (i.e., the temperature gradient across this ice layer) and the rate of heat transfer at that time point must be taken into account when calculating T v , o . Both of the aforementioned parameters influence the temperature gradient across the ice, explaining why h and T g a s have an influence on uncertainty during the crystallization phase. However, the relative importance of the equilibrium temperature did seem lower with an increasing gas flow rate. A rising gas flow rate also resulted in an increase in the ( T v , o T g a s ) term of Equation (9), as higher gas flow rates are associated with lower values of T g a s (see above). ( T v , o T g a s ) was subsequently multiplied by h, including its uncertainty. Therefore, higher gas flow rates result in a higher impact of the uncertainty in h on the overall model uncertainty relative to the contribution of T e q . Similar to the liquid cooling phase, at least 96% of all variability was explained by the uncertainty in h, T e q , and T g a s for every gas flow rate.
Finally, during the solid cooling phase, the observations of the liquid cooling phase were still valid. Again, close to all model output uncertainty (i.e., >96%) was explained by the uncertainty in the regression model parameters, T e q and T g a s .
In essence, the GSA showed that the uncertainty in h was by far the most impactful on overall model uncertainty, followed by the uncertainty in T g a s , and, during the crystallization phase, the uncertainty in T e q . The uncertainty in the remaining four parameters (i.e., Ø v i a l , m v i a l , m w a t e r , and V ˙ ) had a negligible impact on the overall model uncertainty. Therefore, to reduce the model uncertainty, the focus should lie in the uncertainty of h. As mentioned above, using a more elaborate calibration procedure may help in this way to alleviate some of the model uncertainty.

4. Conclusions

Using in-line thermal imaging, it is possible to predict and control the product and vial temperature profiles during spin freezing using the mechanistic model developed in this work. The model was experimentally verified using an uncertainty analysis and was further characterized using a global sensitivity analysis. The uncertainty analysis was used to establish a prediction interval that contained the experimental data in every experimental verification run. The global sensitivity analysis was performed to determine the contribution of each model input parameter’s uncertainty to the uncertainty of the model output (i.e., the predicted temperature). The uncertainty in the heat transfer coefficient, the uncertainty in the gas temperature T g a s , and the uncertainty in the equilibrium temperature T e q were the most important parameters influencing the overall temperature prediction uncertainty. While the relative contribution of these three parameters differed depending on the phase of freezing and the gas flow rate, together they always explained at least 96% of the uncertainty in the model prediction. The remaining four parameters (i.e., Ø v i a l , m v i a l , m w a t e r , and V ˙ ) were therefore considered to be of lesser importance regarding the model prediction uncertainty.
It was found that the calibration of the heat transfer coefficient has a large influence on the model uncertainty. Additionally, the applicability of the model was limited to the conditions where calibration was performed (e.g., the distance between the vial and the gas nozzle and the rotation speed of the vial). Therefore, to deal with these disadvantages, alternatives such as computational fluid dynamics or mechanistic calculation of the heat transfer coefficient should be investigated. To increase the applicability of the calibration approach, it will be investigated how different parameters of the spin freezing setup (e.g., the vial rotation speed) influence the regression model parameters. When the influence of the main parameters can be calculated, calibration data may be transferable to different spin freezing conditions, provided that the necessary adjustments are calculated.
In summary, a mechanistic model predicting the vial temperature during the spin freezing phase of a continuous freeze-drying process was developed and evaluated thoroughly. The results of the verification and subsequent analyses indicate a good process understanding. This model can be applied to control the spin freezing phase during continuous freeze-drying to overcome some of the disadvantages related to batch freeze-drying.

Author Contributions

Conceptualization: S.R.R., A.K., B.V., P.-J.V.B., T.D.B. and J.C.; methodology: G.N., A.K. and P.-J.V.B.; formal analysis: G.N.; software: G.N., P.-J.V.B. and B.V.; investigation:. G.N., S.R.R., A.K., P.-J.V.B., J.L., L.L., B.V., J.C. and T.D.B.; resources: J.C., C.V. and T.D.B.; writing—original draft preparation: G.N. and P.-J.V.B.; writing—review and editing: G.N. and T.D.B. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge financial support from the International Society of Lyophilization—Freeze Drying (ISLFD).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available upon request from the corresponding authors.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

T g Glass transition temperature of the maximum freeze-concentrated drug solution (K)
T e Eutectic temperature (K)
ρ Density of the solution (kg/m3)
c p Specific heat capacity of the solution (J/(kg K)) (Nakagawa et al.)
kThermal conductivity (W/(m K)) (Nakagawa et al.)
TTemperature distribution of the system (K) (Nakagawa et al.)
Q n Total latent heat released due to ice nucleation (W) (Nakagawa et al.)
Q c Total latent heat released due to ice crystallisation (W) (Nakagawa et al.)
Δ H f Latent heat of crystallization (J/kg) (Nakagawa et al.)
k i Nucleation rate (kg/(m3 s K)) (Nakagawa et al.)
T s Temperature in the supercooled liquid (K) (Nakagawa et al.)
ν Freezing front velocity (m/s) (Nakagawa et al.)
sThickness of the undercooled zone (m) (Nakagawa et al.)
T Homogeneous undercooled temperature (K) (Nakagawa et al.)
χ i c e Fraction of ice suspended in the undercooled solution (Nakagawa et al.)
L * Mean ice crystal size (m)
aan empirical constant based on experimental data (ms/K) in Equation (6)
RFreezing front rate (m/s)
GTemperature gradient in the frozen zone (K/m)
T e , l a g Thermocouple lag error (K)
N t Time constant (s)
N t , a Time constant regression model coefficient (1/s)
N t , b Time constant regression model coefficient (m 1 2 s 1 2 )
wVelocity of the air flow (m/s)
Q c r y s t Total heat released due to crystallization of ice (J)
Q ¯ Mean rate of heat transfer (W)
m w a t e r Mass of water in the vial (kg)
hHeat transfer coefficient (W/(m 2 K))
Ø v i a l Diameter of the glass vial (m)
h v i a l Height of the glass vial (m)
T v , o Temperature at the outside of the glass vial wall (K)
T g a s Temperature of the gas (K)
Q ˙ Rate of heat transfer (W)
t c r y s t Duration of the crystal growth phase (s)
V ˙ Volumetric gas flow (m 3 /s)
h a Regression coeffient (J/(m 5 K))
h b Regression error term (W/(m 2 K))
C t o t , w Total heat capacity of the system before freezing (J/K)
d t Length of the modelling time step (s)
c p w a t e r Specific heat capacity of water (J/(kg K))
c p v i a l Specific heat capacity of the vial glass (J/(kg K))
m v i a l Mass of the glass vial (kg)
r s c Radius of the nucleation ice crystal formation (m)
r v , o The outer vial radius (m)
k g l a s s Thermal conductivity of the vial glass (W/(m K))
T v , i Temperature at the inside of the glass vial wall (K)
t h i c e Thickness of the ice layer (m)
r v , i Inner vial radius (m)
k i c e Thermal conductivity of ice (W/(m K))
T e q Equilibrium freezing temperature (K)
χ n u c Fraction of ice crystallized after nucleation
C t o t n u c l Total heat capacity of the volume below the equilibrium freezing temperature (J/K)
T n u c l Nucleation temperature (K)
V n u c l Volume below the equilibrium freezing temperature (m 3 )
ρ w a t e r Density of water (kg/m 3 )
ρ i c e Density of ice (kg/m 3 )
C t o t , i Total heat capacity of the system after freezing (J/K)
c p i c e Specific heat capacity of ice (J/(kg K))
m i c e Mass of ice (kg)
GSAGlobal sensitivity analysis
UAUncertainty analysis
S T i Total order effect
C r Cooling rate (K/s)

References

  1. De Meyer, L.; Van Bockstal, P.J.; Corver, J.; Vervaet, C.; Remon, J.P.; De Beer, T. Evaluation of spin freezing versus conventional freezing as part of a continuous pharmaceutical freeze-drying concept for unit doses. Int. J. Pharm. 2015, 496, 75–85. [Google Scholar] [CrossRef] [Green Version]
  2. Kasper, J.C.; Friess, W. The freezing step in lyophilization: Physico-chemical fundamentals, freezing methods and consequences on process performance and quality attributes of biopharmaceuticals. Eur. J. Pharm. Biopharm. 2011, 78, 248–263. [Google Scholar] [CrossRef]
  3. Nakagawa, K.; Hottot, A.; Vessot, S.; Andrieu, J. Modeling of freezing step during freeze-drying of drugs in vials. AIChE J. 2007, 53, 1362–1372. [Google Scholar] [CrossRef]
  4. Pisano, R. Lyophilization of Pharmaceuticals and Biologicals; Methods in Pharmacology and Toxicology; Springer: New York, NY, USA, 2019; pp. 79–111. [Google Scholar] [CrossRef]
  5. Patel, S.M.; Bhugra, C.; Pikal, M.J. Reduced pressure ice fog technique for controlled ice nucleation during freeze-drying. AAPS PharmSciTech 2009, 10, 1406–1411. [Google Scholar] [CrossRef] [PubMed]
  6. Petersen, A.; Schneider, H.; Rau, G.; Glasmacher, B. A new approach for freezing of aqueous solutions under active control of the nucleation temperature. Cryobiology 2006, 53, 248–257. [Google Scholar] [CrossRef] [PubMed]
  7. Inada, T.; Zhang, X.; Yabe, A.; Kozawa, Y. Active control of phase change from supercooled water to ice by ultrasonic vibration 1. Control of freezing temperature. Int. J. Heat Mass Transf. 2001, 44, 4523–4531. [Google Scholar] [CrossRef]
  8. Nakagawa, K.; Hottot, A.; Vessot, S.; Andrieu, J. Influence of controlled nucleation by ultrasounds on ice morphology of frozen formulations for pharmaceutical proteins freeze-drying. Chem. Eng. Process. Process Intensif. 2006, 45, 783–791. [Google Scholar] [CrossRef]
  9. Hottot, A.; Nakagawa, K.; Andrieu, J. Effect of ultrasound-controlled nucleation on structural and morphological properties of freeze-dried mannitol solutions. Chem. Eng. Res. Des. 2008, 86, 193–200. [Google Scholar] [CrossRef]
  10. Passot, S.; Tréléa, I.C.; Marin, M.; Galan, M.; Morris, G.J.; Fonseca, F. Effect of controlled ice nucleation on primary drying stage and protein recovery in vials cooled in a modified freeze-dryer. J. Biomech. Eng. 2009, 131, 74511. [Google Scholar] [CrossRef] [PubMed]
  11. Kramer, M.; Sennhenn, B.; Lee, G. Freeze-drying using vacuum-induced surface freezing. J. Pharm. Sci. 2002, 91, 433–443. [Google Scholar] [CrossRef] [PubMed]
  12. Cochet, N.; Widehem, P. Ice crystallization by Pseudomonas syringae. Appl. Microbiol. Biotechnol. 2000, 54, 153–161. [Google Scholar] [CrossRef]
  13. Assegehegn, G.; Brito-de la Fuente, E.; Franco, J.M.; Gallegos, C. The Importance of Understanding the Freezing Step and Its Impact on Freeze-Drying Process Performance. J. Pharm. Sci. 2019, 108, 1378–1395. [Google Scholar] [CrossRef] [PubMed]
  14. Searles, J.A.; Carpenter, J.F.; Randolph, T.W. The ice nucleation temperature determines the primary drying rate of lyophilization for samples frozen on a temperature-controlled shelf. J. Pharm. Sci. 2001, 90, 860–871. [Google Scholar] [CrossRef]
  15. Nakagawa, K.; Hottot, A.; Vessot, S.; Andrieu, J. Modeling of freezing step during vial freeze-drying of pharmaceuticals—Influence of nucleation temperature on primary drying rate. Asia-Pac. J. Chem. Eng. 2011, 6, 288–293. [Google Scholar] [CrossRef]
  16. Mortier, S.T.F.C.; Gernaey, K.V.; De Beer, T.; Nopens, I. Global sensitivity analysis applied to drying models for one or a population of granules. AIChE J. 2014, 60, 1700–1717. [Google Scholar] [CrossRef]
  17. Mortier, S.T.F.C.; Van Bockstal, P.J.; Corver, J.; Nopens, I.; Gernaey, K.V.; De Beer, T. Uncertainty analysis as essential step in the establishment of the dynamic Design Space of primary drying during freeze-drying. Eur. J. Pharm. Biopharm. 2016, 103, 71–83. [Google Scholar] [CrossRef] [PubMed]
  18. Pisano, R.; Capozzi, L.C. Prediction of product morphology of lyophilized drugs in the case of Vacuum Induced Surface Freezing. Chem. Eng. Res. Des. 2017, 125, 119–129. [Google Scholar] [CrossRef]
  19. Van Bockstal, P.J.; Mortier, S.T.F.C.; Corver, J.; Nopens, I.; Gernaey, K.V.; De Beer, T. Global Sensitivity Analysis as Good Modelling Practices tool for the identification of the most influential process parameters of the primary drying step during freeze-drying. Eur. J. Pharm. Biopharm. 2018, 123, 108–116. [Google Scholar] [CrossRef] [Green Version]
  20. Van Bockstal, P.J.; Corver, J.; Mortier, S.T.F.C.; De Meyer, L.; Nopens, I.; Gernaey, K.V.; De Beer, T. Developing a framework to model the primary drying step of a continuous freeze-drying process based on infrared radiation. Eur. J. Pharm. Biopharm. 2018, 127, 159–170. [Google Scholar] [CrossRef] [PubMed]
  21. Van Bockstal, P.J.; Corver, J.; De Meyer, L.; Vervaet, C.; De Beer, T. Thermal Imaging as a Noncontact Inline Process Analytical Tool for Product Temperature Monitoring during Continuous Freeze-Drying of Unit Doses. Anal. Chem. 2018, 90, 13591–13599. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Nicholas, J.V.; White, D.R. Use of Thermometers. Traceable Temperatures: An Introduction to Temperature Measurement and Calibration, 2nd ed.; Wiley: Hoboken, NJ, USA, 2001; pp. 125–158. [Google Scholar] [CrossRef]
  23. Majdak, M.; Jaremkiewicz, M. The analysis of thermocouple time constants as a function of fluid velocity. Meas. Autom. Monit. 2016, 62, 284–287. [Google Scholar]
  24. Rohsenow, W.M.; Hartnett, J.P.; Cho, Y. Handbook of Heat Transfer, 3rd ed.; McGraw Hill: New York, NY, USA, 1997; pp. 1–36. [Google Scholar]
  25. Jeng, T.M.; Tzeng, S.C.; Xu, R. Heat transfer characteristics of a rotating cylinder with a lateral air impinging jet. Int. J. Heat Mass Transf. 2014, 70, 235–249. [Google Scholar] [CrossRef]
  26. Sammut, C.; Webb, G.I. (Eds.) Cross-Validation. In Encyclopedia of Machine Learning; Springer: Boston, MA, USA, 2010; pp. 600–601. [Google Scholar] [CrossRef]
  27. Ortuño, M.; Márquez, A.; Gallego, S.; Neipp, C.; Beléndez, A. An experiment in heat conduction using hollow cylinders. Eur. J. Phys. 2011, 32, 1065–1075. [Google Scholar] [CrossRef]
  28. Searles, J.A. Freezing and annealing phenomena in lyophilization. In Freeze-Drying/Lyophilization of Pharmaceutical and Biological Products; CRC Press: Boca Raton, FL, USA, 2004; pp. 109–145. [Google Scholar]
  29. Saltelli, A. The Critique of Modelling and Sensitivity Analysis in the Scientific Discourse; European Commission, Joint Research Centre: Luxembourg, 2006; ISBN 92-79-03130-9. [Google Scholar]
  30. Van Bockstal, P.J.; Mortier, S.T.F.C.; Corver, J.; Nopens, I.; Gernaey, K.V.; De Beer, T. Quantitative risk assessment via uncertainty analysis in combination with error propagation for the determination of the dynamic Design Space of the primary drying step during freeze-drying. Eur. J. Pharm. Biopharm. 2017, 121, 32–41. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Colucci, D.; Fissore, D.; Barresi, A.A.; Braatz, R.D. A new mathematical model for monitoring the temporal evolution of the ice crystal size distribution during freezing in pharmaceutical solutions. Eur. J. Pharm. Biopharm. 2020, 148, 148–159. [Google Scholar] [CrossRef] [PubMed]
  32. Rowe, T. A technique for the nucleation of ice. In Proceedings of the International Symposium on Biological Product Freeze-Drying and Formulation, Bethesda, MD, USA, 24–26 October 1990. [Google Scholar]
  33. Konstantinidis, A.K.; Kuu, W.; Otten, L.; Nail, S.L.; Sever, R.R. Controlled nucleation in freeze-drying: Effects on pore size in the dried product layer, mass transfer resistance, and primary drying rate. J. Pharm. Sci. 2011, 100, 3453–3470. [Google Scholar] [CrossRef]
  34. AboulFotouh, K.; Cui, Z.; Williams, R.O., 3rd. Next-Generation COVID-19 Vaccines Should Take Efficiency of Distribution into Consideration. AAPS Pharm. Sci. Tech. 2021, 22, 126. [Google Scholar] [CrossRef]
  35. Miller, M.A.; Rodrigues, M.A.; Glass, M.A.; Singh, S.K.; Johnston, K.P.; Maynard, J.A. Frozen-state storage stability of a monoclonal antibody: Aggregation is impacted by freezing rate and solute distribution. J. Pharm. Sci. 2013, 102, 1194–1208. [Google Scholar] [CrossRef]
  36. Hansen, L.J.J.; Daoussi, R.; Vervaet, C.; Remon, J.P.; De Beer, T.R.M. Freeze-drying of live virus vaccines: A review. Vaccine 2015, 33, 5507–5519. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Hauptmann, A.; Podgoršek, K.; Kuzman, D.; Srčič, S.; Hoelzl, G.; Loerting, T. Impact of Buffer, Protein Concentration and Sucrose Addition on the Aggregation and Particle Formation during Freezing and Thawing. Pharm. Res. 2018, 35, 101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Roy, I.; Gupta, M.N. Freeze-drying of proteins: Some emerging concerns. Biotechnol. Appl. Biochem. 2004, 39, 165–177. [Google Scholar] [CrossRef] [PubMed]
  39. Aly, R.M. Current state of stem cell-based therapies: An overview. Stem Cell Investig. 2020, 7, 8. [Google Scholar] [CrossRef] [PubMed]
  40. Baboo, J.; Kilbride, P.; Delahaye, M.; Milne, S.; Fonseca, F.; Blanco, M.; Meneghel, J.; Nancekievill, A.; Gaddum, N.; Morris, G.J. The Impact of Varying Cooling and Thawing Rates on the Quality of Cryopreserved Human Peripheral Blood T Cells. Sci. Rep. 2019, 9, 3417. [Google Scholar] [CrossRef] [PubMed]
  41. Natan, D.; Nagler, A.; Arav, A. Freeze-drying of mononuclear cells derived from umbilical cord blood followed by colony formation. PLoS ONE 2009, 4, e5240. [Google Scholar] [CrossRef] [PubMed]
  42. Meneghel, J.; Kilbride, P.; Morris, G.J. Cryopreservation as a Key Element in the Successful Delivery of Cell-Based Therapies—A Review. Front. Med. 2020, 7, 592242. [Google Scholar] [CrossRef]
  43. Rockinger, U.; Funk, M.; Winter, G. Current Approaches of Preservation of Cells During (freeze-) Drying. J. Pharm. Sci. 2021, 110, 2873–2893. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Diagram of vial temperature over time during the separate phases of cooling and freezing. The process starts with the cooling of the vial containing liquid product, until the nucleation temperature is reached. At this point, the temperature of the product within the vial increases up to the equilibrium temperature, and the ice crystallization phase starts. This phase ends when all water has crystallized into ice, whereafter the vial containing the now-solid product cools down until the final temperature. A cross-section of the vial during spin freezing is shown.
Figure 1. Diagram of vial temperature over time during the separate phases of cooling and freezing. The process starts with the cooling of the vial containing liquid product, until the nucleation temperature is reached. At this point, the temperature of the product within the vial increases up to the equilibrium temperature, and the ice crystallization phase starts. This phase ends when all water has crystallized into ice, whereafter the vial containing the now-solid product cools down until the final temperature. A cross-section of the vial during spin freezing is shown.
Pharmaceutics 13 02076 g001
Figure 2. Single-vial continuous freeze-drying system.
Figure 2. Single-vial continuous freeze-drying system.
Pharmaceutics 13 02076 g002
Figure 3. Schematic representation of the cooling system during spin freezing. (A): Scheme showing how pressurized air travels from the gas outlet through a pressure regulator, after which the gas flow rate is subsequently determined by a mass flow controller. The pressurized air then passes through the liquid nitrogen heat exchanger before flowing through the gas diffuser and cooling the rotating vial within the single-vial processing chamber. (B): Top view of the experimental setup, showing the rotating vial, the outlet of cold gas from the gas nozzle, the gas temperature measurement thermocouple, and the IR camera located outside the chamber measuring through a germanium window.
Figure 3. Schematic representation of the cooling system during spin freezing. (A): Scheme showing how pressurized air travels from the gas outlet through a pressure regulator, after which the gas flow rate is subsequently determined by a mass flow controller. The pressurized air then passes through the liquid nitrogen heat exchanger before flowing through the gas diffuser and cooling the rotating vial within the single-vial processing chamber. (B): Top view of the experimental setup, showing the rotating vial, the outlet of cold gas from the gas nozzle, the gas temperature measurement thermocouple, and the IR camera located outside the chamber measuring through a germanium window.
Pharmaceutics 13 02076 g003
Figure 4. Axial view of spin freezing with ice layer thickness progression from left to right. The ice layer grows over time as more water crystallizes into ice, until all water has solidified and the crystal growth phase is complete. The inner vial radius r v , i and the radius of the nucleated zone r s c are shown.
Figure 4. Axial view of spin freezing with ice layer thickness progression from left to right. The ice layer grows over time as more water crystallizes into ice, until all water has solidified and the crystal growth phase is complete. The inner vial radius r v , i and the radius of the nucleated zone r s c are shown.
Pharmaceutics 13 02076 g004
Figure 5. Linear regression model of volumetric gas flow rate in function of the heat transfer coefficient. The blue dots indicate experimental data, and the red solid line indicates the linear regression model. The linear regression had a slope of 71.11 E3 J/(m 5 K) and an intercept of 32.05 W/(m 2 K).
Figure 5. Linear regression model of volumetric gas flow rate in function of the heat transfer coefficient. The blue dots indicate experimental data, and the red solid line indicates the linear regression model. The linear regression had a slope of 71.11 E3 J/(m 5 K) and an intercept of 32.05 W/(m 2 K).
Pharmaceutics 13 02076 g005
Figure 6. Uncertainty analysis of outer vial wall temperature profiles during cooling and freezing for 6 constant gas flow rates. The red solid lines represent outer vial wall temperature data from experiments as measured by the infrared camera. The black solid lines represent the model prediction of the outer vial wall temperature based on the model input parameters, assuming that no uncertainty is present on the model input parameters. The gray areas indicate the prediction interval of the model prediction of the outer vial wall temperature, taking the uncertainty in the model input parameters into account. Data of the following constant gas flow rates are displayed: A: 32 L/min, B: 44 L/min, C: 50 L/min, D: 56 L/min, E: 68 L/min, and F: 80 L/min.
Figure 6. Uncertainty analysis of outer vial wall temperature profiles during cooling and freezing for 6 constant gas flow rates. The red solid lines represent outer vial wall temperature data from experiments as measured by the infrared camera. The black solid lines represent the model prediction of the outer vial wall temperature based on the model input parameters, assuming that no uncertainty is present on the model input parameters. The gray areas indicate the prediction interval of the model prediction of the outer vial wall temperature, taking the uncertainty in the model input parameters into account. Data of the following constant gas flow rates are displayed: A: 32 L/min, B: 44 L/min, C: 50 L/min, D: 56 L/min, E: 68 L/min, and F: 80 L/min.
Pharmaceutics 13 02076 g006
Figure 7. Uncertainty analysis of one representative outer vial wall temperature profile of the 20 imposed cooling profiles (i.e., which had liquid- and solid-phase cooling rates of 20 ° C/min and a crystallization phase duration of 150 s). The red solid line represents the outer vial wall temperature data from the experiment as measured by the infrared camera. The black solid line represents the model prediction of the outer vial wall temperature based on the model input parameters, assuming that no uncertainty was present in the model input parameters. The gray area indicates the prediction interval of the model prediction of the outer vial wall temperature, taking the uncertainty in the model input parameters into account.
Figure 7. Uncertainty analysis of one representative outer vial wall temperature profile of the 20 imposed cooling profiles (i.e., which had liquid- and solid-phase cooling rates of 20 ° C/min and a crystallization phase duration of 150 s). The red solid line represents the outer vial wall temperature data from the experiment as measured by the infrared camera. The black solid line represents the model prediction of the outer vial wall temperature based on the model input parameters, assuming that no uncertainty was present in the model input parameters. The gray area indicates the prediction interval of the model prediction of the outer vial wall temperature, taking the uncertainty in the model input parameters into account.
Pharmaceutics 13 02076 g007
Figure 8. Global sensitivity analysis of the mechanistic model. Total order effect values S t are shown for every input parameter of the model during solid cooling, crystallization, and liquid cooling at different gas flow rates. Parameters with scores under 0.05 for every gas flow rate and phase of spin freezing are not shown. A: 20 L/min, B: 50 L/min, and C: 80 L/min.
Figure 8. Global sensitivity analysis of the mechanistic model. Total order effect values S t are shown for every input parameter of the model during solid cooling, crystallization, and liquid cooling at different gas flow rates. Parameters with scores under 0.05 for every gas flow rate and phase of spin freezing are not shown. A: 20 L/min, B: 50 L/min, and C: 80 L/min.
Pharmaceutics 13 02076 g008
Table 1. Input parameters and their uncertainty level for the UA and GSA.
Table 1. Input parameters and their uncertainty level for the UA and GSA.
FactorUncertainty ValueReason of Inclusion
h (W/(m 2 K)) 4.3058 RMSE on regression model
Ø v i a l (m) 1 × 10 4 Uncertainty of vial properties
m v i a l (kg) 5 × 10 4 Uncertainty of vial properties
m w a t e r (kg) 3 × 10 5 Measurement error
V ˙ (m 3 /s) 6.67 × 10 6 Uncertainty of the gas flow rate measurement
T e q ( ° C)2Uncertainty of the infrared temperature measurement
T g a s ( ° C)2Uncertainty of the gas temperature measurement
Table 2. Relevant experimental data from verification experiments.
Table 2. Relevant experimental data from verification experiments.
Gas
Flow Rate
(L/min)
Cooling Rate
Liquid Cooling Phase
° C/min)
Cooling Rate
Solid Cooling Phase
° C/min)
Length Crystallization
Phase (s)
T nuc
Outer Vial Wall
° C)
T nuc
Inner Vial Wall
° C)
2018.514.3161−10.6−1.8
2624.625.7134−10.8−0.4
3225.028.6123−9.0−0.6
3827.426.3115−9.6−0.3
4425.826.7116−8.7−0.9
5033.654.587−11.4−2.0
5638.767.180−12.1−0.5
6253.879.365−14.4−0.9
6847.091.366−13.7−0.5
8046.582.960−13.2−0.9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nuytten, G.; Revatta, S.R.; Van Bockstal, P.-J.; Kumar, A.; Lammens, J.; Leys, L.; Vanbillemont, B.; Corver, J.; Vervaet, C.; De Beer, T. Development and Application of a Mechanistic Cooling and Freezing Model of the Spin Freezing Step within the Framework of Continuous Freeze-Drying. Pharmaceutics 2021, 13, 2076. https://doi.org/10.3390/pharmaceutics13122076

AMA Style

Nuytten G, Revatta SR, Van Bockstal P-J, Kumar A, Lammens J, Leys L, Vanbillemont B, Corver J, Vervaet C, De Beer T. Development and Application of a Mechanistic Cooling and Freezing Model of the Spin Freezing Step within the Framework of Continuous Freeze-Drying. Pharmaceutics. 2021; 13(12):2076. https://doi.org/10.3390/pharmaceutics13122076

Chicago/Turabian Style

Nuytten, Gust, Susan Ríos Revatta, Pieter-Jan Van Bockstal, Ashish Kumar, Joris Lammens, Laurens Leys, Brecht Vanbillemont, Jos Corver, Chris Vervaet, and Thomas De Beer. 2021. "Development and Application of a Mechanistic Cooling and Freezing Model of the Spin Freezing Step within the Framework of Continuous Freeze-Drying" Pharmaceutics 13, no. 12: 2076. https://doi.org/10.3390/pharmaceutics13122076

APA Style

Nuytten, G., Revatta, S. R., Van Bockstal, P. -J., Kumar, A., Lammens, J., Leys, L., Vanbillemont, B., Corver, J., Vervaet, C., & De Beer, T. (2021). Development and Application of a Mechanistic Cooling and Freezing Model of the Spin Freezing Step within the Framework of Continuous Freeze-Drying. Pharmaceutics, 13(12), 2076. https://doi.org/10.3390/pharmaceutics13122076

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop