Survival of the Fittest: Time-To-Event Modeling of Crystallization of Amorphous Poorly Soluble Drugs

The objective of this study was to gain a quantitative understanding of the link between physicochemical properties and long-term and time-censored amorphous stability of poorly water-soluble drugs using parametric time-to-event modeling. Previously published data on amorphous stability and physico-chemical properties of 25 structurally diverse neutral, poorly soluble compounds were used. To describe the general shape of the survival curve (probability of event at time > t ), Constant, Gompertz, and Weibull hazard functions and their linear combinations were tested. For a selected Weibull hazard base model, the effect of each physicochemical covariate was investigated, with combined in ﬂ uence of enthalpy of fusion ( H f ) and molecular weight ( M r ) showing the highest statistical signi ﬁ cance. The covariate model was used to simulate survival curves and calculate the median survival time for different values of H f and M r . It was found that a decrease in H f or an increase in M r contribute to longer survival times. The derived model equation was validated against external data sets consisting of 11 compounds. It showed better predictive ability than a previously published multiple linear regression model incorporating H f and M r . The proposed Weibull covariate model may assist in faster and more cost-effective decision making in the pre-formulation phase of drug development, where compound properties and appropriate drug formulation strategies are investigated. © 2016 The Authors. Published by Elsevier Inc. on behalf of American Pharmacists Association ® . This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Small-molecule drug candidates with good pharmacological properties frequently suffer from low aqueous solubility.For oral drug products, this can lead to insufficient and erratic bioavailability of the active pharmaceutical ingredient (API), and for that reason, lipophilic APIs with limiting aqueous solubility are often formulated in the amorphous state, typically as solid dispersions in polymer matrices. 1This approach is favorable because the amorphous solid state of a molecule is at a higher energy level compared with its crystalline state, giving rise to higher dissolution rates and solution concentrations of poorly soluble APIs from amorphous versus crystalline formulations. 2However, on account of the comparatively disordered and mobile nature of amorphous solids, APIs can crystallize from amorphous formulations on storage or on contact with gastrointestinal fluids after administration.Such instability is a significant hurdle in the development of new medicines, and it would be highly desirable to be able to predict the inherent amorphous stability of poorly soluble APIs and the stability of such molecules in amorphous solid dispersions.
One physical stability classification methodology assesses the glass-forming ability of APIs after rapid solvent evaporation from a solution or cooling from a melt, 3,4 whereas in a different method the phase transition from the amorphous to a crystalline state is observed directly as a function of time. 5,6In the former method, the assumption is that some compounds have an inherent ability to form amorphous solids from liquids, whereas others more readily assume a crystalline state on transition from the liquid to the solid state.It has been observed that the crystallization classification of model APIs is similar, independent, however, of methodology.From the limited data published to date, it would appear that the crystallization behavior of APIs depends predominantly on the properties of the API and other formulation components, rather than the amorphization methodology itself. 4,7However, environmental conditions, for example, humidity and temperature, can also affect amorphous drug stability 8 and can, therefore, lead to different results reported by different laboratories.Although the intrinsic factors affecting ease of compound crystallization from the amorphous state, such as molecular mobility and thermodynamic properties of the amorphous state, 9 remain poorly understood, efforts have been underway to correlate observed glass-forming ability or physical stability of amorphous APIs with physicochemical and other properties of the APIs.][12][13] The framework developed on pure APIs in this work can be used as a preliminary risk assessment for the development of their solid dispersion formulations.It is expected that the amorphous stability of solid dispersions will be improved in comparison with pure drugs due to the dilution effect and the presence of specific interactions between a drug and a polymer.Amorphous compounds that are shown to have good amorphous stability should also be stable in their solid dispersions.Conversely, compounds that have poor stability may still be formulated as viable solid dispersions, but the products are likely to have higher risk of crystallization.
In a previous study, 6 the amorphous stability of 25 diverse neutral poorly soluble drug compounds were investigated.Principal component analysis and clustering methods were used to select 25 compounds with diverse physicochemical properties and chemical structures from the database of 533 marketed poorly soluble drugs.The selected sample set was shown to be representative for the calculated and predicted variables as side-by-side histograms for both distributions showed the same mean location and variance. 6Several multiple linear regression (MLR) models were proposed to predict long-term amorphous drug stability using only easily accessible physicochemical drug properties as covariates.Due to practical limitations, continuous crystallization records during the approximately 6-month period of the experiment were not obtained.Instead, amorphous drug stability was measured at defined time points, with some compounds remaining stable at the end of the experiment.Time-to-event (TTE) modeling, also known as survival analysis, 14,15 is particularly suited to this type of data and was applied in the present study to the previously published data on 25 diverse compounds.Here we present the new application of TTE modeling to more accurately determine the influence of physicochemical parameters on the long-term amorphous stability of poorly soluble compounds.
A central concept in TTE modeling is the survival function S(t), which describes the probability that an event will occur at a time greater than t.The survival function is related to the hazard function h(t), which can be understood as the instantaneous failure rate (onset of amorphous to crystalline transition in this study) given that the compound has survived to that point in time.The survival and hazard functions are related through Equation 1.
In parametric TTE modeling, a base model is first derived using a process in which various hazard functions h(t) are proposed, and estimates of the parameters of the hazard functions are determined using maximum likelihood estimation.This involves integration of the hazard function from zero to the time of each crystallization event, giving the contribution of each event to the total likelihood, where the parameter estimates are derived by maximizing the sum of all likelihood contributions with respect to the parameter values.This estimation process is applied to each of the candidate hazard functions, and the best function is selected based on model selection criteria and graphical comparison of observed and predicted data.Once a hazard function has been selected, a covariate analysis is performed where the physicochemical properties of the compounds are allowed to influence the parameters within the hazard function.Model selection criteria are then used to select a final covariate model where the included covariates produce a statistically significant improvement compared with the base model.The model can then be used to simulate the expected behavior of new compounds with different values of the included covariates.

Software
The TTE models were developed using NONMEM software, 16 version 7.3 (Icon Development Solutions, Ellicott City, MD).All other analyses and visualization of data were implemented using R software 17 and use of the "survival" library, version 2.38, 18 and the "deSolve" library, version 1.11. 19

Database Preparation
Previously published data of physicochemical properties and amorphous stability of 25 compounds 6 were used to derive a TTE model.For 17 of the compounds, the precise time of detectable transition from amorphous to crystalline structure is unknown.It was only recorded that the transition took place between 2 observation times.Such data are called interval-censored and were flagged on both sides of the interval within the database used for modeling (Supplementary Material A).For the remaining 8 compounds, observations were flagged as right-censored because crystallization was not observed during the 168 days of experiment.This vector of flag values was used as a dependent variable.The elapsed time of the interval and right-censored observations were included as main independent variable.The following measured, calculated, and predicted physicochemical properties of the 25 compounds 6 were also tested during development of the covariate model: enthalpy of fusion (H f ), glass transition temperature (T g ), melting temperature (T m ), configurational entropy (S c ), enthalpy (H c ) and free energy (G c ), relaxation time (t), molecular weight (M r ), hydrogen bond donors and acceptors, rotatable bonds (rotB), number of rings, aromatic rings, aliphatic rings, heavy atom count, ratio of carbon to heteroatoms, polar surface area, lipophilicity (clogP), and water solubility (logS w ) predicted both with CLab 6,20 and ALOGPS 21 (www.vcclab.org).These descriptors have been shown in the literature to have an impact on the glass-forming ability and amorphous stability, 3,11 drug bioavailability, 22 and the stability of compounds formulated as solid dispersions. 23

Model Building and Selection Criteria
For selection of an appropriate base model, 8 different hazard functions were evaluated (Table 1). 24The hazard functions are defined with parameters l and b which determine the shape of the hazard function (Supplementary Material B).Because hazard functions must always be positive, the lambda parameters (l c , l 1 , and l 2 ) were constrained to positive values during parameter estimation, whereas positive or negative estimates of the beta parameters (b c and b 1 ) were permitted.The parameter values in the TTE models were estimated using the first-order conditional estimation method, which uses maximum likelihood estimate to minimize the objective function value (OFV) and, thus, obtain the best model fit to the data. 25he direction in which the OFV decreases to find the global minimum was determined in an iterative process using the gradient method. 26o ensure convergence to the global minima, models were run repetitively, each time starting from different initial parameter estimates.The global minimum was considered to be reached when all repetitive runs resulted in the same final parameter estimates and OFVs.Estimates of parameter uncertainty (standard error) were generated using the NONMEM covariance procedure. 16odels were first required to pass the following acceptance criteria 27 : 1.At the last iteration, the values of the gradients have to be between 10 À6 and 10 to indicate that the OFV successfully converged to its minimum.2. The covariance step, where uncertainty in parameter estimates is determined, has to complete successfully.3. The number of significant digits of parameter estimates has to be !3. 4. The correlation between estimated parameters, cp, has to be À0.95 cp 0.95.This criterion will reject models that are overparametrized.5.The ratio of maximum to minimum eigenvalues of the correlation matrix of parameter estimates (condition number) has to be 1000.This criterion will reject models that are overparametrized.6.The standard error of each parameter estimate is required to be <50% of the estimated parameter value.This ensures that the 95% confidence interval for the parameter estimate excludes zero.
For those models that passed the acceptance criteria, the base model (Table 1) was then selected using Akaike information criterion (AIC).This criterion is used to compare nonnested models, where fixing one parameter in the first model to its null value is not leading to the second model.According to this criterion (Eq.2), a model (A) is statistically superior to a model (B) if AIC <0.
where n is the number of model parameters.
The influence of covariates on the parameters of the selected base model was then tested using the forward inclusion method.Covariate effects were assessed one by one followed by evaluation of the combined effect of significant covariates.The effect of each covariate, COV, was modeled relative to the median value of that covariate, COV median , in the sample set 28 according to Equation 3.
Covariate models of 2 different types were tested.When testing the influence of a covariate on a l parameter, l as given by the equations in Table 1 was replaced by l cov as defined by Equation 4l where q is a new parameter, estimated by the modeling software, which quantitatively describes the influence of the covariate on the hazard function.The functional form of Equation 4enables the possibility that an increase in the covariate can increase or decrease the magnitude of the hazard function (positive or negative estimate of q) while satisfying the constraint that the hazard function should always be positive.When testing the influence of a covariate on a b parameter, a different functional form of the covariate model was used as defined by Equation 5.
This different functional form allows the possibility of changes in the covariate influencing the magnitude and sign of the original b parameter because both positive and negative values of b cov are permitted (they both lead to a positive hazard function).The forms of Equations 4 and 5 also ensure that if the covariate parameter q is zero, then the hazard function reduces to that of the base model.The significance of the covariate models was assessed using the likelihood ratio test (LRT). 14The covariate was considered to have a significant effect if the decrease in the OFV of the covariate model from that of the base model was >6.63 (p value < 0.01).This stringent p value was used due to the inflated type I error that occurs during multiple testing (testing of many possible covariate relationships).For the final covariate model, the estimates of parameter uncertainty produced by the NONMEM covariance procedure were supplemented by a more thorough bootstrap procedure implemented using PsN software. 29This involved taking the structure of the final covariate model and repeatedly fitting it to 1000 data sets, each produced by sampling with replacement from the original data set of 25 compounds.This provided 1000 estimates of each model parameter.Then 95% nonparametric confidence intervals were generated from the 2.5 and 97.5 percentiles of the 1000 parameter estimates.

Model Qualification
1][32] To generate the VPC, the selected covariate model was used to simulate 1000 sets of event times for the 25 compounds.Then, using the survival library in R, KaplaneMeier survival curves with 95% confidence intervals were plotted for the observed events and for the 2.5, 50, and 97.5 percentiles.A good model should show extensive overlap of the confidence intervals derived from the observed and simulated event data with little evidence of bias.The simulated event times were generated through numerical integration (using R library deSolve) of h(t) to give S(t).This was followed by repeated sampling of random numbers from a uniform distribution in the interval 0 to 1 (because S(t) is uniformly distributed on the interval 0 to 1) 24 and finding corresponding times at which the survival function became less than or equal to the sampled random number.Adjustment of the simulated event times according to the scheduled experimental measurement times (1 hour, 3 hours, 1 day, 7 days, 1 month, 2 months, 4 months, and 6 months) 6 was necessary to allow correct comparison with the interval-censored nature of the observed data.For example, if a particular simulated event time was 4.63 days, it would be adjusted to 7 days because that is the first point in time where the crystallization event could have been observed according to the experimental schedule.

Selection of the Base Model
The 8 hazard functions tested for selection as base model were Constant, Gompertz, and Weibull functions and their linear combinations (Table 1).Parameter estimates, OFV, and selection criteria status for the 8 candidate base models are given in Table 2.The last 4 models were immediately rejected due to failure of acceptance criteria.Failure was either due to an unacceptably wide standard error in one or more parameter estimates or to unsuccessful completion of the NONMEM covariance procedure that generates the standard errors.The latter occurred due to numerical instability with the Constant þ Weibull model where the parameter l c optimized to a vanishingly small value.Figure 1 shows the 4 models that passed the acceptance criteria, where the estimated survival curves derived from the models are shown superimposed on the KaplaneMeier plot of the observed crystallization events.The Weibull and Constant þ Gompertz models clearly match the observed KaplaneMeier curve better than the Constant and Gompertz models.The 2 former models have OFVs that are lower than the latter 2; therefore, a selection between the Weibull and Constant þ Gompertz model was the only decision required.
Selection of the base model from the 2 remaining models was initially performed using a consideration of AIC alone.The Constant þ Gompertz model has an OFV that is 2.8 units lower than the Weibull model (Table 2) but has 1 additional parameter (Table 1).This leads to an AIC of À0.8 in favor of the Constant þ Gompertz model.It should be noted from Table 2 that the standard errors of all 3 parameter estimates in the Constant þ Gompertz model are close to 50% of the parameter estimates, whereas the precision of parameter estimates in the Weibull model is considerably better.During early attempts to build covariate models using the initially selected Constant þ Gompertz base model, it was soon found that the additional complexity arising from incorporation of covariate models caused all parameters to have standard errors that were >50% of the associated parameter estimate.Due to its additional complexity, the Constant þ Gompertz model appears to be less stable than the Weibull model and is likely to be overparameterized when covariates are added.Hence, the selection of the most appropriate base model was reassessed, and the Weibull model was chosen on the basis that it has an OFV only 2.8 units higher than the Constant þ Gompertz model.It leads to a similar visual concordance with the observed KaplaneMeier curve (Fig. 1) but with the advantage of one fewer parameter and more precise parameter estimates.

Selection of the Covariate Model
The individual and combined effects of calculated, predicted, and measured compound properties 6 on the parameters l 1 and b 1 in the Weibull hazard function were investigated.The LRT was used to decide if incorporation of a particular covariate led to a significant improvement in the model, with the requirement of a decrease in the OFV of the covariate model with respect to the base model of >6.63.Based on this criterion, none of the covariates were found to have a significant effect on b 1 and 4 covariates were found to have a significant effect on l 1 when tested independently: H f (DOFV ¼ À13.46), M r (DOFV ¼ À10.33), heavy atom count (DOFV ¼ À9.79), and number of rings (DOFV ¼ À7.84).Furthermore, a model involving the influence of both H f and M r on l 1 gave a further significant decrease in OFV (DOFV ¼ À23.91).Additional inclusion of a multiplicative interaction term between H f and M r to the latter model did not further improve this model significantly given the use of an additional estimated parameter (DOFV ¼ À26.14).Attempts to include either heavy atom count or number of rings as an additional covariate in this model did not lead to a further significant reduction in OFV according to the LRT.This is likely due to heavy atom count and number of rings being significantly correlated with M r ; hence, they are unable to contribute significantly to further explanation of variability in the data.The correlation matrix for all variables used in model selection is shown in Supplementary Material C. The final selected covariate  model, therefore, included the influence of H f and M r on l 1 .The hazard function for this model is given in Equation 6(where 72.92 is the median H f and 406.56 is the median M r ), the parameter estimates are given in Table 3, and the NONMEM control file is shown in Supplementary Material D.
The estimate of parameter q 1 is positive which indicates that an increase in H f leads to an increase in h(t) and, hence, a decrease in stability of amorphous compounds.This directional influence is entirely consistent with a larger enthalpy of fusion leading to a greater thermodynamic driving force toward crystallization.
Contrarily, q 2 has a negative estimate such that an increase in M r leads to an increase in amorphous drug stability.The positive effects of high M r and low H f on amorphous stability determined here are consistent with the earlier MLR analysis of the same data set 6 and also with other literature reports. 3,9,10The MLR study found that both M r and H f were individually correlated with amorphous stability (R ¼ 0.59 and R ¼ À0.73, respectively) and the model equation involving both covariates is given by Equation 7.
It was previously reported that molecules with high M r often have a complex structure that impedes orientation in a crystal lattice, which leads to higher stability in a disordered, amorphous state. 3Molecules with high H f require more energy to disrupt the crystalline lattice during melting.As melting precedes the formation of the amorphous material, the energy supplied during melting increases the internal energy of the system and, thus, lowers its physical stability. 6he 95% confidence intervals around the parameter estimates in Table 3 were derived from a nonparametric bootstrap procedure which involved fitting of the model to 1000 different data sets of 25 compounds, each sampled (with replacement) from the original data set. 29The distributions of parameter estimates for the 4 model parameters, derived from the bootstrap procedure, are shown in Figure 2, and the 95% confidence intervals in Table 3 were derived from the 2.5 and 97.5 percentiles of each distribution.The 95% confidence intervals around the estimates of l 1 , b 1 , and q 1 all exclude the null value of zero, but the confidence interval around q 2 does include zero by a small margin.However, the inclusion of q 2 in the model does satisfy the requirement of the LRT because the OFV reduces by 10.45 on addition of the influence of M r to a simpler model that only includes the influence of H f .It is likely that adding more compounds to the sample set for a future analysis would lead to a more precise estimate of q 2 .
It is important to recognize that the generation of events using a TTE model is a stochastic process and repeated simulations of event times from the model will be different and will encompass a range of values. 24A suitable method for qualification of such a model is the VPC, which involves repeated simulations of the model followed by a visualization (using KaplaneMeier curves) of the observed events and the range of simulated events. 30This diagnostic plot helps to ensure that repeated simulations of event times from the model are consistent with the observed event times, without evidence of significant bias.A VPC of the selected Weibull covariate model was generated using 1000 repeated simulations of event times for the set of 25 compounds.The results are shown in Figure 3, which indicates a very good concordance between the KaplaneMeier curve (and associated 95% confidence interval) of the observed crystallization events and the median of the 1000 KaplaneMeier curves (and associated 95% confidence interval) derived from simulated data.There is substantial overlap between the 2 sets of confidence intervals, and there is no indication of any significant bias, which could be observed in different relative positions of observed and simulated KaplaneMeier curves over long periods of time.This confirms that the selected Weibull covariate model gives a good description of the experimental data.

Sensitivity Analysis
A sensitivity analysis was performed to investigate and visualize the effect of the selected covariates on the expected median event time for crystallization.The Weibull covariate model hazard function (Eq.6) was analytically integrated to generate the associated survival function (Eq.1). 24The resulting survival function is given by Equation 8SðtÞ where K is given by Equation 9 The expression for S(t) given by Equation 8 can now be plotted for different values of H f and M r to explore the sensitivity of S(t) to ±20% changes in H f and M r around values of 100 J/g and 500 g/mol, respectively.Given that the model contains significant uncertainty in the estimates of the 4 parameters in Equation 6, it was decided that parameter uncertainty should be incorporated into the sensitivity analysis.Hence, for each fixed pair of M r and H f values, S(t) was calculated 1000 times using the 1000 sets of parameter estimates derived from the bootstrap procedure.At each time point, the median of the 1000 S(t) values was calculated along with 80% confidence intervals leading to the plots shown in Figure 4.These plots indicate that 20% changes in both H f and M r lead to highly significant changes in the survival curves for amorphous stability.The sensitivity to changes in H f is greater than changes in M r such that the example ±20% changes in H f lead to 7.9-fold changes in median crystallization time, whereas ±20% changes in M r lead to 4.1-fold changes in median crystallization time.Note that Figure 4b contains significantly wider confidence bands than Figure 4a, and this is related to the greater uncertainty in estimation of q 2 compared with q 1 (Table 3).
If parameter uncertainty is neglected, the sensitivity analysis can be reduced to a simpler format by taking the expression for S(t) in Equation 8, setting it equal to 0.5, and then solving the equation for t, which is equivalent to finding the median crystallization time rather than the entire survival curve including parameter uncertainty for a particular combination of both H f and M r .This leads to Equation 10 The median crystallization time for various combinations of H f and M r was calculated using Equation 10and is shown in Table 4.The explored ranges of H f and M r led to a 17,000-fold range in It should be remembered that the generation of predictions using a TTE model is a stochastic process (even when parameter estimates are very precise), and once a model has been derived, the process of predicting event times involves the drawing of samples from a probability distribution.Hence, the crystallization times listed in Table 4 represent the median expected crystallization times that might be observed for large populations of compounds representing each of the combinations of H f and M r .There will be considerable variation in the predicted crystallization times of particular compounds around each of the median values.This is consistent with reports in the literature that compounds with similar physicochemical properties and chemical structure can demonstrate very different amorphous stability. 8This is also true for the data set studied here, for example, celecoxib (H f ¼ 72.2 J/g, M r ¼ 358.8 g/mol, stable for 6 days) and etoricoxib (H f ¼ 84.1 J/g, M r ¼ 381.4 g/mol, stable for 84 days).

Validation on External Data Sets
A test set of 11 compounds, obtained from the literature 3 and from AstraZeneca, 6 was used to validate the derived Weibull covariate model.The observed crystallization times of the 11 test set compounds are given in Table 5.Because only the time of first detection of crystallization was known for these compounds, without knowledge of the time of the previous observation where crystallization was not detected, the crystallization times are less informative than the interval-censored data used for building the model.The Weibull covariate model was used to predict 1000 sets of crystallization times for the 11 compounds, and the observed and simulated events are displayed in Figure 5 using the same display format as that used for the VPC in Figure 3.This shows considerable overlap of the 95% confidence intervals derived from the observed and predicted data and exhibits a very similar shape of the observed and median predicted KaplaneMeier curves.Therefore, the range of crystallization event times predicted by the Weibull covariate model are consistent with the observed event times, which indicates that the model has the ability to generate useful quantitative predictions  Further details of the distributions of predicted crystallization event times are given in Supplementary Material E. We have previously correlated our experimentally observed amorphous stability values with those observed elsewhere using a range of different methodologies. 6Furthermore, the relative amorphous stability ranges for the compounds in our validation data set (Table 5) are in line with the high (indoprofen 3,20,26 ), intermediate (droperidol 3,26,33 and nifedipine 3,26,34 ), and low (clotrimazole 3,26,34,35 and felodipine 3,26,34 ) crystallization tendencies reported in the literature.
Although the representation of the model's predictive power given in Figure 5 indicates that the overall predictive performance is about as good as it could be expected, some of the individual predictions in Table 5 have large discrepancies.These discrepancies may be related to differences in preparation methods and storage conditions for the test set compounds compared with the compounds used in model building.For instance, felodipine was a member of the 25 compound training set 6 and classified as an unstable compound (measured crystallization within 5 days of storage) and was also a member of the test set where it had been assessed as a fairly stable compound (crystallization within 84 days). 3It may be possible to further improve the model by including additional covariates in the future which have not yet been explored in the present study.
We should note that each of the 11 test set molecules represents a single sample from a corresponding large population of possible molecules, with each population having the appropriate pair of H f and M r values.Each of those populations will exhibit a range of amorphous stability.Furthermore, given the earlier discussion of the stochastic nature of TTE models, we can understand that the predicted crystallization times arise from distributions of values.Therefore, both, the observed and predicted crystallization times, should be viewed as having significant variability.The presence of some test set compounds with large discrepancy between the 2 times is to be expected.The variability in the predictions is tolerable as the derived model is not intended to replace the need for empirical drug stability studies and it would not be used in isolation to make critical decisions related to product development.The model should be rather viewed as a "risk assessment" tool to rank the relative amorphous stability of compounds, for example, in a discovery or early development setting, where there is a requirement from a biopharmaceutical or solid-state consideration for an amorphous drug formulation.
To compare the predictive abilities of the TTE model (Eq.6) and the published MLR model (Eq.7), the values of geometric meanfold error (GMFE) and bias 36 were calculated for both models.These model predictivity metrics are listed in Table 6.The GMFE of the MLR model was 33% higher than the GMFE of the Weibull covariate model, and the bias was lower for the TTE model.This indicates that the predictive performance of the TTE model is better than the MLR model.This is consistent with the fact that a TTE analysis is, for several theoretical reasons, better suited to analysis of event versus time data (particularly given the interval-censored and right-censored nature of the data) than MLR.

Concluding Remarks
In this work, we discussed a novel application of TTE modeling to better understand the influence of physicochemical parameters on measured long-term amorphous stability with censored observations.This study used a previously published set of 25 representative poorly soluble compounds.The best description of the shape of the survival curve for the measured data was obtained with a Weibull hazard model consisting of 2 structural parameters, l 1 and b 1 .After investigating the effect of different physicochemical properties on l 1 and b 1 , the enthalpy of fusion (H f ) and molecular weight (M r ) significantly improved the model statistics.Using sensitivity analysis, it was shown that a decrease in H f and an increase in M r contribute to longer survival times of amorphous compounds and that amorphous stability depends more strongly on H f than on M r .The Weibull covariate model was used to calculate the median survival times for different values of H f and M r .The resulting survival times can serve as a useful indication of how amorphous drug stability typically responds to changes in H f and M r but with the realization that significant discrepancies are to be expected.Finally, the Weibull covariate model was tested on an external data set of 11 compounds and showed superior predictive power compared with a previously published MLR model, which also has a dependence on H f and M r .Because the 2 used covariates are readily accessible, through calculation and by means of differential scanning calorimetry, they may be of interest to the pharmaceutical industry in time-and cost-effective assessment of compound stability on storage and development of amorphous products.However, it will be beneficial to study additional compounds in the future to increase the data set and consequently reduce the uncertainty in the model parameter estimates.

Figure 1 .
Figure 1.Visual assessment of candidate base models.KaplaneMeier curve of observed data is shown as continuous black line with 95% confidence intervals shown as black dashed lines.KaplaneMeier curves of estimated survival functions derived from 4 candidate base models are shown as colored lines.

Figure 3 .
Figure3.VPC of the selected covariate model showing KaplaneMeier curves of observed events (black line) and median of simulated events (red line).The dashed lines indicate the 95% confidence intervals of observed events and the gray shading show the 95% confidence interval of the simulated events.

Figure 2 .
Figure 2. Distributions of parameter estimates derived from bootstrap analysis of Weibull covariate model.

Figure 4 .
Figure 4. Sensitivity analysis showing differences in survival curves on changes in the model covariates.(a) M r is constant and H f changes by ±20%.(b) H f is constant and M r changes by ±20%.The colored areas show 80% confidence intervals based on percentiles.

Figure 5 .
Figure 5. Assessment of predictive performance of the Weibull covariate model on 11 test set compounds.

Table 1
Hazard Functions Tested During Selection of the Base Model

Table 2
Parameter Estimates and OFV From Base Models a Standard errors of parameter estimates, derived using NONMEM covariance method.b Covariance procedure did not complete successfully.

Table 3
Parameter Estimates of the Selected Weibull Covariate Model These results further demonstrate the large influence that H f and M r are expected to have on the amorphous stability of compounds with different properties.

Table 4
Median Crystallization Times (days) for Different Combinations of H f (J/g) and M r (g/mol) Calculated Using Equation10and Parameter Estimates From Table3 Median crystallization times longer than 168 days are shaded in green (very stable compounds) and shorter than 1 week are shaded in red (very unstable compounds).

Table 5
Amorphous Stability Data for 11 Compounds Along With Prediction From MLR and TTE Models.Amorphous Stability Predicted With TTE Is Represented as the Median of the Simulated Event Times.Event Times Were Generated Up to 168 Days

Table 6
Metrics Comparing Predictive Performance of TTE and MLR Models When Tested on 11 External Compounds