Model to predict inhomogeneous protein-sugar distribution in powders prepared by spray drying

A protein can be stabilized by spray drying an aqueous solution of the protein and a sugar, thereby incorporating the protein into a glassy sugar matrix. For optimal stability, the protein should be homogeneously distributed inside the sugar matrix. The aim of this study was to develop a model that can predict the distribution of protein and sugar in an evaporating droplet using bovine serum albumin (BSA) and trehalose as model protein and sugar, respectively. This was achieved by expanding a previously developed model that was able to predict the growth or shrinkage of inhaled droplets in the airways (Grasmeijer, Frijlink & Hinrichs, 2016). The droplet was considered to consist of a finite number of concentric spherical shells in which the change in concentration of components was calculated in time, enabling the prediction of concentration gradients inside the droplet. It was found that during evaporation of the droplet, an inhomogeneous protein– sugar distribution was formed even when surface active properties were not considered. The relatively large protein molecule was predicted to accumulate much faster at the surface of the droplet than the sugar due to slower diffusion, resulting in a lower sugar/ protein ratio at the surface of the particle than in the center. For a mixture of BSA and trehalose, not considering surface active properties, it was predicted that 60% of protein was incorporated in the powder at a lower sugar/protein ratio than when protein and sugar would have been homogeneously distributed, which may hamper efficient protein stabilization. These predictions can be used to more accurately adapt the initial composition of the solution to ensure proper stabilization of the protein, for example by simply increasing the amount of sugar to increase the sugar/protein ratio at the surface. Or, if this is limited by the desired loading, changing the drying conditions to slow down the drying rate. & 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Depending on their application, spray dried products should often meet specific physical properties. These physical properties can be a specific particle size distribution and morphology, whether the particles are crystalline, glassy or rubbery, and presence or absence of concentration gradients inside the particle, in case the starting solution contains more than one solute. A typical field where such information is particularly useful, is protein stabilization, accomplished with the aid of a sugar. For such products, it is important that the protein is fully encapsulated in a sugar matrix which should be in its glassy state. Therefore, understanding and the ability to predict the final product properties may determine whether a product and process can be developed successfully or not.
There are two main stages during the drying of a droplet in a spray dryer that can be distinguished. The first stage begins at the moment that water starts evaporating from the droplet. The evaporation is rapid and depends fully on the relative humidity and temperature of the drying air around the droplet, and the presence of dissolved components in the droplet (assuming a stagnant air layer, otherwise relative air velocity is important as well). During evaporation, the solute concentration in the droplet will increase until a critical threshold is reached after which it starts to solidify. From that moment on, the evaporation rate slows down due to the slower diffusion of moisture through the growing shell of solidified solute at the surface of the droplet, which is the second drying stage.
The solidification of solute at the transition from first to second stage will usually happen at the surface of the droplet due to several reasons. In general, solute molecules will diffuse slower in solution than water molecules. During evaporation, the surface gradually recedes towards the center of the droplet. To maintain a homogenous mixture, the solute will have to diffuse towards the center, away from the receding surface. However, this diffusion will be slower than the surface recession resulting in accumulation at the surface until the concentration is high enough to form a solid. Furthermore, when a solute is surface active, it will have a surface excess concentration at the liquid-air interface, which may also influence how fast a solid shell will be formed at the surface of a drying droplet. For protein stabilization this issue is important to consider. First of all, many proteins are surface active and thus have a tendency to be at the interface. This in turn promotes unfolding of the protein because it is thermodynamically favorable when the more lipophilic inner parts of the protein are in contact with relatively hydrophobic air surrounding the droplet. Secondly, a solution will usually consist of multiple solutes. Specifically in the case of protein stabilization, a standard formulation will at least contain a protein, a sugar, and perhaps a buffer and salts or other excipients. Because proteins are larger molecules, they diffuse much slower than other components, such as sugars. During drying of the Protein-sugar distribution in powders prepared by spray drying | 115 droplet, this can result in a separation of the solutes, where the protein will accumulate much faster at the surface than the sugar. Therefore, when the drying is complete, there can be a large fraction of protein at the surface that is not fully encapsulated in a sugar matrix. Consequently, protein instability due to accumulation at the surface of a drying droplet can occur because of both its surface active properties and its relatively slow diffusion.
When the solid shell starts to form at the surface of the droplet at the end of the first stage, evaporation is slowed down considerably due to a gradual viscosity increase of the solid shell. Furthermore, during the second stage, the morphology of the dried particles is determined as well. Although generally the morphology of spray dried particles is very important from a product engineering perspective, the focus in this chapter will be more on protein stabilization. Therefore, current efforts are more focused on building the basis for a model that includes the diffusion characteristics of the solutes, i.e. protein and sugar, and the glass transition temperature. Although models exist that include the glass transition temperature, they are not used for pharmaceutical application such as protein stabilization [1]. In addition, the glass transition temperature is in these cases secondary to the formation of the dried particle, whereas we believe that it is the primary property that determines the transition from the first to the second stage.
The main goal will be to obtain a model that can predict the drying of the droplet both during the first as well as the second stage in one continuous model, without creating two separate models for the first and second stage. For this model, an aqueous solution of trehalose and bovine serum albumin (BSA) is considered. Trehalose was used as an example, since it is a well-known protein stabilizer, the diffusion coefficient as a function of the concentration and temperature is known from literature, and the Gordon-Taylor constants of water/trehalose and trehalose/BSA mixtures are known [2][3][4]. The model will be used to follow the distribution of solutes inside the droplet during evaporation using empirically determined diffusion coefficients for trehalose and water, and an estimated diffusion coefficient for BSA [2].

Open source software
Similar to chapter 5, the model was further developed in GNU Octave, which is freely available [5]. Furthermore, the developed model itself is available online as supplementary data as well as in Appendix B of this thesis.

General model structure and assumptions
116 | Chapter 6 The droplet model that was introduced in chapter 5 was further expanded to include concentric spherical shells inside the droplet, as was also done by Gac and Gradorí [6]. In each shell, the diffusion of a solute, either only trehalose or both trehalose and bovine serum albumin (BSA), and water is calculated, and together with the evaporation, the concentration gradients inside the droplet are predicted. Here, the exchange of mass and heat is calculated from the surface of a subshell to the surface of a neighboring subshell ( Figure 1). Furthermore, each subshell is ideally mixed, implying that any heat or mass diffusing into a subshell is instantly homogeneously mixed. Finally, the system is considered closed, meaning that there is no heat or mass transfer to or from the outside. As a consequence, changes in temperature and humidity changes in the air can occur. Droplet and air layer are divided into concentric subshells, and heat and mass transfer occurs between the surfaces, as depicted by the arrows. Initial subshell volumes of the droplet are the same. Note that the ratio of air to droplet is not to scale, and the depicted number of subshells in the droplet and air layer are not the same as used in further calculations.
Protein-sugar distribution in powders prepared by spray drying | 117 Due to the evaporation of water, the droplet will shrink. During the simulation, the volumes of the individual subshells of the droplet will be kept the same, instead of letting the outer subshell volume decrease more than the other subshells. In order to achieve this, the stepwise calculation is performed as follows. First, the mass of water that evaporates from the outer subshell and the diffused mass of solute and liquid throughout the droplet in a time step (Δt) are calculated. Next, the new solute and liquid mass in each subshell can be calculated. Assuming that the mixtures are ideal, the intermediate subshell volumes (sum of component masses divided by their respective densities) are then also known. With this information, the required volume changes that are needed to equalize the subshell volumes are calculated. The result is that subshells that otherwise would have been smaller than other subshells have received a small amount of liquid and solute mass from neighboring subshells, whereas subshells that otherwise would have been larger than other subshells have lost a small amount of liquid and solute mass to neighboring subshells. This will result in an error, as there is some non-diffusional mass transfer introduced which in reality does not take place. The significance of this error was tested by comparing the results using the aforementioned method with those generated using a second method, with which the non-diffusional mass transfer is smaller or absent. The first method will be referred to as the 'shrink method', whereas the latter will be referred to as the 'merge method'. Although the 'merge method' may be more accurate, the 'shrink method' is preferred due to its shorter calculation time. Of importance is that with both methods the droplet will keep shrinking even after a solid shell has been formed. This will be discussed further below.

Calculation of glass transition temperature
With the calculated component mass fractions and the aid of the secondary or ternary Gordon-Taylor equation, the glass transition temperature inside the droplet can be calculated to determine when the surface of the drying droplet becomes glassy (Equation 1).
Here, subscripts s, w, and p refer to trehalose, water, and BSA, respectively. Furthermore, w is the weight fraction, T g is the glass transition temperature (394, 164, and 150 K for trehalose, water, and BSA, respectively [3,[7][8][9][10]), and k is the Gordon-Taylor constant of a trehalose/water mixture (k sw , 7.90 [3]) or a trehalose/BSA mixture (k sp , 0.08 [4]). Finally, for trehalose only, w p becomes zero by which the equation reduces to the secondary Gordon-Taylor equation (Equation 1 without the the bracketed sections).

Calculation of diffusion coefficients
During the calculations, the diffusion coefficient of the different components will change constantly with the temperature and concentration. To take these changes into account, we designed a tailor made equation to fit the data of a study by Ekdawi-Sever et al. with which the diffusion coefficient for trehalose and water can be calculated as a function of the solute mass fraction and the temperature (Equation 2 and 3, Figure 2) [2]. Here, x solute is the mass fraction of solute and T is the temperature of the droplet (K). Note that no exact diffusion coefficient is known for a ternary mixture that includes a protein, such as BSA. In this scenario, the ratio between the radius of gyration of trehalose (3.4 Å, [11]) and BSA (31.5 Å, [12]) was calculated (31.5 / 3.4 = 9.26), and the diffusion coefficient of trehalose was divided by this ratio to obtain a rough estimate of the diffusion coefficient of BSA. In addition, the mass fraction of solute (x solute in equations 2 and 3) was then taken as the sum of the mass fraction of trehalose and BSA. Protein-sugar distribution in powders prepared by spray drying | 119 Both at 303 K and 358 K the function appears to correlate well with the literature values of the diffusion coefficient over the entire concentration range. It is therefore believed that the calculated values will allow us to correctly estimate the diffusion of trehalose and water during drying.

Calculation of vapor pressure and moisture content
Because the model describes the moisture exchange between air and the droplet/particle, it is important to take into account the effect of dissolved components on the vapor pressure of water. Assuming an ideal solution, the influence of dissolved components on the vapor pressure can be calculated using Raoult's law, in which the saturation vapor pressure is multiplied with the molar fraction of water in the outer subshell of the droplet (x m,w ). However, this will only slow down evaporation due to lowering of the saturation vapor pressure, but not take into account the hygroscopic properties of the solute at higher concentrations. To correct for this, the saturation vapor pressure is further multiplied with the Margules function, which is a measure of the deviation from ideality (Equation 4).
The empirical constants A and B were determined with the aid of dynamic vapor sorption (DVS) data for trehalose, and were found to be -5.0 and 4.6, respectively [3]. This was achieved by calculating the moisture uptake of a dry particle in air with various relative humidities, while changing the constants A and B. The resulting moisture contents according to DVS and our model are shown in Figure 3.

Starting conditions
The volume ratio of air / water was set at 1.5*10 5 and the starting droplet radius was set at 3.75 µm, based on properties of a typical lab-scale spray dryer. Furthermore, the starting concentration of trehalose was set at 5 mg/mL, the starting concentration of BSA was set at either 0.01 or 1 mg/mL, the starting temperature of the droplet at 293 K, and the starting temperature and relative humidity of air at 333 K and 0 %, respectively.

Results & discussion
Temperature gradient inside the droplet A simulation was run in which a temperature gradient inside the droplet was taken into account. This was done to ensure that the choice to omit this temperature gradient in further simulations was justified (i.e. has a negligible effect on the result). The omission of the temperature gradient results in a large decrease of the calculation time, and is therefore desirable. To increase the temperature gradient for visibility, the initial air temperature was increased to 393 K, whereas the initial droplet temperature was decreased to 274 K. A simulation of the first 1 ms is shown in Figure 4. An aqueous 5 mg/mL trehalose solution without BSA was considered for these simulations.
Protein-sugar distribution in powders prepared by spray drying | 121 It was found that indeed the temperature difference between the subshells is hard to distinguish even at these exacerbated conditions. Only at a very small timescale, a difference between the surface and the core of the droplet of about 2.5 K can be seen. Considering the much milder conditions in further simulations and the much longer timescales, the assumption that no temperature gradient inside the droplet exists will not significantly affect the result and therefore appears justified.

Optimum number of subshells and choice for either the 'shrink method' or the 'merge method'
Further simulations were performed to choose the optimal number of droplet subshells. The droplet radius at the formation of the solid shell was calculated at different numbers of subshells for an aqueous 5 mg/ml trehalose solution (no BSA).
Obviously, an infinite number of subshells would yield optimal results, but also requires infinitely long calculation times. With a restricted number of subshells, however, an error is introduced. Due to the assumption that each individual subshell is homogeneously mixed, a larger outer subshell should result in the formation of the solid shell at a later stage (at a smaller droplet radius) due to a larger error in the concentration gradients inside the droplet. This is most pronounced when the droplet is modeled as a single subshell, as no concentration gradients inside the droplet is considered and thus no diffusion. Therefore, the whole droplet will solidify when it is close to the maximum density of trehalose (~1580 kg/m 3 ). In between a single subshell (fast but resulting in the smallest particle size) and infinite subshells (slow but resulting in the largest and most accurate particle size) is an optimal number of subshells with a negligible deviation from the result with infinite subshells while still fast. Furthermore, two methods to handle the receding surface of the droplet were compared: one where the subshell volumes were equalized after each iteration ('shrink method'), and another where the outer subshell was merged with the neighboring subshell only when the thickness was less than 10 nm and the mass fraction difference of water was less than 0.001 ('merge method'). These limits were chosen arbitrarily, although they had to be low enough to prevent constant merging. To make the effect of the chosen number of subshells more pronounced, the drying rate was increase by increasing the initial air temperature to 373 K (see Figure 5).
It was found that with the 'merge method' less subshells were required to reach the maximum radius at which the solid shell is formed than with the 'shrink method'. Due to the generally smaller outer subshell during the simulation with the 'merge method', the error introduced with the assumption of homogeneously mixed subshells is reduced. This reduction will of course depend on the chosen minimum thickness and concentration difference at which a merge of the two outer subshells is allowed. If these limits would be reduced further, the maximum radius would be reached at a lower number of subshells, but the calculation time would increase as well. However, due to the smaller volume of the outer subshell, the time steps during the simulation become smaller to ensure the temperature and mass changes in each time step stay small enough to prevent diverging or oscillating results. Therefore, the total calculation times for simulations with the 'merge method' were actually longer than with the 'shrink method' for the number of subshells that would result in equivalent radii (for example, 20 subshells with the 'shrink method' took about 165 s, whereas 10 subshells with the 'merge method' took about 200 s even though they resulted in the same radius at solid shell formation). Furthermore, the same maximum radius at which a solid shell is formed is obtained at 40 subshells. If indeed the 'shrink method' would introduce a significant amount of mass transfer not attributed to diffusion but to the method used to deal with the receding droplet surface, the outer Protein-sugar distribution in powders prepared by spray drying | 123 subshell would be diluted more and thus form a solid shell at a smaller radius. This appears not to be the case. Due to the fact that the 'merge method' is slower, the 'shrink method' was preferred and used for subsequent simulations. Furthermore, above 40 subshells, the increase in radius at solid shell formation was found to be negligible. Therefore, 40 subshells is considered to be the optimum number of subshells as this appears to give the best balance between accuracy and calculation time.

Pure trehalose
To get a better understanding of what is happening inside the droplet, a simulation with standard conditions was run and the changes in mass fraction of the solute in each subshell of the droplet were tracked in time. and the results of this simulation are shown in Figure 6. The outer radius of each of the subshells is shown as well, indicating the size of the droplet during the drying process. It was found that the model indeed predicted the accumulation of the solute near the surface of the droplet. After about 17.9 ms, the mass fraction of the solute in the outer subshell was high enough for the glass transition temperature to surpass the temperature of the droplet at the surface, changing from a rubbery to a glassy state. For reference, the temperature of the droplet and the glass transition temperature in each shell are shown in Figure 7. It is shown that the increase in glass transition temperature is fast and happens less than one ms before the first solid shell is formed. Furthermore, an increase in the droplet temperature is found already before the outer subshell turned into a glass. Due to the rapidly increasing concentration of solute, the saturation vapor pressure at the surface of the droplet also decreases (Raoult's law). This lower saturation vapor pressure results in a decreased evaporation rate, which Protein-sugar distribution in powders prepared by spray drying | 125 thus reduces the evaporative cooling effect that kept the droplet at a lower temperature in the first place (the wet bulb temperature). What is most striking is the large separation between the glass transition temperature of the outer subshell and the other subshells. When the mass fraction of solute in the outer subshell starts to increase rapidly, the evaporation of water is also inhibited due to the much slower diffusion through the concentrated solute solution. Although once in the glassy state the outer subshell will remain as such due to evaporation and slower diffusion, the concentration differences in the other subshells will start to decrease, as the diffusion there is still much faster. This can be seen by a decrease in the differences between the glass transition temperatures of these subshells. However, during further evaporation, the glass transition temperature of these subshells will diverge again, due to the difference in distance over which water has to diffuse which results in increasing differences of the solute concentrations in the subshells . If the number of subshells would be increased from a finite to an infinite number, the difference in glass transition temperature between the subshells would become smaller and eventually gradual. However, since we consider a finite number of subshells, the artifact of the large difference in glass transition temperature between the outer and inner subshells remains. This artifact is caused by the low diffusion rate of water through the solidified outer shell. However, the question remains whether this low diffusion rate will be found in reality as well. There are at least three scenario's where the diffusion rate would be higher: 1) when the solid shell will be formed much later due to other wetting mechanics, such as capillary forces, 2) when the solid shell will still deform and break up due to low mechanical strength, and 3) when diffusion through the solid shell is faster due to porosity. This implies that the drying time with the current model is a worst case scenario, which will immediately be shortened when any other method is used to estimate the diffusion rate of water through the subshells that would result in an effectively faster diffusion rate through the solid shell. Whereas the developed model considers only homogeneous subshells, all three scenarios are based on inhomogeneity's in the subshells caused by different mechanisms.
For the first scenario one could imagine that the increased viscosity of highly concentrated solution at the surface of the drop in combination with inhomogeneity results in possible channels where mass transfer proceeds faster. These channels would allow water to pass through the shell much faster and perhaps also wet it much faster than is possible when only the diffusion through the concentrated solution is considered. However, no attempt was made to model this scenario as its complexity would require a separate study that is outside the scope of this paper.
The second scenario was also discussed in a study by Ekdawi-Sever et al. [2]. Upon formation of the solid shell, it is uncertain whether the shell is thick enough to stay in the position where it was formed. For example, with 40 modeled shells the outer layer shell is only 6 nm thick when the first glassy shell is formed (Figure 8). Due to the possibly low mechanical strength of such a layer, it could still shrink further until the layer has a sufficient thickness to withstand deformation. It could also simply break forming ruptures, through which water can diffuse more rapidly than through the intact shell. From a modeling perspective this would mean that at some critical solid shell thickness, determined by the mechanical characteristics of the considered solutes, shrinking of the droplet stops and further evaporation results in a vacuole in the core of which the under pressure is withstood by the solid shell, resulting in a hollow particle. This will have to be investigated further experimentally with the aid of mechanical/structural engineering.
Protein-sugar distribution in powders prepared by spray drying | 127 The third scenario, can be roughly modeled by assuming that the maximum porosity of the solidified subshell is equal to the water mass fraction upon solidification. When the water evaporates after the solid shell has formed, the resulting pores will be used for water to diffuse through as vapor. This assumption is based on a model used by Stephen, et al. [13], but should not be considered an accurate prediction of the diffusion through the solid shell. Instead it merely serves to illustrate the influence of different scenario's on the drying characteristics ( Figure 9). The diffusion rate through the porous subshell is calculated according to equation 5, where C is a constrictivity/tortuosity factor (3·10 -5 ), ε is the porosity, D w,vapor is the diffusion coefficient of water vapor in air, and D water is the diffusion coefficient of water in the droplet as calculated with equation 3.  With the faster diffusion of water through the solid shell, the predicted glass transition temperature gradient was indeed found to be much smoother. This also resulted in a more gradual heating of the droplet. Whereas the increase in droplet temperature with the slow diffusion was divided in one to two steps due to the large difference in glass transition temperature of the two outer subshells (see Figure 7), here the increase in temperature is divided into multiple steps. The increase of the droplet temperature in multiple steps is an artifact of the rather rapid increase in water diffusion and thereby increased evaporation rate once the shell transitions from a rubber to a glass. When more subshells would be modeled, the steps become smaller and the temperature increase smoother. We expect that this prediction overestimates the rate at which water diffuses through the solid shell, as the porosity that is used here is most likely an overestimate of the actual porosity that would be available for water diffusion through the subshell. Instead, the solution that correlates best with Protein-sugar distribution in powders prepared by spray drying | 129 experimental values is expected to be between the scenario depicted in Figures 7 and  9.
Further simulations of the ternary water, trehalose, BSA mixture will be calculated using both methods to determine whether the rate of diffusion through the solid shell will be relevant from a protein stabilization perspective, although by default the unmodified diffusion model will be presented unless stated otherwise.

Protein stabilization
To predict the distribution of components in a droplet during protein stabilization by spray drying, BSA was added as a third component. The diffusion coefficient of BSA was estimated with the radius of gyration as described in the materials and methods. To get a good indication of the mass distribution of BSA and trehalose inside the dried particle, a simulation was run with a mass ratio of trehalose to protein of both 5 and 500. The mass fraction of protein and trehalose of the total protein and trehalose content, respectively, is shown in each subshell ( Figure 10). Depicted is the mass fraction of solute that resides in each subshell of the total amount of that specific solute in the particle (i.e. the mass fractions of each solute add up to 1). The total mass ratio trehalose/protein was either 5 (left) or 500 (right).
A clear separation of trehalose and BSA was found with both trehalose/BSA mass ratio's. Due to the slower diffusion of the protein, accumulation of protein is indeed predicted near the surface. Furthermore, the increased protein content near the surface replaces the trehalose, which would otherwise be equally distributed throughout the particle (in case of pure trehalose each subshell of equal volume would contain 1 / 40 = 0.025 of the total trehalose content). With a trehalose/BSA mass ratio of 500, this reduction in trehalose content near the surface is not observed due to the low absolute mass of protein in the particle. Most striking is that the accumulation of BSA near the surface is already very pronounced, despite the fact that the surface excess concentration of surface active components (such as proteins) is not taken into consideration. Doing so would result in an even more pronounced accumulation near or at the surface of the dried particle. Of further note is that, due to the accumulation of BSA at the surface relative to trehalose, there will be a difference in the trehalose/BSA mass ratio throughout the dried particle. Because the stability of the protein is in part dependent on this ratio, it could influence the effectiveness of the protein stabilization process if only the initial trehalose/BSA mass ratio would be taken into account during product development. The protein accumulation near the surface was more pronounced with the trehalose/BSA mass ratio of 500 than with the lower trehalose/BSA mass ratio of 5. However, with the starting trehalose/BSA mass ratio of 500 the protein may still be sufficiently protected as the minimum trehalose/BSA mass ratio is still about 250, whereas with a low starting trehalose/BSA mass ratio of 5 the minimum is only about 3.5. To get a better overview of distribution of the trehalose/BSA mass ratio throughout the particle, the cumulative distribution is shown in Figure 11. With both scenario's, close to 60 % of the total protein mass content is encapsulated in less trehalose than when both components would have been homogeneously distributed. Furthermore, results obtained with simulations where a porous solid shell was considered did not deviate much from this result (data not shown). This can be expected, since the distribution of the components will, for the most part, already have occurred before the solid shell is formed. The main difference was that the distribution was slightly wider, leading to a lower minimum and a higher Protein-sugar distribution in powders prepared by spray drying | 131 maximum trehalose/BSA ratio at the surface and the center of the particle, respectively. This is most likely caused by faster evaporation during the second stage, resulting in less redistribution of components inside the droplet after the initial solid shell was formed.
Whether or not and how much the lower trehalose/BSA mass ratio near the surface will affect the stability cannot be said, however, it is likely that it will be reduced. This will also depend largely on the protein that is used, as some protein are much more sensitive than others. Furthermore, the extent to which the sugar/protein mass ratio will change is dependent on the spray drying conditions. For example, faster drying would result in more accumulation at the surface, and thus possible reduced stabilization. In addition, the molecular weight and shape of the sugar and protein used will also play a role. The smaller the difference in the diffusion rate between the sugar and the protein, the smaller the change in sugar/protein mass ratio.
In case a small protein is chosen that diffuses faster than the chosen sugar, one might even find that sugar accumulates at the surface and protein in the center of the drying droplet. Finally, when the surface excess concentration would be considered, the amount of protein that will not be properly encapsulated will increase even further. Therefore, it should be kept in mind when developing stabilized protein formulations by spray drying, that adding surfactants to prevent accumulation of protein at the surface due to its surface active nature, may actually not be sufficient. Instead, simulations such as performed here can be used to more accurately adapt the initial composition of the solution to ensure proper stabilization of the protein, for example by simply increasing the amount of sugar or, if this is limited by the desired loading, changing the drying conditions to slow down the drying rate.

Conclusions
The presented model gives an insight into the distribution of components with respect to protein stabilization during spray drying. Furthermore, the model can easily be adapted to suit other needs and change components. However, the method strongly depends on the use of reliable diffusion coefficient data and Gordon-Taylor coefficients, which are vital for accurate predictions of the component distribution inside the dried particle. Especially for multi-component mixtures, this can become quite complex. Therefore, at a later stage, a welcome addition would be the use of molecular dynamics models to predict the interaction between different solutes and water, and couple these to the model of a drying droplet.
Furthermore, a better understanding of the morphology changes and the mass transfer during the second stage is welcomed, to further improve the accuracy of predictions in this region. However, it was shown that for protein stabilization, the most important factor is the distribution of the components inside the droplet, which