Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Modelling phase separation in amorphous solid dispersions

Modelling phase separation in amorphous solid dispersions Much work has been devoted to analysing thermodynamic models for solid dispersions with a view to identifying regions in the phase diagram where amorphous phase separation or drug recrystallization can occur. However, detailed partial differential equation non-equilibrium mod- els that track the evolution of solid dispersions in time and space are lacking. Hence theoretical predictions for the timescale over which phase separation occurs in a solid dispersion are not available. In this paper, we address some of these deficiencies by (i) constructing a general multi- component diffusion model for a dissolving solid dispersion; (ii) specializing the model to a binary drug/polymer system in storage; (iii) deriving an effective concentration dependent drug diffusion coefficient for the binary system, thereby obtaining a theoretical prediction for the timescale over which phase separation occurs; (iv) calculating the phase diagram for the Felodipine/HPMCAS system; and (iv) presenting a detailed numerical investigation of the Felodipine/HPMCAS sys- tem assuming a Flory-Huggins activity coefficient. The numerical simulations exhibit numerous interesting phenomena, such as the formation of polymer droplets and strings, Ostwald ripen- ing/coarsening, phase inversion, and droplet-to-string transitions. A numerical simulation of the fabrication process for a solid dispersion in a hot melt extruder was also presented. Keywords: amorphous solid dispersion, phase separation, mathematical model, drug diffusion 1 Introduction Drugs that are delivered orally via a tablet should ideally be readily soluble in water. Drugs that are poorly water-soluble tend to pass through the gastrointestinal tract before they can fully dissolve, and this typically leads to poor bioavailability of the drug. Unfortunately, many drugs currently on the market or in development are poorly water-soluble, and this presents a serious challenge to the pharmaceutical industry. Many strategies have been developed to improve the solubility of drugs, such as the use of surfactants, cocrystals, lipid-based formulations, and particle size reduction. The literature on this topic is extensive, and recent reviews can be found in [1–3]. Corresponding author. Email: martin.meere@nuigalway.ie arXiv:1902.05410v2 [cond-mat.soft] 25 Jun 2020 One particularly effective strategy to improve drug solubility is to use a solid dispersion [4–6]. A solid dispersion typically consists of a hydrophobic drug embedded in a hydrophilic polymer [7, 8] matrix, where the matrix can be either in the amorphous or crystalline state. The drug is preferably in a molecularly dispersed state, but may also be present in amorphous particles or even in the crystalline form (though this is usually undesirable); see Figure 1. The drug release concept for most solid dis- persions is based on the so-called spring and parachute effect [9]. When the drug and the hydrophilic polymer dissolve in solution, a supersaturated drug solution is quickly created (the spring). Although the drug concentration then subsequently decreases, the rate of decrease is slowed by drug-polymer interactions in the dispersion, so that the drug can be present at supersaturated levels in the solution for a period of some hours (the parachute). This results in improved bioavailability of the drug when the solid dispersion dosage form is taken orally. Molecularly dispersed drug in the polymer matrix (desirable) Drug molecule Polymer matrix Crystal drug formation Contains amorphous has occurred (undesirable) drug-rich domains Figure 1: Adapted from [5]. In this figure, we show three possible structures for a polymer/drug dispersion. Top: Here the drug is in the molecularly dispersed state, which is usually desirable for a solid dispersion. Bottom left: Here the dispersion contains drug in the crystalline form. Bottom right: Here the dispersion contains amorphous drug-rich domains. Drug loading in most dispersions greatly exceeds the equilibrium solubility in the polymer matrix for typical storage temperatures. Hence these systems are usually unstable, with phase separation eventually occurring [5]. In such cases, the drug will eventually crystallise out or form an amorphous phase separation. However, if the dispersion is stored well below the glass transition temperature [10] for the polymer, and is kept dry, this can happen extremely slowly. The system is then for all practical purposes stable, and is said to be metastable. The humidity of the storage environment can be an issue because even small amounts of moisture can significantly affect the glass transition temperature. Hence polymers that have high glass transition temperatures and that are resistant to water absorption have become popular. An example of one such polymer is Hydroxypropyl Methylcellulose Acetate Succinate (HPMCAS). Phase separation of solid dispersions in storage is clearly undesirable from the point of manufac- turers. Hence much work has been devoted to constructing phase diagrams for solid dispersions with a view to identifying regimes where drug recrystallization or amorphous phase separation can occur. These phase diagrams are constructed with the aid of thermodynamic models. The most widely used thermodynamic model in this context is the Flory-Huggins model [11–13] for polymer solutions. Flory-Huggins theory is a lattice-based model in which the drug and polymer are confined to live 2 on a regular lattice. Flory-Huggins theory is an extension of regular solution theory, as explained in Chapter 7 of [13]. In the context of a drug/polymer system, each drug molecule is taken to occupy one lattice site and each polymer segment is taken to occupy m  1 sites. Under a number of further simplifying assumptions [13], the change in entropy and enthalpy associated with the mixing of the mix polymer and drug are calculated. With these in hand, the change in bulk Gibbs free energy (g ) per mole associated with mixing is readily calculated, and is found to be mix = X ln( ) + X ln( ) +  X  ; (1) d d p p dp d p RT where R is the gas constant, T is the temperature, X ; X = 1X are the mole fractions of the drug d p d and polymer, respectively, and  ;  are the volume fractions of the drug and polymer, respectively. d p The quantity  is referred to as the Flory-Huggins interaction parameter, and it is discussed further dp below. The mole fractions and volume fractions are related via the formulae X mX d p = ;  = : (2) d p X + mX X + mX d p d p When the model is applied to real binary systems, m can be calculated using the formula m = (3) where V , V (molar ) are the molar volumes of the polymer and drug, respectively. p d The mixing of the polymer and drug is spontaneous if G < 0. The Flory-Huggins parameter mix takes the form dp w + w dd pp =  w (4) dp dp where  is a positive parameter, and w , w , w give measures of the drug-polymer, drug-drug and dp dd pp polymer-polymer interaction energy, respectively. If  < 0 then w < (w + w )=2 indicating dp dp dd pp that that the mixed state has lower energy than the separated pure states, so that mixing is favoured. Conversely,  > 0 is indicative of demixing being favoured. However, these statements are indica- dp tive rather than precise, as will be explained in Section 3. We should also note that  is temperature dp dependent, and is usually given the empirical form (T ) = + (5) dp where ; are constants. Flory-Huggins theory has frequently been used to analyse the stability of binary solid dispersion systems in storage; see, for example, [14–23]. In many of these studies, the Flory-Huggins interaction parameter is first estimated using the melting point depression method [24], or using the Hildebrand and Scott method [25], which involves the estimation of solubility parameters. Once estimates for (T ) have been obtained, the Gibbs free energy of mixing G can be calculated, which in turn dp mix enables the construction of phase diagrams for the systems. Phase diagrams assist with the identifi- cation of regions in composition-temperature space where the system is prone to recrystallization or amorphous phase separation. The models we shall develop in the current study are generic and are not tied to making a specific choice of statistical model. However, given the particular importance of Flory-Huggins theory in applications, we shall derive detailed results for this case. Also, all of our numerical illustrations are calculated within the context of Flory-Huggins theory. It should be emphasized that Flory-Huggins theory does involve quite a number of simplifying assumptions which are not appropriate for some systems; see [26] for a recent critique of the model. 3 2 Theoretical formulation 2.1 A multicomponent diffusion model for solid dispersions We develop a multicomponent diffusion model for the evolution of the concentrations of the com- ponents constituting a solid dispersion. We suppose for the moment that there are p components. However, in the analysis we shall consider in the current study, we will in fact have p = 2, with one of the components being the polymer, and the other being the drug. For a dissolving solid dispersion, there are three components p = 3: the polymer, the drug, and the solvent. The chemical potential  (J/mole) of species i (i = 1; 2; :::; p) gives the Gibbs free energy per mole of species i, and is given here by ( [27]) b 2 2 =   r X (6) i i i i where b 0 =  + RT ln(a ) (7) i i and where  is the bulk chemical potential of species i,  is the chemical potential of species i i0 2 2 in the pure state, a is the activity of species i, and the term involving  > 0 (m J/mole) penalises the formation of phase boundaries ( [28], [29], [30], [31], [32]). The parameters  are referred to as gradient energy coefficients ( [33], [34]). Here X is the molar fraction of species i (i = 1; 2; :::; p), and the activities can depend on these molar fractions, so that a = a (X ; X ; :::; X ): i i 1 2 p The molar fraction is related to the molar concentration via X = V c (8) i i i where V (molar ) is the molar volume of species i. The flux of species i (molarm/s) is given by J = c v (9) i i i where c (molar), v (m/s) give the molar concentration and drift velocity, respectively, of species i. i i The drift velocity v gives the average velocity a particle of species i attains due to the diffusion force acting on it, and is given here by v = M F = M r (10) i i i i i where M (moles/kg),F (J/[mmole]) give the mobility and diffusion force, respectively, for species i i i. Equations (9) and (10) give J = M c r : (11) i i i i Conservation of mass for species i implies that @c +r J = 0 (12) @t and using (11) now gives @c = r (M c r ) i i i @t or equivalently @X i i = r D X r (13) i i @t RT with 2 2 = ln(a )  r X (14) i i RT 4 2 2 2 for i = 1; 2; :::; p, and where  =  =RT > 0 (m /molar), and i i D = M RT (Einstein relation) i i is the self-diffusion coefficient for species i. The model formulation given by (13) and (14) based on chemical potentials will be used for the numerical scheme described in Section 4. However, it is also of value to develop a formulation involv- ing diffusion coefficients since these yield immediate information regarding timescales for transport processes, and will also the enable the development of analytical results via a linearization process. Diffusion Coefficients Using (6), (7) and (11) gives RT 2 2 J = M c r = M c ra  r(r X ) i i i i i i i i and then using the fact that the activities depend on the molar fractions gives RT @a 2 2 J = M c rX  r(r X ) : (15) i i i j i a @X i j j=1 Using (8), we can now write (15) as 2 2 J = D rc + D " c r(r c ) (16) i ij j i i i j=1 2 2 2 where " = V  and where the diffusion coefficients D (m /s) are given by i ij i i V X @a j i i D = D : i; j = 1; 2; :::; p (17) ij i V a @X i i j Conservation of mass (12) then implies that (reverting to molar fractions) @X V i i 2 2 = r D (X)rX D  X r(r X ) i = 1; 2; :::; p (18) ij j i i i @t V j=1 where X = (X ; X ; :::; X ), and where we have included the concentration dependence of the dif- 1 2 p fusion coefficients D here to emphasise that this system is in general a coupled system of nonlinear ij diffusion equations. It should benoted that the equations (18) are not independent since X = 1, i=1 and so it is sufficient to solve for p 1 concentrations only. 2.2 Activity coefficients The activities a are usually written as a = X i i i where the = (X ; X ; :::; X ) are referred to as activity coefficients. Equations (17) now give i i 1 2 p V X @ j i i D = D  + i; j = 1; 2; :::; p (19) ij i ij V @X i i j where  is the Kronecker delta. ij The details of the interactions between the species in solution are captured in the modelling by choosing appropriate forms for the activity coefficients = (X ; X ; :::; X ). The construction of i i 1 2 p appropriate forms for the for various solutions is a large subject with a large literature; see, for example, the books [35] and [36]. 5 2.3 The storage problem for a binary mixture In the current study, we shall be modelling the behaviour of solid dispersions in storage. In this case, we have p = 2, with the label 1 referring to the drug and the label 2 referring to the polymer. However, for transparency, we choose here to use the labels d; p rather than 1; 2, where d stands for drug, and p for polymer. Then using (18) and the fact that X = 1 X , we have p d @X 2 2 = r D (X )rX D  X r r X : (20) eff d d d d d @t where the effective concentration-dependent diffusion coefficient for the drug in the solid dispersion is X @ @ d d d D (X ) = D (X ) V D (X )=V = D 1 + : (21) eff d dd d d dp d p d @X @X d d p For the particular case of a binary Flory-Huggins theory (see Section 1), the activity coefficients are given by d d ln( ) = ln + 1 +   ; (22) d dp X X d d p p ln( ) = ln + 1 + m  ; (23) p dp X X p p where the volume fractions  ;  are given by (2). Substituting (22) in (21) gives d p 2 2 2 (m (m 1)X )(m (m 1)X ) 2 m X (1 X ) d d dp d d D (X ) = D : (24) eff d d (m (m 1)X ) It is more instructive to write this expression in terms of the volume fraction of drug. Writing D (X ) = D ( ), we obtain eff d eff d D ( ) = D (1 + (m 1) ) 1 + 1  2  (1  ) : (25) eff d d d d dp d d This expression is particularly useful because it yields insight into how the mobility of the drug in the dispersion depends on the length of the polymer chains (m), the dispersion composition ( ), and the character of the drug-polymer interaction ( ). We shall analyze this expression further in Section 3, dp and also show how it can be used to calculate the timescale over which phase separation may occur. An equivalent formulation for the Flory-Huggins model involving the chemical potential for the drug  is given by (see (13) and (14) above): @X = r (D X r ) (26) d d @t where d d0 = (27) RT and with X (m 1)(1 X ) 1 X d d d 2 2 2 = ln + +  m  r X : (28) dp d m (m 1)X m (m 1)X m (m 1)X d d d We suppose that the solid dispersion occupies a two-dimensional region . The governing equa- tion for the drug concentration in may be written in the conservation form @X +r J = 0; @t 6 where the drug flux J is given by 2 2 J = D (X )rX + D  X r(r X ): (29) d eff d d d d d We need to supplement the governing equation in with boundary conditions on @ , and we choose these here to be J  n = 0; rX  n = 0 on @ : (30) d d The first of these conditions J  n = 0 implies that the drug cannot penetrate the boundary of the domain. The other condition rX  n = 0 is the natural boundary condition for the variational formulation of the problem, and it implies that the interfaces between the polymer rich and the drug rich domains meet the boundary at right angles. Finally, to obtain a well-posed problem, we need to impose an initial condition and we choose this here to take the form X (x; y; t = 0) = X (x; y) for (x; y) 2 ; (31) where X (x; y) is a given function. Gathering together the governing equation, the boundary conditions and the initial condition, we obtain the following initial boundary value problem: @X 2 2 = r D (X )rX D  X r r X in eff d d d d d @t rX  n = 0; r(r X ) n = 0 on @ ; (32) d d X (x; y; t = 0) = X (x; y) for (x; y) 2 2.4 Phase separation in a Flory-Huggins binary mixture The bulk free energy and spinodal decomposition Spinodal decomposition for binary systems has been long understood using thermodynamic rea- soning, and is well described elsewhere; see, for example, Chapter 5 of [37], Chapter 7 of [13], or [38]. Hence our description of the background theory here will be quite brief, and we will empha- sise instead the particular details for the Flory-Huggins system. The bulk free energy density g for the binary mixture constituting the solid dispersion is given by ( [39]) b b g =  X +  X ; (33) b d p d p b b where  ;  give the bulk chemical potential of the drug and polymer, respectively, and where d p b 0 =  + RT ln(a ) for i = d; p: i i This leads to 0 0 g =  X +  X + RT (X ln(X ) + X ln(X )) + RT (X ln( ) + X ln( )): b d p d d p p d d p p d p mix The Gibbs free energy of mixing g is given by mix 0 0 g = g  X  X = RT (X ln(X ) + X ln(X )) + RT (X ln( ) + X ln( )): b d p d d p p d d p p b d p If we now use (22) and (23) and the fact that X = 1 X , we arrive at p d mix g X m(1 X )  mX (1 X ) d d dp d d = X ln + (1 X ) ln + : (34) d d RT m (m 1)X m (m 1)X m (m 1)X d d d 7 mix Figure 2: (a) Plot of the bulk free energy of mixing g as a function of a drug molar fraction 1s 2s 2 mix 2 X . The spinodal points X , X are the solutions to d (g )=dX = 0. In the spinodal region d d b d 1s 2s 2 mix 2 (X ; X ), we have d (g )=dX < 0 and D (X ) < 0. (b) Phase diagram for the binary eff d d d b d mixture. Here is the coexistence curve, is the spinodal curve, T is the temperature for the free energy density diagram in (a), and T is the critical temperature above which the dispersion is homogeneous. mix In Figure 2 (a), we plot a free energy of mixing diagram g as a function of drug molar fraction 1s 2s X . In this diagram, the points X , X are the solutions to d d 2 mix d (g ) = 0; dX 1s 2s and are referred to as the spinodal points. The region (X ; X ) is referred as the spinodal region, d d and for points X in this region, we have 2 mix d (g ) < 0: dX Compositions X in the spinodal region are unstable, and will split into two phases characterized 1u 2u by the compositions X and X as shown in Figure 2 (a); see [37] for more details. The points d d 1u 2u X , X are referred to as the binodal points, and are defined by the common tangent construction d d shown in Figure 2 (a). The binodal and spinodal points define the coexistence and spinodal curves, respectively, and these are plotted in the phase diagram shown in Figure 2 (b). Using equation (34), we obtain 2 mix d (g ) q(X ) = RT (35) 2 3 dX (1 (1 1=m)X ) X (1 X ) d d d where q(X ) = AX + BX + 1 (36) d d and where 1 1 1 1 1 A = (1 2 ) + 1; B = + (1 2 ) 2: (37) dp dp 3 2 2 m m m m m 2 mix 2 Hence there is a spinodal region with d (g )=dX < 0 if q(X ) < 0 in this region. Inspecting b d (36), we see that q(X ) can be negative if q(X ) = 0 has real roots, that is, if d d B 4A > 0; 8 and using (37), this leads to (2 (1 + 1=m)) 4=m > 0 dp which holds true if 2 q 2 1 1 1 > 1 + p = 1 + V =V : dp d p 2 m 2 Hence, we have a spinodal interval if >  (m) (38) dp dp where 1 1 (m)  1 + ; (39) dp 2 m and where  (m) is a critical value for the Flory-Huggins parameter. If (38) holds true, then there is dp 1s 1s a spinodal interval (X ; X )  [0; 1] where d d p p 2 2 B B 4A B + B 4A 1s 2s X = ; X = ; (40) d d 2A 2A and where A, B are given in (37). The diffusion coefficient and spinodal decomposition Using (24) and (35), elementary calculations show that 2 mix d (g ) y b D (X ) = M (X ) (41) eff d d dX where M (X ) = M X (1 X ) d d d d with M = D =RT , and where M (X ) is a concentration-dependent drug mobility; see, for example, d d d 2 mix 2 equation (3.6) of the paper [40]. Hence, for 0 < X < 1, it is clear from (41) that d (g )=dX < 0 b d implies that D (X ) < 0: (42) eff d Hence, an equivalent criterion for spinodal decomposition to occur is that there exist a region in 0  X  1 where D (X ) < 0, that is, that there exist a region where drug diffusion is against the d eff d concentration gradient (uphill diffusion). 3 Qualitative results and discussion Although the model we have derived in the current study is quite general, and is not tied to any specific statistical model for a solid dispersion, the detailed results we shall present in this section are for the Flory-Huggins case. 3.1 The effective diffusion coefficient for the drug in the dispersion From (25), the scaled effective diffusion coefficient for the drug in the dispersion is given by D ( ) 1 eff d = (1 + (m 1) ) 1 + 1  2  (1  ) (43) d d dp d d D m where we recall that D is the temperature-dependent self-diffusion coefficient for the drug. Equation (43) is of particular value since it yields information on how the mobility of the drug in the dispersion 9 Figure 3: Plots of the scaled effective diffusion coefficient for the drug in the polymer dispersion as a function of the drug volume fraction. Here positive values of the diffusion coefficient correspond to standard drug diffusion down the concentration gradient, while negative values correspond to phase separation of the drug and the polymer, with larger negative values (in absolute terms) corresponding to more rapid phase separation. We have plotted the scaled drug diffusion coefficient for (a) the Flory- Huggins interaction parameter  = 3 and various values of the polymer chain length m, and, (b) dp polymer chain length m = 50 and various values of the Flory-Huggins interaction parameter  . See dp the main body of the text for further discussion. depends on the polymer chain length, the dispersion composition, and the character of the drug- polymer interaction. In Figure 3 (a), we have plotted (43) for the Flory-Huggins interaction parameter  = 3 (which dp is in the unstable regime) and various values of the polymer chain length m. It should be emphasized that positive values for D correspond to standard drug diffusion down the concentration gradient, eff while negative values correspond to unstable regimes where phase separation of the drug and polymer can occur. In Figure 3 (a), it is clear that if the drug loading  is sufficiently low, then D > 0 and d eff the solid dispersion is stable. However, for larger (and more realistic) drug loadings, D < 0, and eff the system is unstable. It is interesting to note that the system becomes more unstable as the length of the polymer chains increase. It is also clear from the curves in Figure 3 (a) that the relationship between the initial drug loading in the dispersion and the initial rate of phase separation is not altogether obvious. It is not necessarily the case that increasing drug loading corresponds to increasing initial dispersion instability. Rather, there is in fact a well defined worst choice for the initial drug loading from the point of view of stability in the initial stages. This worst choice corresponds to the minima of the curves displayed in Figure 3 (a), since these minima correspond to the fastest rates of phase separation. For m  1, the minimum of D ( ) occurs at eff d 1 + 2 + (1 + 2 ) 6 dp dp dp min : (44) dp min These theoretical results predict that choosing initial drug loadings  above or below  should min lead to improved dispersion stability in the initial stages. For m;   1, we have   0:67. dp In Figure 3 (b), we plot (43) for the fixed polymer length m = 50, and various values of the Flory- Huggins interaction parameter  . For m = 50, the critical value for  is given by   0:5707 dp dp dp 10 c (see equation (39)). Recall that for  <  , the system is stable for all drug loadings  , and that dp d dp for  >  , there is a regime of unstable drug loadings. This is borne out by the curves displayed dp dp in Figure 3 (b). These curves predict that the system becomes more unstable with increasing values of  , and this is as expected given the dependence of  on the interaction energies - see equation dp dp (4). Figure 4: Schematic of a phase separating solid dispersion where polymer-rich regions with charac- teristic lengthscale l have formed. A formula for the timescale of evolution of such a dispersion is given in the main body of the text; see equation (45). 3.2 Timescale for phase separation in a solid dispersion In Figure 4, we give a schematic of a phase separating solid dispersion where polymer-rich regions have formed. The characteristic lengthscale of these regions is denoted by l. In order for such regions to form, the drug must have diffused away over a lengthscale of order l, and the timescale over which this diffusion occurs is estimated by (see (25)) 2 2 l l 1 = = (45) 0 0 0 0 D (T )j(1 + (m 1) )[1 + (1=m 1) 2 (T ) (1  )]j jD ( )j d dp eff d d d d where  is the initial uniform volume fraction of the drug in the dispersion, and T is a representative storage temperature. It should be emphasized that this formula is just an estimate since, in reality, the drug volume fraction evolves in space and time. Hence, (45) should only be used as a rough rule of thumb. In Section 4, we evaluate this formula by comparing it with detailed numerical results, and satisfactory agreement is generally found. Equation (45) may, in appropriate circumstances, be used to estimate the shelf life of a solid dispersion product. To see this, suppose that l denotes the largest acceptable size for polymer-rich domains (or drug-rich domains) in the product. Then, since  estimates the timescale for these regions to form, it also estimates the timescale for the shelf life of the product. However, care should be taken when using (45) since, apart from the fact that is based on a fixed value of  , it also incorporates a number of significant assumptions - for example, it assumes that the dispersion is perfectly dry, and that Flory-Huggins theory is an appropriate statistical model for the system. 11 3.3 Criteria for a stable solid dispersion Although the drug loading in real solid dispersions is typically high and in the unstable regime, it is nevertheless worthwhile specifying conditions under which the stability of the dispersion is guaran- teed. The results we display here are based on the discussion given in Section 2.4. For  < dp dp where  = (1 + 1= m) , the system is stable irrespective of the choice of the uniform initial dp 0 c 0 drug load  . For  >  , the dispersion is unstable if the initial drug loading  is chosen in the dp d dp d + + interval ( ;  ), but stable if chosen in either of the intervals (0;  ) or ( ; 1), where d d d d 8 9 < 2 = 1 1 1 1 1 2 = 1 + 1  1 + 1 : 2 : 2 m 2 m  ; dp dp dp These results are based on the bulk free energy only, and do not take account of interfacial energy. However, the interfacial energy can be readily incorporated into the analysis, and this is discussed in Appendix A. 4 Numerical results and discussion 4.1 The numerical method For the purposes of numerical calculations, we take the integration domain to be the square region = f(x; y)j 0 < x < L; 0 < y < Lg with boundary @ . The governing equation to be solved is defined by the equations (26), (27) and (28). The boundary conditions are given by r  n = 0 and rX  n = 0 on @ ; (46) and the initial condition takes the form (31). The boundary conditions (46) are equivalent to those given in (30). The governing equation was numerically integrated using the finite element package COMSOL Multiphysics. A mesh sensitivity analysis was performed to investigate the influence of the size of the mesh on the results. The solution was assumed to be mesh independent when there was less than 1% difference in the mole fraction of drug between successive refinements. The final mesh used in the simulations was triangular and consisted of 7553 vertices and 14796 triangles. The numerical solutions all conserved the total mass of drug in the system to within 1%. 4.2 Parameter values We consider parameter values that are appropriate for a solid dispersion consisting of the drug Felodip- ine (FD) and the polymeric excipient HPMCAS. Felodipine is a calcium channel blocker that is com- monly used to treat blood pressure. For this system, the Flory-Huggins interaction parameter is given as a function of temperature by (see [16]) 7830:4 (T ) = 18:767 + : (47) dp Using data taken from [16], the molar volume for FD is V = 300:19 cm /mol and the molar volume of HPMCAS is V = 14007:78 cm /mol, so that V 14007:78 m = =  46:6630: V 300:19 From (39), the critical value for the interaction parameter below which phase separation cannot occur is given by 1 1 (m) = 1 + p = 0:6571: dp 2 m 12 o 2 1 2 1 T ( C)  (T ) D (T ) (m s ) D (X ) (m s ) dp d eff a 18 17 40 6.2383 1:1661 10 8:8605 10 17 15 50 5.4645 1:5494 10 1:0113 10 16 15 60 4.7371 1:3787 10 7:6107 10 15 14 75 3.7245 2:0297 10 8:3587 10 14 13 90 2.7954 1:7356 10 4:9151 10 14 12 100 2.2176 5:7436 10 1:1670 10 13 12 110 1.6699 1:6336 10 2:0806 10 13 12 120 1.1501 4:0902 10 2:2657 10 Table 1: Illustrative values for some of the parameters of the FD/HPMCAS system at various tem- peratures. Here the initial weight fraction of drug is 70%, which corresponds to an initial drug molar fraction X  0:9909. The self-diffusion coefficient for Felodipine was estimated in [41] (Chapter 4, page 133) to be A A 2 3 2 1 D (T ) = exp(A ) exp exp m s (48) d 1 T T where A = 18:03, A = 445:84 K, A = 874:81 K. Some illustrative values for the diffusion 1 2 3 coefficients and the Flory-Huggins interaction parameter are displayed in Table 1. For the numerical simulations displayed in the current study, we take the size of the square domain to be given by L = 2 mm. The thickness of the interfacial regions is dictated by the parameter  , and here we chose the value  = L=50 = 4 10 m. We illustrate how the initial conditions were specified by considering a particular case. We consider the case where the initial weight fraction of drug is 80%. This means that the initial weight of FD divided by the weight of FD plus the weight of HPMCAS is 0.8. This corresponds to an initial molar drug fraction of X = 0:9947. More precisely, we choose the initial molar fraction of drug to be a small random perturbation about this level given by X (x; y; t = 0) = 0:9947(1 + rnd(x; y)) where rnd(x; y) is a normally distributed random function with a mean value of zero and a standard 5 5 deviation of 10 . The standard deviation for all of the initial conditions was taken to be 10 , with one exception - the numerical results displayed in Figure 6 took the larger value 0:007 to simulate a coarse initial mixture in a hot melt extruder. In Figure 5, we plot the phase diagram for the Felodipine/HPMCAS system. All that is required to calculate the phase diagram here is a knowledge of the Flory-Huggins interaction parameter, and this has been given in (47). The spinodal curve T ( ) is obtained by setting D ( ) = 0 in (43) to s d eff d obtain 1 + 1  2  (1  ) = 0; (49) d dp d d where  is given by (5). Solving (49) for T gives the spinodal curve dp 2  (1  ) d d T ( ) = s d 1 + 1  2  (1  ) d d d where = 18:767, = 7830:4, m = 46:6630 for the Felodipine/HPMCAS system. The binodal curve was estimated numerically. The calculation involved simultaneously solving the pair of equa- b 1u b 2u b 1u b 2u 1u 2u 1u 2u tions  (X ) =  (X ) and  (X ) =  (X ) for the binodal points X , X with X < X . p p d d d d d d d d d d This was implemented using the the fsolve command in MAPLE. 13 Figure 5: The phase diagram for the Felodipine/HPMCAS system. 4.3 Numerical results The first numerical calculations we display simulate the hot melt mixing process for a Felodip- ine/HPMCAS system. To implement this, we used a time dependent temperature profile that treats o o the case where the mixture begins at 25 C and rises in temperature at a rate of 10 C per minute until o o it reaches 145 C. We then assumed for simplicity that the melt cooled linearly back down to 25 C over a period of 30 minutes. The results of the calculation are shown in Figure 6. We see in Figure 6 (a) that the initial drug/polymer mixture is quite coarse (badly mixed). As the temperature rises, we see from Figure 6 (b) that the mixture becomes increasingly homogenous. By the time the mixture has achieved its maximum temperature, Figure 6 (c), it is quite well-mixed. Figure 6 (d) shows that the amount of phase separation that occurs during the cooling process is insignificant. Hence, the numerical results shown here predict a successful hot melt extrusion process for the manufacture of a solid dispersion. We note that different heating and cooling regimes are easily simulated using the model. In Figure 7, we superimpose numerical solutions on a phase diagram for the Felodipine/HPMCAS system. Each of these solutions corresponds to an evolution time of 6 months, with the dispersion mixture beginning from an initially approximate uniform state. We see that there is no significant phase separation for the 30 C degree cases, and for the cases in the metastable and stable regions. This predicts that Felodipine/HPMCAS systems should not suffer considerable phase separation un- der normal storage temperatures. However, we should caution that we are modelling the case of zero relative humidity here. For the cases that do exhibit significant phase separation, we note a coarser separation morphology for higher temperatures. In these figures, dark red corresponds to drug-enriched domains (relative to the initial concentration) while dark blue correspond to polymer- enriched domains. We also note the occurrence of polymer droplets and strings - we return to this issue below. Further numerical simulations are displayed in Figures 8, 9, 10, 11, and these correspond to weight fractions of drug of 80%, 60%, 40% and 20%, respectively. Recall that decreasing weight fractions of drug correspond to increasing weight fractions of polymer since the system is binary. In a given figure, each column corresponds to a given temperature as labelled, and reading a column from top to bottom corresponds to increasing time for the dispersion for the given temperature. We have chosen here not to use the same times for the different temperatures since the rate at which a dispersion evolves depends on temperature. We now highlight some notable features of these numerical simulations. 14 Figure 6: Numerical simulation of the behaviour of a solid dispersion in a hot melt extruder. The simulations here are for a FD/HPMCAS solid dispersion and were obtained by numerically integrating the initial boudary value problem defined in Section 4.1. The colours correspond to different mole fractions of the drug as defined by the colour bar. The weight fraction of drug here is 50%, and the other parameter values can be found in Section 4.2 and Section 4.3. These simulations represent a succesful extrusion where an initially coarse mixture is heated and then cooled to form a well-mixed dispersion. Two phases eventually emerge. The numerical results show that the systems eventually evolve into two distinct phases, characterized by deep blue domains (polymer-rich) and deep red do- mains (drug-rich). Ostwald ripening/coarsening. Another notable feature in many of the numerical illustrations is the formation of polymer droplets (blue discs) in the dispersion, followed by a subsequent growth in their size; see, for example, the third column in Figure 8. This is a well-known and common phenomenon in multicomponent solid systems, and is often referred to as Ostwald ripening or coarsening [42]. We also note the general trend that dispersions at higher tempera- ture tend to be coarser. Phase inversion. The system exhibits the phase inversion phenomenon [31] as the polymer content increases. To see this, consider the panels in Figure 8. These correspond to the case where the polymer content is low (20% by weight), and we see the emergence of polymer droplets in drug-dominated domains. Compare these with the panels in the third column of Figure 11. These correspond to the case where where the polymer content is high (80% by weight), and we see the emergence of drug droplets in polymer-rich domains, the reverse of the low polymer content case. Polymer strings and droplet-to-string transitions. We note the formation of polymer strings in some of the panels; see the first and second columns of Figure 11 for examples. The central column in Figure 11 is of particular interest since the behaviour exhibited here is an example of a droplet-to-string transition [43]. In this droplet-to-string transition, drug droplets coalesce to 15 Figure 7: Numerical solutions superimposed on a phase diagram for the Felodipine/HPMCAS system. The numerical solutions shown here are all for a time of six months. The parameter values used can be found in Section 4.2. The dark blue regions are polymer-enriched and the dark red regions are drug-enriched. The green panels correspond to cases where phase separation is not significant. form long drug-rich strings. In the panel for 23 days, we observe that drug droplets are in the process of chaining [43]. Another droplet-to-string transition is shown in Figure 12. The formula (45) for the timescale for phase separation. The detailed numerical results here enable us to test the utility of our simple formula (45) for the timescale for phase separation. Consider, for example, the panel corresponding to 1 day in the third column of Figure 8. Here we see that polymer droplets with characteristic lengthscale of l  0:3 mm have formed. Our formula (45) predicts that such droplets should form over a timescale dictated by 2 2 (0:3) mm 11 hours jD ( = 0:8006)j eff d which is consistent with the time t = 1 day for the panel since 1 day  2 . It should be emphasized that  does not predict the time for the droplets to form, but rather estimates the timescale over which such droplets form. 5 Conclusions Solid dispersions have been the subject of intensive research in recent years because of their potential to improve the solubility of drugs, and numerous excellent studies have been published. However, de- tailed theoretical studies considering the non-equilibrium behaviour of solid dispersions are lacking. Hence, in this study we have developed a general diffusion model for a dissolving solid dispersion. We then considered the particular case of a binary system modelling a solid dispersion in storage, and developed a formula for the effective diffusion coefficient of the drug. We then specialized further to the case of a Flory-Huggins statistical model. Within the context of this theory, we make the following predictions, some of which should be testable experimentally: T = 60 C T = 90 C T = 120 C 1 day 1 hour 0.5 hours 2 days 3 hours 1 hour 1 week 4 hours 1 day 1 month 1 day 1 week 6 months 1 week 2 months Figure 8: Simulations of a FD/HPMCAS solid dispersion obtained by numerically integrating the initial boudary value problem defined in Section 4.1. The colours correspond to different mole frac- tions of the drug as defined by the colour bar. The weight fraction of drug here is 80%, and the other parameter values can be found in Section 4.2. In the above frame of figures, each column corresponds to a different temperature, and reading a column from top to bottom corresponds to increasing time for the dispersion for the given temperature. T = 60 C T = 90 C T = 120 C 1 day 0.5 hours 1.5 hours 1 week 1 hour 1 day 1 month 1 day 1 week 6 months 1 week 1 month 1 year 1 month 9 months Figure 9: See the caption for Figure 8. The weight fraction of drug here is 60%. T = 60 C T = 90 C T = 110 C 1 week 3 hours 3 hours 2 weeks 6 hours 6 hours 1 month 1 week 1 day 6 months 1 month 1 week 1 year 6 months 4 months Figure 10: See the caption for Figure 8. The weight fraction of drug here is 40%. T = 60 C T = 75 C T = 90 C 1 month 2 weeks 2 weeks 2 months 23 days 23 days 4 months 1 month 1 month 9 months 2 months 2 months 1 year 9 months 9 months Figure 11: See the caption for Figure 8. The weight fraction of drug here is 20%. 20 initial time 16 days 20 days 30 days 100 days 4 months 10 months 36 months 80 months Figure 12: Numerical results illustrating a droplet-to-chain transition. In these panels, the mass frac- tion of drug is 20% and the temperature is T = 75 C . The panels should be read from left to right, starting at the top row. In the panel for 16 days, we see the formation of drug droplets. The panels for 20 and 30 days show the drug droplets in the process of chaining to form strings. Subesequent panels show the evolution of the drug strings. 1. A solid dispersion can always be made stable by choosing a sufficiently low drug loading; see Figure 3 (a). 2. For unstable regimes, the relationship between the local drug volume fraction  and the rate of phase separation is not obvious; see Figure 3 (a). There is in fact a well-defined value of that corresponds to the most rapid rate of phase separation, with the rate decreasing for values of  either side of this value. 3. For unstable regimes, the rate of phase separation increases with increasing polymer chain length m; see Figure 3 (a). 4. Dispersions become more unstable with increasing value of the Flory-Huggins interaction pa- rameter  ; see Figure 3 (b). dp 5. Binary drug/polymer systems are capable of exhibiting a rich set of dynamical behaviours. In the numerical simulations performed in the current study, we observed the formation of polymer droplets and strings, the phase inversion phenomenon, Ostwald ripening, and droplet-to-string transitions. The model can be evaluated empirically using microscopy by comparing the theoretical simula- tions with corresponding images seen in the microscope. Hot-stage polarized light microscopy is one notable possibility - see [44] for a discussion of relevant experimental techniques. 21 There is ample scope for extending the modelling work presented in the current study. One limita- tion of the binary model considered here is that it assumes that the polymer is perfectly dry. However, if the dispersions are stored in humid conditions, this is not a good assumption since even small amounts of moisture in the dispersion may significantly affect the mobility of the drug. Another av- enue for extending the modelling work developed here is to use statistical models that capture more of the detail of the drug-polymer interaction in the dispersion; see, for example, SAFT models [35]. Vis- coelastic effects may also play a significant role in the separation process since the polymer molecules are much larger than the drug molecules in a solid dispersion, giving rise to dynamic asymmetry be- tween the components. Such models are significantly more complex than the model we have consid- ered in the current study; see [45] for some discussion of such models. Another valid critique of the current modelling is that it is incapable of distinguishing between crystalline and amorphous drug. Finally, the we have only considered the storage problem here, and have not addressed the dissolution behaviour at all. The dissolution of solid dispersions is at best partially understood, and there are many open issues that mathematical modelling may help resolve. It is noteworthy that the current study is the first (that we are aware of) that models in detail the spatiotemporal evolution of solid dispersions. Another novel feature of the current study is the development of an effective diffusion coefficient for the drug in the dispersion, the utility of which has been demonstrated in the results sections. An unusual feature of our modelling is that in equation (11) for the flux of species i, we have included the concentration c of species i. This concentration term is frequently omitted in other studies, and a compositionally dependent mobility is assumed instead – see, for example, [32] or [31]. Acknowledgments We thank the reviewers for their numerous helpful suggestions to improve the paper. The authors are grateful to S. Succi for fruitful discussions and acknowledge funding from the European Research Council under the European Unions Horizon 2020 Framework Programme (No. FP/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT). M. Meere thanks NUI Galway for the award of a travel grant. Appendix A. Linearization and stability analysis This appendix can be found in the supplementary material. References [1] K. T. Savjani, A. K. Gajjar, J. K. Savjani, Drug solubility: importance and enhancement tech- niques, ISRN Pharmaceutics 2012 (2012) 1–10. [2] H. D. Williams, N. L. Trevaskis, S. A. Charman, R. M. Shanker, W. N. Charman, C. W. Pou- ton, C. J. H. Porter, Strategies to address low drug solubility in discovery and development, Pharmacological Reviews 65 (2013) 315–499. [3] S. Kalepua, V. Nekkantib, Insoluble drug delivery strategies: review of recent advances and business prospects, Acta Pharmaceutica Sinica B 5 (5) (2015) 442–453. [4] C. Brough, R. O. W. III, Amorphous solid dispersions and nano-crystal technologies for poorly water-soluble drug delivery, International Journal of Pharmaceutics 453 (2013) 157–166. 22 [5] Y. Huang, W. Daib, Fundamental aspects of solid dispersion technology for poorly soluble drugs, Acta Pharmaceutica Sinica B 4 (1) (2014) 18–25. [6] K. Dhirendra, S. Lewis, N. Udupa, K. Atin, Solid dispersions: a review, Pakistan Journal of Pharmaceutical Sciences 22 (2) (2009) 234–246. [7] W. N. Souery, C. J. Bishop, Clinically advancing and promising polymer-based therapeutics, Acta Biomaterialia 67 (2018) 1–20. [8] C. G. Garcia, K. L. Kiick, Methods for producing microstructured hydrogels for targeted appli- cations in biology, Acta Biomaterialia 84 (2019) 34–48. [9] J. Brouwers, M. E. Brewster, P. Augustijns, Supersaturating drug delivery systems: the answer to solubility-limited oral bioavailability?, Journal of Pharmaceutical Sciences 98 (8) (2009) 2549– [10] M. Doi, Introduction to polymer physics, Oxford University Press, Oxford, UK, 1996. [11] P. J. Flory, Thermodynamics of high polymer solutions, Journal of Chemical Physics 10 (1942) 51–61. [12] M. L. Huggins, Solutions of long chain compounds, Journal of Chemical Physics 9. [13] P. C. Hiemenz, T. P. Lodge, Polymer Chemistry, 2nd Edition, CRC Press, Boca Raton, Florida, [14] J. Djuris, I. Nikolakakis, S. Ibric, Z. Djuric, K. Kachrimanis, Preparation of carba- mazepine–soluplus solid dispersions by hot-melt extrusion, and prediction of drug–polymer miscibility by thermodynamic model fitting, European Journal of Pharmaceutics and Biophar- maceutics 84 (2013) 228–237. [15] Y. Tian, V. Caron, D. S. Jones, A. M. Healy, G. P. Andrews, Using Flory–Huggins phase dia- grams as a pre-formulation tool for the production of amorphous solid dispersions: a comparison between hot-melt extrusion and spray drying, Journal of Pharmacy and Pharmacology 66 (2013) 256–274. [16] Y. Tian, J. Booth, E. Meehan, D. Jones, S. Li, G. Andrews, Construction of drug-polymer ther- modynamic phase diagram using Flory-Huggins interaction theory: identifying the relevance of temperature and drug weight fraction to phase separation within solid dispersions, Molecular Pharmaceutics 10 (2013) 236–248. [17] K. Bansal, U. S. Baghel, S. Thakral, Construction and validation of binary phase diagram for amorphous solid dispersion using Flory–Huggins theory, AAPS PharmSciTech 17 (2) (2016) 318–327. [18] S. Y. Chan, S. Qia, D. Q. M. Craig, An investigation into the influence of drug–polymer inter- actions on the miscibility, processability and structure of polyvinylpyrrolidone-based hot melt extrusion formulations, International Journal of Pharmaceutics 496 (2015) 95–106. [19] M. A. Altamimi, S. H. Neau, Use of the Flory–Huggins theory to predict the solubility of nifedipine and sulfamethoxazole in the triblock, graft copolymer soluplus, Drug Development and Industrial Pharmacy 42 (3) (2016) 446–455. [20] P. Chakravarty, J. W. Lubach, J. Hau, K. Nagapudi, A rational approach towards development of amorphous solid dispersions: Experimental and computational techniques, International Journal of Pharmaceutics 519 (2017) 44–57. 23 [21] D. Lin, Y. Huang, A thermal analysis method to predict the complete phase diagram of drug–polymer solid dispersions, International Journal of Pharmaceutics 399 (2010) 109–115. [22] K. Wlodarski, W. Sawicki, A. Kozyra, L. Tajber, Physical stability of solid dispersions with respect to thermodynamic solubility of tadalafil in pvp-va, European Journal of Pharmaceutics and Biopharmaceutics 96 (2015) 237–246. [23] Y. Zhao, P. Inbar, H. P. Chokshi, A. W. Malick, D. S. Choi, Prediction of the thermal phase diagram of amorphous solid dispersions by Flory–Huggins theory, Journal of Pharmaceutical Sciences 100 (2011) 3196–3207. [24] P. J. Marsac, S. L. Shamblin, L. S. Taylor, Theoretical and practical approaches for prediction of drug-polymer miscibility and solubility, Pharmaceutical Research 23 (10) (2006) 2417–2426. [25] J. H. Hildebrand, R. L. Scott, The solubility of nonelectrolytes, 3rd Edition, Dover Publications, New York, 1964. [26] B. D. Anderson, Predicting solubility/miscibility in amorphous dispersions: it is time to move beyond regular solution theories, Journal of Pharmaceutical Sciences 107 (2018) 24–33. [27] E. B. Smith, Basic chemical thermodynamics, sixth Edition, Imperial College Press, London, UK, 2014. [28] K. Binder, P. Fratzl, Spinodal decomposition, in: G. Kostorz (Ed.), Phase Transformations in Materials, WILEY-VCH Verlag, Weinheim, 2001, Ch. 6, pp. 411–480. [29] D. M. Saylor, C. Kim, D. V. Patwardhan, J. A. Warren, Diffuse-interface theory for structure formation and release behavior in controlled drug release systems, Acta Biomaterialia 3 (2007) 851–864. [30] D. M. Saylor, J. E. Guyer, D. Wheeler, J. A. Warren, Predicting microstructure development during casting of drug-eluting coatings, Acta Biomaterialia 7 (2011) 604–613. [31] J. Zhu, X. Lu, R. Balieu, N. Kringos, Modelling and numerical simulation of phase separation in polymer modified bitumen by phase-field method, Materials and Design 107 (2016) 322–332. [32] J. Zhu, R. Balieu, X. Lu, N. Kringos, Numerical investigation on phase separation in polymer- modified bitumen: effect of thermal condition, Journal of Materials Science 52 (11) (2017) 6525–6541. [33] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, The Journal of Chemical Physics 28 (1958) 258–267. [34] N. P. . K. Elder, Phase-field methods in material science and engineering, 1st Edition, John Wiley & Sons, Weinheim, Germany, 2010. [35] G. M. Kontogeorgis, G. K. Folas, Thermodynamic models for industrial applications, 1st Edi- tion, John Wiley & Sons, Chichester, West Sussex, UK, 2010. [36] J. M. Prausnitz, R. N. Lichtenthaler, E. G. de Azevedo, Molecular thermodynamics of fluid- phase equilibria, 3rd Edition, Prentice Hall, New Jersey, 1999. [37] D. A. Porter, K. E. Easterling, M. Sherif, Phase transformations in metals and alloys, 3rd Edition, CRC Press, Boca Raton, Florida, 2009. 24 [38] D. L. Elbert, Liquid–liquid two-phase systems for the production of porous hydrogels and hy- drogel microspheres for biomedical applications: A tutorial review, Acta Biomaterialia 7 (2011) 31–56. [39] P. Atkins, J. de Paula, Elements of physical chemistry, 5th Edition, Oxford University Press, Oxford, England, 2009. [40] E. B. Naumana, D. Q. Heb, Nonlinear diffusion and phase separation, Chemical Engineering Science 56 (2001) 1999–2018. [41] J. Gerges, Numerical studies of the physical factors responsible for the ability to vit- rify/crystallize of model materials of pharmaceutical interest, Ph.D. thesis, Universite de Lille (2015). [42] L. Ratke, P. W. Voorhees, Growth and Coarsening: Ostwald Ripening in Material Processing, 1st Edition, Springer, Berlin, 2002. [43] K. B. Migler, String formation in sheared polymer blends: Coalescence, breakup, and finite size effects, Physical Review Letters 86 (6) (2001) 1023–1026. [44] X. Liu, X. Feng, R. O. W. III, F. Zhang, Characterization of amorphous solid dispersions, Journal of Pharmaceutical Investigation 48 (1) (2018) 19–41. [45] H. Tanaka, Viscoelastic phase separation, Journal of Physics:Condensed Matter 12 (15) (2000) 207–264. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Condensed Matter arXiv (Cornell University)

Modelling phase separation in amorphous solid dispersions

Meere, Martin; Pontrelli, Giuseppe; McGinty, Sean
Condensed Matter , Volume 2019 (1902) – Feb 13, 2019
Download PDF
25 pages

Loading next page...
 
/lp/arxiv-cornell-university/modelling-phase-separation-in-amorphous-solid-dispersions-KT0B30VZvI

References (52)

ISSN
1742-7061
eISSN
ARCH-3331
DOI
10.1016/j.actbio.2019.06.009
Publisher site
See Article on Publisher Site

Abstract

Much work has been devoted to analysing thermodynamic models for solid dispersions with a view to identifying regions in the phase diagram where amorphous phase separation or drug recrystallization can occur. However, detailed partial differential equation non-equilibrium mod- els that track the evolution of solid dispersions in time and space are lacking. Hence theoretical predictions for the timescale over which phase separation occurs in a solid dispersion are not available. In this paper, we address some of these deficiencies by (i) constructing a general multi- component diffusion model for a dissolving solid dispersion; (ii) specializing the model to a binary drug/polymer system in storage; (iii) deriving an effective concentration dependent drug diffusion coefficient for the binary system, thereby obtaining a theoretical prediction for the timescale over which phase separation occurs; (iv) calculating the phase diagram for the Felodipine/HPMCAS system; and (iv) presenting a detailed numerical investigation of the Felodipine/HPMCAS sys- tem assuming a Flory-Huggins activity coefficient. The numerical simulations exhibit numerous interesting phenomena, such as the formation of polymer droplets and strings, Ostwald ripen- ing/coarsening, phase inversion, and droplet-to-string transitions. A numerical simulation of the fabrication process for a solid dispersion in a hot melt extruder was also presented. Keywords: amorphous solid dispersion, phase separation, mathematical model, drug diffusion 1 Introduction Drugs that are delivered orally via a tablet should ideally be readily soluble in water. Drugs that are poorly water-soluble tend to pass through the gastrointestinal tract before they can fully dissolve, and this typically leads to poor bioavailability of the drug. Unfortunately, many drugs currently on the market or in development are poorly water-soluble, and this presents a serious challenge to the pharmaceutical industry. Many strategies have been developed to improve the solubility of drugs, such as the use of surfactants, cocrystals, lipid-based formulations, and particle size reduction. The literature on this topic is extensive, and recent reviews can be found in [1–3]. Corresponding author. Email: martin.meere@nuigalway.ie arXiv:1902.05410v2 [cond-mat.soft] 25 Jun 2020 One particularly effective strategy to improve drug solubility is to use a solid dispersion [4–6]. A solid dispersion typically consists of a hydrophobic drug embedded in a hydrophilic polymer [7, 8] matrix, where the matrix can be either in the amorphous or crystalline state. The drug is preferably in a molecularly dispersed state, but may also be present in amorphous particles or even in the crystalline form (though this is usually undesirable); see Figure 1. The drug release concept for most solid dis- persions is based on the so-called spring and parachute effect [9]. When the drug and the hydrophilic polymer dissolve in solution, a supersaturated drug solution is quickly created (the spring). Although the drug concentration then subsequently decreases, the rate of decrease is slowed by drug-polymer interactions in the dispersion, so that the drug can be present at supersaturated levels in the solution for a period of some hours (the parachute). This results in improved bioavailability of the drug when the solid dispersion dosage form is taken orally. Molecularly dispersed drug in the polymer matrix (desirable) Drug molecule Polymer matrix Crystal drug formation Contains amorphous has occurred (undesirable) drug-rich domains Figure 1: Adapted from [5]. In this figure, we show three possible structures for a polymer/drug dispersion. Top: Here the drug is in the molecularly dispersed state, which is usually desirable for a solid dispersion. Bottom left: Here the dispersion contains drug in the crystalline form. Bottom right: Here the dispersion contains amorphous drug-rich domains. Drug loading in most dispersions greatly exceeds the equilibrium solubility in the polymer matrix for typical storage temperatures. Hence these systems are usually unstable, with phase separation eventually occurring [5]. In such cases, the drug will eventually crystallise out or form an amorphous phase separation. However, if the dispersion is stored well below the glass transition temperature [10] for the polymer, and is kept dry, this can happen extremely slowly. The system is then for all practical purposes stable, and is said to be metastable. The humidity of the storage environment can be an issue because even small amounts of moisture can significantly affect the glass transition temperature. Hence polymers that have high glass transition temperatures and that are resistant to water absorption have become popular. An example of one such polymer is Hydroxypropyl Methylcellulose Acetate Succinate (HPMCAS). Phase separation of solid dispersions in storage is clearly undesirable from the point of manufac- turers. Hence much work has been devoted to constructing phase diagrams for solid dispersions with a view to identifying regimes where drug recrystallization or amorphous phase separation can occur. These phase diagrams are constructed with the aid of thermodynamic models. The most widely used thermodynamic model in this context is the Flory-Huggins model [11–13] for polymer solutions. Flory-Huggins theory is a lattice-based model in which the drug and polymer are confined to live 2 on a regular lattice. Flory-Huggins theory is an extension of regular solution theory, as explained in Chapter 7 of [13]. In the context of a drug/polymer system, each drug molecule is taken to occupy one lattice site and each polymer segment is taken to occupy m  1 sites. Under a number of further simplifying assumptions [13], the change in entropy and enthalpy associated with the mixing of the mix polymer and drug are calculated. With these in hand, the change in bulk Gibbs free energy (g ) per mole associated with mixing is readily calculated, and is found to be mix = X ln( ) + X ln( ) +  X  ; (1) d d p p dp d p RT where R is the gas constant, T is the temperature, X ; X = 1X are the mole fractions of the drug d p d and polymer, respectively, and  ;  are the volume fractions of the drug and polymer, respectively. d p The quantity  is referred to as the Flory-Huggins interaction parameter, and it is discussed further dp below. The mole fractions and volume fractions are related via the formulae X mX d p = ;  = : (2) d p X + mX X + mX d p d p When the model is applied to real binary systems, m can be calculated using the formula m = (3) where V , V (molar ) are the molar volumes of the polymer and drug, respectively. p d The mixing of the polymer and drug is spontaneous if G < 0. The Flory-Huggins parameter mix takes the form dp w + w dd pp =  w (4) dp dp where  is a positive parameter, and w , w , w give measures of the drug-polymer, drug-drug and dp dd pp polymer-polymer interaction energy, respectively. If  < 0 then w < (w + w )=2 indicating dp dp dd pp that that the mixed state has lower energy than the separated pure states, so that mixing is favoured. Conversely,  > 0 is indicative of demixing being favoured. However, these statements are indica- dp tive rather than precise, as will be explained in Section 3. We should also note that  is temperature dp dependent, and is usually given the empirical form (T ) = + (5) dp where ; are constants. Flory-Huggins theory has frequently been used to analyse the stability of binary solid dispersion systems in storage; see, for example, [14–23]. In many of these studies, the Flory-Huggins interaction parameter is first estimated using the melting point depression method [24], or using the Hildebrand and Scott method [25], which involves the estimation of solubility parameters. Once estimates for (T ) have been obtained, the Gibbs free energy of mixing G can be calculated, which in turn dp mix enables the construction of phase diagrams for the systems. Phase diagrams assist with the identifi- cation of regions in composition-temperature space where the system is prone to recrystallization or amorphous phase separation. The models we shall develop in the current study are generic and are not tied to making a specific choice of statistical model. However, given the particular importance of Flory-Huggins theory in applications, we shall derive detailed results for this case. Also, all of our numerical illustrations are calculated within the context of Flory-Huggins theory. It should be emphasized that Flory-Huggins theory does involve quite a number of simplifying assumptions which are not appropriate for some systems; see [26] for a recent critique of the model. 3 2 Theoretical formulation 2.1 A multicomponent diffusion model for solid dispersions We develop a multicomponent diffusion model for the evolution of the concentrations of the com- ponents constituting a solid dispersion. We suppose for the moment that there are p components. However, in the analysis we shall consider in the current study, we will in fact have p = 2, with one of the components being the polymer, and the other being the drug. For a dissolving solid dispersion, there are three components p = 3: the polymer, the drug, and the solvent. The chemical potential  (J/mole) of species i (i = 1; 2; :::; p) gives the Gibbs free energy per mole of species i, and is given here by ( [27]) b 2 2 =   r X (6) i i i i where b 0 =  + RT ln(a ) (7) i i and where  is the bulk chemical potential of species i,  is the chemical potential of species i i0 2 2 in the pure state, a is the activity of species i, and the term involving  > 0 (m J/mole) penalises the formation of phase boundaries ( [28], [29], [30], [31], [32]). The parameters  are referred to as gradient energy coefficients ( [33], [34]). Here X is the molar fraction of species i (i = 1; 2; :::; p), and the activities can depend on these molar fractions, so that a = a (X ; X ; :::; X ): i i 1 2 p The molar fraction is related to the molar concentration via X = V c (8) i i i where V (molar ) is the molar volume of species i. The flux of species i (molarm/s) is given by J = c v (9) i i i where c (molar), v (m/s) give the molar concentration and drift velocity, respectively, of species i. i i The drift velocity v gives the average velocity a particle of species i attains due to the diffusion force acting on it, and is given here by v = M F = M r (10) i i i i i where M (moles/kg),F (J/[mmole]) give the mobility and diffusion force, respectively, for species i i i. Equations (9) and (10) give J = M c r : (11) i i i i Conservation of mass for species i implies that @c +r J = 0 (12) @t and using (11) now gives @c = r (M c r ) i i i @t or equivalently @X i i = r D X r (13) i i @t RT with 2 2 = ln(a )  r X (14) i i RT 4 2 2 2 for i = 1; 2; :::; p, and where  =  =RT > 0 (m /molar), and i i D = M RT (Einstein relation) i i is the self-diffusion coefficient for species i. The model formulation given by (13) and (14) based on chemical potentials will be used for the numerical scheme described in Section 4. However, it is also of value to develop a formulation involv- ing diffusion coefficients since these yield immediate information regarding timescales for transport processes, and will also the enable the development of analytical results via a linearization process. Diffusion Coefficients Using (6), (7) and (11) gives RT 2 2 J = M c r = M c ra  r(r X ) i i i i i i i i and then using the fact that the activities depend on the molar fractions gives RT @a 2 2 J = M c rX  r(r X ) : (15) i i i j i a @X i j j=1 Using (8), we can now write (15) as 2 2 J = D rc + D " c r(r c ) (16) i ij j i i i j=1 2 2 2 where " = V  and where the diffusion coefficients D (m /s) are given by i ij i i V X @a j i i D = D : i; j = 1; 2; :::; p (17) ij i V a @X i i j Conservation of mass (12) then implies that (reverting to molar fractions) @X V i i 2 2 = r D (X)rX D  X r(r X ) i = 1; 2; :::; p (18) ij j i i i @t V j=1 where X = (X ; X ; :::; X ), and where we have included the concentration dependence of the dif- 1 2 p fusion coefficients D here to emphasise that this system is in general a coupled system of nonlinear ij diffusion equations. It should benoted that the equations (18) are not independent since X = 1, i=1 and so it is sufficient to solve for p 1 concentrations only. 2.2 Activity coefficients The activities a are usually written as a = X i i i where the = (X ; X ; :::; X ) are referred to as activity coefficients. Equations (17) now give i i 1 2 p V X @ j i i D = D  + i; j = 1; 2; :::; p (19) ij i ij V @X i i j where  is the Kronecker delta. ij The details of the interactions between the species in solution are captured in the modelling by choosing appropriate forms for the activity coefficients = (X ; X ; :::; X ). The construction of i i 1 2 p appropriate forms for the for various solutions is a large subject with a large literature; see, for example, the books [35] and [36]. 5 2.3 The storage problem for a binary mixture In the current study, we shall be modelling the behaviour of solid dispersions in storage. In this case, we have p = 2, with the label 1 referring to the drug and the label 2 referring to the polymer. However, for transparency, we choose here to use the labels d; p rather than 1; 2, where d stands for drug, and p for polymer. Then using (18) and the fact that X = 1 X , we have p d @X 2 2 = r D (X )rX D  X r r X : (20) eff d d d d d @t where the effective concentration-dependent diffusion coefficient for the drug in the solid dispersion is X @ @ d d d D (X ) = D (X ) V D (X )=V = D 1 + : (21) eff d dd d d dp d p d @X @X d d p For the particular case of a binary Flory-Huggins theory (see Section 1), the activity coefficients are given by d d ln( ) = ln + 1 +   ; (22) d dp X X d d p p ln( ) = ln + 1 + m  ; (23) p dp X X p p where the volume fractions  ;  are given by (2). Substituting (22) in (21) gives d p 2 2 2 (m (m 1)X )(m (m 1)X ) 2 m X (1 X ) d d dp d d D (X ) = D : (24) eff d d (m (m 1)X ) It is more instructive to write this expression in terms of the volume fraction of drug. Writing D (X ) = D ( ), we obtain eff d eff d D ( ) = D (1 + (m 1) ) 1 + 1  2  (1  ) : (25) eff d d d d dp d d This expression is particularly useful because it yields insight into how the mobility of the drug in the dispersion depends on the length of the polymer chains (m), the dispersion composition ( ), and the character of the drug-polymer interaction ( ). We shall analyze this expression further in Section 3, dp and also show how it can be used to calculate the timescale over which phase separation may occur. An equivalent formulation for the Flory-Huggins model involving the chemical potential for the drug  is given by (see (13) and (14) above): @X = r (D X r ) (26) d d @t where d d0 = (27) RT and with X (m 1)(1 X ) 1 X d d d 2 2 2 = ln + +  m  r X : (28) dp d m (m 1)X m (m 1)X m (m 1)X d d d We suppose that the solid dispersion occupies a two-dimensional region . The governing equa- tion for the drug concentration in may be written in the conservation form @X +r J = 0; @t 6 where the drug flux J is given by 2 2 J = D (X )rX + D  X r(r X ): (29) d eff d d d d d We need to supplement the governing equation in with boundary conditions on @ , and we choose these here to be J  n = 0; rX  n = 0 on @ : (30) d d The first of these conditions J  n = 0 implies that the drug cannot penetrate the boundary of the domain. The other condition rX  n = 0 is the natural boundary condition for the variational formulation of the problem, and it implies that the interfaces between the polymer rich and the drug rich domains meet the boundary at right angles. Finally, to obtain a well-posed problem, we need to impose an initial condition and we choose this here to take the form X (x; y; t = 0) = X (x; y) for (x; y) 2 ; (31) where X (x; y) is a given function. Gathering together the governing equation, the boundary conditions and the initial condition, we obtain the following initial boundary value problem: @X 2 2 = r D (X )rX D  X r r X in eff d d d d d @t rX  n = 0; r(r X ) n = 0 on @ ; (32) d d X (x; y; t = 0) = X (x; y) for (x; y) 2 2.4 Phase separation in a Flory-Huggins binary mixture The bulk free energy and spinodal decomposition Spinodal decomposition for binary systems has been long understood using thermodynamic rea- soning, and is well described elsewhere; see, for example, Chapter 5 of [37], Chapter 7 of [13], or [38]. Hence our description of the background theory here will be quite brief, and we will empha- sise instead the particular details for the Flory-Huggins system. The bulk free energy density g for the binary mixture constituting the solid dispersion is given by ( [39]) b b g =  X +  X ; (33) b d p d p b b where  ;  give the bulk chemical potential of the drug and polymer, respectively, and where d p b 0 =  + RT ln(a ) for i = d; p: i i This leads to 0 0 g =  X +  X + RT (X ln(X ) + X ln(X )) + RT (X ln( ) + X ln( )): b d p d d p p d d p p d p mix The Gibbs free energy of mixing g is given by mix 0 0 g = g  X  X = RT (X ln(X ) + X ln(X )) + RT (X ln( ) + X ln( )): b d p d d p p d d p p b d p If we now use (22) and (23) and the fact that X = 1 X , we arrive at p d mix g X m(1 X )  mX (1 X ) d d dp d d = X ln + (1 X ) ln + : (34) d d RT m (m 1)X m (m 1)X m (m 1)X d d d 7 mix Figure 2: (a) Plot of the bulk free energy of mixing g as a function of a drug molar fraction 1s 2s 2 mix 2 X . The spinodal points X , X are the solutions to d (g )=dX = 0. In the spinodal region d d b d 1s 2s 2 mix 2 (X ; X ), we have d (g )=dX < 0 and D (X ) < 0. (b) Phase diagram for the binary eff d d d b d mixture. Here is the coexistence curve, is the spinodal curve, T is the temperature for the free energy density diagram in (a), and T is the critical temperature above which the dispersion is homogeneous. mix In Figure 2 (a), we plot a free energy of mixing diagram g as a function of drug molar fraction 1s 2s X . In this diagram, the points X , X are the solutions to d d 2 mix d (g ) = 0; dX 1s 2s and are referred to as the spinodal points. The region (X ; X ) is referred as the spinodal region, d d and for points X in this region, we have 2 mix d (g ) < 0: dX Compositions X in the spinodal region are unstable, and will split into two phases characterized 1u 2u by the compositions X and X as shown in Figure 2 (a); see [37] for more details. The points d d 1u 2u X , X are referred to as the binodal points, and are defined by the common tangent construction d d shown in Figure 2 (a). The binodal and spinodal points define the coexistence and spinodal curves, respectively, and these are plotted in the phase diagram shown in Figure 2 (b). Using equation (34), we obtain 2 mix d (g ) q(X ) = RT (35) 2 3 dX (1 (1 1=m)X ) X (1 X ) d d d where q(X ) = AX + BX + 1 (36) d d and where 1 1 1 1 1 A = (1 2 ) + 1; B = + (1 2 ) 2: (37) dp dp 3 2 2 m m m m m 2 mix 2 Hence there is a spinodal region with d (g )=dX < 0 if q(X ) < 0 in this region. Inspecting b d (36), we see that q(X ) can be negative if q(X ) = 0 has real roots, that is, if d d B 4A > 0; 8 and using (37), this leads to (2 (1 + 1=m)) 4=m > 0 dp which holds true if 2 q 2 1 1 1 > 1 + p = 1 + V =V : dp d p 2 m 2 Hence, we have a spinodal interval if >  (m) (38) dp dp where 1 1 (m)  1 + ; (39) dp 2 m and where  (m) is a critical value for the Flory-Huggins parameter. If (38) holds true, then there is dp 1s 1s a spinodal interval (X ; X )  [0; 1] where d d p p 2 2 B B 4A B + B 4A 1s 2s X = ; X = ; (40) d d 2A 2A and where A, B are given in (37). The diffusion coefficient and spinodal decomposition Using (24) and (35), elementary calculations show that 2 mix d (g ) y b D (X ) = M (X ) (41) eff d d dX where M (X ) = M X (1 X ) d d d d with M = D =RT , and where M (X ) is a concentration-dependent drug mobility; see, for example, d d d 2 mix 2 equation (3.6) of the paper [40]. Hence, for 0 < X < 1, it is clear from (41) that d (g )=dX < 0 b d implies that D (X ) < 0: (42) eff d Hence, an equivalent criterion for spinodal decomposition to occur is that there exist a region in 0  X  1 where D (X ) < 0, that is, that there exist a region where drug diffusion is against the d eff d concentration gradient (uphill diffusion). 3 Qualitative results and discussion Although the model we have derived in the current study is quite general, and is not tied to any specific statistical model for a solid dispersion, the detailed results we shall present in this section are for the Flory-Huggins case. 3.1 The effective diffusion coefficient for the drug in the dispersion From (25), the scaled effective diffusion coefficient for the drug in the dispersion is given by D ( ) 1 eff d = (1 + (m 1) ) 1 + 1  2  (1  ) (43) d d dp d d D m where we recall that D is the temperature-dependent self-diffusion coefficient for the drug. Equation (43) is of particular value since it yields information on how the mobility of the drug in the dispersion 9 Figure 3: Plots of the scaled effective diffusion coefficient for the drug in the polymer dispersion as a function of the drug volume fraction. Here positive values of the diffusion coefficient correspond to standard drug diffusion down the concentration gradient, while negative values correspond to phase separation of the drug and the polymer, with larger negative values (in absolute terms) corresponding to more rapid phase separation. We have plotted the scaled drug diffusion coefficient for (a) the Flory- Huggins interaction parameter  = 3 and various values of the polymer chain length m, and, (b) dp polymer chain length m = 50 and various values of the Flory-Huggins interaction parameter  . See dp the main body of the text for further discussion. depends on the polymer chain length, the dispersion composition, and the character of the drug- polymer interaction. In Figure 3 (a), we have plotted (43) for the Flory-Huggins interaction parameter  = 3 (which dp is in the unstable regime) and various values of the polymer chain length m. It should be emphasized that positive values for D correspond to standard drug diffusion down the concentration gradient, eff while negative values correspond to unstable regimes where phase separation of the drug and polymer can occur. In Figure 3 (a), it is clear that if the drug loading  is sufficiently low, then D > 0 and d eff the solid dispersion is stable. However, for larger (and more realistic) drug loadings, D < 0, and eff the system is unstable. It is interesting to note that the system becomes more unstable as the length of the polymer chains increase. It is also clear from the curves in Figure 3 (a) that the relationship between the initial drug loading in the dispersion and the initial rate of phase separation is not altogether obvious. It is not necessarily the case that increasing drug loading corresponds to increasing initial dispersion instability. Rather, there is in fact a well defined worst choice for the initial drug loading from the point of view of stability in the initial stages. This worst choice corresponds to the minima of the curves displayed in Figure 3 (a), since these minima correspond to the fastest rates of phase separation. For m  1, the minimum of D ( ) occurs at eff d 1 + 2 + (1 + 2 ) 6 dp dp dp min : (44) dp min These theoretical results predict that choosing initial drug loadings  above or below  should min lead to improved dispersion stability in the initial stages. For m;   1, we have   0:67. dp In Figure 3 (b), we plot (43) for the fixed polymer length m = 50, and various values of the Flory- Huggins interaction parameter  . For m = 50, the critical value for  is given by   0:5707 dp dp dp 10 c (see equation (39)). Recall that for  <  , the system is stable for all drug loadings  , and that dp d dp for  >  , there is a regime of unstable drug loadings. This is borne out by the curves displayed dp dp in Figure 3 (b). These curves predict that the system becomes more unstable with increasing values of  , and this is as expected given the dependence of  on the interaction energies - see equation dp dp (4). Figure 4: Schematic of a phase separating solid dispersion where polymer-rich regions with charac- teristic lengthscale l have formed. A formula for the timescale of evolution of such a dispersion is given in the main body of the text; see equation (45). 3.2 Timescale for phase separation in a solid dispersion In Figure 4, we give a schematic of a phase separating solid dispersion where polymer-rich regions have formed. The characteristic lengthscale of these regions is denoted by l. In order for such regions to form, the drug must have diffused away over a lengthscale of order l, and the timescale over which this diffusion occurs is estimated by (see (25)) 2 2 l l 1 = = (45) 0 0 0 0 D (T )j(1 + (m 1) )[1 + (1=m 1) 2 (T ) (1  )]j jD ( )j d dp eff d d d d where  is the initial uniform volume fraction of the drug in the dispersion, and T is a representative storage temperature. It should be emphasized that this formula is just an estimate since, in reality, the drug volume fraction evolves in space and time. Hence, (45) should only be used as a rough rule of thumb. In Section 4, we evaluate this formula by comparing it with detailed numerical results, and satisfactory agreement is generally found. Equation (45) may, in appropriate circumstances, be used to estimate the shelf life of a solid dispersion product. To see this, suppose that l denotes the largest acceptable size for polymer-rich domains (or drug-rich domains) in the product. Then, since  estimates the timescale for these regions to form, it also estimates the timescale for the shelf life of the product. However, care should be taken when using (45) since, apart from the fact that is based on a fixed value of  , it also incorporates a number of significant assumptions - for example, it assumes that the dispersion is perfectly dry, and that Flory-Huggins theory is an appropriate statistical model for the system. 11 3.3 Criteria for a stable solid dispersion Although the drug loading in real solid dispersions is typically high and in the unstable regime, it is nevertheless worthwhile specifying conditions under which the stability of the dispersion is guaran- teed. The results we display here are based on the discussion given in Section 2.4. For  < dp dp where  = (1 + 1= m) , the system is stable irrespective of the choice of the uniform initial dp 0 c 0 drug load  . For  >  , the dispersion is unstable if the initial drug loading  is chosen in the dp d dp d + + interval ( ;  ), but stable if chosen in either of the intervals (0;  ) or ( ; 1), where d d d d 8 9 < 2 = 1 1 1 1 1 2 = 1 + 1  1 + 1 : 2 : 2 m 2 m  ; dp dp dp These results are based on the bulk free energy only, and do not take account of interfacial energy. However, the interfacial energy can be readily incorporated into the analysis, and this is discussed in Appendix A. 4 Numerical results and discussion 4.1 The numerical method For the purposes of numerical calculations, we take the integration domain to be the square region = f(x; y)j 0 < x < L; 0 < y < Lg with boundary @ . The governing equation to be solved is defined by the equations (26), (27) and (28). The boundary conditions are given by r  n = 0 and rX  n = 0 on @ ; (46) and the initial condition takes the form (31). The boundary conditions (46) are equivalent to those given in (30). The governing equation was numerically integrated using the finite element package COMSOL Multiphysics. A mesh sensitivity analysis was performed to investigate the influence of the size of the mesh on the results. The solution was assumed to be mesh independent when there was less than 1% difference in the mole fraction of drug between successive refinements. The final mesh used in the simulations was triangular and consisted of 7553 vertices and 14796 triangles. The numerical solutions all conserved the total mass of drug in the system to within 1%. 4.2 Parameter values We consider parameter values that are appropriate for a solid dispersion consisting of the drug Felodip- ine (FD) and the polymeric excipient HPMCAS. Felodipine is a calcium channel blocker that is com- monly used to treat blood pressure. For this system, the Flory-Huggins interaction parameter is given as a function of temperature by (see [16]) 7830:4 (T ) = 18:767 + : (47) dp Using data taken from [16], the molar volume for FD is V = 300:19 cm /mol and the molar volume of HPMCAS is V = 14007:78 cm /mol, so that V 14007:78 m = =  46:6630: V 300:19 From (39), the critical value for the interaction parameter below which phase separation cannot occur is given by 1 1 (m) = 1 + p = 0:6571: dp 2 m 12 o 2 1 2 1 T ( C)  (T ) D (T ) (m s ) D (X ) (m s ) dp d eff a 18 17 40 6.2383 1:1661 10 8:8605 10 17 15 50 5.4645 1:5494 10 1:0113 10 16 15 60 4.7371 1:3787 10 7:6107 10 15 14 75 3.7245 2:0297 10 8:3587 10 14 13 90 2.7954 1:7356 10 4:9151 10 14 12 100 2.2176 5:7436 10 1:1670 10 13 12 110 1.6699 1:6336 10 2:0806 10 13 12 120 1.1501 4:0902 10 2:2657 10 Table 1: Illustrative values for some of the parameters of the FD/HPMCAS system at various tem- peratures. Here the initial weight fraction of drug is 70%, which corresponds to an initial drug molar fraction X  0:9909. The self-diffusion coefficient for Felodipine was estimated in [41] (Chapter 4, page 133) to be A A 2 3 2 1 D (T ) = exp(A ) exp exp m s (48) d 1 T T where A = 18:03, A = 445:84 K, A = 874:81 K. Some illustrative values for the diffusion 1 2 3 coefficients and the Flory-Huggins interaction parameter are displayed in Table 1. For the numerical simulations displayed in the current study, we take the size of the square domain to be given by L = 2 mm. The thickness of the interfacial regions is dictated by the parameter  , and here we chose the value  = L=50 = 4 10 m. We illustrate how the initial conditions were specified by considering a particular case. We consider the case where the initial weight fraction of drug is 80%. This means that the initial weight of FD divided by the weight of FD plus the weight of HPMCAS is 0.8. This corresponds to an initial molar drug fraction of X = 0:9947. More precisely, we choose the initial molar fraction of drug to be a small random perturbation about this level given by X (x; y; t = 0) = 0:9947(1 + rnd(x; y)) where rnd(x; y) is a normally distributed random function with a mean value of zero and a standard 5 5 deviation of 10 . The standard deviation for all of the initial conditions was taken to be 10 , with one exception - the numerical results displayed in Figure 6 took the larger value 0:007 to simulate a coarse initial mixture in a hot melt extruder. In Figure 5, we plot the phase diagram for the Felodipine/HPMCAS system. All that is required to calculate the phase diagram here is a knowledge of the Flory-Huggins interaction parameter, and this has been given in (47). The spinodal curve T ( ) is obtained by setting D ( ) = 0 in (43) to s d eff d obtain 1 + 1  2  (1  ) = 0; (49) d dp d d where  is given by (5). Solving (49) for T gives the spinodal curve dp 2  (1  ) d d T ( ) = s d 1 + 1  2  (1  ) d d d where = 18:767, = 7830:4, m = 46:6630 for the Felodipine/HPMCAS system. The binodal curve was estimated numerically. The calculation involved simultaneously solving the pair of equa- b 1u b 2u b 1u b 2u 1u 2u 1u 2u tions  (X ) =  (X ) and  (X ) =  (X ) for the binodal points X , X with X < X . p p d d d d d d d d d d This was implemented using the the fsolve command in MAPLE. 13 Figure 5: The phase diagram for the Felodipine/HPMCAS system. 4.3 Numerical results The first numerical calculations we display simulate the hot melt mixing process for a Felodip- ine/HPMCAS system. To implement this, we used a time dependent temperature profile that treats o o the case where the mixture begins at 25 C and rises in temperature at a rate of 10 C per minute until o o it reaches 145 C. We then assumed for simplicity that the melt cooled linearly back down to 25 C over a period of 30 minutes. The results of the calculation are shown in Figure 6. We see in Figure 6 (a) that the initial drug/polymer mixture is quite coarse (badly mixed). As the temperature rises, we see from Figure 6 (b) that the mixture becomes increasingly homogenous. By the time the mixture has achieved its maximum temperature, Figure 6 (c), it is quite well-mixed. Figure 6 (d) shows that the amount of phase separation that occurs during the cooling process is insignificant. Hence, the numerical results shown here predict a successful hot melt extrusion process for the manufacture of a solid dispersion. We note that different heating and cooling regimes are easily simulated using the model. In Figure 7, we superimpose numerical solutions on a phase diagram for the Felodipine/HPMCAS system. Each of these solutions corresponds to an evolution time of 6 months, with the dispersion mixture beginning from an initially approximate uniform state. We see that there is no significant phase separation for the 30 C degree cases, and for the cases in the metastable and stable regions. This predicts that Felodipine/HPMCAS systems should not suffer considerable phase separation un- der normal storage temperatures. However, we should caution that we are modelling the case of zero relative humidity here. For the cases that do exhibit significant phase separation, we note a coarser separation morphology for higher temperatures. In these figures, dark red corresponds to drug-enriched domains (relative to the initial concentration) while dark blue correspond to polymer- enriched domains. We also note the occurrence of polymer droplets and strings - we return to this issue below. Further numerical simulations are displayed in Figures 8, 9, 10, 11, and these correspond to weight fractions of drug of 80%, 60%, 40% and 20%, respectively. Recall that decreasing weight fractions of drug correspond to increasing weight fractions of polymer since the system is binary. In a given figure, each column corresponds to a given temperature as labelled, and reading a column from top to bottom corresponds to increasing time for the dispersion for the given temperature. We have chosen here not to use the same times for the different temperatures since the rate at which a dispersion evolves depends on temperature. We now highlight some notable features of these numerical simulations. 14 Figure 6: Numerical simulation of the behaviour of a solid dispersion in a hot melt extruder. The simulations here are for a FD/HPMCAS solid dispersion and were obtained by numerically integrating the initial boudary value problem defined in Section 4.1. The colours correspond to different mole fractions of the drug as defined by the colour bar. The weight fraction of drug here is 50%, and the other parameter values can be found in Section 4.2 and Section 4.3. These simulations represent a succesful extrusion where an initially coarse mixture is heated and then cooled to form a well-mixed dispersion. Two phases eventually emerge. The numerical results show that the systems eventually evolve into two distinct phases, characterized by deep blue domains (polymer-rich) and deep red do- mains (drug-rich). Ostwald ripening/coarsening. Another notable feature in many of the numerical illustrations is the formation of polymer droplets (blue discs) in the dispersion, followed by a subsequent growth in their size; see, for example, the third column in Figure 8. This is a well-known and common phenomenon in multicomponent solid systems, and is often referred to as Ostwald ripening or coarsening [42]. We also note the general trend that dispersions at higher tempera- ture tend to be coarser. Phase inversion. The system exhibits the phase inversion phenomenon [31] as the polymer content increases. To see this, consider the panels in Figure 8. These correspond to the case where the polymer content is low (20% by weight), and we see the emergence of polymer droplets in drug-dominated domains. Compare these with the panels in the third column of Figure 11. These correspond to the case where where the polymer content is high (80% by weight), and we see the emergence of drug droplets in polymer-rich domains, the reverse of the low polymer content case. Polymer strings and droplet-to-string transitions. We note the formation of polymer strings in some of the panels; see the first and second columns of Figure 11 for examples. The central column in Figure 11 is of particular interest since the behaviour exhibited here is an example of a droplet-to-string transition [43]. In this droplet-to-string transition, drug droplets coalesce to 15 Figure 7: Numerical solutions superimposed on a phase diagram for the Felodipine/HPMCAS system. The numerical solutions shown here are all for a time of six months. The parameter values used can be found in Section 4.2. The dark blue regions are polymer-enriched and the dark red regions are drug-enriched. The green panels correspond to cases where phase separation is not significant. form long drug-rich strings. In the panel for 23 days, we observe that drug droplets are in the process of chaining [43]. Another droplet-to-string transition is shown in Figure 12. The formula (45) for the timescale for phase separation. The detailed numerical results here enable us to test the utility of our simple formula (45) for the timescale for phase separation. Consider, for example, the panel corresponding to 1 day in the third column of Figure 8. Here we see that polymer droplets with characteristic lengthscale of l  0:3 mm have formed. Our formula (45) predicts that such droplets should form over a timescale dictated by 2 2 (0:3) mm 11 hours jD ( = 0:8006)j eff d which is consistent with the time t = 1 day for the panel since 1 day  2 . It should be emphasized that  does not predict the time for the droplets to form, but rather estimates the timescale over which such droplets form. 5 Conclusions Solid dispersions have been the subject of intensive research in recent years because of their potential to improve the solubility of drugs, and numerous excellent studies have been published. However, de- tailed theoretical studies considering the non-equilibrium behaviour of solid dispersions are lacking. Hence, in this study we have developed a general diffusion model for a dissolving solid dispersion. We then considered the particular case of a binary system modelling a solid dispersion in storage, and developed a formula for the effective diffusion coefficient of the drug. We then specialized further to the case of a Flory-Huggins statistical model. Within the context of this theory, we make the following predictions, some of which should be testable experimentally: T = 60 C T = 90 C T = 120 C 1 day 1 hour 0.5 hours 2 days 3 hours 1 hour 1 week 4 hours 1 day 1 month 1 day 1 week 6 months 1 week 2 months Figure 8: Simulations of a FD/HPMCAS solid dispersion obtained by numerically integrating the initial boudary value problem defined in Section 4.1. The colours correspond to different mole frac- tions of the drug as defined by the colour bar. The weight fraction of drug here is 80%, and the other parameter values can be found in Section 4.2. In the above frame of figures, each column corresponds to a different temperature, and reading a column from top to bottom corresponds to increasing time for the dispersion for the given temperature. T = 60 C T = 90 C T = 120 C 1 day 0.5 hours 1.5 hours 1 week 1 hour 1 day 1 month 1 day 1 week 6 months 1 week 1 month 1 year 1 month 9 months Figure 9: See the caption for Figure 8. The weight fraction of drug here is 60%. T = 60 C T = 90 C T = 110 C 1 week 3 hours 3 hours 2 weeks 6 hours 6 hours 1 month 1 week 1 day 6 months 1 month 1 week 1 year 6 months 4 months Figure 10: See the caption for Figure 8. The weight fraction of drug here is 40%. T = 60 C T = 75 C T = 90 C 1 month 2 weeks 2 weeks 2 months 23 days 23 days 4 months 1 month 1 month 9 months 2 months 2 months 1 year 9 months 9 months Figure 11: See the caption for Figure 8. The weight fraction of drug here is 20%. 20 initial time 16 days 20 days 30 days 100 days 4 months 10 months 36 months 80 months Figure 12: Numerical results illustrating a droplet-to-chain transition. In these panels, the mass frac- tion of drug is 20% and the temperature is T = 75 C . The panels should be read from left to right, starting at the top row. In the panel for 16 days, we see the formation of drug droplets. The panels for 20 and 30 days show the drug droplets in the process of chaining to form strings. Subesequent panels show the evolution of the drug strings. 1. A solid dispersion can always be made stable by choosing a sufficiently low drug loading; see Figure 3 (a). 2. For unstable regimes, the relationship between the local drug volume fraction  and the rate of phase separation is not obvious; see Figure 3 (a). There is in fact a well-defined value of that corresponds to the most rapid rate of phase separation, with the rate decreasing for values of  either side of this value. 3. For unstable regimes, the rate of phase separation increases with increasing polymer chain length m; see Figure 3 (a). 4. Dispersions become more unstable with increasing value of the Flory-Huggins interaction pa- rameter  ; see Figure 3 (b). dp 5. Binary drug/polymer systems are capable of exhibiting a rich set of dynamical behaviours. In the numerical simulations performed in the current study, we observed the formation of polymer droplets and strings, the phase inversion phenomenon, Ostwald ripening, and droplet-to-string transitions. The model can be evaluated empirically using microscopy by comparing the theoretical simula- tions with corresponding images seen in the microscope. Hot-stage polarized light microscopy is one notable possibility - see [44] for a discussion of relevant experimental techniques. 21 There is ample scope for extending the modelling work presented in the current study. One limita- tion of the binary model considered here is that it assumes that the polymer is perfectly dry. However, if the dispersions are stored in humid conditions, this is not a good assumption since even small amounts of moisture in the dispersion may significantly affect the mobility of the drug. Another av- enue for extending the modelling work developed here is to use statistical models that capture more of the detail of the drug-polymer interaction in the dispersion; see, for example, SAFT models [35]. Vis- coelastic effects may also play a significant role in the separation process since the polymer molecules are much larger than the drug molecules in a solid dispersion, giving rise to dynamic asymmetry be- tween the components. Such models are significantly more complex than the model we have consid- ered in the current study; see [45] for some discussion of such models. Another valid critique of the current modelling is that it is incapable of distinguishing between crystalline and amorphous drug. Finally, the we have only considered the storage problem here, and have not addressed the dissolution behaviour at all. The dissolution of solid dispersions is at best partially understood, and there are many open issues that mathematical modelling may help resolve. It is noteworthy that the current study is the first (that we are aware of) that models in detail the spatiotemporal evolution of solid dispersions. Another novel feature of the current study is the development of an effective diffusion coefficient for the drug in the dispersion, the utility of which has been demonstrated in the results sections. An unusual feature of our modelling is that in equation (11) for the flux of species i, we have included the concentration c of species i. This concentration term is frequently omitted in other studies, and a compositionally dependent mobility is assumed instead – see, for example, [32] or [31]. Acknowledgments We thank the reviewers for their numerous helpful suggestions to improve the paper. The authors are grateful to S. Succi for fruitful discussions and acknowledge funding from the European Research Council under the European Unions Horizon 2020 Framework Programme (No. FP/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT). M. Meere thanks NUI Galway for the award of a travel grant. Appendix A. Linearization and stability analysis This appendix can be found in the supplementary material. References [1] K. T. Savjani, A. K. Gajjar, J. K. Savjani, Drug solubility: importance and enhancement tech- niques, ISRN Pharmaceutics 2012 (2012) 1–10. [2] H. D. Williams, N. L. Trevaskis, S. A. Charman, R. M. Shanker, W. N. Charman, C. W. Pou- ton, C. J. H. Porter, Strategies to address low drug solubility in discovery and development, Pharmacological Reviews 65 (2013) 315–499. [3] S. Kalepua, V. Nekkantib, Insoluble drug delivery strategies: review of recent advances and business prospects, Acta Pharmaceutica Sinica B 5 (5) (2015) 442–453. [4] C. Brough, R. O. W. III, Amorphous solid dispersions and nano-crystal technologies for poorly water-soluble drug delivery, International Journal of Pharmaceutics 453 (2013) 157–166. 22 [5] Y. Huang, W. Daib, Fundamental aspects of solid dispersion technology for poorly soluble drugs, Acta Pharmaceutica Sinica B 4 (1) (2014) 18–25. [6] K. Dhirendra, S. Lewis, N. Udupa, K. Atin, Solid dispersions: a review, Pakistan Journal of Pharmaceutical Sciences 22 (2) (2009) 234–246. [7] W. N. Souery, C. J. Bishop, Clinically advancing and promising polymer-based therapeutics, Acta Biomaterialia 67 (2018) 1–20. [8] C. G. Garcia, K. L. Kiick, Methods for producing microstructured hydrogels for targeted appli- cations in biology, Acta Biomaterialia 84 (2019) 34–48. [9] J. Brouwers, M. E. Brewster, P. Augustijns, Supersaturating drug delivery systems: the answer to solubility-limited oral bioavailability?, Journal of Pharmaceutical Sciences 98 (8) (2009) 2549– [10] M. Doi, Introduction to polymer physics, Oxford University Press, Oxford, UK, 1996. [11] P. J. Flory, Thermodynamics of high polymer solutions, Journal of Chemical Physics 10 (1942) 51–61. [12] M. L. Huggins, Solutions of long chain compounds, Journal of Chemical Physics 9. [13] P. C. Hiemenz, T. P. Lodge, Polymer Chemistry, 2nd Edition, CRC Press, Boca Raton, Florida, [14] J. Djuris, I. Nikolakakis, S. Ibric, Z. Djuric, K. Kachrimanis, Preparation of carba- mazepine–soluplus solid dispersions by hot-melt extrusion, and prediction of drug–polymer miscibility by thermodynamic model fitting, European Journal of Pharmaceutics and Biophar- maceutics 84 (2013) 228–237. [15] Y. Tian, V. Caron, D. S. Jones, A. M. Healy, G. P. Andrews, Using Flory–Huggins phase dia- grams as a pre-formulation tool for the production of amorphous solid dispersions: a comparison between hot-melt extrusion and spray drying, Journal of Pharmacy and Pharmacology 66 (2013) 256–274. [16] Y. Tian, J. Booth, E. Meehan, D. Jones, S. Li, G. Andrews, Construction of drug-polymer ther- modynamic phase diagram using Flory-Huggins interaction theory: identifying the relevance of temperature and drug weight fraction to phase separation within solid dispersions, Molecular Pharmaceutics 10 (2013) 236–248. [17] K. Bansal, U. S. Baghel, S. Thakral, Construction and validation of binary phase diagram for amorphous solid dispersion using Flory–Huggins theory, AAPS PharmSciTech 17 (2) (2016) 318–327. [18] S. Y. Chan, S. Qia, D. Q. M. Craig, An investigation into the influence of drug–polymer inter- actions on the miscibility, processability and structure of polyvinylpyrrolidone-based hot melt extrusion formulations, International Journal of Pharmaceutics 496 (2015) 95–106. [19] M. A. Altamimi, S. H. Neau, Use of the Flory–Huggins theory to predict the solubility of nifedipine and sulfamethoxazole in the triblock, graft copolymer soluplus, Drug Development and Industrial Pharmacy 42 (3) (2016) 446–455. [20] P. Chakravarty, J. W. Lubach, J. Hau, K. Nagapudi, A rational approach towards development of amorphous solid dispersions: Experimental and computational techniques, International Journal of Pharmaceutics 519 (2017) 44–57. 23 [21] D. Lin, Y. Huang, A thermal analysis method to predict the complete phase diagram of drug–polymer solid dispersions, International Journal of Pharmaceutics 399 (2010) 109–115. [22] K. Wlodarski, W. Sawicki, A. Kozyra, L. Tajber, Physical stability of solid dispersions with respect to thermodynamic solubility of tadalafil in pvp-va, European Journal of Pharmaceutics and Biopharmaceutics 96 (2015) 237–246. [23] Y. Zhao, P. Inbar, H. P. Chokshi, A. W. Malick, D. S. Choi, Prediction of the thermal phase diagram of amorphous solid dispersions by Flory–Huggins theory, Journal of Pharmaceutical Sciences 100 (2011) 3196–3207. [24] P. J. Marsac, S. L. Shamblin, L. S. Taylor, Theoretical and practical approaches for prediction of drug-polymer miscibility and solubility, Pharmaceutical Research 23 (10) (2006) 2417–2426. [25] J. H. Hildebrand, R. L. Scott, The solubility of nonelectrolytes, 3rd Edition, Dover Publications, New York, 1964. [26] B. D. Anderson, Predicting solubility/miscibility in amorphous dispersions: it is time to move beyond regular solution theories, Journal of Pharmaceutical Sciences 107 (2018) 24–33. [27] E. B. Smith, Basic chemical thermodynamics, sixth Edition, Imperial College Press, London, UK, 2014. [28] K. Binder, P. Fratzl, Spinodal decomposition, in: G. Kostorz (Ed.), Phase Transformations in Materials, WILEY-VCH Verlag, Weinheim, 2001, Ch. 6, pp. 411–480. [29] D. M. Saylor, C. Kim, D. V. Patwardhan, J. A. Warren, Diffuse-interface theory for structure formation and release behavior in controlled drug release systems, Acta Biomaterialia 3 (2007) 851–864. [30] D. M. Saylor, J. E. Guyer, D. Wheeler, J. A. Warren, Predicting microstructure development during casting of drug-eluting coatings, Acta Biomaterialia 7 (2011) 604–613. [31] J. Zhu, X. Lu, R. Balieu, N. Kringos, Modelling and numerical simulation of phase separation in polymer modified bitumen by phase-field method, Materials and Design 107 (2016) 322–332. [32] J. Zhu, R. Balieu, X. Lu, N. Kringos, Numerical investigation on phase separation in polymer- modified bitumen: effect of thermal condition, Journal of Materials Science 52 (11) (2017) 6525–6541. [33] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. i. interfacial free energy, The Journal of Chemical Physics 28 (1958) 258–267. [34] N. P. . K. Elder, Phase-field methods in material science and engineering, 1st Edition, John Wiley & Sons, Weinheim, Germany, 2010. [35] G. M. Kontogeorgis, G. K. Folas, Thermodynamic models for industrial applications, 1st Edi- tion, John Wiley & Sons, Chichester, West Sussex, UK, 2010. [36] J. M. Prausnitz, R. N. Lichtenthaler, E. G. de Azevedo, Molecular thermodynamics of fluid- phase equilibria, 3rd Edition, Prentice Hall, New Jersey, 1999. [37] D. A. Porter, K. E. Easterling, M. Sherif, Phase transformations in metals and alloys, 3rd Edition, CRC Press, Boca Raton, Florida, 2009. 24 [38] D. L. Elbert, Liquid–liquid two-phase systems for the production of porous hydrogels and hy- drogel microspheres for biomedical applications: A tutorial review, Acta Biomaterialia 7 (2011) 31–56. [39] P. Atkins, J. de Paula, Elements of physical chemistry, 5th Edition, Oxford University Press, Oxford, England, 2009. [40] E. B. Naumana, D. Q. Heb, Nonlinear diffusion and phase separation, Chemical Engineering Science 56 (2001) 1999–2018. [41] J. Gerges, Numerical studies of the physical factors responsible for the ability to vit- rify/crystallize of model materials of pharmaceutical interest, Ph.D. thesis, Universite de Lille (2015). [42] L. Ratke, P. W. Voorhees, Growth and Coarsening: Ostwald Ripening in Material Processing, 1st Edition, Springer, Berlin, 2002. [43] K. B. Migler, String formation in sheared polymer blends: Coalescence, breakup, and finite size effects, Physical Review Letters 86 (6) (2001) 1023–1026. [44] X. Liu, X. Feng, R. O. W. III, F. Zhang, Characterization of amorphous solid dispersions, Journal of Pharmaceutical Investigation 48 (1) (2018) 19–41. [45] H. Tanaka, Viscoelastic phase separation, Journal of Physics:Condensed Matter 12 (15) (2000) 207–264.

Journal

Condensed Matter – arXiv (Cornell University)

Published: Feb 13, 2019

There are no references for this article.