Anh D. Phan*abc,
Katsunori Wakabayashic,
Marian Paluchd and
Vu D. Lamef
aFaculty of Materials Science and Engineering, Phenikaa Institute for Advanced Study, Phenikaa University, Hanoi 12116, Vietnam. E-mail: anh.phanduc@phenikaa-uni.edu.vn
bFaculty of Information Technology, Artificial Intelligence Laboratory, Phenikaa University, Hanoi 12116, Vietnam
cDepartment of Nanotechnology for Sustainable Energy, School of Science and Technology, Kwansei Gakuin University, Sanda, Hyogo 669-1337, Japan
dInstitute of Physics, University of Silesia, SMCEBI, 75 Puku Piechoty 1a, 41-500 Chorzow, Poland
eInstitute of Materials Science, Vietnam Academy of Science and Technology, Hanoi, Vietnam
fGraduate University of Science and Technology, Vietnam Academy of Science and Technology, Hanoi, Vietnam
First published on 4th December 2019
Theoretical approaches are formulated to investigate the molecular mobility under various cooling rates of amorphous drugs. We describe the structural relaxation of a tagged molecule as a coupled process of cage-scale dynamics and collective molecular rearrangement beyond the first coordination shell. The coupling between local and non-local dynamics behaves distinctly in different substances. Theoretical calculations for the structural relaxation time, glass transition temperature, and dynamic fragility are carried out over twenty-two amorphous drugs and polymers. Numerical results have a quantitatively good accordance with experimental data and the extracted physical quantities using the Vogel–Fulcher–Tammann fit function and machine learning. The machine learning method reveals the linear relation between the glass transition temperature and the melting point, which is a key factor for pharmaceutical solubility. Our predictive approaches are reliable tools for developing drug formulations.
The relaxation processes of amorphous materials can be experimentally measured using broadband dielectric spectroscopy (BDS) and differential scanning calorimetry (DSC).1,2 BDS technique determines the structural relaxation time corresponding to thermal variation. DSC can measure the glass transition temperature and analyze phase separation in experimental samples at different cooling rates. The relevant timescale of molecular motions measured by BDS spans from picosecond above melting temperature to hundreds seconds in vicinity of the glass transition temperature. This technique can be used to investigate both structural (primary or long-time) and transient (secondary or short-time) relaxation processes.
Recently, the Elastically Collective Nonlinear Langevin Equation (ECNLE) theory has been developed to understand structural relaxation time of amorphous systems,6–11 in which an amorphous material is modeled as a fluid of molecular particles. The ECNLE theory considers a single molecular motion affected by the nearest-neighbor interactions and cooperative motions of molecules outside a particle cage formed by the neighboring molecules.6–10 This physical picture leads to a local barrier of a dynamic free energy and a collective barrier for each molecule caused by its nearest-neighbor interactions and effects of cooperative molecular rearrangements, respectively. These barriers are density-dependent and mutually correlated.
Plugging the two barriers calculated using the ECNLE theory into Kramer's theory gives the structural relaxation time.6–11 To find the temperature dependence of alpha relaxation time, one proposed an analytical conversion (thermal mapping) from averaged particle density to temperature based on experimental dimensionless compressibility data associated with equation of state.6–10 The theoretical calculations have also provided quantitative good predictions for the glass transition temperature and dynamic fragility of colloidal systems,8,9 supercooled liquids,7,9 and polymer melts.6,10 However, the experimental equation-of-state data needed for the original thermal mapping has rarely measured in amorphous drugs. Thus, we have recently proposed another density-to-temperature conversion using the thermal expansion and experimental glass transition temperature values.11 Our approach has successfully described temperature-dependent molecular dynamics in one- and multi-component amorphous drugs.
Most ECNLE calculations have assumed a universal correlation between local and collective barrier for all substances when inserting into Kramer's theory. The assumption simplifies roles of chemical and biological complexities on the glass transition. Consequently, the quantitative agreement between theoretical calculations and experiments is imperfect. Moreover, since the dynamic fragility of a material depends strongly on the form of τα(T), a small theory-experiment deviation in τα(T) leads to an inaccurate prediction of the fragility.
In this work, we introduce an adjustable parameter characterizing for a non-universal local-collective correlation in pharmaceuticals. The new version of the ECNLE theory accurately and simultaneously predicts the glass transition temperature and dynamic fragility of amorphous drugs. We employ machine learning technique to reveal a undiscovered relation between melting point and glass transition. We also predict the glass transition temperature based on BDS data. Our ECNLE numerical results are quantitatively compared to experimental data and machine learning calculations.
Theoretical understanding of how glassy dynamics of amorphous materials varies with the cooling rate is limited. Previous theoretical studies12–14 have given phenological/qualitative descriptions for the Tg change as a function of cooling rate. Based on the theoretical development of this paper, we propose, for the first time, an analysis to provide quantitative determinations for this phenomenon and insightful discussions for experiments.
(1) |
The thermal noise force satisfies 〈δf(0)δf(t)〉 = 2kBTζsδ(t), where kB is the Boltzmann constant and T is temperature. According to ECNLE theory, the dynamic free energy is6,7,15,16
(2) |
(3) |
C(r) = 0 for r > d. | (4) |
Recall that the ECNLE theory ignores effects of rotational motions and only consider translational motions, which are angularly-averaged. The first term of eqn (2), which depends strongly on the fluid structure and density, is the dynamic trapping potential and favors the particle localization. While the second term independent of the system structure represents the ideal fluid state.
Fig. 1 shows an illustration of the dynamic free energy as a function of r and key physical quantities in the local caging constraint. In dilute suspension (Φ ≤ 0.43),15–17 Fdyn(r) monotonically decreases with an increase of r and particles move without constraint. When Φ > 0.43, the tagged particle is dynamically arrested within a particle cage formed by its neighbors and one observes an emergence of a free-energy barrier. The particle cage radius, rcage, as depicted in Fig. 1 is determined as the first minimum position in the radial distribution function, g(r). Since S(q) and g(r) are a Fourier-transform pair, one has Thus, the radius of the particle cage is about 1.3 − 1.5d. When the dynamic free energy reaches the local minimum and maximum, one obtains the localization length (rL) and the barrier position (rB). The energy difference between these two positions is the local hopping barrier FB = Fdyn(rB) − Fdyn(rL). The jump distance from the localized position to the barrier position is defined as Δr = rB − rL. K0 = |∂2Fdyn(r)/∂r2|r=rL and KB = |∂2Fdyn(r)/∂r2|r=rB are absolute curvatures at the localization length and barrier position.
The rearrangement of particles in the first shell causes a small expansion on the surface of the particle cage and generates a harmonic displacement field u(r) in surrounding medium via collective motions of other particles. In bulk systems, one can obtain analytical form of the distortion field by Lifshitz's continuum mechanics analysis18
(5) |
(6) |
Solving eqn (5) with the boundary condition at the cage surface gives
(7) |
The spatial harmonic displacement energy of a particle at separation distance r from the center of its arrested cage is K0u2(r)/2. Since the local time-averaged density is ρg(r), the number of particles found at a distance between r and r + dr is ρg(r)4πr2dr. One can integrate the elastic energies of particles outside the cage to quantify effects of cooperative motions. The collective elastic barrier, Fe, is
(8) |
For r ≥ rcage, one can approximate g(r) ≈ 1. The collective motions of molecules play a more important role than the local dynamics in the glass transition at high densities or low temperatures6–8 as depicted in Fig. 1.
The activated relaxation is governed by both local and non-local processes. If a universal correlation between a local and non-local process is assumed, the total barrier is simply Ftotal = FB + Fe. The alpha relaxation time for a particle to diffuse from its particle cage is quantified by Kramer's theory:
(9) |
(10) |
To convert our hard-sphere calculations into the temperature dependence of the structural relaxation time, we proposed11 a thermal mapping, which is based on the thermal expansion process during temperature variation, to convert from a volume fraction of hard-sphere fluid to temperature of experimental material. The mapping is
(11) |
The parameter T0 depends on molar mass and particle size. Our numerical calculations indicate the structural relaxation time τα ≈ 100 s at Φ ≈ 0.611. Thus, one can approximately obtain T0 = Tg − (0.611 − Φ0)/βΦ0. The experimental value of Tg can be found in many literatures. Based on the calculation, we investigated the temperature dependence of the structural relaxation time of many amorphous drugs and their mixtures (binary and ternary composites).11 The theoretical calculations without any adjustable parameters quantitatively agree with various experimental data over 14 decades in time. While simulations can determine relaxation times only over first 3–6 decades and do not access experimental timescale as illustrated in Fig. 1.
Fig. 2 shows log10τα of five representative pure amorphous drugs21–24 as a function of 1000/T calculated using eqn (9)–(11). Overall, the theoretical curves are relatively close to the experimental counterpart, except for calculations of vitamin A. A deviation of experimental data of the vitamin-A drug from theoretical calculations is expected. This is because the theory ignores many chemical and structural complexities such as hydrogen-bonding, network formers, and flexible molecular docking. Moreover, the approach seems to provide less quantitatively accurate predictions for the dynamic fragility of amorphous materials.
Fig. 2 The temperature dependence of structural relaxation time of chloramphenicol,21 indapamide,22 ezetimibe,22 biclotymol,23 and vitamin-A acetate.24 Open points are experimental data in literatures and solid curves correspond to our ECNLE calculations. |
(12) |
One adopts the physical quantity to classify into two main categories: “strong” or “fragile” for glass-forming materials. For m ≤ 30, the glass formers are strong. The glass formers having m ≥ 100 are fragile. The remaining materials are called intermediate glass-forming materials.
The dynamic fragility is very sensitive to the slope of τα(T) at T = Tg. Thus, a good agreement between theory and experiment in τα versus T does not mean that another consistency in calculations of m occurs. In ref. 10, authors introduced an adjustable parameter ac to scale the collective elastic barrier as Fe → ac2Fe. The parameter ac captures chemical and biological complexities, conformational configuration, and chain connectivity in different thermal liquids and polymers. The parameter assesses the relative importance of the collective elastic distortion by assuming a non-universal coupling of the cage-scale hopping and collective rearrangements of fluid particles. The previous work10 obtained contemporaneously the quantitative accordance between theoretical ECNLE calculations and experimental values of both dynamic fragility and glass transition temperature for 17 polymers.
Motivated by the idea in ref. 10, in our calculations, we adjust parameters T0 and ac to achieve the best fit to the experimental temperature dependence of structural relaxation times. Fig. 3 shows experimental data and our theoretical calculations for τα(T) of 22 pure amorphous materials. We carry out the same procedure as calculations in Fig. 2 except for now Ftotal = FB + ac2Fe. Our numerical results agree quantitatively well with a wide range of experimental data. Remarkably, the activated events of carvedilol, celecoxib, chloramphenicol, and polystyrene below Tg, where τα ranges from 100 s to 104 s, are well-described using the ECNLE theory. Many previous works26,27,43,44 have observed a distinctive deviation, so-called a dynamic structural decoupling, in the relaxation process at low temperatures. For example, in Fig. 3a, the growth of τα(T) of carvedilol drug below Tg is abruptly deviated from what it is supposed to be. Currently, there is poor theoretical understanding for the interesting but challenging feature. In the framework of the ECNLE theory, the decoupling could be related to a temperature dependence of thermal expansion coefficient β in our thermal mapping.
Fig. 3 The temperature dependence of structural relaxation time of twenty-two various amorphous drugs and polymers listed in Table 1.11,18–39 Open points are experimental data in literatures and solid curves correspond to our ECNLE calculations. PVP is an abbreviation of polyvinyl pyrrolidone K30. |
From calculations in Fig. 3, we can estimate the glass transition temperature and dynamic fragility. The local-nonlocal coupling parameter ac, the characteristic temperature T0, the melting temperature Tm, and theoretical and experimental values for Tg and m of the studied materials are listed in Table 1. Clearly, the theoretical Tg is in perfect accordance with the experimental counterpart. The different accuracy of our calculations for the fragility is somehow expected and reflects a complicated Tg − m correlation.
Materials | Tg(th) (K) | Tg(expt) (K) | m(th) | m(expt) | Tm(expt) (K) | ac | T0 (K) |
---|---|---|---|---|---|---|---|
Carvedilol25 | 308 | 310 (ref. 25) | 91.5 | 387.5 (ref. 38) | 2.1 | 450 | |
Celecoxib26,27 | 328 | 328 (ref. 26) | 97.8 | 110 (ref. 26) | 432 (ref. 38) | 2.1 | 470 |
Chloramphenicol21 | 301.1 | 301 (ref. 21) | 89 | 116 (ref. 21) | 423.5 (ref. 38) | 2.1 | 443 |
Griseofulvin28 | 358 | 359 (ref. 28) | 88.7 | 84.6 (ref. 28) | 489 (ref. 39) | 1.2 | 533 |
Indomethacin29 | 314.1 | 315 (ref. 29) | 86.4 | 77, 67, 64 (ref. 29) | 432 (ref. 39) | 1.5 | 476 |
Ketoconazole30 | 308 | 316.3 (ref. 30) | 71.37 | 419 (ref. 38) | 1.0 | 493 | |
Probucol31 | 294 | 294.7 (ref. 31) | 79.4 | 85 (ref. 31) | 398 (ref. 38) | 1.5 | 456 |
Polystyrene32 | 373 | 373 (ref. 10) | 105.53 | 116, 143, 97, 121 (ref. 10) | 513 (ref. 40) | 1.7 | 528 |
Bicalutamide33 | 325.3 | 325.4 (ref. 33) | 78.03 | 84 (ref. 33) | 464 (ref. 39) | 1.1 | 505 |
Biclotymol23 | 288 | 288 (ref. 23) | 85.42 | 85 (ref. 23) | 2.1 | 430 | |
Polyvinylpyrrolidone K30 (ref. 26) | 431.2 | 431 (ref. 26) | 58.47 | 70 (ref. 26) | 523 (ref. 41) | 0.3 | 672 |
Tripropyl phosphate34 | 132.1 | 134 (ref. 34) | 40.3 | 194 (ref. 38) | 2.4 | 266 | |
Vitamin-A acetate24 | 236.3 | 244.3 (ref. 24) | 76.6 | 83 (ref. 24) | 330 (ref. 42) | 3.4 | 349 |
Nisoldipine35 | 303 | 305 (ref. 35) | 70.2 | 81 (ref. 35) | 425 (ref. 42) | 1.0 | 488 |
Nifedipine35 | 315 | 315 (ref. 35) | 74.9 | 84 (ref. 35) | 444 (ref. 39) | 1.1 | 494 |
Nimodipine35 | 284.1 | 285 (ref. 35) | 62 | 82 (ref. 35) | 398 (ref. 38) | 0.9 | 474 |
Indapamide22 | 373.4 | 373.5 (ref. 22) | 75.3 | 73 (ref. 22) | 433 (ref. 38) | 0.7 | 578 |
Ezetimibe22 | 333.3 | 333.1 (ref. 22) | 91.1 | 93 (ref. 22) | 436 (ref. 38) | 1.5 | 495 |
Kollidon VA64 (ref. 11) | 376.2 | 378 (ref. 11) | 79 | 0.8 | 573 | ||
Simvastatin11 | 301 | 303 (ref. 11) | 74.8 | 73 | 411 (ref. 38) | 1.3 | 471 |
Flutamide36 | 272 | 271 (ref. 36) | 71.6 | 385 (ref. 38) | 1.4 | 438 | |
Ibuprofen37 | 222 | 225 (ref. 37) | 68 | 87 (ref. 37) | 346 (ref. 39) | 2.4 | 356 |
Here we introduce, for the first time, another approach based on machine learning techniques. After obtaining a predictive model by applying the support vector regression (SVR) in the Scikit-learn Python library54 to experimental data of the temperature dependence of τα at high temperatures, we can predict new relaxation times at lower temperatures than the coolest temperature in the training dataset. The SVR with the radial basis function (RBF) kernel has two controlled parameters including a regularization parameter CRBF and a RBF kernel parameter γ. While γ is considered as the inverse of the standard deviation of the RBF kernel, CRBF determines the penalty of large slack variables. In our calculations, we chose γ = 0.1 and obtained equal performance when CRBF ≥ 100. As shown in Fig. 4, τα(T) versus 1000/T predicted by the SVR with CRBF = 103 grows smoothly and are close to experimental data for some representative amorphous drugs. The number of BDS data points for simvastatin, ketoconazole, bicalutamide, griseofulvin, vitamin-A acetate, and nisoldipine used for training are 24, 8, 17, 17, 22, and 28, respectively. The SVR calculations give Tgs for simvastatin, ketoconazole, bicalutamide, griseofulvin, vitamin-A acetate, and nisoldipine are 300.2, 308.64, 323, 358.2, 237.8, and 304 K, respectively. These values quantitatively agree experimental data and our ECNLE results in Table 1. This is a new reliable approach for investigating molecular dynamics of amorphous materials near Tg, particularly when the structural decoupling appears.26,27,33,43,44 The linear regression cannot be used in this protocol since it enforces a Arrhenius nature on any experimental data and keeps it unchanged in the predictive process.
Fig. 4 The temperature dependence of structural relaxation time of simvastatin11, ketoconazole,30 bicalutamide,33 griseofulvin28, vitamin-A acetate,24 and nisoldipine.35 Open points are experimental data in literatures and solid curves correspond to machine learning-based calculations. |
To find new minimalist correlations or physical insights among the quantities in Table 1, we employ a linear regression model in the scikit-learn library.54 This regression algorithm provides the simplest relation to describe a target variable from a set of descriptor variables. By adopting Tg and Tm experimental values of 71 glassy drugs listed in ref. 46 as a training dataset, one obtains Tm ≈ 1.362 Tg. It is important to note that these drugs are different from amorphous materials in this work. Then, we apply the predicted relation to our twenty substances to evaluate the validity of the model. Numerical results in Fig. 5 indicate that the Tg − Tm correlation works well. This finding suggests the ECNLE theory can be exploited to estimate the melting temperature with a reasonable deviation.
Fig. 5 Cross-plot of predicted and experimental values of the melting temperature for all 20 drugs and polymers from Table 1. The blue line indicates perfect agreement. |
In many prior works,47–49 the melting point of amorphous drugs exhibits an essential role in the solubility determination of the drugs. The linear relation between Tg and Tm suggests that one can employ ECNLE theory to evaluate the solubility of amorphous materials and their mixtures.
According to an assumption introduced by Cooper and Gupta in 1982,12 the molecular relaxation time is approximately equal to the experimental observation time at Tg. This assumption leads to a new definition of a cooling rate, h,
(13) |
The minus sign in eqn (13) represents an inverse variation between mobility and temperature.
Near glass transition temperature, eFtotal/kBT ≫ 1 and the alpha structural relaxation time in eqn (9) approximately becomes
(14) |
Recall that the total barrier here is Ftotal = FB + ac2Fe as discussed in Section III. After straightforward transformations, one obtains
(15) |
Since τs is order of picoseconds (10−12 s), the first term is much smaller than the second term near Tg. Thus,
(16) |
When the temperature dependence of structural relaxation time obeys the Arrhenius behavior, the total barrier Ftotal is a constant and our eqn (16) can be deduced to be the same mathematical form as previous studies13,14
(17) |
Moreover, combining eqn (14) and (16) gives
(18) |
The above equation reveals explicitly a correlation between the cooling rate and the glass transition temperature. Recall that various experimental studies have empirically indicated that lnh is linearly proportional to −1/Tg in a specific range of Tg. If the relation is universal, it suggests the first two terms in eqn (18) may play a minor role compared others and then our analysis can clearly explain the experimental observation.
Fig. 6 shows how the cooling rate influences the glass transition temperature of polystyrene,50 PVP,51 nifedipine,52 and indomethacin.51 Our theoretical calculations are relatively consistent with experimental data. Although both theoretical curves and data points are not perfect straight lines, the rough linearity is possibly observed in a certain range of the glass transition temperature. One can also crudely view these smooth curves as a combination of high- and low-Tg linear branches. Furthermore, as discussed in Section IV, Tm ≈ 1.362Tg, we expect plotting lnh versus 1000/Tm (not shown here) does not change the variation trends compared to Fig. 6. The finding is in accordance with other experimental works.53
Fig. 6 The logarithm of cooling rate as a function of inverse glass transition temperatures of various amorphous drugs and polymers. Points are experimental data and solid curves correspond to our theoretical calculations using eqn (16). |
The effects of cooling rate on Tg can be also analyzed using the VFT-type relaxation time, which is
(19) |
(20) |
Interestingly, both ECNLE theory and VFT-function analysis give us the same form of the nontrivial correlation among the cooling rate, glass transition temperature, and fragility
(21) |
For small lnq (≤0) or low cooling rates, Tg of polystyrene is nearly unchanged, approximately 373 K. Thus, one can utilize the definition of τα(Tg) = 100 s to estimate Tg in the range of the cooling rate. It suggests that if τα of all substances is measured at the same low cooling rate, hτα(Tg) ≈ 100h and Tg is linearly proportional to m. However, since different experiments have been carried out at different cooling rates, it is hard to obtain a Tg − m trivial correlation.
This journal is © The Royal Society of Chemistry 2019 |