Digital Twin for Lyophilization by Process Modeling in Manufacturing of Biologics

: Lyophilization stabilizes formulated biologics for storage, transport and application to patients. In process design and operation it is the link between downstream processing and with ﬁnal formulation to ﬁll and ﬁnish. Recent activities in Quality by Design (QbD) have resulted in approaches by regulatory authorities and the need to include Process Analytical Technology (PAT) tools. An approach is outlined to validate a predictive physical-chemical (rigorous) lyophilization process model to act quantitatively as a digital twin in order to allow accelerated process design by modeling and to further-on develop autonomous process optimization and control towards real time release testing. Antibody manufacturing is chosen as a typical example for actual biologics needs. Literature is reviewed and the presented procedure is exempliﬁed to quantitatively and consistently validate the physical-chemical process model with aid of an experimental statistical DOE (design of experiments) in pilot scale.


Introduction
Lyophilization is widely used to increase the shelf-life for a variety of substances [1][2][3][4][5][6][7][8]. The process consists of three steps. After freezing, the primary drying step is used to remove the solid ice by sublimation. Secondary drying serves to reduce the amount of water, which is bound to the remaining porous product matrix. Especially pharmaceutical biologics in vials are to a great extent only commercially available due to freeze-drying. Oncology is one of the fastest growing markets for biologics, Table A1 (in Appendix A) shows several recipes of freeze-dried IgG formulations [9]. For these expensive and sensitive products it is important to design a safe and reliable process. During lyophilization the protein is exposed to a variety of stresses like pH changes, freeze-concentration, denaturation and aggregation caused by dehydration [5,10,11]. Other phenomena that can lead to protein denaturation during freeze drying are cold denaturation [12] and denaturation induced by interactions with ice-water interface [13,14]. Therefore an adequate formulation has to be used to protect the protein and ensure product integrity and drug activity [15]. On the one hand different components can be used to optimize the formulation but on the other hand the choice of excipients has an impact on the drying time [10,11]. Various reviews deal with the topic of excipients in lyophilization processes [11,[16][17][18]. Lyophilization development therefore has to take the dependent formulation and process optimization into account to establish a safe process for the lyophilization of proteins [19,20].
At first the protein formulation is frozen. In the freezing step water is converted to ice. Here an annealing step can be incorporated to alter the ice crystal morphology [21]. This step is critical for published [22]. During the next step the so called primary drying the frozen ice is removed by sublimation. The shelf temperature is raised to provide the necessary energy for the sublimation and the chamber pressure is set to a point below the triple point. Primary drying is limited by the collapse temperature at which the frozen solution loses structure and is usually the longest process step. The last step, the secondary drying, removes the remaining water by desorption. Here Tg determines the critical temperature. Since its value is residual moisture dependent, this value is not constant during secondary drying [23]. A qualitative freeze-drying process is exemplified in Figure  1 in the p-T phase diagram of water.
Conventional process development for the determination of shelf temperature and chamber pressure is normally performed as listed below: -determination of collapse temperature , -specification of shelf temperature during primary drying , , -calculation of sublimation pressure/ice pressure , -specification of chamber pressure , -determination of scorch temperature , -specification of shelf temperature during secondary drying , , Figure 1. Freeze drying process shown in p-T phase diagram. Adapted from [24].
Collapse temperature of an amorphous frozen solution is related to the glass transition temperature .
can be determined by Differential Scanning Calorimetry (DSC) and by light transmission freeze-drying microscopy (LT-FDM) [25,26]. The collapse temperature defines the maximum product temperature during primary drying which should not be exceeded to guarantee intact cake morphology. Consequently the shelf temperature has to be set to a value that allows fast drying while maintaining the product temperature under . Here the shelf temperature can be set to higher values than since the endothermic sublimation inside the vial leads to a slow increase of product temperature during lyophilization. Analytical methods to determine the collapse temperature have recently been reviewed [27]. Vapour pressure of ice as a function of product temperature can be read from the phase diagram of the solvent, in most cases water. To provide the necessary pressure difference for sublimation, the chamber pressure during primary drying is set to ~35% of the ice pressure [28]. The shelf temperature during secondary drying , is limited by the temperature, at which the product is thermally damaged. For example, proteins denature at a certain Conventional process development for the determination of shelf temperature and chamber pressure is normally performed as listed below: determination of collapse temperature T collapse , -specification of shelf temperature during primary drying T S,PD , -calculation of sublimation pressure/ice pressure p subl , -specification of chamber pressure p C,PD determination of scorch temperature T scorch , -specification of shelf temperature during secondary drying T S,SD .
Collapse temperature T collapse of an amorphous frozen solution is related to the glass transition temperature T g . T g can be determined by Differential Scanning Calorimetry (DSC) and T collapse by light transmission freeze-drying microscopy (LT-FDM) [25,26]. The collapse temperature defines the maximum product temperature during primary drying which should not be exceeded to guarantee intact cake morphology. Consequently the shelf temperature has to be set to a value that allows fast drying while maintaining the product temperature under T collapse . Here the shelf temperature can be set to higher values than T collapse since the endothermic sublimation inside the vial leads to a slow increase of product temperature during lyophilization. Analytical methods to determine the collapse temperature have recently been reviewed [27]. Vapour pressure of ice as a function of product temperature can be read from the phase diagram of the solvent, in most cases water. To provide the necessary pressure difference for sublimation, the chamber pressure during primary drying is set tõ 35% of the ice pressure [28]. The shelf temperature during secondary drying T S,SD is limited by the temperature [29]. Due to the fact that no ice is left after primary drying, , can be higher than . The chamber pressure during secondary drying has little effect on drying rates [30]. It is well known from literature, that corner vials and vials which face the tray sides or freezedryer walls experience a radial heat input [2,[31][32][33]. This was already tested thoroughly, e.g. by experiments with heated dryer walls or isolation from radial effects [34]. With this in mind it can be concluded that two critical vials can be identified as representatives for the two relevant boundary cases.
The first critical vial receives the lowest heat flow. It needs the longest time to achieve thermal equilibrium and determines the minimum process time both for the primary and secondary drying phase. The second critical vial gets the highest heat input and determines the highest possible shelf temperature to ensure intact cake structure. Here, heat input takes place not only by heat conduction with the shelf but also by thermal radiation. Principle product temperature curves for both cases during freezing, annealing and primary drying are shown in Figure 2. It is of utmost importance to identify the critical vials because they define the maximum shelf temperature or respectively the minimum primary drying time to ensure cake structure and to achieve the aimed residual moisture. Critical vials can be identified by analyzing the stationary product temperature as a function of the vial position of the shelf. With the help of these critical vials, the shelf temperature can be adjusted so that the drying process is completed in the shortest possible time and the quality criteria for the lyophilized product are met.
Long process times, which can reach up to several days, are mostly caused by the slow process of sublimation. Therefore poor process design can lead to an inefficient operation with considerable optimization potential [9,31,35,36]. The process duration makes it to a bottleneck for whole process design. Freeze-drying represents one of the last step in the process chain, the interface between purification and formulation. Process design has to find the compromise between fast drying and product safety. No product losses are allowed during the final manufacturing process. Additionally lyophilization has a low energy efficiency with around 5% [9,31].
The US Food and Drug Administration (FDA) published the Guidance for Industry PAT [37] in order to establish a regulatory framework that ensures and supports efficiency in development and manufacturing of pharmaceutical products. For this purpose process modeling offers a powerful tool to deepen process knowledge and identify critical process parameters but also to ensure process control by adjustment of process parameters to guarantee product safety.
The potential of process modeling for lyophilization was acknowledged by industry. An industry-led consortium called LyoHUB was founded in 2014. To ensure high quality and lower process cost, several approaches are pursued in parallel. These include advanced lyophilization technologies and techniques by research equipment, products and process [9]. The latter includes the T P,max t prim.drying,min It is of utmost importance to identify the critical vials because they define the maximum shelf temperature or respectively the minimum primary drying time to ensure cake structure and to achieve the aimed residual moisture. Critical vials can be identified by analyzing the stationary product temperature as a function of the vial position of the shelf. With the help of these critical vials, the shelf temperature can be adjusted so that the drying process is completed in the shortest possible time and the quality criteria for the lyophilized product are met.
Long process times, which can reach up to several days, are mostly caused by the slow process of sublimation. Therefore poor process design can lead to an inefficient operation with considerable optimization potential [9,31,35,36]. The process duration makes it to a bottleneck for whole process design. Freeze-drying represents one of the last step in the process chain, the interface between purification and formulation. Process design has to find the compromise between fast drying and product safety. No product losses are allowed during the final manufacturing process. Additionally lyophilization has a low energy efficiency with around 5% [9,31].
The US Food and Drug Administration (FDA) published the Guidance for Industry PAT [37] in order to establish a regulatory framework that ensures and supports efficiency in development and manufacturing of pharmaceutical products. For this purpose process modeling offers a powerful tool to deepen process knowledge and identify critical process parameters but also to ensure process control by adjustment of process parameters to guarantee product safety.
The potential of process modeling for lyophilization was acknowledged by industry. An industry-led consortium called LyoHUB was founded in 2014. To ensure high quality and lower process cost, several approaches are pursued in parallel. These include advanced lyophilization technologies and techniques by research equipment, products and process [9]. The latter includes the development and use of process analytical technologies (PAT) as well as process modeling and simulation, which is a major topic of this work. Combined, these two approaches can achieve enhanced process control and automation.
The primary goal of this work is to develop a validated physico-chemical process model for prediction of primary and secondary drying phase of the lyophilization of a protein formulation. In order to act as a digital twin to allow faster process development and further-on advanced process control. Therefore an established protein formulation that has been shown to guarantee product stability has been used. The process model is developed through a coupled heat and mass balance and can describe the time-dependent product temperature and the residual moisture. Novel is the efficient model parameter determination method in laboratory scale and the general validation concept by an experimental DoE in pilot scale. By such a distinct validation the process model is enabled to file decisions within a regulated environment: The simulated results are in good agreement with the experiments concerning accuracy and precision. Significant variables are determined with the help of correlation loading plots.

Fundamentals of Process Modeling for Lyophilization
Process models offer the possibility to strongly improve and speed up process development, if they are validated. Then, they can be implemented in the Quality-by-Design (QbD) approach. QbD is more than DOE application, it is a consistent concept demanded by regulatory authorities in order to enable data-driven constant process improvement where PAT is included as one supporting technology [49][50][51][52][53]. PAT is more than inline analytics, it is a consistent technology approach. This strategy is shown in Figure 3. Any validated process model developed could be utilized for autonomous process control implementation in combination with any sensor concept within the PAT approach. Real time release testing (RTRT) could be the final benefit. First inline studies have been published with Raman (drug conformity) and near infrared spectroscopy (residual moisture) [54]. The red dotted outline marks the model validation workflow. This part is covered in this paper. Model assisted process design in combination with Design of Experiments (DoE) offers a promising approach to reduce the amount of experimental workload. For this purpose the model has to be validated, which means accuracy and precision must be proven by experiments. Both is shown by comparing simulations with experiments. From this data correlation-loading diagrams are obtained which show the influence of each parameter on the set operation point.
The risk assessment helps to identify the critical parameters with the highest impact on the outcome of the process. It can be visualized by the Ishikawa diagram or the Occurrence-Impact diagram, shown in Figure 4. Risk assessment is not a part of the model development, but crucial to identify the parameters to be included in the model. The critical quality attributes (CQA) for a lyophilized substance are product integrity and stability to prevent an immune-response caused by aggregation, drug potency, appearance and reconstitution time. Process parameters that influence the CQAs are called critical process parameters (CPP). CPPs in lyophilization are: autonomous process control implementation in combination with any sensor concept within the PAT approach. Real time release testing (RTRT) could be the final benefit. First inline studies have been published with Raman (drug conformity) and near infrared spectroscopy (residual moisture) [54]. The red dotted outline marks the model validation workflow. This part is covered in this paper. Model assisted process design in combination with Design of Experiments (DoE) offers a promising approach to reduce the amount of experimental workload. For this purpose the model has to be validated, which means accuracy and precision must be proven by experiments. Both is shown by comparing simulations with experiments. From this data correlation-loading diagrams are obtained which show the influence of each parameter on the set operation point. Figure 3. Quality-by-Design based process development strategy [55]. Figure 3. Quality-by-Design based process development strategy [55]. The risk assessment helps to identify the critical parameters with the highest impact on the outcome of the process. It can be visualized by the Ishikawa diagram or the Occurrence-Impact diagram, shown in Figure 4. Risk assessment is not a part of the model development, but crucial to identify the parameters to be included in the model. The critical quality attributes (CQA) for a lyophilized substance are product integrity and stability to prevent an immune-response caused by aggregation, drug potency, appearance and reconstitution time. Process parameters that influence the CQAs are called critical process parameters (CPP). CPPs in lyophilization are: -Freezing: temperature, time, annealing temperature/time -Primary drying: temperature, pressure, time -Secondary drying: temperature, pressure, time The Ishikawa diagram visualizes the possible impacts arranged by categories. The Occurrence-Impact diagram helps to make a quantitative assessment of error scenarios. The main influences during lyophilization based on occurrence are chamber pressure and based on impact shelf temperature.  For freeze-drying in vials, different types of models have been presented. Preceding model approaches in literature have dealt with the subdivision of the balance volume into areas: The upper, dried part, the frozen part at the lower end of the vial and the sublimation front as indicated in Figure  5 left. The challenge is addressed in literature as a two-phase moving boundary problem. The two phases are the frozen solution and gas mixture of already existing atmosphere and sublimated water vapor. Modeling has been developed with the help of the Finite Element Method (FEM). Models using the so-called moving-mesh method divide the balance volume in a dried and a frozen part and calculate the node, where phase change is considered [40,[56][57][58]. For freeze-drying in vials, different types of models have been presented. Preceding model approaches in literature have dealt with the subdivision of the balance volume into areas: The upper, dried part, the frozen part at the lower end of the vial and the sublimation front as indicated in Figure 5 left. The challenge is addressed in literature as a two-phase moving boundary problem. The two phases are the frozen solution and gas mixture of already existing atmosphere and sublimated water vapor. Modeling has been developed with the help of the Finite Element Method (FEM). Models using the so-called moving-mesh method divide the balance volume in a dried and a frozen part and calculate the node, where phase change is considered [40,[56][57][58].
The availability of FEM and the knowledge about the influence of radial effects led to a further development of the models, whereby the model was considered two-dimensional. Assuming radial symmetry, a 3-dimensional description of the vial can be achieved. The new models were not only able to describe the position of the sublimation interface, but also to provide information of the presumed form. A schematic illustration showing the modes of heat transfer is given in Figure 5 right. The availability of FEM and the knowledge about the influence of radial effects led to a further development of the models, whereby the model was considered two-dimensional. Assuming radial symmetry, a 3-dimensional description of the vial can be achieved. The new models were not only able to describe the position of the sublimation interface, but also to provide information of the presumed form. A schematic illustration showing the modes of heat transfer is given in Figure 5 right.
To describe multi-dimensional heat and mass transfer analytically derived equations were fitted with coefficients [38]. FEM requires a lot of computational power, which prolongs calculation and hinders the use of such models in common everyday process control practice [60,61]. Ultimately, none of these models could be validated due to the fact, that spatial distribution of residual moisture and temperature during process cannot be measured.
Newer model approaches describe lyophilization in vials with the diffuse interface model [61]. The previous models for the description of multi-dimensional mass and heat transport are used for the description of advanced freeze-drying methods. Here, the drying of spin-freezed vials by heat radiation or the drying of spray-freezed particles in vials should be mentioned [33,62]. Another focus for model based process design is the modeling of scale-up and batch uniformity [36,[63][64][65][66].
The most common approach to describe primary as well as secondary drying is the sorptionsublimation model, which is used in this process modeling approach. These models are able to calculate the removal of frozen and bound water during lyophilization [67]. The presented case study contains the drying of a solution of 2.5 w-% sucrose in purified water. With the aid of modeling experimental effort and use of product solution can be minimized. Additionally an experimental model parameter determination concept was developed to establish values for these constants.
With regard to develop a model with sufficient accuracy and precision a reproducible validation process has been developed. The workflow shown in Figure 6 is used. To describe multi-dimensional heat and mass transfer analytically derived equations were fitted with coefficients [38]. FEM requires a lot of computational power, which prolongs calculation and hinders the use of such models in common everyday process control practice [60,61]. Ultimately, none of these models could be validated due to the fact, that spatial distribution of residual moisture and temperature during process cannot be measured.
Newer model approaches describe lyophilization in vials with the diffuse interface model [61]. The previous models for the description of multi-dimensional mass and heat transport are used for the description of advanced freeze-drying methods. Here, the drying of spin-freezed vials by heat radiation or the drying of spray-freezed particles in vials should be mentioned [33,62]. Another focus for model based process design is the modeling of scale-up and batch uniformity [36,[63][64][65][66].
The most common approach to describe primary as well as secondary drying is the sorption-sublimation model, which is used in this process modeling approach. These models are able to calculate the removal of frozen and bound water during lyophilization [67]. The presented case study contains the drying of a solution of 2.5 w-% sucrose in purified water. With the aid of modeling experimental effort and use of product solution can be minimized. Additionally an experimental model parameter determination concept was developed to establish values for these constants.
With regard to develop a model with sufficient accuracy and precision a reproducible validation process has been developed. The workflow shown in Figure 6 is used.
The workflow is divided into four steps or criteria, which have to be completed to develop a validated model. The workflow is divided into four steps or criteria, which have to be completed to develop a validated model.

Step 1
With known model tasks, the first step includes model derivation and implementation. Correct model implementation can be verified by calculating characteristic numbers or by comparison to literature data. This shows that the model describes the basic physical relationships in an expectable way.

Step 2
The second step checks the model for sensitivity. It is used to identify the critical parameters, whose variations lead to significantly different results. These parameters are examined in a given design space. This can correspond, for example, to the design space of the QbD approach, in which the model describes the equipment. The sensitivity of the model is checked first by single parameter studies. During these studies, only one parameter at a time is varied while all other parameters are held at their respective mean values. The results of this study not only serve to verify the sensitivity but also provide evidence that the model does not produce obviously false results. Afterwards, multiparameter studies are conducted. To reduce the amount of simulations, Design of Experiments (DoE) can be used. A full factorial DoE design would have too many simulations. Therefore a screening Plackett-Burmann plan is established. The model and process parameters for the simulated DoE can be obtained by different sources (prior knowledge, literature or previous process developments) and have to be in a reasonable range for the given study. At first the expected influence of the parameters on the target value is considered. In this study the target value is the residual moisture of the lyophilized cake.
After planning the DoE is simulated and the results are fed back into a statistical tool for visualization. Pareto plots of the standardized effects can lead to a fast identification of significant parameters. These make it possible to distinguish between the parameters that have a significant influence on the result and those whose variance within the chosen design space is negligible.

Step 1
With known model tasks, the first step includes model derivation and implementation. Correct model implementation can be verified by calculating characteristic numbers or by comparison to literature data. This shows that the model describes the basic physical relationships in an expectable way.

Step 2
The second step checks the model for sensitivity. It is used to identify the critical parameters, whose variations lead to significantly different results. These parameters are examined in a given design space. This can correspond, for example, to the design space of the QbD approach, in which the model describes the equipment. The sensitivity of the model is checked first by single parameter studies. During these studies, only one parameter at a time is varied while all other parameters are held at their respective mean values. The results of this study not only serve to verify the sensitivity but also provide evidence that the model does not produce obviously false results. Afterwards, multi-parameter studies are conducted. To reduce the amount of simulations, Design of Experiments (DoE) can be used. A full factorial DoE design would have too many simulations. Therefore a screening Plackett-Burmann plan is established. The model and process parameters for the simulated DoE can be obtained by different sources (prior knowledge, literature or previous process developments) and have to be in a reasonable range for the given study. At first the expected influence of the parameters on the target value is considered. In this study the target value is the residual moisture of the lyophilized cake.
After planning the DoE is simulated and the results are fed back into a statistical tool for visualization. Pareto plots of the standardized effects can lead to a fast identification of significant parameters. These make it possible to distinguish between the parameters that have a significant influence on the result and those whose variance within the chosen design space is negligible.

Step 3
With the results from Step 2, a model parameter determination concept can be established. Since the significant parameters are determined in the previous step it is known which parameters have to be measured with high accuracy. Furthermore, experimental and modeling results are compared for model validation. The experimental error is obtained by repetition of an operating point. The model also shows errors. This is determined from the known measurement deviations and inaccuracies during model parameter determination. To calculate the model error simulations are conducted where all input variables are hold at average, single parameter studies and variation of all parameters to maximum respectively minimum value. Furthermore the parameters are randomized between these boundaries. This is done by Monte Carlo simulations. The resulting model error should be smaller than the experimental one.

Step 4
The fourth and last step uses the experimental results from a DoE. The experimental input is simulated and the results are fed into a statistical analysis tool. This is used to find the correlations between input variables and target values. The correlations are determined by the Partial Least Square (PLS) method and visualized by Correlation-Loading diagrams. The knowledge of correlating effects helps in process development to specify the stability of design and operation space chosen. After all steps have been processed, the model should be validated for the previously defined design space.
The workflow for model assisted process development was already shown for different unit operations, such as solid-liquid (phyto-) extraction [68], aqueous two-phase extraction of monoclonal antibodies [55], upstream fermentation of monoclonal antibodies [69] and chromatographic separation [49]. Implementing this approach, a validated model is able to reduce experimental workload significantly, which expedites process design [49,70,71].

Process Model Task
The aim of the process model presented in this work is to describe residual moisture both for ice and bound water as well as the simulation of product temperature. The latter serves to ensure product safety regarding the collapse temperature. For this purpose, a model parameter determination concept is to be developed in which the required parameters are collected with only a few experiments. Subsequently, the target parameters are to be calculated by entering known parameters such as vial geometry, solution properties as well as the process parameters shelf temperature and chamber pressure. This enables the prediction of the drying end point for both drying phases and the identification of the optimal operation point.

Process Model Depth
The depth of the process modeling must be chosen by taking several aspects into account:  Figure 7. The developed process model can describe the dynamic behavior of lyophilization processes based on physico-chemical considerations. Therefore it has a higher model depth than steady-state models. Computational fluid dynamics (CFD) models have a higher model depth because they give further insight to gas dynamics during the process. Computational fluid dynamics (CFD) models have a higher model depth because they give further insight to gas dynamics during the process. In this work a one-dimensional, discretized sorption-sublimation model with uniform sublimation front is used. This kind of model provides a physico-chemical description of lyophilization processes. The considered balance volume for modeling of lyophilization processes inside a vial is shown in Figure 8. The approach uses an energy balance for simulation of the product temperature, a mass balance for description of vapor and solid phases during primary as well as a mass balance of the solid phase for desorption during the secondary drying. The freezing step is not considered in this work in a first approach.

Energy Balance
The energy balance consists of the different contributions of the three phases. In the frozen layer heat is transported through conduction. This can be described by the one-dimensional heat Equation: In this work a one-dimensional, discretized sorption-sublimation model with uniform sublimation front is used. This kind of model provides a physico-chemical description of lyophilization processes. The considered balance volume for modeling of lyophilization processes inside a vial is shown in Figure 8. The approach uses an energy balance for simulation of the product temperature, a mass balance for description of vapor and solid phases during primary as well as a mass balance of the solid phase for desorption during the secondary drying. The freezing step is not considered in this work in a first approach. Computational fluid dynamics (CFD) models have a higher model depth because they give further insight to gas dynamics during the process. In this work a one-dimensional, discretized sorption-sublimation model with uniform sublimation front is used. This kind of model provides a physico-chemical description of lyophilization processes. The considered balance volume for modeling of lyophilization processes inside a vial is shown in Figure 8. The approach uses an energy balance for simulation of the product temperature, a mass balance for description of vapor and solid phases during primary as well as a mass balance of the solid phase for desorption during the secondary drying. The freezing step is not considered in this work in a first approach.

Energy Balance
The energy balance consists of the different contributions of the three phases. In the frozen layer heat is transported through conduction. This can be described by the one-dimensional heat Equation: ⋅ , ⋅ = ⋅ (1) Figure 8. Balance volumes-energy and mass.

Energy Balance
The energy balance consists of the different contributions of the three phases. In the frozen layer heat is transported through conduction. This can be described by the one-dimensional heat Equation: with ρ as density, c p as specific heat capacity, T as product temperature and λ as thermal conductivity. The heat is conducted to the sublimation front. Here the frozen water is sublimated and the water vapour flows through the dried matrix above the sublimation interface: .
with . m subl as mass flow by sublimation and ∆h subl as sublimation enthalpy. The convective heat transport is neglected because the contributions in high vacuum are low in relation to the sublimation. The dried layer is described in analogy to the frozen layer by the one-dimensional heat equation.
This three Equations can be summarized by: In this equation the properties density and thermal conductivity are mixed variables for the vial content, depending on the mass fractions of frozen and ice-free product of the respective discrete. These are derived by mass balance. The heat capacity of the product is substituted by the so-called 'apparent heat capacity'c p,apparent to consider the phase change [73]. This ensures that the latent heat in each discrete is taken into account only when the front is positioned there. Here, the heat capacity is calculated by the mass fractions of the solid phase as sensible part and the heat of sublimation as latent part. These kind of calculation is already used in commercial simulation software, e.g., in the COMSOL Multiphysics ® Software. This software has already been used in modeling of lyophilization [74].

Overall Mass Balance for Water
The water mass is described by two fractions: the frozen part and the ice-free part. During primary drying, only the frozen part varies. To calculate this part in the balance volume, both vapor and solid phase of the water have to be described: with m W as overall mass of water, m W,g as mass of water in the vapor phase and m W,s as mass of water in the solid phase. Each mass change is calculated by an own mass balance and combined in this overall balance of the water content.

Mass Balance of Vapor Phase
The mass balance of the vapor phase is influenced by vapor convection through the dried part of the cake and phase change by sublimation. Convection is described by the continuity equation (first term on the right side of the equal sign of Equation (5)) and sublimation by the Hertz-Knudsen formula. The latter is derived from the kinetic gas theory and used to describe phase change on free phase boundary surfaces. Since only the mass of the ice changes during the primary drying, only its mass has to be considered: ρ W,g as density of water vapor, u g as vapor velocity, A vial as cross-sectional area of the vial, M W as molecular weight of water, R as universal gas constant, p subl as sublimation pressure as function of temperature, p f ront as pressure at the sublimation front and A subl as area where phase change takes place. The gas volume flow through a porous medium considering the material properties can be described by the Carman-Kozeny equation. The equation is adapted to describe the unknown gas velocity u g : with ∆p as pressure difference between discretes, η W as dynamic vapor viscosity of water and K as hydraulic flow resistance. K describes the product resistance of the lyophilisate. This parameter is the overall flow resistance which takes all material properties into account that hinder the vapor flow. K can be calculated by Equation (7) [75]: represents the porosity, H the height of the porous medium, C is a constant that is usually taken to be 5 [75] sometimes assumed as tortuosity and SSA stands for the specific surface area. Measurement techniques to obtain values for the porosity and specific surface area are reviewed in [76]. Any non-destructive characterization of such highly porous cakes is quite sophisticated and of high efforts [77]. Besides, all the single parameter values are still quite inaccurate, therefore in practice the overall coefficient K is applied.

Mass Balance of Solid Phase
Change in mass of solid ice only takes place by sublimation, which is described by the Hertz-Knudsen formula. Since the decrease in solid ice is equal to the increase in the vapor phase by sublimation, the pressure difference is the same, however reversed:

Assumptions and Conclusion for Primary Drying Phase
For calculation of the mass balance of ice, several assumptions have to be taken into account. The Hertz-Knudsen formula for describing the sublimation includes p f ront and A subl . The pressure at the sublimation front p f ront increases due to sublimation of ice and decreases by removal of water vapor, described by the continuity equation. The area of phase change A subl equals the cross-sectional area of the vial at the start of the primary drying, but increases during the process due to radial effects [2] This effect leads to an increase in sublimating water.
Both variables are dynamic during primary drying and cannot be calculated. Since sublimation with the same values for pressure difference and area is much faster than convection, latter is the transport phenomena which governs the overall water transport rate. This kind of quasi-stationary principle is used to neglect the sublimation term and therefore eliminate p f ront and A subl from the mass balance. Hence, p f ront can be calculated to be equal to p subl as long as there is still frozen ice present. Therefore the overall mass balance for water can now be described by only the continuity equation. This term describesṁ subl and therefore the energy and mass balance are coupled.

Mass Balance During Secondary Drying
After the solid ice is removed, the bound water in the solid ice-free product decreases on the account of desorption: with w bw as mass share of bound water in the solid product, k bw as desorption constant and w bw,eq as mass share of bound water at equilibrium. Desorption of water from solid materials is strongly temperature dependent. Therefore, an Arrhenius' approach was chosen to describe the rate constant k: with E A as activation energy and T as temperature. The exponent x was included to take the pressure dependence into account. For calculation, the temperature-dependent sublimation enthalpy was chosen as E A , product temperature for T and water activity for x. The sublimation enthalpy was used in a first approach because it is higher than the desorption enthalpy [59] and therefore seemed more reasonable because this would rather lead to an under-than an over-prediction of the desorption constant. The water activity is defined to be calculated as ratio of partial pressure of water in the surrounding system to saturation vapor pressure. The composition of the vapor phase within the product cake was assumed to be nearly completely water, so the pressure was taken to be equal to the partial pressure. The saturation vapor pressure was calculated to be equal to the sublimation pressure. Therefore, the rate constant for desorption of water was calculated as shown in Equation (11): with a W as water activity. The mass share of bound water on sucrose at equilibrium w bw,eq for Equation (10) is adopted from literature [78].
The physico-chemical process model based on coupled heat and mass balance allows to predict the time-dependent product temperature and residual moisture content during primary and secondary drying. In order to solve the model equations the heat flow between the glass vial bottom and the product inside the vial has been equated and the temperature at the vial bottom is taken as shelf temperature. Furthermore, the evolution of the sublimation front during primary drying can be described. The resulting algebraic partial differential equation (PDE) system is solved via orthogonal collocation with Jacobi polynom approximation and Gear integration algorithm within standard flowsheeting software like Aspen Custom Modeller™. The objective is to link lyophilization to all other unit operation as total process simulations [49].
The following chapters describe the material and methods for experimental model parameter determination as well as a consistent procedure for distinct model validation.

Product Mixture and Instruments
For preparation of sucrose solution 25 g/L d(+)-Sucrose (>99.5%, p.a., Carl Roth GmbH + Co. KG, Karlsruhe, Germany) was dissolved in purified water (arium™pro, Sartorius AG, Göttingen, Germany). Both masses were determined by a laboratory scale LC 1200 S (Sartorius AG, Göttingen, Germany). For weighting of the empty, filled and lyophilized vials a precision scale LA 310 S obtained by Sartorius was used. A fill volume of 1 mL is used and 135 vials are loaded onto the middle shelf.

Experimental Runs
The freezing step for all experiments used annealing and was adapted from literature [1]. First, the shelf temperature was lowered to −45 °C and held for 2 h. Subsequently, annealing started and the shelf temperature was raised to −20 °C and held for 1 h, before it was lowered again to −45 °C. After another hold for 2 h, primary drying started with varying values for chamber pressure and shelf temperature based on the DoE shown in Table 1. The values for the secondary drying temperature, chamber pressure and duration are varied as depicted in the DoE. All temperature ramps are set to 1 K/min, which includes the temperature ramps during primary and secondary drying.
The product temperature during the runs is measured by 8 WTMplus sensors. The position of the sensors and the weighed vials is shown in Table A2 (in Appendix A).

Experimental Runs
The freezing step for all experiments used annealing and was adapted from literature [1]. First, the shelf temperature was lowered to −45 • C and held for 2 h. Subsequently, annealing started and the shelf temperature was raised to −20 • C and held for 1 h, before it was lowered again to −45 • C. After another hold for 2 h, primary drying started with varying values for chamber pressure and shelf temperature based on the DoE shown in Table 1. The values for the secondary drying temperature, chamber pressure and duration are varied as depicted in the DoE. All temperature ramps are set to 1 K/min, which includes the temperature ramps during primary and secondary drying. The product temperature during the runs is measured by 8 WTMplus sensors. The position of the sensors and the weighed vials is shown in Table A2 (in Appendix A).

Analytics
To determine residual moisture, the vial is weighted empty, filled with solution and after the process. From the known mass of sucrose, the excess weight can be determined as water.

Software Tools
Both DoEs were generated with JMP (JMP Inc., SAS Institute, Cary, NC, USA). It was also used for statistical evaluation and creating Pareto plots. Correlation-Loading plots as well as PLS regression were done with UNSCRAMBLE (Camo Analytics, Oslo, Norway).

Model Parameter Determination
Parameters such as the overall heat transfer coefficient, product resistance and stationary product temperature must be determined in one previous characterizing experiment. The procedure to obtain these values is described in this chapter in detail. Other necessary input values and their determination method are shown in Figure 10.

Thermal Conductivity of the Vial
The vial heat transfer coefficient is determined by an experiment in analogy to literature which is stopped right before the end of primary drying [59]. The mass of sublimated water is known by weighting, which enables calculation of total heat necessary to sublime. This value ∆ in combination with process duration ∆ gives the heat flow. Before weighing the vials are thawed at room temperature. The temperature values for the measured vials were taken from the closest vial that was probed with a WTMplus. Now can be calculated by heat conduction: with as cross-sectional area of the vial, , as shelf temperature during primary drying and , as average product temperature during primary drying. The latter is measured by wireless temperature measurement sensors. The position of the WTMplus sensors is shown in Table A2 (in Appendix A). To obtain the thermal conductivity of the vial , has to be multiplied by the thickness of the vial bottom. The thermal conductivity of the vial has three main contributions [80]: quantifies the contribution from direct conduction of the shelf to the vial, whereas and

Thermal Conductivity of the Vial
The vial heat transfer coefficient k vial is determined by an experiment in analogy to literature which is stopped right before the end of primary drying [59]. The mass of sublimated water is known by weighting, which enables calculation of total heat necessary to sublime. This value ∆Q in combination with process duration ∆t gives the heat flow. Before weighing the vials are thawed at room temperature. The temperature values for the measured vials were taken from the closest vial that was probed with a WTMplus. Now k vial can be calculated by heat conduction: with A vial as cross-sectional area of the vial, T shel f ,PD as shelf temperature during primary drying and T product,av as average product temperature during primary drying. The latter is measured by wireless temperature measurement sensors. The position of the WTMplus sensors is shown in Table A2 (in Appendix A). To obtain the thermal conductivity of the vial λ vial , k vial has to be multiplied by the thickness of the vial bottom. The thermal conductivity of the vial has three main contributions [80]: (14) k cond quantifies the contribution from direct conduction of the shelf to the vial, whereas k r and k gc describe the contributions of heat transfer by radiation respectively gas conduction [80]. Values of k vial depend on chamber pressure, vial location, container typeand possibly on the used freeze dryer equipment [81].

Product Resistance
Product resistance K is determined by the same experiment as the vial heat transfer coefficient. The calculated mass flow . m PD is inserted in Formula (6), which is extended by the average vapor density ρ w,g,av and the cross-sectional area of the vial A vial : m PD,av = p C,PD − p subl,av η W,av ·K ·ρ W,g,av ·A vial (15) with u W,av as average gas velocity during primary drying, p C,PD as primary drying chamber pressure, p subl,av as sublimation pressure as function of the average product temperature during primary drying and η W,av as average dynamic viscosity. The workflow for the determination of the parameters is shown in Figure 11. The workflow for the determination of the parameters is shown in Figure 11.

Product Temperature
The maximal product temperature is a function of vial position in the freeze-dryer. This effect is well-known and described in literature where it is mainly attributed to radiation effects. This work uses the highest measured product temperature as position-dependent shelf temperature.

Water Properties
The water properties density and dynamic viscosity are calculated as a function of pressure and temperature by UNIFAC (Universal Quasichemical Functional Group Activity Coefficients). is determined by DSC and has been determined by freeze-dry microscopy. The used formulation has a of −34 °C and a of −29 °C [1].

Process Model Development and Validation
As the workflow for process-and model development as depicted in Figure 8 has been applied successfully for other unit operations [55,69,[82][83][84][85] it is now consequently applied for lyophilization as well. Only in Step 1 the model is verified by using literature data sets which have a different sucrose concentration than the one used in the experimental DoE. The model verification shows

Product Temperature
The maximal product temperature is a function of vial position in the freeze-dryer. This effect is well-known and described in literature where it is mainly attributed to radiation effects. This work uses the highest measured product temperature as position-dependent shelf temperature.

Water Properties
The water properties density and dynamic viscosity are calculated as a function of pressure and temperature by UNIFAC (Universal Quasichemical Functional Group Activity Coefficients). T g is determined by DSC and T collapse has been determined by freeze-dry microscopy. The used formulation has a T g of −34 • C and a T collapse of −29 • C [1].

Process Model Development and Validation
As the workflow for process-and model development as depicted in Figure 8 has been applied successfully for other unit operations [55,69,[82][83][84][85] it is now consequently applied for lyophilization as well. Only in Step 1 the model is verified by using literature data sets which have a different sucrose concentration than the one used in the experimental DoE. The model verification shows whether the developed process model is able to adequately describe the underlying physico-chemical effects and their interaction. It is necessary to evaluate if the model is implemented correctly and if the results are expectable. In all the next steps data from own DoE experiments are used and compared with simulations in order to validate the process model.

Step 1
For a general verification of the computerized model, the results are compared to a study from literature. For this, a study from Pikal et al. is chosen [59] since all necessary input parameters for the process model are published. The given input parameters are shown in Table 2. The simulated product temperature is plotted against time and compared to the literature temperature. The results are shown in Figure 12 and show a good conformity.  The simulated product temperature is plotted against time and compared to the literature temperature. The results are shown in Figure 12 and show a good conformity. The published residual moisture during experiments were ~13.8 w-% at the end of primary drying and ~0.8 w-% at the end of secondary drying. While the simulated moisture at the end of primary drying with 17.7 w-% showed good agreement, the content of bound water is calculated to 4.93 w-%. However, because of the good agreement of temperatures and residual moisture after primary drying the model can be considered verified because it can describe the basic mass and heat transfer phenomena in an appropriate way. The published residual moisture during experiments were~13.8 w-% at the end of primary drying and~0.8 w-% at the end of secondary drying. While the simulated moisture at the end of primary drying with 17.7 w-% showed good agreement, the content of bound water is calculated to 4.93 w-%. However, because of the good agreement of temperatures and residual moisture after primary drying the model can be considered verified because it can describe the basic mass and heat transfer phenomena in an appropriate way.

Step 2
Before creating the DoE, all parameters to vary in experiment and model have to be listed. These parameters, their considered design range and expected impact on the target value are shown in Table 3. The freezing protocol is the same for all experiments as described above to ensure the same starting point for drying. The used method is considered to be safe for the lyophilization of IgG solutions since the product integrity has been shown for processes that used shelf temperatures higher than the used design range and the published average product temperatures are not exceeded in the experiments made [1]. The expected effect is estimated by the study designer based on prior knowledge or literature. The last column explains, why the design range is chosen.
Even if freeze-drying is a bottleneck due to its long process duration, cake appearance, product stability and potency are the decisive quality criteria for the lyophilization process step. Especially the residual moisture after lyophilization is a quantitative characteristic that allows conclusions about product stability and potency.
For this work, a Plackett-Burman DoE was conducted. This plan includes 15 different experiments that were conducted inside a pilot freeze dryer. Furthermore 6 experiments were done to obtain the model parameters. With the help of the model parameters simulations were carried out and later compared to the experimental data. The center point is repeated two times for statistical evaluation. The statistical evaluation showed that the requested p-value of less than 0.01 is achieved. The visualization via Pareto plot is shown in Figure 13. Table 3 shows only the main effects. Due to the low resolution, interactions cannot be considered by means of a Plackett-Burman plan. It can be seen that only the shelf temperature and duration of secondary drying have a significant effect. Any temperature and pressure effect of the primary drying parameters are below the significance limit because all operating points were run under successful primary drying and therefore no more solid ice was present in the system that deviates the residual moisture or the appearance of the final cake. Temperature and duration of secondary drying are then significant because the mechanism of desorption is strongly temperature dependent. As expected, residual moisture depends most on energy input by heat conduction and lesser on the drying time.
For this work, a Plackett-Burman DoE was conducted. This plan includes 15 different experiments that were conducted inside a pilot freeze dryer. Furthermore 6 experiments were done to obtain the model parameters. With the help of the model parameters simulations were carried out and later compared to the experimental data. The center point is repeated two times for statistical evaluation. The statistical evaluation showed that the requested p-value of less than 0.01 is achieved. The visualization via Pareto plot is shown in Figure 13.  Table 3 shows only the main effects. Due to the low resolution, interactions cannot be considered by means of a Plackett-Burman plan. It can be seen that only the shelf temperature and duration of secondary drying have a significant effect. Any temperature and pressure effect of the primary drying

Step 3
The third step in the model validation workflow consists of the comparison of experimental and model results. Monte Carlo simulation is used to conduct this multi-parameter study. In this simulation the model parameters are randomly distributed in the given range that has been obtained by experiments. The modeled results should have a smaller deviation than the experimental results in order to obtain a validated model. For this work, the center point of the Plackett-Burrman DoE is chosen because it is already conducted three times. The model input as well as the respective deviations are given in Table 4. The results of product temperature during primary and secondary drying of experiments and simulation runs are summarized in Figure 14. An overall of 76 simulations were conducted. For a better insight, only the average, maximum and minimum values are shown in Figure 14. The other curves of the simulated temperatures lie between these values.
In general, simulations and experiments show good agreement. It can be seen that the simulated temperatures towards the end of secondary drying are below the measured product temperatures. This deviation also explains the higher residual moisture, which decreases with increasing temperature. The experiments showed a residual moisture of 4.75 w-% (+47.15/−42.58), while the simulation runs resulted in 9.21 w-% (+1.74/−2.78). Nevertheless, the results are in good agreement with experiments, and the experimental error is far greater than that of the model. With that in mind, the model has been proven valid. The results of product temperature during primary and secondary drying of experiments and simulation runs are summarized in Figure 14. An overall of 76 simulations were conducted. For a better insight, only the average, maximum and minimum values are shown in Figure 14. The other curves of the simulated temperatures lie between these values.

Step 4
The results of the Monte Carlo simulations are fed into a statistical tool. To find correlations between the input variables, PLS is used. The outcome is visualized in a Correlation-Loading plot, as shown in Figure 15. Here, the predictors are shown in blue, significant ones in black. The red circle marks the boundary, where 50% of variance can be explained by the partial components 1 and 2, which are displayed on both axes. PLS helps to identify so-called "principal components". When two quantities face each other in the diagram, there is a negative correlation between them. If they are at the same point, there is a positive correlation and they change covariantly. If variables are perpendicular to each other, no correlation exists. The further outside a quantity is located in the diagram, the greater the variance, which can be explained by the variance of the opposite quantities. The sum of PC 1 and 2 displays the total variance, which can be explained by the considered parameters. Here, a sum of 80% can be achieved by only two principal components. In summary, the process model for freeze drying is validated for both primary and secondary drying in appropriate accuracy and precision. It can be used for both process control and process development. The significant influencing variables have been identified and a model parameter determination concept has been developed. Model parameters must be re-determined for each type of freeze dryer as well as different compositions of the feed.

Discussion and Conclusions
This paper proposes development and distinct validation of a predictive physico-chemical process model for freeze drying as a digital twin. The first step towards autonomous process operation and control within any PAT concept. To start, the freezing process is assumed to be constant and the primary and secondary drying are described in the model. Furthermore a safe formulation system for the lyophilization of proteins was used. The quality by design approach was considered during the four step quantitative model validation. The first step dealt with the derivation of the model and the assumptions made. The applicability of the model was shown by comparing it with literature data.
Step 2 dealt with the sensitivity of the model. A design space was defined that ensures product integrity. Within these limits it could be shown that the residual moisture does not react sensitively to the process variables of primary drying. As a result it can be stated that in future process developments the focus of primary drying can be placed on the sublimation of the solid ice without exerting an influence on secondary drying.
The required parameters can be determined within one piloting-scale experiment. The work required is thus a period of one drying process with preparation and post-processing and a feed input from a single shelf. In the case investigated, one week of work and 150 mL feed solution were required.
The third step was used to validate the predictive physico-chemical model. Experimental errors were compared to simulation errors using Monte Carlo simulation studies. The error of the simulation was significantly smaller than that observed in the experiment, thus accurate and precise enough. The envelope curve of experiments included the envelope curve of simulations. Thus the model is quantitative successfully validated and can be used for further process development.
The last step was conducted to identify the correlations between operation parameters of influence during freeze-drying. The results support risk assessment in process development and emphasize which variables must be determined with great care in the experimental model parameter determination concept.
The presented approach for the physico-chemical process model is proven data-driven to be successful. The model is verified and distinctly quantitatively validated. Besides the process model as a digital twin, a deeper understanding of lyophilization was gained, which is valuable for further work and process development as well as sensor and process control implementation. In addition, the freezing step should be taken into account as well in future.
Author Contributions: L.S.K. performed experiments, modeling and simulation and wrote this paper. A.J. supported the experimental part, helped in the writing and reviewed this paper. P.K. and F.H. supported with laboratory and piloting equipment as well as expert knowledge regarding operation and objectives. J.S. is responsible for conception and supervision. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors want to thank the Bundesministerium für Wirtschaft und Energie (BMWi) especially M. Gahr (Projektträger FZ Jülich) for funding this scientific work.