Developing a framework to model the primary drying step of a continuous freeze-drying process based on infrared radiation

Accepted Manuscript Developing a framework to model the primary drying step of a continuous freeze-drying process based on infrared radiation Pieter-Jan Van Bockstal, Jos Corver, Séverine Thérèse F.C. Mortier, Laurens De Meyer, Ingmar Nopens, Krist V. Gernaey, Thomas De Beer PII: DOI: Reference:

S0939-6411(17)31164-5 https://doi.org/10.1016/j.ejpb.2018.02.025 EJPB 12705

To appear in:

European Journal of Pharmaceutics and Biopharmaceutics

Received Date: Revised Date: Accepted Date:

10 October 2017 9 February 2018 19 February 2018

Please cite this article as: P.V. Bockstal, J. Corver, S.T.F.C. Mortier, L.D. Meyer, I. Nopens, K.V. Gernaey, T.D. Beer, Developing a framework to model the primary drying step of a continuous freeze-drying process based on infrared radiation, European Journal of Pharmaceutics and Biopharmaceutics (2018), doi: https://doi.org/10.1016/ j.ejpb.2018.02.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Developing a framework to model the primary drying step of a continuous freeze-drying process based on infrared radiation Pieter-Jan Van Bockstala,1,∗, Jos Corvera,∗∗, S´everine Th´er`ese F.C. Mortiera,b , Laurens De Meyera , Ingmar Nopensb , Krist V. Gernaeyc , Thomas De Beera a Laboratory

of Pharmaceutical Process Analytical Technology, Department of Pharmaceutical Analysis, Faculty of Pharmaceutical Sciences, Ghent University, Ottergemsesteenweg 460, 9000 Ghent, Belgium b BIOMATH, Department of Mathematical Modelling, Statistics and Bioinformatics, Faculty of Bioscience Engineering, Ghent University, Coupure Links 653, 9000 Ghent, Belgium c Process and Systems Engineering Center (PROSYS), Department of Chemical and Biochemical Engineering, Technical University of Denmark, Building 229, 2800 Kgs. Lyngby, Denmark

Abstract The continuous freeze-drying concept based on spinning the vials during freezing and on non-contact energy transfer via infrared (IR) radiation during drying, improves process efficiency and product quality (uniformity) compared to conventional batch freeze-drying. Automated control of this process requires the fundamental mechanistic modelling of each individual process step. Therefore, a framework is presented for the modelling and control of the continuous primary drying step based on non-contact IR radiation. The IR radiation emitted by the radiator filaments passes through various materials before finally reaching the spin frozen vial. The energy transfer was computed by combining physical laws with Monte Carlo simulations and was verified with experimental data. The influence of the transmission properties of various materials on the emitted IR radiation profile was evaluated. These results assist in the selection of proper materials which could serve as IR window in the continuous freeze-drying prototype. The modelling framework presented in this paper fits the model-based design approach used for the development of this prototype and shows the potential benefits of this design strategy by establishing the desired engineering parameters and by enabling the engineer to assess mechanical tolerances and material options. Keywords: Continuous freeze-drying, Infrared radiation, View factor, Monte Carlo simulation

∗ Corresponding

author first authorship Email addresses: [email protected] (Pieter-Jan Van Bockstal), [email protected] (Jos Corver), [email protected] (S´ everine Th´ er` ese F.C. Mortier), [email protected] (Laurens De Meyer), [email protected] (Ingmar Nopens), [email protected] (Krist V. Gernaey), [email protected] (Thomas De Beer) 1 Phone number: +32(0)9 264 83 55; Fax number: +32(0)9 264 82 36 ∗∗ Shared

Preprint submitted to European Journal of Pharmaceutics and Biopharmaceutics

February 13, 2018

1. Introduction Freeze-drying is a low temperature drying process commonly used to improve the stability of biopharmaceutical drug products during storage and distribution [1]. Conventional pharmaceutical freeze-drying of unit doses is a batch-wise process which involves three consecutive steps: freezing, primary drying and secondary drying [2]. The complete batch of glass vials containing the aqueous drug solution is loaded onto the shelves in the drying chamber. The freezing step is initiated by lowering the temperature of these shelves, leading to ice nucleation in the formulation and crystal growth upon further cooling. The formed ice crystals are removed through sublimation during the primary drying step by lowering the pressure inside the drying chamber and raising the temperature of the shelves to provide the sublimation heat. Any remaining unfrozen water dissolved in the amorphous solid is removed during the secondary drying step by further increasing the temperature of the shelves, also under vacuum, resulting in the dried end product. This conventional batch approach is inherently associated with several disadvantages, which causes various issues in terms of quality assurance of drug products and inefficiencies of the production process [3, 4]. The freezing step is uncontrolled and the energy transfer to each individual vial of the batch during drying is inhomogeneous. Both phenomena could result in a possible uncontrolled vial-to-vial end product variability, within a batch and between batches. The end product quality is only tested for a limited number of vials, before releasing the complete batch. This quality approach does not meet the most recent guidelines issued by the regulatory authorities regarding Quality-by-Design (QbD) and Process Analytical Technology (PAT), which state that quality should be guaranteed by building it in the product [5, 6]. To overcome this poor quality assurance approach and other disadvantages related to the batch process, a continuous concept for the freeze-drying of unit doses is proposed, where each single process step is integrated in a continuous production flow (Figure 1) [3, 4, 7]. The continuous freeze-drying process is initiated by rapidly rotating glass vials filled with the aqueous drug formulation along their longitudinal axis (Figure 1) [3, 4]. The flow of a cold, inert and sterile gas cools the solution, eventually inducing ice nucleation [7]. At the end of this spin freezing step, the product is solidified over the entire inner vial wall resulting in a much thinner product layer compared to the conventional frozen plug at the bottom of the vial, reducing the total process time for pharmaceutical compositions with a factor 10 to 40, depending on the vial dimensions and product specifications [3, 4]. The spin frozen vials are continuously transferred to a long belt in a temperature-controlled chamber, which allows the crystallization and solidification of the solutes under standardized conditions (annealing) (Not included in Figure 1). When the desired morphological structure is obtained, the spin frozen vials are further processed to the 2

Figure 1: Schematic illustration of the continuous freeze-drying system consisting of multiple modules connected via load-locks

primary drying unit. Here, the pressure is kept at a constant value between 10 to 30 Pa. An appropriate load-lock system separates both modules to facilitate the vial transfer while maintaining the specific conditions of pressure and temperature in each chamber [8]. Continuous primary drying requires an adequate and uniform energy transfer towards the entire vial wall surface to achieve an efficient and homogeneous drying behaviour. Therefore, non-contact IR radiation was evaluated and was shown to be advantageous compared to conduction in providing the energy for the drying of spin frozen vials [4]. Each spin frozen vial rotates along its longitudinal axis in front of an individual IR heater, allowing individual energy input regulation. Rotation of the vials is a prerequisite to obtain a homogeneous heat transfer over the entire vial surface. Movement of the vials takes place in discrete steps so each vial is almost constantly centred in front of an individual IR heater. In case secondary drying requires a pressure level different from primary drying, a second continuous drying module is necessary, connected via an appropriate load-lock system for vial transfer [7]. At the end of the continuous lyophilisation process, vials are removed from the vacuum conditions in the drying module via another load-lock system and transferred to the final unit, which allows the stoppering of the processed vials, in a controlled and inert nitrogen environment. The vial throughput (i.e., scale-up) can be increased by multiplying the continuous production lines, as illustrated in Figure 2. Hence, process and formulation optimization can start at the very beginning of product development as scale-up will not require re-optimization and re-validation of the process. 3

Figure 2: Scale-up by multiplication approach of a continuous freeze-drying system

The energy transfer to sustain sublimation and desorption during the continuous processing of spin frozen vials is effectuated through non-contact IR radiation [4]. Each vial rotates freely in front of an individual IR heater, which allows 100% monitoring of Critical Quality Attributes (CQAs) and control of the process settings at the individual vial level. Therefore, the vial-to-vial and batch-to-batch variability can be significantly reduced. Process optimization and control are supported by the application of mechanistic models which assess the radiative heat transfer [9]. With the current prototype system, the IR radiators are positioned inside the vacuum drying chamber. In this configuration, the IR heaters are considered as a uniformly heated flat surface of 45 mm by 45 mm, as specified by the supplier. For the industrial system, for Good Manufacturing Practice (GMP) reasons, the IR radiators will be placed outside the vacuum chamber behind an IR transparent window, as illustrated in Figure 3, in order to be compatible with the Cleaning-in-Place (CIP) and Sterilization-in-Place (SIP) procedures. The IR radiation is generated by heating an electrical filament, which is covered by a tube of fused silica. One side of this tube is coated with gold to direct the radiation in the desired direction. The fused silica tube absorbs a fraction of the IR radiation emitted by the heated filament. Also the window positioned between the IR heater and the (borosilicate) glass vial absorbs a fraction of the radiation, depending on the transmission properties of the material. Specific materials have to be selected for the IR windows: transmission of the IR radiation must be optimal, while the vacuum in the chamber has to be maintained. The objective of this publication is to provide a framework for the modelling and control of the continuous primary drying step based upon IR radiation. These models can assist in the prototype system design via a model-based design strategy to determine the influence of various engineering parameters, e.g., the geometry of the radiator, the distance between the radiator and the glass vial, the glass vial dimensions, the IR window material and the thickness of the IR window. Initially, the

4

configuration where the vial is situated inside the drying chamber is considered. Here, the IR radiator is assumed to be a uniformly heated flat surface with the cylindrical vial centred in front of this plane. This approach was used by Van Bockstal et al. to compute the radiative energy transfer based on the current prototype system [9]. In addition, the situation with the radiator situated outside the vacuum chamber, radiating through an IR transparent window, and the situation where the vial is not centred in front of the radiating plane were considered since these are the configurations that are anticipated for later industrial applications.

Figure 3: Illustration of the path for radiant heat

In order to assess the effect of radiant heat, the laws of Stefan-Boltzmann, Planck and Lambert-Beer have to be taken into account. The Stefan-Boltzmann law (Equation 1) is used to determine the radiant power that is transmitted between two surfaces of unequal temperature, Planck’s law (Equation 4) indicates the spectral power distribution with the temperature as a parameter and the Lambert-Beer law (Equation 5) specifies the transmittance/absorptivity of a certain material. For the continuous freeze-drying set-up, the Stefan-Boltzmann law is given by [9]:

4 P = Arad F σ(Trad − aTv4 )

(1)

with P the power provided by the emitting surface (i.e., the IR heater or the chamber walls) to the target surface (i.e., the spin frozen vial) (W), Arad the area of the emitting surface (m2 ), F the view factor (-), σ the Stefan-Boltzmann constant (5.67 10−8 W/m2 /K4 )),  the emission coefficient of the emitting surface (-), Trad the temperature of the emitting surface (K), a the absorptivity of the target surface (-) and Tv is the temperature of the target surface (K). In general, a is estimated as the value of  for the target 5

surface [10]. F is defined as the percentage of total radiation which leaves a specific surface (i.e., the IR heater or the chamber walls) and goes directly to the target surface (i.e., the spin frozen vial) [11]. The IR heater is considered to be a diffuse emitter, meaning, this surface emits radiation uniformly in all directions. Consequently, the view factor solely depends on the relative geometric orientation of the emitting surface to the target surface. The view factor Fi,j is given by [11]:

Fi,j =

radiation leaving surface i that goes directly to surface j total radiation leaving surface i

(2)

This definition excludes the possibility of radiation emitted by surface i reflecting of a third surface before finally reaching surface j, i.e., only the direct radiation is considered. The view factor will always be a value between 0 (i.e., two surfaces that do not see each other) and 1 (i.e., two surfaces facing each other with a very small gap in between) [11]. The determination of the view factor F is essential for the calculation of the total amount of energy that is transferred from one body to another. Due to the fact that the view factor is not obvious in most cases, it must be calculated for the specific application. Therefore, mathematical integration over the respective surfaces is required [11]: Z

Z

Ai Fi,j = Aj Fj,i = Aj

Ai

cos(θi ) cos(θj ) dAi dAj πr2

(3)

with Ai and Aj the respective surfaces, r the distance between differential areas dAi and dAj on the surfaces, θi the angle between the normal to dAi and a vector connecting areas dAi and dAj and θj the angle between the normal to dAj and a vector connecting areas dAi and dAj . For the configuration under study, i.e., a cylindrical vial in front of a planar radiator, analytical solutions are available which provide a rough estimate of F by assuming a configuration which approximates the actual set-up [12–14]. These analytical solutions are developed for configurations where the plane and the cylinder have either infinite or the same length. In addition, the cylinder is assumed to be at a fixed position. Most of the analytical methods assume that the cylinder is symmetrically placed towards the plane, while for one other solution the centre of the cylinder is aligned with the edge of the plane. As in the configuration under study the vial is smaller compared to the radiator, these analytical solutions will strongly deviate from the actual value for F leading to an under- or overestimation of the energy transfer. Additionally, the analytical methods do not allow to calculate F for the vial at different positions relative to the radiating plane. Therefore, alternative methods can be applied for the more accurate determination of F , such as the Monte Carlo-based method, the Cross string method [15], the hemi-cube method [16], the Monte Carlo method combined with the finite element technique [16], 6

the finite volume method [17], etc. The Monte Carlo method, which is used in this study, is based on the statistical characteristics of physical processes, or analogous models that mimic physical processes [15]. In optical engineering, ray tracing is the most commonly applied method for determining the impact of light radiation interacting with lenses and mirrors. The rays are considered to follow a straight path between optical obstacles, which can easily be modelled. Although other methods exists to solve the double integral, this leads to more complexity without necessarily more accurate results. Using ray tracing in the fashion of Monte Carlo simulations, leads to fast results, with sufficient means to assess the accuracy. By performing statistical sampling, virtual experiments are performed (i.e., application of ray tracing on a computer) and the solution is approached. This method has been demonstrated to be successful for complex geometries [15, 18]. However, the computational load and the statistical fluctuation are disadvantages. The Monte Carlo method for the determination of the view factor is described by Brenner [11, 19]. In order to assess the impact of the amount of radiant heat which is emitted, first the spectral information must be determined based on Planck’s law. The spectral radiance Bλ (W/sr/m3 )), in function of the wavelength λ (m) and the absolute temperature T (K), is based on the formula of Planck [10]:

Bλ (λ, T ) =

2hc2 1 λ5 e λkhc BT − 1

(4)

with h the Planck constant (6.63 x 10−34 J s), c the speed of light (3.00 x 108 m/s) and kB the Boltzmann constant (1.38 x 10−23 J/K). Knowing the temperature of the radiator, the spectral information of the radiant heat can be determined. Once the transmittance TB (λ) (%) is known, the effective spectral radiance Bλ,ef f (W/sr/m3 ) can be calculated as TB Bλ The influence of the transmission of the different materials on the path of the emitted IR radiation should be included in the model (Figure 3). The radiator filament is encapsulated in fused silica, for which the transmitted portion of radiant heat needs to be assessed. The spectral transmittance of the borosilicate glass vial must be considered as well. However, all radiant heat reaching the vial contributes to the actual heat transfer to the (frozen) product, either through direct radiation (the portion that is transmitted), or through conductive heat transfer from the heated glass (the portion that is absorbed). In the current prototype layout, the radiator is placed directly inside the vacuum chamber [4]. In future systems, it is anticipated that the radiator will be placed outside the vacuum unit in order to achieve a system that can be adapted to CIP and SIP procedures. For that matter, it is necessary to include the assessment of various window materials that can be applied for optical and vacuum applications. Exemplary materials are CaF2 ,

7

silicon, germanium, ZnSe, etc. The effectiveness of these materials for this application is also assessed in this study. The transmission spectrum for several materials is often given by the supplier for a certain thickness of the material. To determine the specific transmittance of the material, Lambert-Beer’s law has to be applied: I = 10− c d I0

(5)

with I and I0 the intensity of the transmitted and incident light beam, respectively (W/m2 ),  the extinction coefficient (m3 /mol/m), c the molar concentration of absorbing medium (mol/m3 ) and d the path length (m). Given a certain transmittance for a given material thickness, the transmittance for other values of thickness can be derived. It has to be noted that in this study it is assumed that the extinction coefficient  is a constant, i.e., independent from the wavelength. To derive a mechanistic model from temperature of the radiator filament till effective transmitted IR radiation to the substance inside the vial, all the above aspects have to be included.

2. Objectives This contribution will treat the following steps: (1) determine the view factor of various radiator-vial situations, (2) determine the influence of the radiation profile of the radiator, and (3) model the absorption characteristics of different media. A Global Sensitivity Analysis (GSA) on the view factor is included to support the necessary input to mechanical tolerance calculations to assure a robust process. The results of this modelling work are then summarized, leading to conclusions about the application of the modelling framework including a proposed control scenario for the process.

3. Materials and methods 3.1. Model to calculate the view factor for a flat radiator The computation of the view factor F for this specific situation where the IR heater and the chamber walls are treated as a flat plate and the spin frozen vial is considered a cylinder, is based on the Monte Carlo method described by Brenner [19]. The Monte Carlo method is a simulation approach in which a defined amount of rays N is propagated from random positions on the radiator surface at randomly chosen angles [11]. For each generated ray, it is evaluated whether it will directly hit the target surface or not. F is estimated by the ratio of the number of rays that strike the target surface, M , to the total amount of emitted rays, N . 8

First, the principles to determine the view factor for a plane surface as radiator with a cylindrical target are presented. For this configuration, it is assumed that no radiant heat reaches the internal surface of the vial. There are three situations for a parallel orientation of the cylindrical target to the radiator plane, illustrated in Figure 4. On the left-hand side, the situation with a cylindrical target smaller than the radiator is given, and therefore, the width of the radiator w is larger than the diameter of the cylinder d. On the right-hand side, the other extreme case is illustrated, where w < d. A clear distinction between these situations is important during the computation of the view factor. For each configuration, the basic situation is described where the vial is positioned in the centre of the radiating plane. However, by including the shift parameters, i.e., xshif t and zshif t , other configurations where the vial is not in the centre of the radiator can also be assessed. These configurations become relevant when the radiation contribution to neighbouring vials needs to be quantified in case of a full series of radiators or during the movement of vials in front of one or multiple radiators. w>d

w=d

k

Front view

y

x

Top view

h

h

k

d

d

w

w

d

x

k

h d w

d

d

s

s

s z

w
w

w

w

Figure 4: The three situations for a cylindrical target parallel to the plane of the radiator (Target in black)

Two characteristics for each generated ray need to be defined: the starting position from the radiator surface and the direction. The starting position of the generated ray is determined by selecting a random number between 0 and 1 for both dimensions of the surface, which is then multiplied by the length h and the width w of the surface, respectively. Therefore, the origin of each generated ray is described by two coordinates. To select a random direction of the generated ray, two angles are defined: the polar angle θ and the azimuthal angle φ (Figure 5). For each angle, a random number between 0 and 1 is generated, pθ and pφ [19]. The

9

corresponding angles are determined using [19]: √ θ = arcsin pθ

(6)

φ = 2πpφ

(7)

Figure 5: The direction of the generated ray to determine the view factor

The length of the generated ray L is gradually increased by the sequential addition of increments, dL, until it hits the target or till it is impossible to hit the target from the end-position of the ray. The length of the generated ray is determined by [19]: xray = x + L cos φ sin θ

(8)

yray = y + L sin φ sin θ

(9)

zray = L cos θ

(10)

for each dimension of the coordinate system. Here, x, y and z are the starting coordinates and xray , yray and zray are the final coordinates of the generated ray. The length of the ray is increased by adding an increment dL. The length of dL is dependent on the orientation of the cylinder with respect to the radiator. Two cases are distinguished, i.e., the cylinder is oriented parallel or perpendicular to the plane (Figure 6). For the parallel case, dL is determined as the diameter of the cylinder d divided by 10 [19]. In some cases it is recommended to divide d by 100 or more. The value of this parameter is further denoted as αdL . The initial length of the ray L is determined as the distance between the surface and the border of the cylinder minus the length of dL: s − d/2 − dL [19]. In the case where the axis of the cylinder is perpendicular to the 10

plane of the radiator, the initial length of the ray is equal to s − k/2 − dL with dL equal to k divided by 10 or more.

Perpendicular

Parallel

d k

Front view

k

h

s y

d y

w

x

d

Top view

w

x

d

h

s z

z

w

x

x

w

Figure 6: The orientation of the radiator plane versus the cylindrical target: parallel (Left) and perpendicular (Right). The variables which are essential for the model to calculate the view factor, i.e., k, d, h, s and w, are indicated in the Figure (Target in black)

The end-position of the ray should be within the position of the target in order to have a hit, therefore, the following constraints are verified in case of a parallel orientation of the cylinder to the plane:   

h−k 2

≤ yray ≤

h−k 2

+k (11)

  p(xray − xshif t − w/2)2 + (zray − zshif t − s)2 ≤ d/2 with h and w the length and width of the flat radiator plane, respectively, s the distance between the centre of the cylinder and the radiator plane, d the diameter of the cylinder and k the height of the cylinder. If the end-position of the ray is not situated within the cylinder, the length of the ray is incremented. In case it is impossible for the ray to ever hit the target, independent from the number of increments which is added, a next ray is generated, but the counter which tracks the number of hits is not incremented. Therefore, the following equations are evaluated, depending on the size of the target with respect to the size of the

11

radiator:

   0 ≤ xray − xshif t ≤ w    w ≥ d =⇒ 0 ≤ yray ≤ h      0 ≤ zray − zshif t ≤ s + d/2

(12)

   −(d − w)/2 ≤ xray − xshif t ≤ w + (d − w)/2    w < d =⇒ 0 ≤ yray ≤ h      0 ≤ zray − zshif t ≤ s + d/2

(13)

If the end-point of the ray exceeds the box, defined by equation 12 or 13, the ray is not categorized as a hit. As long as the endpoint of the ray is within this box, without fulfilling the constraints defined by equation 11, an increment is added to the ray. For the perpendicular case, the conditions to have a hit are given by:    s − k/2 ≤ yray ≤ s + k/2   p(xray − xshif t − w/2)2 + (zray − zshif t − h/2)2 ≤ d/2

(14)

The conditions for a further increase of the length of the ray are:    0 ≤ xray − xshif t ≤ w    w ≥ d =⇒ 0 ≤ yray ≤ s + k/2      0 ≤ zray − zshif t ≤ h

(15)

   −(d − w)/2 ≤ xray − xshif t ≤ w + (d − w)/2    w < d =⇒ 0 ≤ yray ≤ s + k/2      0 ≤ zray − zshif t ≤ h

(16)

3.2. Specific solution for a single-tube radiator In one of the engineered solutions for the continuous freeze-drying prototype, the radiators must be treated as line-sources of finite dimensions (Figure 7). In terms of the view factor, this is similar to the condition as identified in equation 13. Additionally, the radiation profile of the individual radiator must be taken into account. This profile, as indicated in Figure 8, has been supplied by the manufacturer of the radiator tube with the backside covered with a gold reflector. The angle-dependent efficiency was determined at various angular positions, which was fitted by an empirical function. With this information, each ray as 12

simulated in the Monte Carlo method was corrected by this weighting function, after normalizing. In this way an effective view factor could be determined, which is used in the subsequent simulations.

w h

Figure 7: Point source versus plane radiator (Target in black)

Figure 8: The radiation profile for a specific radiatior

3.3. Definition of parameters The parameters to calculate the view factor for the specific experimental set-up are listed in table 1. The first five parameters identify the stationary geometric situation, where the spin frozen vial is centred in front of the IR heater. In some cases, the vial is not perfectly centred in front of the radiation source, e.g., the position of the vial with respect to the walls of the drying chamber or, for future use of the model, in cases where the vials move with respect to the radiators. Therefore, the shift parameters xshif t and zshif t are included in the calculations. Here, xshif t is defined as the shift of (the centre of) the cylinder in the x-direction relative to the central position in front of the plane surface. Here, zshif t is defined as the shift of (the centre of) the cylinder in the z-direction relative to the central position in front of the plane surface. Furthermore, these parameters are also used to calculate the potential interference between two neighbouring radiators. For the current configuration, xshif t and zshif t are both equal to 0 as the spin frozen vials is assumed to be perfectly centred in front of the IR radiator. 13

Table 1: Position of the vial in the freeze-dryer

Parameter Width of the flat radiator for a large radiator (x-direction) Width of the flat radiator for a small radiator (x-direction) Length of the flat radiator (y-direction) Diameter of the cylinder Length of the cylinder Distance between center of the cylinder and the flat radiator Position of cylinder shifted in x-direction Position of cylinder shifted in z-direction

wL wS h d k s xshif t zshif t

Numerical value 45 mm 5 mm 45 mm 25 mm 30 mm 40 mm 0 mm 0 mm

3.4. Variance-based GSA The variance-based GSA technique, using the design proposed by Saltelli et al. [20], has been conducted to examine the influence of the different parameters involved in the model on the view factor [21]. For each parameter, two indices are calculated for each factor, i.e., the first order index Si and the total order indices ST i . Si is the actual fraction of variance accounted for by each factor, and ST i represents the total effect of factor i, i.e., the sum of the first order effect and all interactions with other factors [22]. The GSA has been done for a sample size of 6,000 (Sobol sampling), i.e., 6,000 combinations of the parameters have been investigated. The ranges for the parameter values are chosen between 95 and 105% of the nominal value (Table 1), except for αdL which is varied between 10 and 1,000 (the nominal value is set to 10). The confidence in the value of the length-based parameters is high, and, therefore, a low value for the uncertainty level can be chosen. 3.5. Experimental model verification The presented model framework was experimentally verified by comparing the computed IR energy transfer with the amount of sublimated ice during the primary drying step of spin frozen vials containing pure water. This verification was conducted for several energy input levels. Based on calibration data provided by the manufacturer, the link between the electric power input and the temperature of the four filaments of the IR heater was established. Using Planck’s law (equation 4), the spectral radiance Bλ (λ, T ) was calculated for each filament. Here, the absorptivity of the fused silica of the radiator and the borosilicate of the vial (thickness of 1 mm) was taken into account. During the experimental verification, the original proof-of-concept set-up was used, i.e., no IR window was positioned between the IR heater and the spin frozen vial [4]. Based on equation 1, the energy transfer to the spin frozen vial can be calculated [9]. Here, Arad was 0.0001414 m2 , while the emissivity of the filament was 0.53. The IR energy transfer was quantified by measuring the amount of sublimated ice during primary drying for 14

several energy input levels. A 10 mL type I glass vial (Schott, M¨ ullheim, Germany) filled with 3.9 mL of pure water was spin frozen according to the procedure described by Van Bockstal et al. [4, 9]. The spin frozen vial was hung at a distance of 40 mm from the IR heater, without any contact with the shelf, in the drying chamber of an Amsco FINN-AQUA GT4 freeze-dryer (GEA, K¨oln, Germany). The vial was continuously rotating at 5 rpm to achieve a homogeneous energy transfer from the IR heater to the vial. Immediately after hanging the vial inside the drying chamber, the pressure was lowered to 13.3 Pa. Within 5 minutes, the pressure in the drying chamber was below the triple point of water, hence avoiding any melting of ice. After exactly 17 minutes, the time necessary to achieve the desired pressure level, the IR heaters were activated and primary drying was conducted for 20 minutes at a specific IR heater filament temperature, before sublimation was interrupted. The amount of sublimated ice was determined gravimetrically, including a correction for the energy contributed by the surroundings and the amount of ice sublimated before activation of the IR heaters [9]. The provided energy to the spin frozen vial for each filament temperature was calculated via:

P =

∆m ∆Hs ∆t M

(17)

where P is the power provided by the IR heater to the spin frozen vial (W), ∆m the amount of sublimated ice (kg), ∆t the primary drying time interval (1200 s), ∆Hs is the latent heat of ice sublimation (51139 J/mol) and M the molecular weight of water (0.018 kg/mol). The system was assumed to be at steady-state, i.e., the radiation energy was only used for ice sublimation.

4. Results and discussion 4.1. View factor for a flat radiator with a cylindrical target To investigate the influence of the number of rays N leaving the radiator plane, the confidence interval on the view factor is determined for several values of N by repeating the calculation 100 times. In Figure 9, the intervals between the mean +/− the standard deviation are outlined for an N -value ranging from 102 to 109 . The width of the confidence interval is decreasing with the increase of the N -value. Based on these results, in further analysis, a value of 106 was used for N to compute the view factor for different situations.

The view factor F for a large flat radiator (wL = 45 mm) with a cylindrical target (i.e., the glass vial) parallel to the plane of the radiator is equal to 0.1248, given the parameter settings as defined in table 1. As a comparison, the analytical solution for a cylinder and plane of infinite length results in a value of 0.163 15

0.132

0.13

F (-)

0.128

0.126

0.124

0.122

0.12

0.118 10 0

10 2

10 4

10 6

10 8

10 10

N (-)

Figure 9: Determination of optimal number of samples for Monte Carlo simulation

for the same parameter settings, indicating the inaccuracy of the analytical method to estimate F for the configuration under study [13]. The results for the different planes, i.e., the x-y-, x-z- and y-z-plane, are presented in the upper part of Figure 10 for this set-up. It is clear that in the x-y-plane the rays which hit the target (black dots) are within the target boundaries (red frame), while the rays which miss the target (gray dots) are located outside the cylinder. The reason why the endpoints of some rays which missed the target are located within the target boundaries is due to the fact that the third dimension is not presented in these two-dimensional representations. The results for the same target in comparison with a small radiator (wL = 5 mm) for the parallel case are given in the lower part of Figure 10. For a small radiator, the view factor is slightly higher compared to a large radiator, i.e., 0.1419 versus 0.1248. For the perpendicular case, the results for the large radiator are presented in Figure 11. The view factor equals 0.1212 and 0.1458 respectively for a large and small radiator (results not shown).

4.2. GSA results for a flat radiator with a cylindrical target By using a GSA, the influence of the different parameters on the view factor was determined. The results for the parallel case with a large radiator are presented in Figure 12. The first order indices Si correspond with the sensitivity of the parameter without taking the interaction with other parameters into account, whereas the total order indices ST i represent the sensitivity by taking these interactions into account. The negative value for the first order index of h is due to numerical issues, because, in theory, all values should be positive. The distance between the radiator and the target s has the largest influence on the view factor. Therefore, it is recommended to know the exact value of s to ensure the radiation energy is not higher or lower than expected, which should be taken into account when engineering the system. If the energy provided 16

Figure 10: The rays which hit (black) and miss (gray) the cylindrical target (red) for a large (Upper) and shallow (Lower) radiator (blue) (parallel case) in different planes

Figure 11: The rays which hit (black) and miss (gray) the cylindrical target (red) for a large radiator (blue) (perpendicular case) in different planes

by the radiator to the spin frozen vial is higher than expected, the product which is radiated can receive too much energy, leading to an increased product temperature and, eventually, loss of cake structure (i.e., cake collapse) [9]. The effect of the estimated deviation on s on the response, i.e., radiation energy transfer to the spin frozen vial, can be quantified by performing an Uncertainty Analysis (UA) on the mechanistic model, similar to batch freeze-drying [23, 24]. The diameter of the target d also has a large influence on the number of rays which will hit the target. The calculated value of the view factor is dependent on the value of αdL . However, it can be concluded that the influence of αdL is only limited. k, the length of the target, is involved in interaction, as the total order index is obviously higher than the first order index. The results for a small radiator in a parallel orientation are identical. For the perpendicular case, the distance between the radiator and the target s is the most important

17

0.5

s

0.6

s

0.55

0.4 0.5

d

0.45

0.3

d

ST

S

0.4

0.2

0.3

w

0.1

k

0.25

,dL

k

0

-0.1

0.35

1

2

3

4

5

0.2

h

w

0.15

6

0.1

1

2

3

4

h 5

dL 6

Figure 12: Parameter significance ranking for the parameters involved in the calculation of the view factor (parallel case, large radiator): First order indices S (Left) and total order indices ST (Right)

parameter, followed by the diameter and the length of the target (results not shown). Again, the influence of αdL is limited. 4.3. Results for single-tube radiator Based on the radiation profile, a weighting factor wrad (-) for the view factor is determined (Figure 8). Therefore, each ray hitting the cylinder is multiplied by the factor related to the emission angle. An empirical curve is fitted to the radiation profile. Here, a multi-exponential relation is used. In Figure 13, the predicted weighting factor in function of the angle is presented with the measured radiation profile. The weighting factor in function of the angle is given by:

wrad = −1.24 e0.0157θ + 2.54 e0.00841θ

(18)

The effect on the view factor is relatively low for the radiation profile under study. When applying a radiator with a smaller diameter than the glass vial, most of the radii hitting the vial will originate at small angles. Therefore, the shape of the radiator profile will not have a large impact in this case. For the parallel case with a small radiator, the view factor with and without taking the radiation profile into account is calculated, respectively given by F and Fef f . It can be concluded that with the single-tube radiator the influence of the radiation profile is rapidly decreasing with distance (Figure 14). This is due to the fact that the hitting rays are mostly originating from an angle with maximum intensity.

18

1.4

1.2

w (-)

1

0.8

0.6

0.4

0.2 0

10

20

30

40

50

60

70

80

90

3 (degrees)

Figure 13: Curve fitting for the radiation profile

0.45

6 F Feff

0.4

5 0.35

(F-F eff )/F eff (%)

F (-)

0.3 0.25 0.2 0.15

4

3

2

0.1 1 0.05 0 0.02

0.04

0.06

0.08

0.1

s (m)

0 0.02

0.04

0.06

0.08

0.1

s (m)

Figure 14: Influence of the radiation profile on the view factor in function of the distance s

4.4. Transmission through different media In this section, the resulting calculations are shown for a representative situation during continuous freeze-drying. Based upon experimental experience, the anticipated temperature of the filament will not exceed 1200 ◦ C. Therefore, the method is illustrated using temperature settings equal to this value and below. The spectra of the emitted radiation are compared to the transmitted spectra. The spectral radiance Bλ,ef f for silicon in function of the wavelength λ is given in the left part of Figure 15 for several temperature values. Both the temperature and the wavelength have an influence on the spectral radiance Bλ,ef f . To determine the transmitted power, the contribution over all wavelengths is integrated 19

(presented in the right part of Figure 15). The influence of the thickness of the material is obvious (Right part of Figure 15). For a thickness of 1 mm approximately 85% of the radiance is transmitted (and 15% is absorbed by silicon), whereas only 20% of the radiance is transmitted for a thickness of 10 mm. The temperature has a limited impact on the ratio Bλ,ef f /Bλ for the temperature range under study. However, when the thickness of the material increases, the relative impact of the temperature on the transmitted radiation is higher. For various other materials, the influence of the material thickness is presented in Figure 16. Based on these results, for the application of windows for vacuum systems (10 mm thickness), CaF2 seems the preferred material as the transmission of IR radiation emitted by the radiator filament is the highest for the temperature range under study. Special attention must be given to borosilicate as the material for the vial. Since the thickness of this glass is only 1 mm, a significant amount of heat is directly impinging on the material to be freeze-dried. #109

0.9 100°C 200°C 300°C 400°C 500°C 600°C 700°C 800°C 900°C 1000°C 1100°C 1200°C

10 8

B

6,eff

(W/sr/m3 )

12

1mm 2mm 3mm 4mm 5mm 6mm 7mm 8mm 9mm 10mm

0.8 0.7 0.6

eff

14

B /B (-)

16

0.5

6

0.4

4

0.3

2

0.2

0

0.1 2

4

6

8

10

6 (m)

12

14 #10-6

200

400

600

800

1000

1200

T (°C)

Figure 15: The spectral radiance in function of the wavelength for several values of the temperature for silicon and a thickness of 5 mm (Left) and the influence of the thickness of the silicon on the ratio Bλ,ef f /Bλ (Right)

4.5. Experimental model verification The modelling framework was verified by comparing the computed IR energy transfer with experimental data. A vial filled with pure water was spin frozen and submitted under vacuum to a variation of settings of the current through the filament [9]. The amount of ice sublimed is a measure of the effective energy transfer to the vial. By using the data from the manufacturer of the radiator, a link between the electrical power settings to the filament temperature has been found. Therefore, the filament temperature could be correlated to the effective sublimation heat, which could be compared to the results of the model. Using Planck’s law the spectral radiance Bλ (λ, T ) of the radiator under study, which consists of 4 filaments, is calculated. The 20

Germanium

1

1

0.5

0.5

0

0 200 400 600 800 1000 1200

200 400 600 800 1000 1200

CaF 2

1

1

0.5

ZnSe

0.5

eff

B /B (-)

Fused silica

0

0 200 400 600 800 1000 1200

200 400 600 800 1000 1200

Borosilicate

1

1mm 2mm 3mm 4mm 5mm

0.5 0 200 400 600 800 1000 1200

6mm 7mm 8mm 9mm 10mm

T (°C) Figure 16: The influence of the thickness of the material on the ratio Bλ,ef f /Bλ for different materials

absorptivity of the fused quartz of the radiator and the borosilicate of the vial (thickness of 1 mm) is taken into account. Moreover, the view factor from the 4 filament radiators to the vial is determined as 0.0368. Therefore, based on equation 1, the energy supplied by the radiator can be calculated. The output of the model was compared with the experimentally determined power. Based on the mass sublimated per time unit (i.e., the sublimation rate), the energy supply is calculated. Both the experimental and the calculated energy transfer, in function of the filament temperature, are shown in Figure 17. A good agreement between both is found. 4.6. Practical use of the model Based on the results it can be concluded that the different steps in heat transfer using radiant heat can be modelled. Moreover, the geometric and engineering parameters can be determined and using a GSA the most influential factors can be determined. The spectral situation may lead to a certain complexity. The aim of this study was to create a framework to study and implement process control for sublimation and desorption in pharmaceutical freeze-drying. Based upon (physicochemical) formulation characteristics, it should be possible to determine the necessary settings of the radiators depending on the phase in the process [9]. The required heat power is derived from 21

5 4.5 4 3.5

P (W)

3 2.5 2 1.5 1 0.5 0 600

700

800

900

1000

1100

1200

T filament (°C)

Figure 17: Comparison of experimentally determined power (dots) and calculated power (line)

the latent heat that has to be supplied to the physical transition from ice to vapour. Using Wien’s law, i.e., to determine the maximum of the spectral radiance in function of the temperature, it is then possible to define the required temperature of the radiator filament (omitting absorption for the moment). This is the starting point for determining the required setting of the radiator. The procedure can be illustrated by a flow chart (Figure 18). In the flow chart, an application of the method is illustrated. The starting point is the input information from the freeze-drying process characteristics, leading to the required thermal energy that has to be supplied in a certain phase of the process. Using the Stefan-Boltzmann law, this translates to the initial estimate for the required temperature for the filament in the radiator. The spectrum of such radiation is determined using Planck’s law and using the properties of window-materials the effective heat reaching the vial is determined. After comparison with the required heat, an iterative cycle is performed until the set specification is reached. This iterative process is executed each time an adjustment regarding the freeze-drying process is required.

5. Conclusion A framework was developed for the modelling and control of the continuous freeze-drying process of spin frozen vials. Different elements of the modelling framework have been assessed in this manuscript. The determination of the view factor is based on Monte Carlo simulation. The view factor has been calculated for different radiator-vial situations, i.e., for both small and large radiators and with the vial parallel or 22

START

Load parameters

Initialisation: Determine required radiant heat & initial radiator temperature (filament)

Determine spectral emission by filament

Determine spectral transmission through windows

Calculate effective radiant heat reaching vial

Update temperature of the radiator

No

Difference less than 1% between required and effective heat

Yes

Set temperature of the radiator

END

Figure 18: Procedure to determine optimal radiator settings during freeze-drying

perpendicular to the radiator plane. For an accurate computational estimation of the result, at least 106 rays should be generated. The GSA showed that the distance between the radiator and the target has the most influence on the view factor. Therefore, the position of the vial with respect to the radiator is essential to ensure the radiant heat is not exceeding the critical value leading to structural loss of the dried cake. The radiation profile, specific for a certain radiator, is taken into account by introducing a weighting factor. The transmission of IR radiation through various IR window materials is evaluated and a comparison is made. Lastly, the modelling framework is verified with experimental data, indicating a good agreement between the computed energy transfer and the experimentally determined value. The suggested framework allows to provide an accurate calculation of the radiant heat during freeze-drying. Moreover, by explicitly incorporating the case of positioning the radiator outside of the drying chamber, it is also feasible to achieve a configuration which is acceptable to the pharmaceutical industry. The framework assists in establishing the desired engineering parameters and enables the engineer to assess mechanical tolerances and material options. Hence, the framework fits in the model-based design approach which is 23

applied for the development of the continuous freeze-drying prototype. The model helps to reduce the number of time-consuming experiments for the assessment of different elements of this prototype. Future studies will focus on the experimental application of the framework and, more specifically, the presented procedure for the determination of the optimal dynamic radiator settings during continuous freeze-drying of spin frozen vials.

Acknowledgment Financial support for this research from the Fund for Scientific Research Flanders (FWO Flanders Ph.D. fellowship S´everine Th´er`ese F.C. Mortier) is gratefully acknowledged.

References [1] M. J. Pikal, Freeze drying, in: Encyclopedia of pharmaceutical technology, Marcel Dekker, New York, 2002, pp. 1299–1326. [2] T. Jennings, Lyophilization: Introduction and basic principles, Informa healthcare, New York, 2008. [3] L. De Meyer, P. J. Van Bockstal, J. Corver, C. Vervaet, J. P. Remon, T. De Beer, Evaluation of spin freezing versus conventional freezing as part of a continuous pharmaceutical freeze-drying concept for unit doses, Int. J. Pharm. 496 (2015) 75–85. [4] P.-J. Van Bockstal, L. De Meyer, J. Corver, C. Vervaet, T. De Beer, Non-contact infrared-mediated heat transfer during continuous freeze-drying of unit doses, J. Pharm. Sci. 106 (2017) 71–82. [5] Food and Drug Administration, International Conference on Harmonisation Q8 (2004). URL http://www.fda.gov/downloads/Drugs/Guidances/ucm073507.pdf [6] Food and Drug Administration, Guidance for Industry PAT: A Framework for Innovative Pharmaceutical Development, Manufacuring, and Quality Assurance (2004). URL https://www.fda.gov/downloads/Drugs/GuidanceComplianceRegulatoryInformation/Guidances/ucm070305.pdf [7] J. Corver, Method and system for freeze-drying injectable compositions, in particular pharmaceutical. WO2013036107 (2013). [8] B. G. Ramsay, Vacuum chamber load lock structure and article transport mechanism. US 6609877 B1 (2003). [9] P.-J. Van Bockstal, S. Mortier, L. De Meyer, J. Corver, C. Vervaet, I. Nopens, T. De Beer, Mechanistic modelling of infrared mediated energy transfer during the primary drying step of continuous freeze-drying process, Eur. J. Pharm. Biopharm. 114 (2017) 11–21. [10] R. Bird, W. Stewart, E. Lightfoot, Transport phenoma, John Wiley & Sons, New York, 2006. [11] G. F. Nellis, S. A. Klein, Radiation, in: Heat Transfer, Cambridge University Press, Cambridge, 2009. [12] J. A. Wiebelt, S. Y. Ruo, Radiant-interchange configuration factors for finite right circular cylinder to rectangular plane, Int. J. Heat Mass Transfer 6 (1963) 143–146. [13] A. Feingold, K. Gupta, New analytical approach to the evaluation of configuration factors in radiation from spheres and infinitely long cylinders, J. Heat Transfer 92 (1970) 69–76.

24

[14] E. Sparrow, R. Cess, Radiation Heat Transfer, augmented Edition, Hemisphere Publishing Corporation, Washington, D.C., 1978. [15] A. R. Hajji, M. Mirhosseini, A. Saboonchi, A. Moosavi, Different Methods for Calculating a View Factor in Radiative Applications : Strip to In-Plane Parallel Semi-Cylinder, J. Eng. Thermophys. 24 (2015) 169–180. [16] M. R. Vujicic, N. P. Lavery, S. G. R. Brown, View factor calculation using the Monte Carlo method and numerical sensitivity, Commun. Numer. Methods Eng. 22 (2006) 197–203. [17] J. C. Chai, H. S. Lee, S. V. Patankar, Finite Volume Method for Radiation Heat Transfer, J. Thermophys. Heat Transfer 8 (1994) 419–425. [18] J. R. Howell, The Monte Carlo Method in Radiative Heat Transfer, J. Heat Transfer. 120 (1998) 547–560. [19] J. Brenner, Design specifications for wet-bulb aspirator apparatus, Ph.D. thesis, University of Wisconsin-Madison (2010). [20] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola, Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Commun. 181 (2010) 259–270. [21] A. Saltelli, S. Tarantola, F. Campolongo, M. Ratto, Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models, John Wiley & Sons, Chichester, 2004. [22] T. Homma, A. Saltelli, Importance measures in global sensitivity analysis of nonlinear models, Reliab. Eng. Syst. Saf. 52 (1996) 1–17. [23] S. Mortier, P. J. Van Bockstal, J. Corver, I. Nopens, K. Gernaey, T. De Beer, Uncertainty analysis as essential step in the establishment of the dynamic Design Space of primary drying during freeze-drying, Eur. J. Pharm. Biopharm. 103 (2016) 71–83. [24] P.-J. Van Bockstal, S. Mortier, J. Corver, I. Nopens, K. Gernaey, T. De Beer, 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. 121 (2017) 32–41.

25

Graphical Abstract