Prediction of Ice Crystal Size Distribution after Freezing of

Subscriber access provided by NEW YORK UNIV

Article

Prediction of ice crystal size distribution after freezing of pharmaceutical solutions Andrea Arsiccio, Antonello A. Barresi, and Roberto Pisano Cryst. Growth Des., Just Accepted Manuscript • DOI: 10.1021/acs.cgd.7b00319 • Publication Date (Web): 20 Jun 2017 Downloaded from http://pubs.acs.org on June 27, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Crystal Growth & Design is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Crystal Growth & Design

Prediction of ice crystal size distribution after freezing of pharmaceutical solutions Andrea Arsiccio, Antonello A. Barresi, and Roberto Pisano

∗

Department of Applied Science and Technology, Politecnico di Torino 24 corso Duca degli Abruzzi, Torino, 10129 Italy E-mail: [email protected]

Abstract In pharmaceutical industry, much attention was given to the problem of prediction of the average size of solvent crystals formed during freezing of pharmaceutical solutions. All the methods proposed in literature to respond to this request are based on empirical laws, while the present work focuses on the development of a mechanistic model. This model was found to give physical theoretical background on the relationship between solvent crystal size, velocity of the freezing front, and temperature gradients within the product being frozen, which was postulated by the empirical laws previously proposed. Model simulations were validated upon experimental observations obtained by Scanning Electron Microscopy. More specically, the average size of solvent crystals was determined from that of pores formed after lyophilization of the frozen solutions, provided that no matrix collapse occurred. The developed model was demonstrated to be valid over a wide range of freezing conditions, and for solutions containing both amorphous and crystalline solutes. As further step, a simplied model for dilute solutions has also been developed and experimentally validated.

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 28

Introduction Prediction of the average size of solvent crystals formed during freezing of a solution is a major issue in various elds and would prove benecial for food, chemicals, and pharmaceutical industries. To give you an example, size of solvent crystals is of utmost importance for those pharmaceutical solutions that have to be freeze-dried, 14 because it strongly inuences both the product dynamics, maximum temperature reached by the frozen product during drying and the processing time itself, 5 and the capability of additives to preserve the active pharmaceutical ingredient against freezing and drying stresses. As a further example, ice crystal size has a strong impact on the organolectic properties of frozen food, 6 while in the case of metal alloys crystal size greatly aects the mechanical properties of the material. 7 As a consequence of the importance of the subject, many researchers tried to nd a quantitative relationship that correlates the thermal history of the product being frozen with the average size of solvent crystals; however, all these models are essentially empirical formulas. For example, early in 1991 Bald 8 proposed an empirical formula that correlates, in the case of dendritic growth, the average size of ice crystals to cooling rate ( ∂T /∂t),

Dp = α

∂T −β ∂t

(1)

where α and β are experimentally determined parameters. In the literature there are a number of empirical laws, having all the same structure, which better detail the eect of cooling rate, correlating the crystal size with the freezing front velocity (ν ) and the temperature gradient within the frozen product ( θ),

Dp = cν −λ1 θ−λ2

(2)

where c, λ1 and λ2 change with the type of product, application and processing conditions. However, attention has to be paid before their application as these parameters were found to

2


Page 3 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


change even within the same application. A summary of the various formulas proposed in literature is given in Table 1. It is therefore desirable to develop a relationship, based on physics and chemistry, that can be applied over a wide range of conditions and independently of the nal application. This paper aims to cover this gap, developing a mechanistic mathematical Table 1: Empirical formulas for crystal sizing Application Metal solidication at low rates Freezing of apples Freezing of pharmaceutical solutions Alloy solidication at high rates Metal solidication at low rates

Correlation Dp ∝ ν −1 θ−1 Dp ∝ ν −0.5 θ−0.5 Dp ∝ ν −0.5 θ−0.5 Dp ∝ ν −0.25 θ−0.5 Dp ∝ ν −0.5

Reference Kurz and Fischer-1 7 Bomben and King 9 Nakagawa et al. 10 Kochs et al. 11 Kurz and Fischer-2 7

model capable to predict the average size of solvent crystals, and its distribution within the product being frozen, knowing the freezing front velocity and temperature gradients. If it is true that the proposed model requires knowledge of a physical parameter, related to the supercial morphology and crystal habit, that at present can hardly be estimated a priori, it must be evidenced that this parameter can be determined by regression of experimental data without loosing any theoretical background of the system and, thus, making it reliable independently of the initial composition of the solution and of the freezing protocol used. Model simulations were compared to experimental observations obtained by Scanning Electron Microscopy (SEM) of lyophilized solutions. The size of pores formed after lyophilization, in fact, corresponds to the size of solvent crystals formed after freezing. The model developed was found to provide physical insight into the process behavior, clarifying the relationship between freezing conditions and the average size of ice crystals. A simplied model was also derived, and experimentally validated, for the case of dilute solutions.

Mathematical formulation The model here presented applies to the crystal growth phase, once nucleation occurred. We aim to develop a lumped parameter model because it can easily be solved; in fact, in 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

lumped models partial dierential equations of the continuous time and space behavior are transformed into ordinary dierential equations with a nite number of parameters, making easier their numerical solution. We further assume that temperature is constant within the liquid layer as well as at the freezing front; in particular, temperature at the freezing front is assumed to be equal to the equilibrium freezing temperature of the solvent. This assumption is even more realistic if we take into account that temperature gradients within the frozen zone are orders of magnitude greater than those within the liquid zone. The model further assumes that the experiment is conducted below the collapse temperature of the formulation, so that no collapse of the sample is observed and the pore size is eectively representative of the average size of ice crystals removed during drying. In model derivation, we divided the control volume into three zones, solid, mushy and liquid, where the term mushy zone refers to a suspension of solid crystals in liquid. We now consider a given volume of product, having cylindrical shape with base diameter

D and height ∆z , in which crystallization takes place. In other words, we can consider ∆z as the fraction of the mushy zone where crystal growth is occurring. So, the control domain considered changes its position during the process, as it follows the mushy zone. In other words, this is a typical Stefan problem, dealing with a moving boundary which separates the liquid zone from the frozen one. As reported elsewhere, 12 the nucleation front moves from the bottom to the top of the product during conventional freezing, while it follows the opposite direction when using some controlled freezing techniques. On the contrary, in most cases the freezing front, dened here as the line dividing the completely frozen zone from the rest of the product, moves from the bottom to the top of the product for both conventional and controlled freezing. However, something dierent could happen for example when facing particular phenomena, such as melting back. A schematic of the control domain is shown in Figure 1, and might be referred to, for example, a portion of product conned within a glass vial that is a typical conguration for the lyophilization of pharmaceuticals.

4


Page 4 of 28

Page 5 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Mushy zone D

Δz (a)

(b)

Figure 1: Scheme of the control domain investigated by the model here presented. As can be seen in Figure 2, we assume that all the crystals are cylinders with diameter

Dp and with an inclination with respect to the axial direction, z , that is described by the angle δ and is linked to tortuosity τ of the porous structure resulting by the freeze-drying of the frozen product, τ = cos δ −1 . As it will be evident in the following, the variable τ will cancel out during mathematical derivation, thus the exact value of the crystal inclination angle does not aect the nal result. τ=

1 cos δ

z

δ

Dp

Δz

τ Δz

δ

Figure 2: Crystals shape as assumed by the model. During the crystal growth process, the freezing front moves along the vial height; thus, as already said, the axial interval ∆z considered, which necessarily follows the freezing front, changes its position during the process. Thus, we will have to discretize the vial height into a nite number n of these axial intervals, each in a dierent position, and the number and 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 28

diameter of the crystals within the ith of these intervals, as well as the other variables which change during the process, will be denoted with subscript i. In addition, it must be pointed out that the total volume, and thus the total height, of the product within the vial changes during the process because of the dierent density of water and ice. For sake of simplicity, a constant ∆z will be assumed for all values of i. From the viewpoint of a lumped parameter model, the enthalpy balance for the given volume is: 

 

 

 



         heat removed by  heat transferred from heat generated by  enthalpy change due to           + + + =0 the frozen product  the lateral surface   crystallization  new interface generation        

−kf θi π

dVice,i dSice,i D2 − πD∆zh∆T + %ice ∆Hf − γ=0 4 dt dt

(3)

where kf and θ are the thermal conductivity and temperature gradient within the frozen zone respectively, ∆T = T − Tair the temperature dierence between air and product being frozen, Vice and Sice the ice crystals volume and surface respectively. Integrating Equation 3 in time from 0 to ∆ti , which is the time interval requested for the mushy zone to move ahead of ∆z , we get: 2

D −kf θ¯i π ∆ti − πD∆zh∆T ∆ti + ∆Vice,i %ice ∆Hf − ∆Sice,i γ = 0 4

(4)

where θ¯i is the mean value of the temperature gradients within the frozen zone during ∆ti . Then, if we assume the initial volume of crystal nuclei to be negligible with respect to the nal dimension of ice crystals, we obtain: 2 Dp,i = Ni π τ ∆z 4

(5)

∆Sice,i = Ni πDp,i τ ∆zaS,i

(6)

∆Vice,i

6


Page 7 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Ni and Dp,i are the nal number and dimension of ice crystals respectively and aS accounts for the real ice crystals habit and surface morphology and is equal to the ratio between the real surface area and the ctitious surface area calculated considering the ice crystals as perfect cylinders of height ∆z :

aS,i =

real surface area Ni πDp,i ∆zτ

(7)

This coecient is actually the product of two terms; the rst, b1 , is related to the real ice crystal habit and can be assumed to be constant while the second, b2,i , takes into account surface irregularities and can be hypothesized to be inversely proportional to the total surface area of crystals. In fact, the greater is the total surface area, the smaller is the contribution given by surface irregularities.

aS,i = b1 b2,i

(8)

If we substitute Equations 5 and 6 in Equation 4 and divide all the terms by ∆ti , we get: 2 Dp,i D2 ¯ − πD∆zh∆T + Ni π %ice τ νi ∆Hf − Ni πDp,i τ νi γaS,i = 0 −kf θi π 4 4

(9)

where ν is the freezing front velocity. We make an additional assumption, the mass of ice that crystallises in 14 πD2 ∆z is directly proportional to θ¯; in fact, if θ¯ increases, the heat released by ice crystallization, and then removed, increases as well:

Ni π

2 Dp,i %ice τ ∆z = aθ¯i 4

(10)

where a is a proportionality coecient. This assumption is much more physically grounded than simply assuming Dp to be proportional to θ−λ2 . Substitution of Equation 10 into Equation 9 leads to the following expression:

D2 νi Ni πDp,i τ νi γaS,i = −kf θ¯i π − πD∆zh∆T + aθ¯i ∆Hf 4 ∆z 7


(11)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 28

which says that the rate of generation of crystals surface area within the control domain increases with the temperature gradients within the frozen zone, and thus with the total mass of water which crystallizes, and with the freezing front velocity. This is again in accordance with experimental observations; in fact, at equal mass of ice formed, the rate of generation of new surface area diminishes if the freezing process is slower. This result again conrms the validity of the assumption by Equation 10. Next, we consider the mass balance for the given volume of liquid being frozen: n X i=1

Ni π

2 Dp,i %ice τ ∆z = φm 4

(12)

where m is the total mass of water of the solution and φ is a coecient close or equal to one, at least for dilute solutions, that accounts for the mass of water which remains in the amorphous phase and does not crystallize. In crystalline materials, nearly all water is frozen and thus we can assume φ = 1. 13 However, in the case of amorphous solutes, the solid solution phase usually contains between 15% and 50% of unfrozen water, but except for very concentrated solutions the value of φ is anyway close to one. 1 In the following, we will assume for semplicity that φ is equal to 1. Combining Equations 10 and 12 it is possible to determine an expression for a,

m a= P n θ¯i

(13)

i=1

and thus Equation 10 can be rewritten as, 2 Dp,i θ¯i Ni π %ice τ ∆z = m P n 4 θ¯i i=1

8


(14)

Page 9 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Substituting Equation 14 into Equation 9 and rearranging gives:

Dp,i =

16mθ¯i γaS,i νi n n P P %ice (4mθ¯i νi ∆Hf − kf θ¯i D2 π∆z θ¯i − 4π∆2 z θ¯i Dh∆T ) i=1

(15)

i=1

where, as anticipated, τ cancels out. Equation 15 describes the eects of the freezing front velocity ν and temperature gradients θ on Dp . The freezing process is extremely complex and several other factors should, rigorously, be taken into account. However, previous works 7,911 have shown that, on constant nucleation temperature, freezing behavior is satisfactorily described by two parameters, the temperature gradients within the frozen zone and the freezing front velocity. However, it is also true that the proposed model indirectly accounts for other factors, namely nucleation temperature, cooling rate and product geometry, during the calculations of θ and ν . Thus,

θ and ν implicitly take into account all the other relevant factors involved. As previously said, the aS coecient, which takes into account the real crystal habit and the surface irregularities, is the product of a constant term and of a term which is assumed to be inversely proportional to the surface area of crystals. Equation 10 states that the volume of ice crystallized in 14 πD2 ∆z is proportional to θ¯. Hence, if the volume is proportional to

θ¯, the surface area of the crystals can be assumed to be proportional to θ¯2/3 . As a result:

aS,i = b1 b2,i =

b 2/3 θ¯i

(16)

where b is determined by regression of experimental observations for pore size vs axial position. Unlike previous approaches, Equation 15 has been obtained from heat and mass balance equations. Therefore, it might be used to provide insight into the freezing process behaviour, and was here demonstrated to be applicable over wide ranges of conditions. Equation 15 has also the advantage of automatically respecting the mass balance. However, the model is not entirely theoretical. In fact, even if it uses balance equations,

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 28

the numerical value of one parameter, namely γb, was calculated from experimental data. In spite of this, all the terms of the equation derived have a physical meaning and a clear explanation, even if some of them can hardly be predicted theoretically.

Simplied model for dilute solutions For dilute solutions, the model can be further simplied. More specically, we can assume that the term N Dp is large and hence heat transferred from the lateral surface by natural convection is negligible in Equation 3. For example, if water constitutes 95% of the whole system mass and air inside the freezing chamber is stationary, heat removed by natural convection is about 0.5% of the heat generated by water crystallization. Moreover, we can write that, 2 Dp,i D2 ≈ Ni π τ επ 4 4

(17)

where ε is the ratio between the volume of ice and the total volume of the system. Thus, combining Equation 9 with Equation 17:

Dp,i ≈

4εγaS,i νi ε%ice ∆Hf νi − kf θ¯i

(18)

where aS is given in Equation 16. The crystals inclination angle, that denes tortuosity, has been introduced for sake of generality in writing balances, but it is not relevant for the calculation, as it cancels out in the nal Equations for ice crystal-size.

10


Page 11 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Materials and methods Materials and Instrumentation The formulations investigated in this work were mannitol solutions at dierent concentrations (5% and 10% w/w) and sucrose or dextran solutions at 5% w/w, prepared from distilled water and mannitol, sucrose or dextran powder; all reagents were purchased from Sigma Aldrich and used as received. Here, we used sugar-based formulations because they are the most important excipients used for those drugs to be freeze-dried. For example, mannitol is used as bulking agent, whereas sucrose and dextran are used as cryo- and lyo-protectants. The vials employed (type 1, 10R, 45x24 mm, Schott AG, Germany), having 24 mm diameter and 45 mm height, were accurately lled with 3 ml of sample solution, while igloo stoppers (type 1319 4432/50/ Westar, West Pharmaceutical Services, Terrassa, Spain) were used as closures. A freeze-dryer LyoBeta 25 (Telstar, Terrassa, Spain) was used for the freeze-drying cycles. In order to monitor the temperature of the shelves and of the product we used a system of T-type copper/constantan miniature thermocouples. The experiments performed for the validation of the model were carried out using two dierent freezing modalities, namely conventional freezing and the controlled freezing technique introduced by Oddone et al. 14 For details on the freezing procedure, see the Supporting Information, section S1. The operating conditions of the experimental tests used for the validation of the mathematical model are listed in Table 2, where Tn and Tm are the holding temperature before and after the induction of nucleation respectively. Operating conditions were chosen so as to avoid product collapse.

Characterization of product morphology The pore dimension of the products obtained after freeze-drying, which corresponds to the dimension of ice crystals formed during freezing, was analyzed using a Scanning Electron 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 28

Table 2: Operating conditions of experimental tests solution A B C D E F G H I

Nucleation type

mannitol 5% w/w sucrose 5% w/w dextran 5% w/w mannitol 10% w/w mannitol 5% w/w mannitol 5% w/w mannitol 5% w/w mannitol 5% w/w sucrose 5% w/w

conventional conventional conventional conventional conventional controlled controlled controlled controlled

∂T , ∂t

K min−1 0.8 0.8 0.8 0.8 0.3 -

Tn , K

Tm , K

263 268 268 268

263 268 263 263

Microscope (SEM, FEI type, Quanta Inspect 200, Eindhoven, the Netherlands). SEM images were taken at various positions along the previously metallized sample in order to obtain a quantitative estimation of within-vial heterogeneity. This quantitative analysis was done in terms of size (Dp ) and distribution of the pore size. For each sample, three images were taken at the top, center and bottom of the product. Approximately 100 pores per image were selected and each of them was approximated to an ellipse. The dimension Dp we attributed to the pore was the diameter of the circle having the same area to perimeter ratio of the approximating ellipse.

Model validation protocol The validation of the model involved several steps. Firstly, the product temperature proles during freezing were measured experimentally and compared with numerical simulations carried out using the software COMSOL. The simulated geometry consisted of the vial and the product, see Figure 1(b). A mapped grid was used for the domain that represents the product and a triangular unstructured grid for the domain that represents the side of the vial. 552 cells were used for the product, while 351 were employed for the vial wall. Simulation results were veried to be grid independent. 12


Page 13 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A short description of the main features of the model implemented in COMSOL, which was originally developed by Nakagawa et al. for conventional freezing 10 and Pisano and Capozzi for controlled freezing, 15 is given in the Supporting Information, section S2. After calculation of temperature proles, the freezing front velocity ν and temperature gradients within the frozen zone θ were estimated as described in the Supporting Information (section S3) and substitued in Equations 15 and 18, thus providing the model estimation of the average solvent crystals dimension and of its distribution. Finally, the model predictions were compared to experimental data obtained with Scanning Electron Microscopy (SEM).

Results and Discussion The number of crystal nuclei basically depends on nucleation temperature. But if the time left for temperature equilibration is too short, large temperature gradients are generally observed within the liquid being frozen, and crystal nuclei are not uniformly distributed. 16 In addition to this, these nuclei, once formed, could diuse, for example as a consequence of Brownian motion. As general guideline, low equilibration temperatures result in signicant crystal size heterogeneity. Moreover, the nal product structure is linked to the crystals aspect ratio which, in its turn, depends on the crystals growth rate. If the growth rate is controlled by diusion, crystals mainly grow in the direction where the mass transfer rate is higher, assuming a hoppered shape. On the contrary, if the rate determining step is the surface kinetics process, that is the process in which atoms are incorporated into the crystal, the aspect ratio is mainly determined by the dierences in the intrinsic growth rate of the dierent crystal surfaces. 17,18 The dierences in number and dimension of crystals along the product height can also be explained if we assume that, in some situations, crystals preferentially grow in the radial direction with respect to the axial one and vice versa. In the case of preferential growth along the radial direction, crystals are mostly conned within the zone where they are formed and

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 28

thus the product can be divided into layers where crystal growth is independent of one another and where the number of crystals is determined exclusively by the number of nuclei formed during nucleation. On the contrary, in the case of preferential growth in the axial direction, these layers are no more independent of one another and the number of ice crystals formed in a layer is also inuenced by the number of crystals coming from the surrounding layers. The parameter b dened in Equation 16 could also implicitly be linked to all these microscopic aspects of cristal growth, even if the model derived is a macroscopic one. An example of SEM images of the samples analyzed is shown in Figure 3.

Center

Top

Bottom

400 μm

400 μm

400 μm

400 μm

400 μm

400 μm

Test A

Test F Figure 3: Scanning electron microscope pictures of mannitol 5% w/w lyophilised samples in the case of conventional freezing (test A in Table 2) and controlled freezing (test F in Table 2).

A list of the numerical values of the parameters used in the model is reported in Tables 3 and 4. As can be seen from Table 4, the value of the parameter γb does not change for the two amorphous solutes tested, namely sucrose and dextran, while it is signicantly dierent for mannitol. Moreover it is not aected by cooling rate, nucleation temperature or solute 14


Page 15 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


concentration. In addition, γb increases for the controlled freezing technique with respect to conventional freezing, but the degree of increase is greater for mannitol solutions than for sucrose solutions. Table 3: Numerical values of the parameters employed

kf h %ice ∆Hf

2.5 5 918 333500

W m−1 K−1 W m−2 K−1 kg m−3 J kg−1

Table 4: Numerical values of parameter γb

solute mannitol sucrose dextran

γb, J K2/3 m−8/3 conventional freezing controlled freezing 7·104 26·104 23·104 32 ·104 23·104 -

As a rst step, the model was validated for dierent formulations, A and B in Table 2, and using the same freezing protocol, i.e. conventional freezing. As shown in Figures 4(a) and 5(a), in the case of the lowest solid content (5% w/w), both the detailed model and the simplied one were found to be appropriate for the description of the freezing of the formulations considered. We can thus deduce that both models can be applied to both crystalline and amorphous systems. An aspect that must be highlighted is the extreme sensibility of these models to small variations in ν and θ. It is therefore of utmost importance to check the adherence of the temperature proles predicted by the model to those measured experimentally. Tests B and C were carried out using two amorphous solids, sucrose and dextran, with the same γb value. Even the cooling rate and freezing protocol were the same for both tests. However, the two samples were characterized by a signicantly dierent nucleation temperature, 267 K for sucrose and 258 K for dextran. This is reected by the much smaller 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

pores for test C if compared to test B, as can be seen comparing Figures 5(a) and 5(b). Therefore, model predictions consistently predicted that crystal size increases as nucleation temperature increases. While tests A-C were carried out on constant solid content, 5% (w/w), test D was carried out using 10% as solid content. This further test conrmed that model predictions were consistent with experimental observations in the range of 5 to 10% (w/w). This was also true for the simplied model, which has been derived for dilute solutions. In this study we focused on pharmaceutical formulations to be freeze-dried, whose solid content is generally of the order of or smaller than 10% (w/w). So far, the same cooling rate of 0.8 K min −1 has been considered. We decided to demonstrate the validity of the models for dierent operating conditions (0.3 K min −1 , case E in Table 2). Even in this case, both models predicted fairly well the dimension of ice crystals, as can be seen in Figure 4(c). The mean ice crystal diameter for test A was 23 µm, while for test E was 26 µm, again in line with the fact that the average crystal diameter increases with decreasing cooling rate. However, it can also be noted that the experimental dispersion of pores in Test A was much larger than that for Test E. This could be explained by the fact that ice crystals grow more uniformly at lower cooling rates, because time left to crystal growth is greater and water molecules have sucient time to arrange themselves in a regular way. Another possible explanation could be related to the phenomenon of Ostwald ripening. Parameters such as cooling rate and nucleation temperature are taken into account by the model by ν and θ. If either cooling rate decreases or nucleation temperature increases, ν and

θ decrease; according to the model, this means that, at equal mass of ice formed, the total surface area decreases, thus meaning that the ice crystal diameter is increased. These results demonstrate the validity of the models derived and their wide range of possible applications. For all tests performed, the percentage variation between the predictions of the detailed model and the simplied one is similar. In fact, although for tests B and C, the dierence between the predictions of the detailed model and the simplied one seems bigger with

16


Page 16 of 28

Page 17 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


respect to the case of tests A, D and E, this dierence disappears if we consider that the pores of products B and C were much bigger with respect to those of products A, D and E. The models were also validated using a dierent freezing protocol, namely controlled freezing. Oddone et al. 14,19 showed that conventional freezing produces mostly small and spherulitic ice crystals, while control of freezing gives large, chimney-like pores. In fact, controlled freezing gives the possibility to induce nucleation at given temperature and indirectly control temperature gradients within the liquid being frozen by adjusting Tn and Tm . So, three dierent tests using mannitol (F, G and H in Table 2) were carried out changing these two values and, thus, modifying the product nucleation temperature. As can be seen in Figures 6(a) and 6(b), both models gave accurate quantitative estimations of the average pore size and its variation along the axial direction of the product, demonstrating that they can eectively be used independently of the freezing protocol. The pore size generally increases as Tn increases; in fact for test H the average pore dimension predicted by the model was greater (100 µm, data not shown) with respect to test F (85 µm), in accordance with experimental data. A further test was carried out using an amorphous solute, sucrose. An example of results is given in Figure 6(c). We can again observe a fairly good agreement between model predictions and experimental observations, both in terms of average pore size and its variation along the axial direction. In contrast with the general eect of Tn , the pore size of samples produced by test G was signicantly smaller than that observed for test F, see Figures 6(a) and 6(b). We hypothesized that this result was due to the fact that after stabilization, the solution was not completely frozen and thus it freezes as soon as the temperature of the heat transfer uid is reduced to -45◦ C. In fact, it must be remarked that test G was carried out at a higher value of Tm (268 K) and, thus, this temperature was likely too high to allow completion of freezing during the stabilization step. As a consequence, a part of the solution was frozen later on and freezing of that part occurred as in a conventional freezing cycle. This phenomenon

17



mannitol 5%, conv.

80 pore size,

m

(a) Test A

60 40 20 0

0

2

4

z,

6

8

10

mm

mannitol 10%, conv.

80

(b) Test D

pore size,

m

60 40 20 0

0

2

4

z,

6

8

mm

mannitol 5%, conv. low cool. rate

80

(c) Test E

m

60 pore size,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 28

40 20 0

0

2

4

6

z, mm

8

10

Figure 4: Comparison between detailed model predictions (solid red line), simplied model predictions (dashed green line) and SEM observations (box plots) for the average pore size in the case of tests A, D and E (mannitol solutions, with conventional freezing). For box plots, top bar is maximum observation, lower bar is minimum observation, top of box is third quartile, bottom of box is rst quartile, middle bar and internal square are median and average value respectively and bullets are original data. For box plots, the x-axis is that of boxes. 18


Page 19 of 28

sucrose 5%, conv.

200

(a) Test B

pore size,

m

150 100 50 0

0

2

4

6

z, mm

8

10

dextran 5%, conv.

90

m

(b) Test C

pore size,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


75 60 45 0

2

4

6

8

10

z, mm

Figure 5: Comparison between detailed model predictions (solid red line), simplied model predictions (dashed green line) and SEM observations (box plots) for the average pore size in the case of tests B and C (sucrose or dextran solutions, with conventional freezing). For box plots, top bar is maximum observation, lower bar is minimum observation, top of box is third quartile, bottom of box is rst quartile, middle bar and internal square are median and average value respectively and bullets are original data. For box plots, the x-axis is that of boxes.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 28

was already observed by Oddone et al. 19 and the simulated temperature proles for test G conrmed that the growing phase could not be completed during the holding time at

Tm , but occurred at least in part during the temperature ramp from Tm to -45◦ C, insets of Figures 6(a) and 6(b). The model here presented could predict crystal size even for this very particular situation. Finally, we decided to compare the predictions of the mechanistic model here presented with those of the empirical laws reported in Table 1. For this comparison we chose two different experimental conditions, namely test E, which was performed employing conventional freezing, and test G, which was conducted using the controlled freezing technique. The result of such comparison is shown in Figure 7, where it is possible to see that the mechanistic model, in both the detailed and simplied version, worked better than the empirical laws. The empirical laws by Bomben and King 9 and Kochs et al. 11 seemed to work fairly well. On the contrary, the equations proposed by Kurz and Fischer 7 were less accurate; moreover they did not have a stable behaviour, as can be seen from the fact that they were quite accurate in one situation (conventional freezing for Kurz and Fischer-2 and controlled freezing for Kurz and Fischer-1), while they were imprecise in the other cases. Even for amorphous solutes, the empirical laws by Bomben and King and Kochs et al. worked better than those proposed by Kurz and Fischer, even if not as well as the mechanistic model here presented (data not shown). The values of the c coecient employed for applying the empirical laws are reported in Table 5. Table 5: Numerical values of parameter c employed

c, m1+λ1 −λ2 s−λ1 Kλ2 Bomben and King 9 Kochs et al. 11 Kurz and Fischer-1 7 Kurz and Fischer-2 7

Test E

1.6·10−6 28·10−6 0.1·10−6 0.09·10−6

20

Test G

7·10−6 120·10−6 0.5·10−6 0.32·10−6


°C

0

-1

=

T

. n

tr

n

o

c

,

%

5

l

o

it

n

n

a F

t

s

e

T

)

(a

5

0

1 0

9

0

2 0 0 0 0

4

6

0

6

8

0

T 2 e m p 2 e r a t 2u r e , 3 K

5

7

3

p o r e s iz e , m

3

2

,h

e

im

T

1

0

0

2

2

5

4

0

1

8

6 °C

-5

= n

T

.

m tr

n

m o

c

,

%

5

l

o

it

n

n

a

m

z,

4

2

0

0

3

G

t

s

e

T

)

(b

5

0

1 0

9

0

2 0 0 0 0

4

6

8

0

T 2 e m p 2 e r a t 2u r e , 3K

3

p o r e 6 s i z e 7,  m 0 5

3

2

,h

e

im

1

T

0

0

2

2

5

4

0

1

8

6 m

m

z, 4

2

0

0

3

sucrose 5%, contr.

m

200

pore size,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


m

Page 21 of 28

(c) Test I

150 100 50 0

2

4

z, mm

6

8

10

Figure 6: Comparison between detailed model predictions (solid red line), simplied model predictions (dashed green line) and SEM observations (box plots) for the average pore size in the case of tests F, G and I (controlled freezing). The insets show the temperature proles of the product (solid line) and of the shelf (dashed line).

21



mannitol 5%, conv.

pore size,

m

40 30 20 10

0

2

4

6

8

z, mm

10

mannitol 5%, contr.

100

m

80 pore size,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 28

60

40

0

2

4

z, mm

6

8

10

Figure 7: Comparison between detailed model predictions (solid red line), simplied model predictions (dashed green line) and empirical laws described in: Bomben and King 9 (4 ), Kochs et al. 11 (), Kurz and Fischer-1 7 (- - - - - -) and Kurz and Fischer-2 7 (♦) for tests E and G. SEM observations are reported in the box plots. Only the boxes are reported for sake of clarity

22


Page 23 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Conclusions In this paper we developed a detailed model and a simplied one, valid for dilute solutions, capable to predict the solvent crystal dimension of frozen products. We also veried the models predictions adherence to experimental values. The models were demonstrated to be eective for (a) various formulations, (b) dierent concentrations (at least up to 10% w/w), and (c) dierent operating conditions (cooling rate and nucleation temperature). Unlike previous approaches, all the terms of our models have a strong physical meaning and a clear explanation. Moreover, the relationship between the solvent crystal diameter, the freezing front velocity and the temperature gradients within the frozen zone, which was postulated but never found a theoric demonstration, has nally found a physical explanation in our models. The detailed model developed has a wider range of validity; however, the simplied model is easier to implement and has proved to work well in the concentration range tested. However, not all the problems have been solved; in fact there are still some terms in the models that must be determined from experimental data. So future work should analyze more deeply the growth of crystals during freezing, in order to make crystal size completely predictable without need of experimental data, thus giving the possibility to taylor in advance the crystal structure of the nal product to a desired morphology.

Supporting Information . Details of freezing procedure, brief description of freezing models, procedure for determination of model parameters, ν and θ.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Nomenclature a

parameter of Equation 10, kg m K −1

aS

empirical coecient that accounts for the real ice crystals surface, -

b

coecient of Equation 16, K 2/3 m−2/3

b1

corrective coecient which takes into account the real crystal habit, -

b2

corrective coecient which takes into account the surface irregularities, -

c

empirical parameter of Equation 2

cˆp

specic heat capacity, J kg −1 K−1

cˆp,`

specic heat capacity of the liquid, J kg −1 K−1

D

vial base diameter, m

Dp

crystals diameter, m

h

heat transfer coecient, W m −2 K−1

∆Hf

latent heat of crystallization, J kg −1

k

thermal conductivity, W m −1 K−1

kf

thermal conductivity of the frozen zone, W m −1 K−1

k`

thermal conductivity of the liquid, W m −1 K−1

m

mass of water, kg

n

number of discretization intervals,-

N

crystals number, -

Qn

source term representing the generation due to nucleation, J s −1 m−3

Qc

source term representing the generation due to crystal growth, J s −1 m−3

Sice

surface of ice crystals, m 2

T

temperature, K

Tair

air temperature, K

Tm

holding temperature after the induction of nucleation, K

24


Page 24 of 28

Page 25 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Tn

holding temperature before the induction of nucleation, K

∆T

temperature dierence between air and product, K

t

time, s

∆t

time interval, s

Vice

volume of ice, m3

z

axial coordinate, m

zf

axial coordinate of the freezing front, m

∆z

axial interval, m

∆zf

variation of the freezing front position, m

Greek letters α

empirical parameter of Equation 1

β

exponent of Equation 1, -

γ

solid-solid interfacial tension, J m −2

δ

crystals inclination angle

ε

ratio between the volume of ice and the total volume of the system,-

θ

temperature gradient within the frozen zone, K m −1

θ¯

time average of the temperature gradient within the frozen zone, K m −1

λ1

exponent of Equation 2

λ2

exponent of Equation 2

ν

freezing front velocity, m s −1

%

density, kg m−3

%ice

density of ice, kg m−3

%`

density of the liquid, kg m −3

τ

tortuosity, -

φ

fraction of water which crystallizes, -

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abbreviations SEM

Scanning Electron Microscopy

References (1) Franks, F. Freeze-Drying of Pharmaceuticals and Biopharmaceuticals ; RSC Publishing: Cambridge, UK, 2007. (2) Rey, L.; May, J. Freeze-Drying/Lyophilization of Pharmaceutical and Biological Prod-

ucts ; Marcel Dekker, Inc: New York, USA, 2004. (3) Hottot, A.; Vessot, S.; Andrieu, J. Chem. Eng. Process. 2007, 46, 666674. (4) Bosca, S.; Barresi, A.; Fissore, D. Drying Technol. 2015, 33, 10391050. (5) Searles, J.; Carpenter, J.; Randolph, T. J. Pharm. Sci. 2001, 90, 860871. (6) Singh, R.; Heldman, D. Introduction to Food Engineering ; Academic Press: San Diego, USA, 2014. (7) Kurz, W.; Fischer, D. Fundamentals of Solidication ; Trans Tech Publications: Switzerland, 1992. (8) Bald, W. Food Freezing: Today and Tomorrow ; Springer-Verlag: New York/London, 1991. (9) Bomben, J.; King, C. J. Food Technol. 1982, 17, 615632. (10) Nakagawa, K.; Hottot, A.; Vessot, S.; Andrieu, J. AIChE J. 2007, 53, 13621372. (11) Kochs, M.; Korber, C.; Heschel, I.; Nunner, B. Int. J. Heat Mass Transfer 1991, 34, 23952408.

26


Page 26 of 28

Page 27 of 28

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12) Oddone, I.; Fulginiti, D.; Barresi, A.; Grassini, S.; Pisano, R. Drying Technol. 2015,

33, 16211630. (13) Nail, S.; Akers, M. Development and Manufacture of Protein Pharmaceuticals ; Kluwer Academic/Plenum Publisher: New York, USA, 2002. (14) Oddone, I.; Pisano, R.; Bullich, R.; Stewart, P. Ind. Eng. Chem. Res. 2014, 53, 18236 18244. (15) Pisano, R.; Capozzi, L. Submitted to Chem. Eng. Res. Des. (16) Kasper, J.; Friess, W. Eur. J. Pharm. Biopharm. 2011, 78, 248263. (17) Nishinaga, T.; Nishioka, K.; Harada, J.; Sasaki, A.; Takei, H. Advances in the Under-

standing of Crystal Growth Mechanisms ; Elsevier Science: Amsterdam, 2012. (18) Hannay, N. Changes of State -Treatise on Solid State Chemistry ; Springer Science & Business Media: New York, USA, 2013. (19) Oddone, I.; Bockstal, P. V.; Beer, T. D.; Pisano, R. Eur. J. Pharm. Biopharm. 2016,

103, 167178.

27



For Table of Contents Use Only Prediction of ice crystal size distribution after freezing of pharmaceutical solutions Andrea Arsiccio, Antonello A. Barresi, and Roberto Pisano*

Department of Applied Science and Technology, Politecnico di Torino 24 corso Duca degli Abruzzi, Torino, 10129 Italy E-mail: [email protected]

crystal size, m

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 28

90 75 60 45 30

400

0

2

4

6

8

z, mm

m

10

A mechanistic model is developed for the prediction of solvent crystal-size distribution after freezing of pharmaceutical solutions. Model simulations were validated upon experimental observations obtained by SEM for various formulations and varying freezing conditions.

28


Prediction of Ice Crystal Size Distribution after Freezing of

Recommend Documents